© C. Ferrari et al. 2021

1 Introduction

Beyond large-scale meteorite impacts or geological activity, there are many processes at work that effectively alter the surfaces of Saturn’s atmosphere-less icy moons. Irradiation by solar ultraviolet photons or bombardment by magneto-spheric plasma particles are among the processes known to alter the regolith properties through chemical alteration, amorphisation, sputtering, or molecular desorption (Johnson et al. 2013). The regolith may also be polluted by exogenous elements originating from Saturn’s E ring or Phoebe ring (Howett et al. 2018). The structural alteration of the surface over time and the intricacy of these processes are difficult to constrain (Baragiola 2003). However, they may be observed through surface features or radiometric asymmetries.

Thanks to numerous multi-wavelength observations, the Cassini mission (2004–2017) has provided a much more detailed picture of the surface of these moons. Asymmetric patterns between their leading and trailing hemispheres have been confirmed or detected, either in terms of colour and Bond albedo at visible and near-infrared (VNIR) wavelengths or in terms of thermal inertia at thermal infrared wavelengths (Hendrix et al. 2012; Howett et al. 2010). The brightness increase observed in the VNIR suggests that the E ring icy grains could impact these surfaces by sandblasting or by deposition (Verbiscer et al. 2007; Kempf et al. 2018; Howett et al. 2018 for a review). The Cassini VIMS (Visible and Infrared Mapping spectrometer) data provide evidence that the leading hemisphere of Dione, Tethys, and Rhea exhibit deeper absorption bands of water ice, suggesting the presence of fresh E ring particles (Scipioni et al. 2013, 2014; Stephan et al. 2010, 2012, 2016). Hendrix et al. (2012) suggest that the regolith grains are smaller on the trailing hemisphere of Mimas than on its leading hemisphere, which is consistent with a bombardment with E ring particles.

The alteration of icy surfaces by high energy electrons produces lens-shaped patterns at the equator of Mimas, Tethys, and Dione (Paranicas et al. 2001, 2014). In Saturn’s system, keV electrons are expected to bombard their trailing hemisphere while MeV electrons should collide with their leading hemisphere. The extent of the impacted areas varies according to the moon (Paranicas et al. 2014; Nordheim et al. 2017; Howett et al. 2018). Lens-shaped patterns have been discovered in the leading hemispheres of icy moons that look darker in the visible and brighter in the IR/UV colour ratio maps (Schenk et al. 2011). They are correlated with thermal anomalies observed by CIRS and interpreted as higher thermal inertia in the leading hemispheres (Howett et al. 2011, 2012). No such anomalies were discovered on Rhea (Howett et al. 2014). It is indeed expected that MeV electrons alter the surface down to a depth of a few centimetres, where the diurnal thermal waves penetrate and the thermal inertia anomalies are detected (Paranicas et al. 2014; Schaible et al. 2017; Nordheim et al. 2017). However, the predicted lens-shaped features in the trailing hemisphere, induced by below-MeV electrons at shallower depths, have not been observed and may be hidden by some competing process such as E ring deposition or cold plasma bombardment (Nordheim et al. 2017; Howett et al. 2018). The darkening and reddening observed on the trailing hemispheres of Tethys, Dione, and Rhea (Schenk et al. 2011) may result from radiolysis induced by colliding protons or ions within this cold plasma (Johnson et al. 2013; Schenk et al. 2018; Howett et al. 2018). The VIMS near-infrared spectra did not reveal any change in size or composition correlated with the lenses (Hendrix et al. 2012).

If the observed lens-shaped patterns can be explained by an alteration of the surface by MeV electrons, the mechanism of grain sintering and the reason for thermal inertia increases in the icy regolith remain poorly understood (Ferrari & Lucas 2016; Schaible et al. 2017). To explore this direction further, we propose a thermophysical model that relates thermal inertia to regolith properties, namely: grain size and composition, layer structure, porosity, and the quality of grain contacts.

The thermal emission of Saturn’s icy moons has been thoroughly observed by the CIRS instrument during Cassini’s flybys. Analysis of this voluminous dataset by Howett et al. (2010, 2011, 2012, 2014, 2016) yielded major advances in the detection of regional thermal patterns. It brought new constraints on the Bond albedo and the thermal inertia, and helped build scenarios on their origin. Thermal inertia were found to be very weak, around 10–20 J m⁻² K⁻¹ s^−1∕2 (hereafter noted as MKS). These constraints were derived assuming smooth surfaces, unit emissivity (ε_obs = 1), and vertical homogeneity of the icy regolith. The hypothesis made on emissivity arises from the fact that most of the analyses were performed on daytime data of the CIRS mid-infrared focal plane FP3 to maximise the spatial resolution, but at the cost of a degeneracy between temperature and emissivity, FP1 measurements being used to determine night-time temperatures. High porosity of the regolith is generally cited to explain such a low thermal inertia. Ferrari & Lucas (2016) related such low values to the presence either of amorphous ice grains, whose thermal conductivity is much lower than that of the crystalline form, or by the presence of rough grains with irregular faces that are likely to reduce the solid conduction through grain contact. They showed that such low values can then be obtained with normal porosity of particle beds. Thermal emission largely depends on regolith albedo and emissivity that balance the absorbed and emitted radiation on the surface and that depend on both grain size and composition. Therefore, the unit emissivity assumption may lead to an underestimation of the measured surface temperature and impact the derivation of thermal inertia.

The topography, mesoscale roughness, and porosity of the upper regolith layers also impact the energy deposition on the surface and the local thermal history, mainly due to the various slopes and shadowing effects (Spencer 1990; Lagerros 1996; Ferrari 2018). Davidsson et al. (2015) investigated the effect of topography on thermal emission with artificial topographies. They showed that it strongly depends on the roughness type, even when the degrees of roughness are identical. They outlined the bias introduced by these topographic models on the determination of thermal inertia and surface roughness. So it is also expected that the temperature observed in the CIRS field-of-view (FOV) depends on the local topography.

