© The Authors 2024

1 Introduction

While the Voyager Infrared Interferometer Spectrometer and Radiometer (IRIS) paved the way for thermal exploration of the outer Solar System in the spectral range 1.4–170 µm, the Cassini Composite Infrared Spectrometer (CIRS) opened up a new field for studying the thermal properties of Saturn’s rings and the surface of icy moons through its still unique broad spectral coverage from 7 to 1000 µm, which made it possible to observe most of their spectrum and wavelength of maximum thermal emission (Flasar et al. 2004).

The brightness temperature of Saturn’s rings has been known for a long time to dramatically decrease by tens of kelvins in between 20 µm and millimetre wavelengths (Esposito et al. 1984). The focal plane FP1 of CIRS was the first to fully cover the spectral region 16.7–1000 µm (10–600 cm⁻¹) and accurately described this decrease for wave numbers w_n < 50 cm⁻¹, identified as a roll-off in emissivity (Flasar et al. 2004; Spilker et al. 2005). It was indeed suspected to result from the increasing transparency of water ice in this range and to depend on the size distribution of ring particles or of regolith grains with which they might be covered. Assuming ring particles are made of pure water ice, Spilker et al. (2005) studied the spectra of the lit side of the three main rings, A, B, and C, acquired in 2004, early in the mission. All exhibited similar roll-offs, and a simple model of emissivity, assuming regolith-free ring particles following a power-law size distribution, was shown to reasonably reproduce the observed spectra, and no exact adjustment was proposed. Morishima et al. (2012) analysed all FP1 spectra obtained with the lowest spectral resolution (15.5 cm⁻¹) before the equinox in 2010 and showing ring temperatures above 70 K in order to reject noisier spectra. Inverting a diffraction-corrected hybrid Mie-Hapke emissivity model of ring particles covered with a regolith of pure water-ice grains with these data, they found the emissivity to depend on temperature, most probably reflecting that grains of various sizes have different temperatures depending on illumination and viewing angles. The size distribution of regolith grains was found to be broad, ranging from 1 µm to 1–10 cm, with a power-law index of p ~ 3, almost similar and independent of temperature in the two main thickest A and B rings.

CIRS observations of Saturn’s icy moons have also mostly been analysed assuming that their surface emissivity is equal or close to unity. The aim was then to determine the space-time distribution of surface temperatures to derive maps of the regolith thermal inertia. The discovery of thermal anomalies on the surface of icy moons has certainly made an essential contribution to understanding the differential effects of space weathering between their leading and trailing hemispheres (Howett et al. 2011, 2012; Schenk et al. 2018; Paranicas et al. 2014; Nordheim et al. 2017). Carvano et al. (2007) provided the only but qualitative study of the emissivity of Phoebe, Iapetus, Enceladus, Tethys, and Hyperion in the range 50–400 cm⁻¹, that is, above the roll-off region. The lack of spectral signatures in the emissivity of these moons in the dataset available at this early time in the mission was interpreted as essentially due to the very high porosity (>95%) of clumps of small grains that may cover the surface. In the case of Tethys, an intimate mixture of water ice with amorphous carbon (AmC) could also suppress spectral features if ice grains are larger than contaminant grains, or if the fraction of AmC were at least 75%. Nor did Howett et al. (2016) detect any spectral signature in Rhea’s spectra captured in 2013 and 2015 flybys in a similar spectral range.

It has long been known that crystalline water ice is the main component of Saturn’s rings and the surfaces of the icy moons. Visible and near-infrared (NIR) ring spectra could be satisfactorily reproduced by a mixture of regolith grains with typical sizes from tens of micrometres to millimetres, the water ice being slightly contaminated (a few percent) by still-debated tholins (Th) and AmC or silicate inclusions to reproduce observed spectral reddening and darkening, respectively (Poulet et al. 2003; Cuzzi et al. 2018). Low absorption by CO₂ has also been observed for some of these moons (Clark et al. 2008; Filacchione et al. 2022).

Finally, the analysis of the huge CIRS dataset on Saturn’s rings and moons is still incomplete, both in terms of seasonal effects over the 13 yr that the mission lasted and on spectral behaviours, in particular, that of emissivities (Spilker et al. 2018; Schenk et al. 2018). It is still not known how much the emissivity varies in the spectral range of FP1 due to contaminants, in particular, in the spectral regions where water ice becomes transparent, and how this may modify the observed thermal emission and the retrieved temperatures. Other instruments have taken over from Earth orbit to observe the thermal emission from icy worlds in spectral domains related to that of CIRS. The Atacama Large Millimiter/submillimeter Array (ALMA) within its receiver bands 7 to 10 can explore the sub-millimeter range down to 320 µm (950 GHz), that is, in between 10 and 32 cm⁻¹, and the Mid-Infrared Instrument (MIRI) of the James Webb Space Telescope (JWST) can probe their thermal emission in between 4.9 and 28 µm (357–2040 cm⁻¹). Upcoming space instruments will also target the Galilean satellites, for instance the Submil-limetre Wave Instrument (SWI) on board the JUpiter ICy moons Explorer (JUICE) around 520 µm (17.7–21 cm⁻¹) and 250 µm (36–42 cm⁻¹), or the Europa Thermal Emission Imaging System (E-THEMIS) on board Europa-Clipper in the range 7-80 µm (125–1429 cm⁻¹; Christensen et al. 2024). These instruments will therefore also explore these spectral regions where water ice can be transparent and the emissivity is sensitive to both the presence of contaminants and the size distribution of the regolith grains.

In the present study, we estimate quantitatively how much the hemispherical emissivity of a regolith covering ring particles or surface of moons is sensitive, in the wavelength range 5–1000 µm (10–2000 cm⁻¹), to the size distribution of grains, the type of contaminant, its fraction, and how it is mixed with water ice. The influence of temperature on the optical properties of water ice is also considered. In section 2, a hybrid Mie-Hapke model for regolith hemispherical emissivity is proposed that includes these factors together with the diffraction correction of Mie calculations and an asymmetry factor of the regolith. A global sensitivity analysis yields quantitative estimates of the importance of these factors, which may be used to reduce their number depending on the spectral range of the data on which the model is inverted. Its results are presented and discussed in Sect. 3.

2 Mie-Hapke hybrid model