In this paper, we put together the thermal model developed by Ferrari & Lucas (2016), which is controlled by regolith grains and layer properties, and the known topography of Saturn’s icy moons at various spatial resolutions ranging from 1 km to 8 km. In Sect. 2, we couple heat transfer in the near-surface with topography to predict the thermal emission as a function of regolith properties, topography, the area covered by the footprints of the FOV on the surface, and the epoch. The effect of these parameters on the directional and temporal variations of Dione thermal emission is investigated in Sect. 3. Finally, the ability of the model to mimic the measurements along a latitudinal scan on the surface of this moon around noon local time is evaluated.

2 Thermal emission of icy regoliths

The thermal emission of Saturn’s icy moons was captured by the CIRS focal plane FP1 in the wavelength range from 17 to 1000 μm (Flasar et al. 2004). The footprints corresponding to the 3.9 mrad-wide FOV of the instrument vary from a few tens to a few thousand kilometres (Flasar et al. 2004). They comprise different local slopes and orientations relative to the sun and the instrument, so that various thermal histories and temperatures are captured at the time of observation. Thermal emission is a mixture of black/gray body thermal emissions. The temperature, T_F, and the apparent emissivity, ε_F, measured in the FOV may vary as well with the viewing geometry. We describe the articulation between surface models and a thermal model to account for the effect of topography and regolith properties on the thermal emission. We take advantage of an observation of Dione to study the impact of the complexity of the shape model on the measured T_F and ε_F.

2.1 Topography of icy moons

Thanks to the lifespan of the Cassini mission and the multi-angular viewing geometries during its flybys, shape models have been derived for some icy moons (Gaskell et al. 2008). Specific formats (containers) known as DSK (Digital Shape Kernel) have been generated to use these models within the NASA’s Navigation and Ancillary Information Facility (NAIF) SPICE system environmentof observation ephemeris and geometry (Acton et al. 2018). These shape models are made of triangular plates and vertices positioned in the inertial body reference frame. The DSK environment provides routines to compute the position, visibility and illumination angles of any plate at a given epoch. The Python wrapper (SpicePy, Annex et al. 2020) of the NAIF/SPICE software package N066 was used to link the shape file with the thermal model and the CIRS pointing vectors. The surfaces have be tessellated at various ground sample distances (GSD) and vectors are stored in an implicitly connected quadrilateral format (ICQ, Gaskell et al. 2008). The GSD can be approximated by with R_M the radius of the moon and q a parameterthat reflects the grid sampling that is usually taken as a power of 2. It corresponds to a latitudinal resolution of deg per plate. According to Gaskell’s nomenclature, the shape model at a resolution of q = 64 has to be considered as a global topography model (GTM) of good spatial resolution. Shape models can be found up to q = 512. This topographic model is referred as the q-topo case (Fig. 1).

The surface of Dione has been also approximated by a tessellated ellipsoid of variable resolution, parametrised by a latitudinal sampling of 180∕n deg per plate, that is . The ellipsoid is generated and tessellated using MATLAB, with a longitudinal sampling twice as large, since the number of samples per axis should be constant. It is then encapsulated in a DSK (Digital Shape Kernel) format to be exploited in the same way as the topographic model. The tessellated ellipsoid model is referred as the n-ell case. This surface model allows for studying the effect of the dispersion of local times of plates in the footprint on the thermal emission, independently of the effect of the topography, which itself introduces an additional dispersion due to the different slopes, orientations, shadows, and visibility conditions (taken into account by the q-topo model).

The viewing and illumination angles are calculated at the barycentre of each plate within the footprints of the observation. The relativedistance from the plate to the center of the footprint is also stored to estimate its relative contribution to the total heat flux in the footprint, taking into account the FOV spatial sensitivity (Sect. 2.2.1).

Finally, we have the ‘flat’ case, where the surface is approximated by a single plate whose normal is that of the underlying ellipsoid at the FOV intercept point. This is the simplest surface model, where neither the effect of dispersion in local times in the footprint nor that of the topography are taken into account. The illumination and viewing geometries are calculated at this point at the time of acquisition and the model calculates the thermal history of this virtual surface element.

Fig. 1 Representation of the Gaskell Global Topography Model (GTM) of Dione. Southern latitudes of the leading hemisphere are visible about longitude 90°W with the Evander basin on the lower left near the limb and the Sabinus crater below the center.

2.2 Thermal emission of the surface elements

The model for the thermal emission of planetary regolith inherits the formalism developed by Ferrari & Lucas (2016). It links the effective conductivity and the thermal inertia of regolith with its compositional and structural properties. Thermal inertia can be highly time and space dependent as the effective thermal conductivity may have a significant temperature and grain size dependent radiative component, even in icy cold regoliths (Ferrari & Lucas 2016). Heat transfer by conduction through contacts may be indeed very limited due to grain roughness and supplanted by raditive heat transfer through the pores. Here, this formalism iscoupled with the heat diffusion equation so as to account for the thermal history of any surface element of the GTM.

2.2.1 Heat diffusion and boundary conditions

Temperature distribution in the regolith, T(z, t), is a function of depth, z, and time, t. It is determined by solving the heat diffusion equation: (1)

with ρ(z) as the effective volume density, k_E as the effective thermal conductivity, and C_P(T) as the specific heat capacity. The value of ρ(z) depends on the regolith porosity, p, which may vary with depth but which we keep constant here (Table 1). The effective thermal conductivity can be either constant (k_E = k_E0) or may depend on temperature, porosity, and grain size (k_E = k_E(T, z)). Radiative heat transfer may indeed occur in parallel with conduction heat transfer, so that k_E (T, z) = k_C + k_R. Ferrari & Lucas (2016) studied the variations of the effective thermal conductivity with grain size, ice phase, the type of grain contact (tight or loose), and the model used for the radiative conductivity k_R. The latter is proportional to the emissivity, the grain size, R, the temperature, T³, and the porosity, p (Table 1, Breitbach & Barthels 1980; Gundlach & Blum 2012). Ferrari & Lucas (2016) also showed that: (1) Thermal inertia values as low as 10–20 MKS at Saturn temperatures (80 ± 30 K) could be obtained with regolith of normal porosity (p = 0.3–0.8) but with loose grain contact or amorphous ice phase; (2) The thermal inertia of the regolith is minimum for grain sizes in the range of 100 μm–1 mm. For lower sizes, it increases as contacts between particles increase and radiative conductivity decreases. For larger grain size, the radiative conductivitydominates as it scales with the grain size; (3) The thermal inertia of amorphous ice at 80 K is minimum for grain size in the range 10–100 μm, while it can reach up to 100 MKS for grains of several centimetres in size.