Mie-Hapke hybrid models for regoliths couple the formalism developed by Hapke (2012) to estimate their reflectance and emittance based on the scattering and absorption efficiencies of spherical grains provided by the Mie theory. These models require knowledge of the optical properties such as the complex refractive index or dielectric constant of the materials involved and must include a correction for the independent scattering hypothesis of Mie’s theory to be suitable for the closely packed medium of regoliths.

2.1 Composition of grains

Within Saturn’s magnetosphere, which is bombarded by plasma and highly energetic particles, the surface chemistry ofrings and icy moons is almost certainly complex. Radiolysis and photolysis may indeed lead to a mixture that is rich in various icy components, such as nitrogen oxides (Delitsky & Lane 2002). We limit ourselves here to a mixture of water ice and contaminants capable of reddening and darkening the spectrum, such as the proposed Th and AmC. A more complete sensitivity analysis would require refractive indices of complex icy mixtures on the wide spectral range of FP1, which are not available.

2.1.1 Water ice

Several authors have produced measurements of the complex refractive index n = n_r + jn_k of water ice at low temperatures between 10 and 140 K, very few of which cover the full spectral range 10–2000 cm⁻¹ (Fig. 1). In mid- and far-infrared domains, the dominant absorption mechanisms are bending mode at 6 µm (1666 cm⁻¹), molecular hindered rotation at 12 µm (833 cm⁻¹), and lattice vibrations at 45 and 65 µm (222 and 154 cm⁻¹), which translate into a high imaginary part n_k of the refractive index (Fig. 1 bottom). In these specific ranges, the real part n_r oscillates between 1 and 2. At sub-millimetre wavelengths, about the roll-off in n_k (below 50 cm⁻¹), n_r is about constant and ice becomes transparent (Warren & Brandt 2008; Warren 2019, Fig. 1 top). The optical properties depend on the temperature and phase of the ice, either crystalline (Ih) or amorphous (Am), the latter irreversibly transforming into the former at temperatures between 125 and 150 K (e.g. Jenniskens & Blake 1996).

The optical constants for crystalline water ice that have long been the most frequently used result from a compilation of measurements by Warren (1984), covering the entire range from UV to mm, and updated for a temperature of 266 K (Warren & Brandt 2008), that is, well above the temperatures of interest here (Fig. 1). Hudgins et al. (1993) produced measurements at 10, 100, and 140 K for water ice from 50 to 4000 cm⁻¹ and reported an apparent phase transition between 100 and 140 K. They made available additional data at intermediate temperatures (40, 80 and 120 K) from 500 to 4000 cm⁻¹ . Some divergence in the range 400–600 cm⁻¹ can be seen between these data that corresponds to a change in detectors and a knitting that operates for 10, 100, and 140K at the junction of spectra recorded both in the mid- and far-infrared. The estimates above 455 cm⁻¹ for amorphous ice by Mastrapa et al. (2009) are consistent with Hudgins et al. (1993) for comparable temperatures, but remain of marginal spectral coverage. In the most absorptive region (80–300 cm⁻¹), the optical constants by Curtis et al. (2005) are comparable to those of Hudgins et al. (1993). Elsewhere, they differ significantly. n_k estimated by Hudgins et al. (1993) appears to be either smaller or larger in between 300 and 600 cm⁻¹ and noisier below 100 cm⁻¹. They are consistent above 500 cm⁻¹ and at temperatures of 40, 80, and 120 K. For both sets of measurements, n_k always appears to increase slightly with decreasing temperature.

Mätzler et al. (2006) summarized experimental and theoretical work on the dielectric properties of Ih ice between 1 MHz and 1 THz (≥300 µm, ≤ 33 cm⁻¹), in particular, that of Hufford (1991). The latter formulated the frequency and temperature variations of the dielectric constant from measurements and the theoretical work by Mishima et al. (1983) for temperatures in between 80 and 202 K. We have used this synthesis to deduce the dependence of the indices n_r and n_k on wave number between 10 and 33 cm⁻¹ and temperature. n_r can be considered as constant with temperature, while n_k (Fig. 1, bottom) is significantly impacted by temperature. This same formalism at millimetre wavelengths was included by de Kleer et al. (2021) in their thermal model to analyse ALMA observations of Ganymede. The emissivity was calculated from the complex dielectric constant, which depends on porosity and on a meteoritic dust fraction of 20 % mixed with ice at the molecularscale, and on temperature. It was found to be 0.78 ± 0.04 at 343 GHz (0.87 mm).

As the crystalline phase is observed at the very surface, we chose, just like Morishima et al. (2012), to merge the Curtis et al. (2005) data at 136 K to the sub-millimetre model as described above, through a cubic-spline extrapolation in between 33 and 50 cm⁻¹. The sub-millimetre model for n_k(w_n, T₀) can be calculated for any ice temperature T₀ (Fig. 2 bottom). Between 566 and 2000 cm⁻¹, Hudgins et al. (1993) data at 120 K for Am ice phase were chosen, after being smoothed to remove some noise in the transparency window near 9.5 µm (1050 cm⁻¹). These data are indeed very close to the Curtis et al. (2005) measurements at 136 K for the Ih phase at 530 cm⁻¹, and they are almost identical above 700 cm⁻¹ to the Ih phase at 140K. This merging constitutes the Ih-T₀ optical properties used in the following (Fig. 2). This illustrates the difficulty of finding a complex refractive index for the Ih ice phase versus temperature over the full range 10–2000 cm⁻¹ in order to conduct a sensitivity analysis of emissivity versus temperature. The sub-millimetre range, about the roll-off (≤33 cm⁻¹), is the only range where this can be studied because it has been modelled. Therefore, no temperature dependence of optical properties was considered above 33 cm⁻¹.

Fig. 1 Complex refractive index of water ice, either crystalline (Ih) or amorphous (Am), as proposed by various authors at different temperatures. (top) Real part nr and (bottom) imaginary part nk. Legends correspond to authors T (K) ice phase with authors C05 (Curtis et al. 2005), M83 (Mishima et al. 1983), H93 (Hudgins et al. 1993), M09 (Mastrapa et al. 2009), and W08 (Warren & Brandt 2008).

Fig. 2 Complex refractive index of Ih water ice and contaminants. (top) Real part nr and (bottom) imaginary part nk of the Ih. Ih T0 point to the refractive index of Ih water ice at temperature T0 (see text). The optical constants (nr, nk) for an intra-particle mixture of Ih 90 ice with AmC (Preibisch et al. 1993) or Th (Khare et al. 1984) at 3% and 1% levels respectively are also plotted.

2.1.2 Contaminants

Poulet et al. (2003) adjusted a B ring spectrum in the NIR domain with 10-to-1300 µm sized grains of water ice with 0.74% intra-particle Th inclusions mixed with 2.3% of 10 µm sized AmC grains. Filacchione et al. (2012) also selected an intraparticle mixture of ice with a Th fraction lower than 0.5% to model NIR spectra of Mimas, Enceladus, Dione, Tethys, and Rhea obtained by the Visual and Infrared Mapping Spectrometer onboard Cassini, which they combined with an additional intimate (volumic) mixture of AmC grains to reproduce the reddening and darkening of their spectrum. The NIR spectra of Mimas yield a satisfactory fit with the Hapke theory with 0.1% of Th and 2% of AmC inclusions (Filacchione et al. 2012). Potential contaminants such as these were therefore chosen for this study.

The optical constants of AmC were taken from Preibisch et al. (1993). They are consistent with many other estimates (Zubko et al. 1996). The optical constants for Th originate from Khare et al. (1984) as they remain those that cover our spectral range of interest (Brassé et al. 2015). Both contaminants exhibit significantly different optical properties, but no significant spectral features in the region of interest except for Th at about 6 µm (1666 cm⁻¹; Fig. 2). They are more absorbent than water ice. The effect of mixing water ice at the molecular scale with 1% of Th increases the absorption, but does not alter the real part of the refractive index. The effect is mostly significant in the region of the roll-off and is slightly lower than the mixing effect with 3% of AmC (Fig. 2). The sensitivity to the contaminant fraction will have to be observed in areas in which water ice is relatively transparent, that is to say, in the ranges 10–50 cm⁻¹, 300–600 cm⁻¹, and 900–1300 cm⁻¹. Other types of mixtures will also be considered to quantify the sensitivity of emissivity to the fraction of contaminants (Sect. 2.4).

2.2 Hemispherical emissivity of a regolith

Hapke (2012) proposed the model for isotropic multiple scattering approximation (IMSA) for the hemispherical emissivity ɛ_h(w_n, a) of a regolith of grains with size a. In case of isotropic scattering, it is written at wave number w_n as

$${\epsilon }_{\text{h}}\left(w_{\text{n}},a\right)\sim \frac{2\gamma \left(w_{\text{n}},a\right)}{1+\gamma \left(w_{\text{n}},a\right)}\left(1+\frac{r_0\left(w_{\text{n}},a\right)}{6}\right),$$(1)

where

$$r_0\left(w_{\text{n}},a\right)=\frac{1-\gamma \left(w_{\text{n}},a\right)}{1+\gamma \left(w_{\text{n}},a\right)}$$(2)

and

$$\gamma \left(w_{\text{n}},a\right)=\sqrt{1-{\omega }_0\left(w_{\text{n}},a\right)}$$(3)

as a function of the single-scattering albedo of the grain ω₀(w_n, a), provided possibly by Mie’s theory. The hemispherical reflectance is approximated as follows:

$$r_{\text{s}}\left(w_{\text{n}},a\right)\sim 1-{\epsilon }_{\text{h}}\left(w_{\text{n}},a\right)\sim r_0\left(w_{\text{n}},a\right)\left(1-\frac{\gamma \left(w_{\text{n}},a\right)}{3\left(1+\gamma \left(w_{\text{n}},a\right)\right)}\right).$$(4)

2.3 Diffraction-corrected Mie scattering albedos

The estimation of single-scattering albedos ω₀(w_n, a) by Mie’s theory includes the diffraction pattern: Grains larger than the wavelength (with a size parameter x = 2πa/λ greater than a few units) have an important asymmetry factor ξ_a =< cosθ >~ 1, a single-scattering albedo ω₀(w_n, a) ~ 0.5 because their scattering and extinction efficiencies are Q_S ~ 1 and Q_ext ~ 2, respectively (Fig. 3, top). Within a close-packed medium, diffraction peaks of grains interfere, so that their scattering efficiency and the asymmetry factor of their phase function should differ from Mie’s predictions. According to Wald (1994), independent scattering in regolith may be approached for transparent large particles when the scattering lobe and diffraction peak coincide. He proposed a diffraction-corrected single-scattering albedo, which is valid for the mid-infrared domain (2–14 µm) and absorbing snow particles much larger than the wavelength, that is, ≥ 50µm in their case. García-Santos et al. (2016) reviewed efforts made in Earth science to validate compactness correction methods in emissivity models of snow and mineral surfaces in the spectral range 8–14 µm. The Mie-Hapke hybrid model corrected for diffraction with various methods (static structure factor; Wald (1994) and delta-Eddington (see below) corrections) was tested on emissivity measurements of quartz and gypsum samples, with grain sizes ranging from 15-to-277 µm, and against another model. The Hapke model appeared to perform best with the delta-Eddington correction for sizes a >75 µm. For smaller grain sizes, the correction according to Wald (1994) gave better results. Given the wavelength range of interest here (5–1000 µm) and the possible range of grain sizes considered (1 µm to 5 cm, Sect. 2.4), a correction according to Wald (1994) appears partially irrelevant here because for a given size, the scattering regime might indeed change from Rayleigh (x<< 1) to geometrical optics (x>>1).

Carvano et al. (2007) chose not to include diffraction correction in their model, following results from Moersch & Christensen (1995) established on quartz regoliths in the spectral range of 400–1400 cm⁻¹ (7–25 µm) with grain sizes smaller than 250 µm. Morishima et al. (2012) chose the delta-Eddington approximation (δ-Edd) proposed by Joseph et al. (1976) to estimate the single-scattering albedo ω_0,δ–Edd(w_n, a) after diffraction correction,

$${\omega }_{0,\delta -\text{Edd}}\left({\omega }_{\text{n}},a\right)=\frac{1-{\xi }_a^2}{1-{\xi }_a^2{\omega }_0},{\omega }_0$$(5)

where ξ_a ≡ ξ(w_n, a) is the asymmetry factor, and ω₀ ≡ ω_Mie(w_n, a) for a grain of size a as calculated with Mie’s theory. They mentioned that their final results might depend on the diffraction correction method.