In such a model, thermal inertia is therefore expected to vary spatially depending on local time, temperature, latitude, and regolith properties (porosity, grain size or ice phase). Any plate radiates into space from its skin-layer (z = 0) with the boundary condition: (2)

where μ₀ stands for the cosine of the illumination angle relative to the local normal of the plate (q-topo or n-ell shape models) or to the normal of the ellipsoid at the center of the FOV intercept (flat surface model). In addition, T_S is the surface temperature, A its bolometric Bond albedo, its average hemispherical emissivity, S the solar constant, and D_AU the heliocentric distance measured in AU at the epoch of acquisition. The bottom layer at depth may receive an internal up-welling thermal flux ϕ₀, which is null here. The diffusion equation is solved with a Crank–Nicolson scheme. Vertical sampling is fixed to ten samples per skin depth δ_p, which is estimated for the average temperature at the epoch of acquisition for a diurnal cycle. Time sampling is adapted to ensure convergence on temperatures within 0.1–0.2 K in a reasonable calculation time, that is, up to 0.3 sec per plate over a full diurnal cycle of Dione and for very low thermal inertia. This is below the error bar of ~0.5 K on observed temperatures.

In cases where the regolith does not conduct heat, either because the thermal conductivity is very low (k_E ~ 0) or the porosity is very large (p ~ 1), the thermal inertia is near zero and the surface temperature is calculated by: (3)

At this point, we can check that the surface temperature, T_S, which is a solution of Eq. (2), converges on T₀ for very low thermal inertia.

2.2.2 Bond albedo and emissivity

The surface of Saturn’s icy moons is covered with water ice, most likely crystalline ice (e.g. Cruikshank et al. 2005). Both crystalline and amorphous phases may coexist at various depths depending on their thermal and space weathering history (Ferrari 2018). Their VNIR reflectance spectra typically exhibit a reddening below 1 μm and a quenching of strong water ice absorption above 3.6 μm. This indicates the presence of tholins as no other darkening components can explain such effects (Cruikshank et al. 2005). Pitman et al. (2010) constrained the bolometric Bond albedo A of the leading and trailing hemispheres of icy moonsusing VIMS data: A ranges from 0.37 ± 0.08 to 0.93 ± 0.11 for Dione and Enceladus trailing hemispheres, respectively. A decrease in A may be due either to an increased grain size or to an increased fraction of water ice contaminants (Emery et al. 2005). The fraction of water ice contaminants varies from a few to up to tens of percent and the superficial grain size ranges from 1 μm to a few hundred μm, depending on the photometric model (Cruikshank et al. 2005; Scipioni et al. 2013).

Modelling the VNIR data is beyond the scope of this paper, however, the bolometric Bond albedo, A, and the average hemispherical emissivity, , which depend on grain size, composition, and temperature, are primordial factors controlling the energy balance (Eq. (2)). Thus, we developed a coherent approach ensuring consistency between the solar reflectance and the infrared emittance. For that purpose, we introduce Hapke’s formalism (Hapke 2012).

Let consider an intimate mixture of water ice and tholins, the optical constants of which are available over a large wavelength domain from the UV to the millimetre waves (Fig. 2). Hapke (2012) expressed the plane albedo, A_h, and the hemispherical emissivity ϵ_h of a regolith made of grains with size, R, as a function of their single scattering albedo, ω₀, their cosine asymmetry factor, ξ, and the regolith porosity. Emissivity may be also directional, ϵ_d(μ), depending on the emission angle (Table 1). For a quasi-continuous medium made of arbitrary single-scatterers, we can use Hapke’s Isotropic Multiple Scattering Approximation (IMSA) model to write the hemispherical reflectance r_h and the hemispherical emissivity ϵ_h: (4)

and (5)

with κ_P the porosity factor set to 1, and r₀ = (1 − γ_λ)∕(1 + γ_λ). The directional hemispherical reflectance r_h(μ₀, λ), also called plane albedo, is the total fraction of incident collimated radiation per unit area scattered into the upward hemisphere. In case of anisotropic scattering within the slabs, the single scattering albedo ω₀ can be replaced by : (6)

and (7)

This approximation works quite well with the hemispherical reflectance, however, the performance is reduced when considering the directional reflectance (Hapke 2012). The spherical albedo A_S (λ) is the total fraction of incident irradiance scattered by a body in all directions. In the case of isotropic scatterers, it is written as: (8)

The bolometric Bond albedo A is the average of the spherical albedo weighted by the solar spectral irradiance B_λ (T_⊙): (9)

The value of B_λ(T_⊙) is estimatedby the Standard Zero Air Mass Solar Spectral Irradiance ASTM E490 over the 0.2–1000 μm range. The corresponding solar constant S = 1366 W m⁻². The radiance emitted in the direction μ from an isothermal layer at temperature T is I(λ, μ) = ϵ_d(μ, λ)B_λ(T), where the directional emissivity ϵ_d(μ, λ) = γ_λH(μ) can be approximated by: (10)

The average hemispherical emissivity over the thermal infrared domain at Saturn distance for a surface element covered with regolith is: (11)