For this study, the equivalent slab approximation (esa) as developed by Hapke (2012) is proposed for comparison. It estimates for high x values the diffraction-corrected single-scattering albedo ω_0,esa by calculating the non-diffracted component of the scattered light q_s = Q_S − Q_d and by assuming the extinction efficiency Q_e = 1. This approximation was implemented here as the esa diffraction correction, so that

$${\omega }_{0,\text{ esa }}\left(w_{\text{n}},x\ge x_0\right)=\frac{q_{\text{s}}}{Q_e}\text{ and }{\omega }_{0,\text{ esa }}\left(w_{\text{n}},x<x_0\right)\equiv {\omega }_{\text{Mie}}\left(w_{\text{n}},a\right),$$(6)

where x₀ = 2πaw_n0 corresponds to the minimum wave number w_n0 at which ω_0,esa < ω_Mie with ω_Mie > 0.45 (Fig. 3, top). Mie calculations were performed with the PyMieScatt Python library (Smulin et al. 2018). The albedos and hemispherical emissivities were calculated for 41 sizes ranging between 1 µm and 5 cm (Fig. 3). Above this size, ε_h(w_n, a) is nearly constant in wave number and size. For micrometre-sized grains, all models provide the same results as long as the size parameter x < 0.2, corresponding to the Rayleigh regime. For a = 10 µm, the diffraction correction methods differ significantly above 300 cm⁻¹ as water ice becomes less absorbing and x becomes large (Fig. 3, bottom and Fig. 2). Above a = 100 µm, the effects of diffraction corrections are clearly visible, with ω₀ well below 0.5 for the high values of x, above 100 cm⁻¹ (x scales with w_n). As the grain size grows, this is even truer at lower wave numbers. The roll-off in emissivity is well reproduced at w_n ~ 50 cm⁻¹ for grain sizes ranging from a = 100 µm to a few millimetres (Fig. 3, bottom left). The diffraction correction method has some impact on its spectral position. These calculations are compatible with the studies of Wald (1994), Moersch & Christensen (1995), or Morishima et al. (2012).

Fig. 3 Single-scattering albedo and hemispherical emissivity as a function of grain size and wave number or wavelength. (top) Single-scattering albedo ω0(wn, a) of a grain with size a (bottom) hemispherical emissivity εh(wn, a) of regolith covered with these grains made of pure crystalline water ice Ih90. Calculations with the Mie theory (dotted line) or after diffraction corrections such as delta-Eddington (dash-dotted line) or equivalent slab approximation (solid line) are displayed. (Left) 10–700 cm−1 range and (right) 5–14.3 µm range (or 700–2000 cm−1).

2.4 Mixing with contaminants and the effect of the size distribution

Filacchione et al. (2012) studied the absorption bands of water ice at 1.5 and 2.0 µm in spectra of the B ring, from which they derived grain sizes ranging between 40 µm and 0.3-to-0.6 cm, assuming pure water ice. Poulet et al. (2003) adjusted a B ring spectrum with 10-to-1300 µm sized grains of water ice with intra-particle inclusions of 0.7% of Th intimately mixed with 2% of 10 µm sized AmC grains. Grain sizes of 1 µm to 5 cm were therefore considered hereinafter together with mono-or poly-dispersed power-law size distributions. As often in NIR reflectance spectral modelling, the effect of any size distribution has been little explored for icy bodies in the thermal infrared, except for Morishima et al. (2012).

Three different types of mixture of water ice with contaminants were considered here: at the molecular level, the intra-particle mixture, secondly as an aerial mixing, and finally as an intimate (volumic) mixture with the same (distribution of the) grain size. The resulting emissivity for an intra-particle mixture was calculated by coupling the optical constants of water ice and the contaminants at a volume fraction f with the Maxwell-Garnett effective medium theory and following Eqs. (1)–(3). The hemispherical emissivity of an aerial mixture of a contaminant covering a fractional area f within the field of view reads

$${\epsilon }_{\text{h}}\left(w_{\text{n}},a_{\text{i}},a_{\text{p}}\right)=(1-f){\epsilon }_{\text{h},\text{i}}\left(w_{\text{n}},a_{\text{i}}\right)+f{\epsilon }_{\text{h},\text{p}}\left(w_{\text{n}},a_{\text{p}}\right),$$(7)

where ε_hi(w_n) and ε_h,p(w_n) are the hemispherical emissivities of the icy and contaminant areas, respectively, either made of the same size or not. In the case of an intimate mixture with a volume fraction f and scattering or extinction efficiencies Q_s,i, Q_e,i, or Q_s,p, and Q_e,p for the icy or the contaminant populations, respectively, the total single-scattering albedo is

$${\omega }_0\left(w_{\text{n}},a_{\text{i}},a_{\text{p}}\right)=\frac{(1-f)\pi a_{\text{i}}^2{Q}_{s,i}\left(w_{\text{n}},a_{\text{i}}\right)+f\pi a_{\text{p}}^2{Q}_{s,p}\left(w_{\text{n}},a_{\text{p}}\right)}{(1-f)\pi a_{\text{i}}^2{Q}_{e,i}\left(w_{\text{n}},a_{\text{i}}\right)+f\pi a_{\text{p}}^2{Q}_{e,p}\left(w_{\text{n}},a_{\text{p}}\right)},$$(8)

where ai and ap are their grain sizes, respectively.

Figure 4 displays the hemispherical emissivities of regolith with grain sizes ranging over three orders of magnitudes as a function of contaminant fraction and mixing type when a_i = a_p = a in the case of esa diffraction correction. The results are comparable with the δ-Edd correction. Grains of 10 micrometres have a high emissivity below 300 cm⁻¹ (Fig. 4, top left). Above this, ε_h decreases and appears to be mainly dependent on the AmC fraction of 5% when intra-mixed or, to a lower extent, if intimately mixed (vol). It is almost insensitive to contaminants for an aerial mixing at this fraction level. The impact of 1% Th mixing is negligible. For grains that are ten times as large (100 µm, Fig. 4, middle left), the roll-off in emissivity below 50 cm⁻¹ shows significant variation with mixing type (intra or vol) in the case of AmC at 5% fraction. Aerial mixing remains of low impact, just like any type with 1% Th. AmC remains influential above 300 cm⁻¹ . For larger grains (Fig. 4, bottom left), the impact of contaminants is restricted to the roll-off region and the 9.5 µm transparency (Fig. 4, bottom right) only for a few percent of intra-mixed contaminants. The hemispherical emissivity is then constant above 40 or 50 cm⁻¹ (geometrical optics regime). The contamination by a maximum fraction of 1% Th appears to be comparatively negligible for the hemispherical emissivity of water ice in the roll-off region.

For a power-law size distribution of grains following n(a) da = n₀ a^−p da with minimum and maximum sizes a_min and a_max, the scattering albedo of the size distribution reads

$${\omega }_{0,\text{sd}}\left(w_{\text{n}}\right)=\frac{{{{\displaystyle \int }}^{\text{}}}_{a_{\min }}^{a_{\max }}\pi a^{2-p}{Q}_{\text{sca}}\left(w_{\text{n}},a\right)\text{d}a}{{{{\displaystyle \int }}^{\text{}}}_{a_{\min }}^{a_{\max }}\pi a^{2-p}{Q}_{\text{ext}}\left(w_{\text{n}},a\right)\text{d}a}=\frac{Q_{\text{sc}a,\text{sd}}\left(w_{\text{n}}\right)}{Q_{\text{ext},\text{sd}}\left(w_{\text{n}}\right)}.$$(9)

Figure 5 illustrates the dependence of ε_h(w_n) on p, a_min, and a_max in the case of an intra-mixture. For a steep slope p=3, ɛ_h(w_n) is more sensitive to the minimum grain size a_min than to a_max in the roll-off region (Fig. 5, left). Other transparency regions (300–600 cm⁻¹ and 7–11 µm) are equally sensitive to both extreme sizes. For p = 2, ɛ_h(w_n, a) is insensitive to sub-millimetre a_min sizes, but strongly affected by the choice of maximum size a_mаx within the roll-off region, where this size is comparable to the wavelength. In the case of δ-Edd diffraction correction, the impact of a_max is higher in the roll-off region as the emissivity variation with large grain sizes is higher in this case (Fig. 3).

Finally, in the case of an intimate mixture, the scattering albedo is defined according to

$${\omega }_{0,\text{sd}}\left(w_{\text{n}}\right)=\frac{(1-f)Q_{\text{sca},\text{sd},\text{i}}\left(w_{\text{n}}\right)+f{Q}_{\text{sca},\text{sd},\text{p}}\left(w_{\text{n}}\right)}{(1-f)Q_{\text{ext},\text{sd},\text{i}}\left(w_{\text{n}}\right)+f{Q}_{\text{ext},\text{sd},\text{p}}\left(w_{\text{n}}\right)}.$$(10)

Fig. 4 Hemispherical emissivity ɛh(wn, a) of regolith grains with various sizes a. (top) 10 µm, (middle) 100 µm, (bottom) 1 mm, and various types of mixing: no contaminants (Ih90), with 5% AmC (black) or 1% Th (red) either mixed at intra-particle scale (dotted line), or with aerial (dashed line) or volumic (dash-dotted line) mixing. esa diffraction correction is assumed. Left: 10–700 cm−1 range and right: 5–14.3 µm range (700–2000 cm−1).

2.5 Average asymmetry factor

Finally, in the case of anisotropic and multiple scattering within the regolith and for hemispherically averaged quantities, the solution of radiative transfer equations can be approximated by applying the similarity principle, which consists of replacing ω₀ by

$${\omega }_0^{\ast }\left(w_{\text{n}}\right)=\frac{1-\xi }{1-\xi {\omega }_0\left(w_{\text{n}}\right)}{\omega }_0\left(w_{\text{n}}\right),$$(11)

where ξ is the regolith average asymmetry factor (Hapke 2012). The hemispherical emissivity ${\epsilon }_{\text{h}}^{\ast }\left(w_{\text{n}}\right)$ can be deduced with equations Eqs. (1)–(3) and Eq. (11). If ξ < 0, the regolith is more retro-diffusive and therefore less emissive, and if $\xi >0,{\epsilon }_{\text{h}}^{\ast }\left(w_{\text{n}}\right)$ will increase compared to the isotropic case (Fig. 6). With ξ varying from −0.6 to 0.6, ${\epsilon }_{\text{h}}^{\ast }\left(w_{\text{n}}\right)$ may vary by about ±0.1−0.2.

3 Global sensitivity analysis and discussion

Given the many factors of the model (Table 1), their a priori uncertainties, and their potential interactions on the value of ${\epsilon }_{\text{h}}^{\ast }\left(w_{\text{n}}\right)$, a global sensitivity analysis was carried out to quantitatively determine their importance at any wave number with a view to reducing their number (and simplify the model) before the data analysis. This analysis was conducted with the Sobol method (Sobol 2001; Saltelli et al. 2010) of the SAlib library¹. It quantifies how the total variance on ${\epsilon }_{\text{h}}^{\ast }\left(w_{\text{n}}\right)$ (the model output) at a given wave number can be apportioned to input factors uncertainties. It provides two coefficients for a factor X_i, the first-order sensitivity S_i(w_n), and the total effect sensitivity ST_i(w_n). S_i(w_n) scales the expected relative reduction of the output variance that would be achieved if the factor were fixed. ST_i(w_n) is the expected relative variance of ${\epsilon }_{\text{h}}^{\ast }\left(w_{\text{n}}\right)$ that would be left if all factors but X_i could be fixed. S T_i(w_n) ~ 0 is a necessary and sufficient condition for factor to be non-influential while S_i(w_n) ~ 0 is not. In the case of non-linear models, the factors may have small main effects and large interactions with other factors, so that S_i(w_n) is small and ST_i(w_n) is large. This was observed in some cases, in particular, for models assuming a single grain size. In order to decide which factors are ultimately relatively insignificant for the model we present, the total sensitivity coefficient S T_i(w_n) was therefore chosen.

Fig. 5 Hemispherical emissivity εh(wn, a) as a function of p=2 (dashed line) or p=3 (solid line) and (top) amin or (bottom) amax. The grain composition is water ice Ih90, intra-mixed with 5% AmC. esa diffraction correction is assumed. Left: 10–700 cm−1 range and right: 5–14.3 µm range (700–2000 cm−1).