with B_λ(T_S) as the blackbody radiation of the plate at surface temperature, T_S. The single scattering albedo of a grain with size R is calculated by the Mie theory using the complex refractive index of crystalline and amorphous water ice. When appropriate, the geometrical optics approximation is used to speed up the calculations. Bertie & Whalley (1967) data were complemented by Hudgins optical constants at longer wavelengths for crystalline water ice around 100 K (Hudgins et al. 1993). In case of pure amorphous water ice, we assumed that k ~ λ^−0.5 in the far infrared region (Fig. 2). The optical constants of tholins are taken from Khare et al. (Khare et al. 1984). We verified that an increase in the tholins fraction reddens and darkens the VNIR spectrum as the imaginary part of the complex refractive index (the extinction coefficient) of tholins decreases and is much larger than that of water ice between 0.2 and 0.7 μm (Fig. 2).

As expected, the bolometric Bond albedo A decreases with the grain size, R, and the tholins fraction up to 20% (Fig. 3a). It is close to unity for pure ice micrometre grains as suspected on the Enceladus trailing hemisphere (Pitman et al. 2010). For a tholins fraction of a few percent, A ~0.6, which is a common value on icy moons, for grain sizes between 10 μm and a few hundred μm. The minimum observed albedo A ~ 0.4 is obtained for a fraction of 20%.

The model shows that the average hemispherical emissivity does not change whatever the tholin fraction (Fig. 3b). It highly depends on grain size: –0.6 for sub-micro-metric grains; it is lowest at 0.3 for micro-metric grains and grains ≤ 0.1 μm; it steeply increases for grain sizes between 10 μm and 40 μm; finally, –0.97 for larger grains. Given a grain size, it exhibits some dispersion when the surface temperature varies from 50 to 110 K (Fig. 3b). Similar curves are obtained with amorphous water ice. While depends on T_S, we approximate (80 K) in the boundary condition (Eq. (2)) to simplify the resolution of the transient thermal cycle. It is therefore assumed that depends on grain size but not on temperature.

As a consequence, the regolith temperature is also expected to vary with grain size, since it is scaling with : grains ≤ 0.2μm are poor emitters and absorbers; grains ~1 μm are poor emitters and good absorbers; and larger ones are good absorbers and emitters. Figure 3c displays the characteristic temperature T_slow = T₀∕2^1∕4 of a slowly rotating moon with low thermal inertia and covered with icy grains polluted by a variable fraction of tholins (0%, 5%, 20%). For pure ice, T_slow increases with R ≥ 0.3 μm. As the albedo decreases with increased tholins fraction, the temperature of μm-to-10-μm-sized grains becomes prominent. The temperature of grains ≥60 μm is almost constant for tholins fractions greater than a few percent.

This indicates that the bolometric Bond albedo, A, is mainly controlled by the grain size and pollutant fraction in ice, while the average hemispherical emissivity, , in the infrared domain is mainly controlled by the grain size, the icy composition and marginally the temperature. We may keep the tholins fraction as a model parameter that fixes A on the basis of the ice phase and the grain size. However, in the following, A will be set up a priori, taken from Pitman et al. (2010). Observations in the visible domain are indeed the most appropriate to constrain the bolometric Bond albedo. It is also useful to limit the number of parameters to perform the inversion of the model within acceptable time frames.

Fig. 2 Optical constants of hexagonal water ice (Ih) and tholins (from Khare et al. 1984). Both real (nr) and imaginary (nk) components of the refractive indices are plotted.

2.2.3 Diurnal thermal cycles

Figure 4 provides examples of the diurnal thermal cycle of a surface element of the GTM (q = 64 and GSD = 8 km) located on the equator of Dione (west longitude 90°). They are part of a simulation aimed at reproducing a CIRS observation which consists in a latitudinal scan of its leading hemisphere acquired on September 30, 2015 and named CIRS_222DI_DIONE001_PIE. Howett et al. (2010, 2014) showed that the thermal inertia Γ of Dione was about 10 MKS. A constant thermal conductivity k_E0= 10⁻³ W m⁻¹ K⁻¹ yields similar thermal inertia with the values of C_P(T), ρ₀, and p provided in Table 1.

For a constant effective thermal conductivity k_E0 (Fig. 4a), the regolith average temperature over the diurnal cycle increases conversely to the average hemispherical emissivity of grains: micro-metric grains (0.2–2 μm) with lower emissivity are hotter (Fig. 3b). Regolith covered with grains larger than a few hundred μm exhibit lower average temperature because of their higher emissivity. The thermal inertia decreases with higher porosity, so that the amplitude of the diurnal thermal gradient expectedly increases.

When the thermal conductivity is both radiative and conductive (Fig. 4b), the thermal cycle depends more on grain size, temperature, and bulk conductivity of the ice phase. It also varies with the type of contacts between grains. Ferrari & Lucas (2016) showed that thermal inertia decreases conversely to grain size in the range of μm-to-mm for a regolith porosity p ≥ 0.3. It increases again for cm-sized grains. In case of amorphous ice with loose contacts, the thermal inertia for grains size ≤ 1 μm reaches 60 MKS as the contact conductivity scales with 1∕R and decreases below 22 MKS for larger ones. The thermal gradient significantly increases. For cm-sized grains the thermal inertia increases again and the gradient is reduced. The thermal gradient and average temperature are all the more dependent on the grain size and bulk conductivity of the ice phase. Because their bulk conductivity indeed differs by almost a factor of ten, the thermal inertia of amorphous ice is lower and the diurnal thermal gradient is larger (Fig. 4b).

Fig. 3 Modelled regolith properties versus grain size R (a) Bolometric Bond albedo A (Eq. (9)) versus grain size, R, for various tholins fractions (0%, 5%, 20%). The range of bolometric Bond albedos for Saturn’s icy moons is marked as horizontal dotted lines (Pitman et al. 2010). (b) Average hemispherical emissivity versus grain size, R, for various tholins fractions (0%, 5%, 20%) and surface temperature TS = 80 K (Eq. (11)). The error bars figure the dispersion of when TS varies from 50 K to 110 K for a tholins fraction of 5% in crystalline water ice grains. The curves are similar in cases of amorphous water ice (Eq. (11)). (c) Typical surface temperature Tslow of a plate covered with regolith of grain size, R, for various tholins fractions (0%, 5%, 20%) in the case of a slowly rotating moon.