Fig. 6 Hemispherical emissivity ${\epsilon }_{\text{h}}^{\ast }\left(w_{\text{n}}\right)$ as a function of average asymmetry factor ξ for a size distribution of gains with amin = 1 µm , amax = 5 mm, and p=2 (dashed line) or p=3 (solid line). The grain composition is water ice Ih90 with intra-mixing of 5% AmC. esa diffraction correction is used.

3.1 Sensitivity versus mixing type and diffraction correction

Figure 7 displays the total effects in the case of intra-mixture with an AmC contaminant. In the case of a mono-dispersed size distribution (Fig. 7, top), the grain size a is the most important factor, except in the region of strong ice absorption, where the average asymmetry factor ξ has a similar impact. Both fraction f and temperature T₀ are negligible, <0.05 in the roll-off region for f and <0.001 for T₀. f might be conserved if noise is limited in the FIR or in transparency windows. The impact of ξ is more extended spectrally if diffraction is not corrected for (Mie). In the case of a power-law size distribution (Fig. 7, bottom), T₀ remains of the lowest and negligible importance. In the roll-off region, all factors a_min, a_max, p, ξ, and f have comparable influence. Above 50 cm⁻¹, the variance is mainly sensitive to ξ in the ice absorption regions and to p and a_min elsewhere. The analysis shows that p and a_min interact above 50 cm⁻¹ as their S_i(w_n) are twice as low as their S T_i(w_n). If the a priori knowledge of the factors are those set out, T₀ can be dropped (i.e. fixed to some value). In the case of a power-law size distribution, all factors are influential in some specific spectral regions and should be conserved if the data inversion is conducted over the whole range. Some may be dropped (fixed), depending on the region of fit. The results do not change in the case of Th contamination, but in the case of a sub-percent fraction, the impact of fwill be even smaller in the roll-off region. The behaviour changes little between the esa or δ-Edd diffraction correction methods, as expected.

These same coefficients are shown versus type of mixing in Fig. 8. In the case of a single grain size (Fig. 8 top), the size a is dominant on the variance regardless of the mixing type. The aerial mixture is equivalent to the intra-mixture with a dominant size of pure icy (a_i) grains; the size of pure AmC grains (a_p) is unimportant, ST_i(a_p) < 0.01. With vol mixing, the size of pure AmC grains (a_p) cannot be dropped. For all mixing types, the asymmetry factor ξ keeps its maximum influence in the most absorbing region (100–300 cm⁻¹, 700–900 cm⁻¹, and about 1660 cm⁻¹). The S T_i for f is largest in the roll-off region, but stays <0.01 for both aerial and intimate mixtures. For T₀, S T_i < 0.001 for intra-mixture and <0.01 for both the aerial and intimate mixtures. In the case of a power-law size distribution (Fig. 8 bottom), the sensitivity to factors does not vary much with mixing type, but it does vary in the roll-off region. We note that the impact of f is relatively stronger within the absorption regions in case of aerial and vol mixture compared to intra-mixture, which is greater in transparency windows. The coefficients for aerial and vol mixing can hardly be distinguished. The overall behaviour remains that a_min, a_max, p, and ξ are of somewhat equally important in the roll-off region and in the 900–1100 cm⁻¹ range, whereas a_min, p, and ξ are most influential above 50 cm⁻¹.

Fig. 7 Total effect sensitivities S Ti in the case of an intra-mixture with AmC contaminant and Mie (dotted line), δ-Edd (dash-dotted line) or esa (solid line) diffraction correction methods and (top) mono-dispersed (bottom) power-law size distributions. The dashed line marks the 3% level.

Fig. 8 Total effect sensitivities STi in the case of contamination with an AmC mixture, either intra (dotted line), aerial (dashed line), or volumic (solid line) mixing for (top) mono-dispersed (bottom) power-law size distributions. The esa diffraction correction is used here.

3.2 Discussion

The global sensitivity analysis shows that given their a priori range of uncertainties (Table 1), the most important factors acting on the hemispherical emissivity ${\epsilon }_{\text{h}}^{\ast }\left(w_{\text{n}}\right)$ remain the size distribution of regolith grains and the way in which regolith scatters around infrared light, translating into the average asymmetry factor ξ and the diffraction correction method (Fig. 8). The factor ξ may be best constrained where the effect of sizes is minimum, that is, in the ice absorption regions. Regardless of the mixing type, the index p or minimum size a_min may be constrained over the spectral range above the roll-off (>50 cm⁻¹), where ice is less absorbing. In the case of a mono-dispersed size distribution, the constrained size will be that of the intra-mixed grains, while the size of contaminants may be constrained from spectral observations above 100 cm⁻¹ if vol mixing is assumed (Fig. 8 top). The hemispherical emissivity ${\epsilon }_{\text{h}}^{\ast }\left(w_{\text{n}}\right)$ in the region of the roll-off is sensitive to other factors for a poly-dispersed size distribution. It may be possible to constrain a_max and perhaps f in the case of intra-mixing if observations at larger wave-numbers (> 50 cm⁻¹) are also available to constrain ξ, p and a_min. This also requires a good spectral resolution and signal-to-noise ratio in this region.

Initial comparisons of this model with Cassini-CIRS or more recent JWST-MIRI observations are promising for constraining the properties of icy regoliths, regardless of whether they are contaminated. From CIRS spectra of the outer B ring, Morishima et al. (2012) derived $p={2.9}_{-0.4}^{+0.2}$ and no strong constraint on $a_{\max }={3.2}_{-2.5}^{+97}$ cm, assuming a_min=1 µm, pure water ice, a δ-Edd diffraction correction, no asymmetry factor, and two populations of grains, cold and warm. For their study, however, they favoured a good S/N at the cost of a good spectral resolution on the roll-off region, that is, 15.5 cm⁻¹ over this 40 cm⁻¹ wide region. We are currently studying a set of CIRS spectra of the B ring with our model. Figure 9 displays as an example a few occurrences of this model, including contaminants that are intra-mixed together with two observations of the emissivity ε_F of the thickest part of the unlit B ring by CIRS-FP1. This emissivity within the field of view results from ε_F = I(w_n)/βB(w_n, T_F), where I(w_n) is the observed spectrum, β = 1, and B(w_n, T_F) the Planck function for the single fitted temperature T_F. The modelled population follows a power-law size distribution with p = 2.8, a_max= 5 cm, and a_min = 1 µm with an average asymmetry factor ξ = 0. It is observed that the contaminant and mixing type both affect the roll-off slope.