2.3 Temperature and apparent emissivity in CIRS FP1 FOV

The temperature estimated from the measured spectrum in the FP1 footprint depends on the thermal history of all the surface elements the footprint contains. Considering N visible elements, the modelled spectrum I_F(λ) can be calculated as the weighted sum of the elementary spectra as (12)

with μ_j the cosine of the local emission angle of plate j, S_j its surface, and I_j(λ, μ_j) its modelledemission spectrum. This latter can be either (13)

for an isotropic hemispherical emissivity or (14)

for a directional emissivity ε_d,j. The emission spectrum I_j(λ, μ_j) of each surface element is weighted by the spatial sensitivity w_j of the FP1 detector that is a function of its distance to the center of the FOV (Sect. 2.1). The sensitivity to the distance is assumed to be a 2D Gaussian with a FWHM of 2.42 mrad and an aperture diameter of 3.9 mrad (Flasar et al. 2004).

The apparent temperature, T_F, and emissivity, ε_F, in the FOV are derived from the best fit of I_F (λ) = ε_FB_λ(T_F) over the wavenumber (w_n) range from 20 to 350 cm⁻¹ (wavelength range from 29 to 500 μm) around the observed emission peak. This is the same range wherein the CIRS FP1 emission spectra I_obs are fitted in order to optimise the signal-to-noise ratio and to minimise the errors on T_obs and ε_obs. The noise equivalent spectral radiance (NESR) may be prominent above this wavenumber range depending on the spectral resolution used and the measured temperature and emissivity (Flasar et al. 2004). Figure 5 shows simulated and observed emission spectra within two footprints selected among the 350 obtained during the observation (Sect. 3, Fig. 6). The temperature and emissivity are lower in southern footprint #5 than in footprint #3 (Sect. 3, Fig. 7) and the data are noisy above 400 cm⁻¹. For the chosen model parameters – which do not result from any inversion procedure – it can be seen that the inclusion of spatial sensitivity in the FP1 FOV has variable impact on the modelled emission spectrum I_F and is most detectable at southern latitudes (footprint #5). Increasing the thermal conductivity K_E0 from 0 to 0.001 W m⁻¹ K⁻¹ reduces the apparent temperature (the peak shifts to lower wavenumbers) in these footprints which have been acquired close to 12:30 local time.

Fig. 4 Thermal cycles of a surface element located at the center of the leading hemisphere of Dione. The bolometric Bond albedo A is assumed 0.63 (Pitman et al. 2010). (a) The thermal conductivity is constant kE0 = 0.001 W m−1 K−1 and the grain size, R, varies between 0.2 μm and 2 cm, assuming a porosity of p = 0.3 (solid line) or 0.8 (dash-dot line). The ice phase is amorphous. The resulting thermal inertia range between 10 and 20 MKS. (b) Both radiative and solid conductions through loose contacts are assumed here. The grain size, R, varies between 0.2 μm and 2 cm and the ice phase is either crystalline (Ih, solid line) or amorphous (Am, dash line) ice phase. The porosity is p = 0.5.

3 Results

The impact of topography and regolith properties on the directional variations of T_F and ε_F is discussed for the observation CIRS_222DI_DIONE001_PIE. Its design follows a south-north latitudinal scan, from about − 50°S to + 16°N, at an average west longitude of 123°W (Fig. 6). The spatial extent of FP1 footprints is about 270 km at the equator. N = 3413 surface elements of the GTM are flown over during that sequence, with an average of 600 elements per footprint. Local time is almost constant during the scan, between 12:30 and 12:50. For the GTM (q = 64), the GSD is about 8 km (given R_M = 561 km), which fits the latitudinal sampling of the n-ell ellipsoidal shape model with n = 222.

The incidence angle at the FOV intercept with the ellipsoidal surface of the moon varies from 80° to 25° along the scan while the emission angle varies from 55° to a few degrees. Theobserved temperature T_obs decreases from 95 K to 78 K with the incidence angle (Figs. 7a and 8b). The apparent emissivity ε_obs highly varies along this scan, from 0.5 at southern latitudes (high emission angles) to 0.92 close to the normal viewing angle (Fig. 7b). At ε_obs ≤ 0.75, the measured temperature T_obs is constant around 80 K. Releasing the constraint of unit emissivity (ε_obs = 1) highlights the directional anisotropy of the apparent emissivity as a function of the emission angle Em.

3.1 Effect of topography

Simulations of apparent temperature and emissivity within the footprint are performed assuming a regolith with a bolometric Bond albedo A = 0.63 and a grain size, R = 2 mm (Fig. 7). Their average hemispherical emissivity is indeed comparable to the observed value ε_obs near nadir (Figs. 3b and 7b). The thermal conductivity is either null (by default) or constant with k_E0 = 0.001 W m⁻¹ K⁻¹ (when specified), so that the thermal inertia is either Γ = 0 or Γ = 15 MKS for a porosity p = 0.5.