Hedman et al. (2024) observed Saturn’s rings in November 2022 with the JWST-MIRI instrument, in particular, the optically thick B ring in the spectral region 5–18 µm, where solar reflected light and thermal emission both contribute to the spectrum. Their best fit to the thermal component in the MIRI channel 3 (15.41–17.98 µm), assuming I(λ) = /βB(λ, T), yields T = 86.3 K and β = 0.663. They point out that this value of the filling factor β is far too low given an expected value for the B ring of 0.95 (Spilker et al. 2018). They also observed a decrease in brightness temperature with increasing wavelength, which they could hardly explain if the emissivity were constant with wavelength. We simulated the brightness temperature of particle regolith T_b from B(λ, T_b) = ε_h(λ)βB(λ, T_R) with ε_h(λ) estimated for two size distributions, with p = 3, a_max = 5 cm, and either a_min= 1 µm or 1 mm, assuming β = 0.95 (Fig. 10 top). The solar elevation above the ring plane is B₀ = 12.66°, and the solar distance is D_UA = 9.84775 AU. The regolith temperature T_R was estimated from the radiative balance between the absorbed solar spectrum in the visible domain (VIS) and the thermal emission in the infrared (IR),

$${\mu }_0{\displaystyle {\int }_{VIS}\left(1-r_{\text{s}}(\lambda )\right)}\frac{C(\lambda )}{D_{\text{UA}}^2}\text{d}\lambda ={\displaystyle {\int }_{\text{IR}}{\epsilon }_{\text{h}}(\lambda )}\beta B_{\lambda }\left(T_R\right)\text{d}\lambda ,$$(12)

where µ₀ = sin(B₀), C(λ) is the air-mass zero solar spectrum ASTM-E490² measured at Earth (Frölich & Lean 1998), and r_s(λ) is the hemispherical reflectance (Eq. (4)). For this specific calculation, the range of optical constants for water ice was extended to the UV NIR domain with Warren & Brandt (2008) data. The brightness temperature is then about 82–85 K for λ > 11 µm and decreases with wavelength when sub-millimetre grains are present in the regolith (Fig. 10 bottom). If this is not the case, the hemispherical emissivity is fairly constant in this range, and so is T_b. The slope in this region is both a function of the regolith thermal history at the epoch of observation, the local hour angle on the rings, and of the fraction of contaminants that modulate the radiative balance. A peak can be observed near λ=9.5µm. It results from the drop in emissivity observed in Fig. 10, top, which corresponds to a local peak in the hemispherical reflectance r_s, and therefore, to a peak in the solar contribution, which is still detectable above the thermal emission (with a nearly constant brightness temperature). The shape of this peak is highly dependent on the size distribution, the contaminant, its fraction, and the diffraction correction applied (not shown). In the case of pure ice, the thermal contribution is much smaller with an equilibrium temperature ~66 K because solar light is only weakly absorbed in the visible domain. The peak is absent in the case of intra-mixed AmC at the 3% level. It has been detected by Hedman et al. (2024). With this simple model, the spectral shape and range of the brightness temperatures observed by Hedman et al. (2024) (reported in the bottom panel of their figure 13) can then be reproduced reasonably well with a filling factor β = 0.95 with the addition of contaminants and current estimates of the size distribution of regolith grains (Morishima et al. 2012). As expected from the global sensitivity analysis, this spectral region may bring complementary information on the regolith properties and composition. This would deserve a quantitative adjustment.

Bockelée-Morvan et al. (2024) targeted the leading and trailing hemispheres of Ganymede with the JWST-MIRI instrument in this same spectral region. They observed a remarkable quasi-linear decrease in the brightness temperature with increasing wavelength by more than 10 K in between 7 and 11 µm (their Fig. 14), for which they proposed three possible origins that may be combined: either a spectrally constant emissivity <1, a decreasing spectral emissivity versus wavelength, or a variety of surface temperature mixing in the PSF. The brightness temperature estimated with our model, including both solar and thermal contributions and assuming the regolith in thermal equilibrium at a solar distance of 4.96 AU (i.e. replacing µ₀ by 1/4 in Eq. (12)), may exhibit a steep slope with increasing wavelength for sub-millimetre grains (Fig. 11). It is shown here for the case of intra-mixed AmC with ƒ = 3%. The size of regolith grains and the nature of the diffraction correction significantly affect the slope of T_b and its linearity in this spectral region. As in the previous case of Saturn’s rings, the slope of T_b depends on the temperature, that is, on the energy balance and thermal history, and on the contaminant mixing in the water ice. Therefore, a convolution of this model with a terrain model and a thermal model remains essential to take into account the surface temperature heterogeneity over the observed hemisphere and to interpret these observations. The influence of the two contaminants (AmC or Th) at the level of a few percent does not yield significant variations in the slope, but their fraction is more influential for λ < 6 µm or λ > 10 µm.

Fig. 9 Modelled hemispherical emissivity εh of a regolith of grains following a power-law size distribution with p = 2.8, amax = 5 cm, and amin = 1 µm. Water ice is crystalline (Ih90) and contaminants are AmC or Th, intra-mixed with fraction f. Two estimates of the emissivity εF of the unlit B ring with the CIRS instrument on April 5, 2017, are displayed (blue dots, see text for details). Both esa (solid line) and δ-Edd (dash-dotted line) diffraction corrections are shown, assuming ξ = 0.

Fig. 10 Model of the mid-infrared brightness of the optically thick B ring of Saturn. (top) Hemispherical emissivity of a regolith of grains following a power-law size distribution with p = 2.8, amax= 5 cm, and amin=1 µm (solid line) or amin=1 mm (dashed line). Water ice is crystalline (Ih90), and contaminants (AmC or Th) are intra-mixed with a fraction ƒ. The esa diffraction correction is used here and ξ = 0. (bottom) Corresponding modelled brightness temperatures Tb(K) of these regoliths on particles. The solar contribution is plotted for both contaminants in the case of amin=1 mm (dotted line).

Fig. 11 Modelled brightness temperature Tb (K) of a regolith of grains with a given size of 10 µm (solid line), 50 µm (dashed line), or 500 µm (dotted line) at equilibrium temperature at a Jupiter distance of 4.96 AU in the wavelength range of the JWST-MIRI instrument. The grains are composed of pure crystalline ice (Ih90) with intra-mixed AmC with a fraction ƒ = 3%. The esa (black) or δ-Edd (red) diffraction corrections are used here, and ξ = 0. The solar contribution (orange) is plotted for all sizes in the case of an esa diffraction correction.