Considering a ‘flat’ surface model, the apparent emissivity is found to be constant along the scan: . The apparent temperature is consistent with T_obs, although slightly too cold for southern latitudes (footprints #4 and #5). The value of ε_F varies from 0.93 to 0.86 (~7.5%) with Em when the hemispherical directional emissivity e_d(μ) – Eq. (14) – is included in the model, which is much lower compared to the observed variation.

In the case of a tessellated ellipsoid (n-ell), ε_F varies from 0.92 to 0.71 (~23%) with Em. At the spatial resolution of ~8 km (n = 222), this model indeed takes into account the dispersion in local time due to the dispersion in longitude of the different surface elements in the field of view. If it is not yet sufficient to fit ε_obs(Em), the directional variation introduced here (with zero conductivity and no spatial sensitivity of the FP1 detector) is already significant. Increasing the spatial resolution n does not change the results.

Taking into account the local topography at the equivalent 8 km scale (q-topo with q = 64) slightly increases the slope down to ε_F = 0.70 at footprint #5, with null thermal inertia and no spatial sensitivity taken into account. Improving the spatial resolution of the GTM down to the kilometre scale (q = 512) does not significantly change T_F(ε_F) and ε_F (Em) in this observation. Even at southern latitudes, the variation of emissivity as a function of q ≥ 64 does not exceed 1%.

On the other hand, the introduction of the spatial sensitivity of FP1 reduces the directional emissivity for Em ≥ 35 deg (footprints #4 and #5) as it appears to reduce the radiometric weight (w_j) of hotter plates, reducing the total specific intensity, I_F, at larger wavenumbers (Fig. 5). The apparent temperature decreases only slightly, while the apparent emissivity increases more significantly. Setting the thermal inertia to 15 MKS naturally decreases the temperature by about ~15 K, as expected for a scan performed around noonlocal time, while ε_F(Em) reaches its highest slope with a 27% variation between normal viewing and Em ~ 55 deg. Increasing the thermal inertia may help better fit ε_obs(Em) but it would even more decrease the temperatures at noon and widen the shift between T_F and T_obs (Fig. 7a).

Accounting for the dispersion in local times in the footprints with a tessellated ellipsoid (n-ell shape model) definitely impacts the directionality of the apparent emissivity. Including the 8 km-scale topography has a second-order influence on this. Improving this q-resolution by a factor of 8 does not have much impact on the results of the model at Em ~ 55 deg. Taking into account the spatial sensitivity of the FP1 detector is important. Finally, we notice that with the chosen parameters corresponding to reasonable a priori in terms of Bond albedo and thermal inertia, we can only approximate the data obtained while assuming uniform regolith properties. However, we can globally explain the behaviour of T_F (ε_F) and ε_F (Em), while introducing a tessellated surface and releasing the condition ε_obs = 1.

Fig. 5 Observed Iobs, fitted εobsI(Tobs) and modelled spectra IF for footprints sampled during the FP1 observation CIRS_222DI_DIONE001_PIE (see Figs. 6, 7). The modelled spectra are calculated for amorphous ice with grain size, R = 2 mm, and regolith properties A = 0.63, p = 0.5, and KE0 = 0 (Γ=0) or KE0= 0.001 W m−1 K−1 (Γ ~ 15). The effect of the spatial sensitivity of FP1 is included (+SS) or not. The Dione GTM at spatial resolution q = 64 is used. Top: footprint #3 located at latitude and west longitude (−30.89°, 123.84°). Bottom: footprint #5 located at latitude and west longitude (–50.40°, 130.02°).

Fig. 6 Map of the leading hemisphere of Dione at 154 m resolution provided by the USGS. The lens-shaped region of expected MeV electrons bombardment appears slightly shaded, centered about the equator with a maximum latitudinal extent of 70°. FOV footprints of the CIRS FP1 during observation CIRS_222DI_DIONE001_PIE are plotted in gradually lighter colour from the start to the end of the 3h35 min sequence. One in twenty is printed over 350 footprints.

3.2 Effect of thermal conductivity and grain size

The grain size controls the energy balance at the surface and the thermal cycle by affecting the average hemispherical emissivity. The less emissive micro-metric grains are the hottest, given the Bond albedo A = 0.63 (Fig. 3b, c). For a constant thermal conductivity k_E0 and a given shape model, when varying R from 0.2 μm to 2 cm (Fig. 8-Top a, b), the modelled temperature T_F is too cold, except for micro-metric grains. A regolith made of grains of a few hundred microns to a few millimetres has an apparent emissivity close to that observed. For micro-metric or sub-micro-metric grains, T_F agrees with T_obs but ε_F remains much lower than ε_obs (Fig. 8-Top c).

If the thermal conductivity is both radiative and conductive (Fig. 8-Bottom), the grains must be composed of amorphous water ice for the model to yield thermal inertia ranging between 2 and 23 MKS, in cases of tight contacts and p = 0.7 (Ferrari & Lucas 2016). Introducing size and temperature dependent thermal conductivity directly impacts T_F (Fig. 8-Bottom a, b) but only slightly changes the apparent emissivity ε_F. This latter is mainly driven by the dispersion in local hour angle (Fig. 8-Bottom c versus Fig. 8-Top c). Such conductivityintroduces a large discrepancy in the thermal history among plates covered with grains larger than a few hundred of μm: the thermalconductivity is indeed dominated by radiative conductivity which scales with R: the larger R the larger k_E ~ k_R. Therefore, cm-sized grains have larger thermal inertia and their temperature decreases near noon local time (Fig. 8-Bottom a, b). The thermal conductivity is the lowest for grain size in the range 1–100 μm, so that their temperature is the highest at local noon. It increases again for very small grains as conduction through loose contacts scales with 1∕R and dominates the radiative conductivity. In this case, the modelled temperature is too cold compared to the measurements.

Finally, the question of regional variations in regolith properties can be also addressed with this model, for example a regional variation in grain size, which may originate from sintering after bombardment by MeV electrons in the lens-shaped region of icy moons leading hemispheres (Fig. 6). We considered here a regolith composed of 20 μm-sized icy grains covering the leading hemisphere outside of the central 70°-wide lens-shaped region, which is itself covered by 2 mm-sized grains made of amorphous ice. The choice of the size is dictated by the temperature and the emissivity observed at the two extremes of the scan, near the equator and at southern latitudes. At southern latitudes outside the oval shape, T_obs ~ 80 K and ε_obs ~ 0.55, whereas at the equator T_obs ~ 95 K and ε_obs ~ 0.92 (Fig. 8-Top, Fig. 3). We performed four simulations assuming either A =0.63 (Pitman et al. 2010; Howett et al. 2010) or A = 0.49 (Howett et al. 2014) with either null or constant thermal inertia (k_E0 = 0.001 W m⁻¹ K⁻¹) and porosityp = 0.8 (Fig. 9). When taking into account the regional variability and the local topography at 8 km-scale (GTM with q = 64) together with the FP1 spatial sensitivity, the modelled ε_F(Em) is close to ε_obs(Em) in all cases (Fig. 9b). This is only slightly sensitive to the Bond albedo and the low thermal inertia. The apparent temperature T_F is best reproduced by a regolith with either albedo A = 0.49 and Γ = 10 MKS or A = 0.63 and Γ = 0 MKS. Considering regional variations in regolith properties around the lens-shaped region of the leading hemisphere, on the Dione moon particularly, seems to help in reproducing the observed apparent temperatures and emissivity, whereas a single population can hardly accomplish this (Fig. 8). Of course, an estimation of the physical properties of these regions with the model should be carried out for each of Saturn’s moons by analysing as many relevant observations as possible.

Fig. 7 Modelled temperatureTF and apparent emissivity εF versus topographic models: (a) T(ε) (b) ε(Em). In the simulations A = 0.63, p = 0.5, R = 2 mm and the ice phase is amorphous. The thermal conductivity is constant, either kE0 = 0 (Γ = 0 when not specified in the legend) or kE0 = 0.001 W m−1 K−1 (Γ ~ 15 MKS). TF and εF are calculated for five footprints (#1 to #5) regularly distributed along the emission angle Em. The observed (Tobs, εobs) in the 350 footprints of observation CIRS_222DI_DIONE001_PIE are plotted (+). The topographic models are: flat for the flat case assuming average hemispherical emissivity εh (Eq. (13)), flat+ed for the flat case assuming directional emissivity ϵd (Eq. (14)), n=222 for the ellipsoidal n-ell model, q=64, 128 or 512 for the q-topo model at various spatial resolutions of the GTM. The effect of the spatial sensitivity (SS) of FP1 is shown for the GTM with q = 64 and kE0 = 0 (q=64+SS) or for kE0 = 0.001 W m−1 K−1 and p = 0.5 (q=64+SS, Γ = 15).

Fig. 8 Modelled (TF, εF) and observed (Tobs, εobs, +) apparent temperature and emissivity in the FP1 FOV during observation CIRS_222DI_DIONE001_PIE vs incidence (Inc) or emission (Em) angles. The bolometric Bond albedo A = 0.63 and the ice phase is amorphous. The topography of Dione GTM (q = 64) and the spatial sensitivity of FP1 are included. The emissivity of surface elements is assumed to be isotropic (Eq. (13)). Top: Modelled (TF, εF) versus grain size, R, for thermal inertia Γ ~ 15 MKS assuming constant thermal conductivity kE = 0.001 W m−1 K−1 and p = 0.5. Bottom: Modelled (TF, εF) versus grain size, R, for an effective conductivity being both radiative and conductive. The contacts are tight and p = 0.7 so that the thermal inertia is in the range 2–23 MKS.

Fig. 9 Modelled (TF, εF) and observed (Tobs, εobs, +) apparent temperature and emissivity as a function of Bond albedo A and thermal inertia Γ in the case of variations of regolith properties about a lens-shaped feature centred on the leading hemisphere, with maximum latitudinal width of 35° (shaded region, Fig. 6). Inside the oval, the grain size is R = 2 mm and outside of it R = 20μm. The Bond albedo is uniform, either A = 0.49 or A = 0.63. The thermal conductivity of amorphous ice is constant, either kE0 = 0 or kE0 = 0.001 W m−1 K−1 with porosity p = 0.8 (Γ = 10 MKS) for both regions. The topography of GTM (q = 64) and the spatial sensitivity of FP1 are included. The emissivity of surface elements is assumed to be isotropic (Eq. (13)).

4 Conclusion

In this paper, we present a thermophysical model with the aim of analysing the temperatures and the apparent emissivities of Saturn’s icy moons observed with CIRS FP1. We include their topography and a parametrisation via the regolith physical properties such as grain size, water ice phase, porosity, and the type of grain contact. We conclude that:

• 1.

  The dispersion in local times of plates within the footprints has a definite impact on the directional anisotropy of the apparent emissivity. Including the 8 km-scale topography has a second-order influence on this. Improving this resolution down to the 1 km-scale has small impact on both apparent temperature and emissivity at Em ~ 55 deg. Taking into account the spatial sensitivity of the FP1 detector is important. The directional anisotropy of the apparent emissivity cannot be explained by the directional emissivity ε_d(R, μ) at the scale of the plate as proposed by Hapke (2012).

• 2.

  The model can account for regional variations in regolith properties and test current hypotheses on its evolution by space weathering, such as a change in grain size between a lens-shaped region and southern and northern latitudes, which may be covered with smaller grains.

• 3.

  Considering the hemispherical emissivity varying with grain size and the topography to analyse CIRS FP1 observations of icy moons with this model will yield constraints on the physical properties of their surface and may lead to some changes in the current estimates of their thermal inertia and their regional variations.

This model is adapted to address the question of the actual properties of atmosphere-less planetary surfaces evolving under diverse weathering processes by analysing their infrared emission. The analysis of CIRS FP1 observations of Mimas, Dione, and Rhea is underway in order to test current hypotheses on regolith processing by space weathering, such as regolith changes in grain size.

Acknowledgements

This work was supported by the Centre National d’Etudes Spatiales, thanks to Planètes et petits corps du Système solaire and TOSCA programs, through proposals CIRS-Anneaux et Surfaces and more recently ADSINTHE (Anisotropie Directionnelle des Sols en INfrarouge THErmique). The polyscope software has been used to generate the GTM rendering with arbitrary illumination and viewing geometry (Nicholas Sharp et al. 2019, www.polyscope.run).

References

All Tables

All Figures

	Fig. 1 Representation of the Gaskell Global Topography Model (GTM) of Dione. Southern latitudes of the leading hemisphere are visible about longitude 90°W with the Evander basin on the lower left near the limb and the Sabinus crater below the center.
In the text

	Fig. 2 Optical constants of hexagonal water ice (Ih) and tholins (from Khare et al. 1984). Both real (nr) and imaginary (nk) components of the refractive indices are plotted.
In the text

Fig. 3 Modelled regolith properties versus grain size R (a) Bolometric Bond albedo A (Eq. (9)) versus grain size, R, for various tholins fractions (0%, 5%, 20%). The range of bolometric Bond albedos for Saturn’s icy moons is marked as horizontal dotted lines (Pitman et al. 2010). (b) Average hemispherical emissivity versus grain size, R, for various tholins fractions (0%, 5%, 20%) and surface temperature TS = 80 K (Eq. (11)). The error bars figure the dispersion of when TS varies from 50 K to 110 K for a tholins fraction of 5% in crystalline water ice grains. The curves are similar in cases of amorphous water ice (Eq. (11)). (c) Typical surface temperature Tslow of a plate covered with regolith of grain size, R, for various tholins fractions (0%, 5%, 20%) in the case of a slowly rotating moon.

In the text

Fig. 4 Thermal cycles of a surface element located at the center of the leading hemisphere of Dione. The bolometric Bond albedo A is assumed 0.63 (Pitman et al. 2010). (a) The thermal conductivity is constant kE0 = 0.001 W m−1 K−1 and the grain size, R, varies between 0.2 μm and 2 cm, assuming a porosity of p = 0.3 (solid line) or 0.8 (dash-dot line). The ice phase is amorphous. The resulting thermal inertia range between 10 and 20 MKS. (b) Both radiative and solid conductions through loose contacts are assumed here. The grain size, R, varies between 0.2 μm and 2 cm and the ice phase is either crystalline (Ih, solid line) or amorphous (Am, dash line) ice phase. The porosity is p = 0.5.

In the text

Fig. 5 Observed Iobs, fitted εobsI(Tobs) and modelled spectra IF for footprints sampled during the FP1 observation CIRS_222DI_DIONE001_PIE (see Figs. 6, 7). The modelled spectra are calculated for amorphous ice with grain size, R = 2 mm, and regolith properties A = 0.63, p = 0.5, and KE0 = 0 (Γ=0) or KE0= 0.001 W m−1 K−1 (Γ ~ 15). The effect of the spatial sensitivity of FP1 is included (+SS) or not. The Dione GTM at spatial resolution q = 64 is used. Top: footprint #3 located at latitude and west longitude (−30.89°, 123.84°). Bottom: footprint #5 located at latitude and west longitude (–50.40°, 130.02°).

In the text

Fig. 6 Map of the leading hemisphere of Dione at 154 m resolution provided by the USGS. The lens-shaped region of expected MeV electrons bombardment appears slightly shaded, centered about the equator with a maximum latitudinal extent of 70°. FOV footprints of the CIRS FP1 during observation CIRS_222DI_DIONE001_PIE are plotted in gradually lighter colour from the start to the end of the 3h35 min sequence. One in twenty is printed over 350 footprints.

In the text

Fig. 7 Modelled temperatureTF and apparent emissivity εF versus topographic models: (a) T(ε) (b) ε(Em). In the simulations A = 0.63, p = 0.5, R = 2 mm and the ice phase is amorphous. The thermal conductivity is constant, either kE0 = 0 (Γ = 0 when not specified in the legend) or kE0 = 0.001 W m−1 K−1 (Γ ~ 15 MKS). TF and εF are calculated for five footprints (#1 to #5) regularly distributed along the emission angle Em. The observed (Tobs, εobs) in the 350 footprints of observation CIRS_222DI_DIONE001_PIE are plotted (+). The topographic models are: flat for the flat case assuming average hemispherical emissivity εh (Eq. (13)), flat+ed for the flat case assuming directional emissivity ϵd (Eq. (14)), n=222 for the ellipsoidal n-ell model, q=64, 128 or 512 for the q-topo model at various spatial resolutions of the GTM. The effect of the spatial sensitivity (SS) of FP1 is shown for the GTM with q = 64 and kE0 = 0 (q=64+SS) or for kE0 = 0.001 W m−1 K−1 and p = 0.5 (q=64+SS, Γ = 15).

In the text

Fig. 8 Modelled (TF, εF) and observed (Tobs, εobs, +) apparent temperature and emissivity in the FP1 FOV during observation CIRS_222DI_DIONE001_PIE vs incidence (Inc) or emission (Em) angles. The bolometric Bond albedo A = 0.63 and the ice phase is amorphous. The topography of Dione GTM (q = 64) and the spatial sensitivity of FP1 are included. The emissivity of surface elements is assumed to be isotropic (Eq. (13)). Top: Modelled (TF, εF) versus grain size, R, for thermal inertia Γ ~ 15 MKS assuming constant thermal conductivity kE = 0.001 W m−1 K−1 and p = 0.5. Bottom: Modelled (TF, εF) versus grain size, R, for an effective conductivity being both radiative and conductive. The contacts are tight and p = 0.7 so that the thermal inertia is in the range 2–23 MKS.

In the text

Fig. 9 Modelled (TF, εF) and observed (Tobs, εobs, +) apparent temperature and emissivity as a function of Bond albedo A and thermal inertia Γ in the case of variations of regolith properties about a lens-shaped feature centred on the leading hemisphere, with maximum latitudinal width of 35° (shaded region, Fig. 6). Inside the oval, the grain size is R = 2 mm and outside of it R = 20μm. The Bond albedo is uniform, either A = 0.49 or A = 0.63. The thermal conductivity of amorphous ice is constant, either kE0 = 0 or kE0 = 0.001 W m−1 K−1 with porosity p = 0.8 (Γ = 10 MKS) for both regions. The topography of GTM (q = 64) and the spatial sensitivity of FP1 are included. The emissivity of surface elements is assumed to be isotropic (Eq. (13)).

In the text