4 Conclusions

A hybrid Mie-Hapke model for the hemispherical emissivity of an icy regolith contaminated with either amorphous carbon or tholins along with various types of mixing was proposed here. It can be coupled with two diffraction-correction methods and various types of size distributions, and it includes the average asymmetry factor of the regolith. A global sensitivity analysis quantified the importance of these factors on the variance of its hemispherical emissivity versus wave number. We plan to analyse a wider dataset of the CIRS instrument with this model, for either rings or moons, to complement existing studies, including the consideration of lightly contaminated water ice or seasonal phenomena. This model provides a self-consistent tool for interpreting multi-modal observations of the thermal emission from icy surfaces. It offers interesting insights into recent mid-infrared observations of Saturn’s rings and the Galilean moon Ganymede by the JWST-MIRI instrument.

Acknowledgements

This work was supported by French Programme National de Planétologie (PNP) and Centre National d’Etude Spatiales (CNES) under projects PISTE and ADSINTHE respectively.

References

All Tables

All Figures

	Fig. 1 Complex refractive index of water ice, either crystalline (Ih) or amorphous (Am), as proposed by various authors at different temperatures. (top) Real part nr and (bottom) imaginary part nk. Legends correspond to authors T (K) ice phase with authors C05 (Curtis et al. 2005), M83 (Mishima et al. 1983), H93 (Hudgins et al. 1993), M09 (Mastrapa et al. 2009), and W08 (Warren & Brandt 2008).
In the text

	Fig. 2 Complex refractive index of Ih water ice and contaminants. (top) Real part nr and (bottom) imaginary part nk of the Ih. Ih T0 point to the refractive index of Ih water ice at temperature T0 (see text). The optical constants (nr, nk) for an intra-particle mixture of Ih 90 ice with AmC (Preibisch et al. 1993) or Th (Khare et al. 1984) at 3% and 1% levels respectively are also plotted.
In the text

Fig. 3 Single-scattering albedo and hemispherical emissivity as a function of grain size and wave number or wavelength. (top) Single-scattering albedo ω0(wn, a) of a grain with size a (bottom) hemispherical emissivity εh(wn, a) of regolith covered with these grains made of pure crystalline water ice Ih90. Calculations with the Mie theory (dotted line) or after diffraction corrections such as delta-Eddington (dash-dotted line) or equivalent slab approximation (solid line) are displayed. (Left) 10–700 cm−1 range and (right) 5–14.3 µm range (or 700–2000 cm−1).

In the text

Fig. 4 Hemispherical emissivity ɛh(wn, a) of regolith grains with various sizes a. (top) 10 µm, (middle) 100 µm, (bottom) 1 mm, and various types of mixing: no contaminants (Ih90), with 5% AmC (black) or 1% Th (red) either mixed at intra-particle scale (dotted line), or with aerial (dashed line) or volumic (dash-dotted line) mixing. esa diffraction correction is assumed. Left: 10–700 cm−1 range and right: 5–14.3 µm range (700–2000 cm−1).

In the text

	Fig. 5 Hemispherical emissivity εh(wn, a) as a function of p=2 (dashed line) or p=3 (solid line) and (top) amin or (bottom) amax. The grain composition is water ice Ih90, intra-mixed with 5% AmC. esa diffraction correction is assumed. Left: 10–700 cm−1 range and right: 5–14.3 µm range (700–2000 cm−1).
In the text

	Fig. 6 Hemispherical emissivity ${\epsilon }_{\text{h}}^{\ast }\left(w_{\text{n}}\right)$ as a function of average asymmetry factor ξ for a size distribution of gains with amin = 1 µm , amax = 5 mm, and p=2 (dashed line) or p=3 (solid line). The grain composition is water ice Ih90 with intra-mixing of 5% AmC. esa diffraction correction is used.
In the text

	Fig. 7 Total effect sensitivities S Ti in the case of an intra-mixture with AmC contaminant and Mie (dotted line), δ-Edd (dash-dotted line) or esa (solid line) diffraction correction methods and (top) mono-dispersed (bottom) power-law size distributions. The dashed line marks the 3% level.
In the text

	Fig. 8 Total effect sensitivities STi in the case of contamination with an AmC mixture, either intra (dotted line), aerial (dashed line), or volumic (solid line) mixing for (top) mono-dispersed (bottom) power-law size distributions. The esa diffraction correction is used here.
In the text

Fig. 9 Modelled hemispherical emissivity εh of a regolith of grains following a power-law size distribution with p = 2.8, amax = 5 cm, and amin = 1 µm. Water ice is crystalline (Ih90) and contaminants are AmC or Th, intra-mixed with fraction f. Two estimates of the emissivity εF of the unlit B ring with the CIRS instrument on April 5, 2017, are displayed (blue dots, see text for details). Both esa (solid line) and δ-Edd (dash-dotted line) diffraction corrections are shown, assuming ξ = 0.

In the text

Fig. 10 Model of the mid-infrared brightness of the optically thick B ring of Saturn. (top) Hemispherical emissivity of a regolith of grains following a power-law size distribution with p = 2.8, amax= 5 cm, and amin=1 µm (solid line) or amin=1 mm (dashed line). Water ice is crystalline (Ih90), and contaminants (AmC or Th) are intra-mixed with a fraction ƒ. The esa diffraction correction is used here and ξ = 0. (bottom) Corresponding modelled brightness temperatures Tb(K) of these regoliths on particles. The solar contribution is plotted for both contaminants in the case of amin=1 mm (dotted line).

In the text

Fig. 11 Modelled brightness temperature Tb (K) of a regolith of grains with a given size of 10 µm (solid line), 50 µm (dashed line), or 500 µm (dotted line) at equilibrium temperature at a Jupiter distance of 4.96 AU in the wavelength range of the JWST-MIRI instrument. The grains are composed of pure crystalline ice (Ih90) with intra-mixed AmC with a fraction ƒ = 3%. The esa (black) or δ-Edd (red) diffraction corrections are used here, and ξ = 0. The solar contribution (orange) is plotted for all sizes in the case of an esa diffraction correction.

In the text