A continuum of icy satellites’ radar properties explained by the coherent backscatter effect

The radar properties of icy satellites of Jupiter and Saturn commonly differ by more than an order of magnitude from those of rocky planets because of the lower absorptivity of water ice than that of rock. However, the specific mechanisms behind these differences are not confidently known yet. Here we show that the global radar albedos and the circular polarization ratios of icy satellites are correlated and vary together along a continuum, which is not satisfactorily explained by existing models. We modified a backscatter model that includes the coherent backscatter opposition effect (CBOE) and found that it fits successfully the observed continuum, indicative of CBOE’s primary role in driving the circular polarization ratios of icy satellites. We predict that the linear polarization ratios and the incidence angle variation properties also vary along the observed continuum, and that, globally, Europa and Enceladus have very heterogeneous surfaces/subsurfaces with very little microwave absorption. The modified CBOE model suggests that some worlds beyond Saturn’s orbit that are rich in non-water ices may have even more extreme radar properties. The unusual radar properties of icy satellites appear to be correlated and distributed along a spectrum of values. Only modelling including the coherent backscatter opposition effect can reproduce this behaviour.

The icy satellites of Jupiter and Saturn are markedly different radar scatterers than rocky worlds (Mercury, Venus, solid Earth, Earth's moon, Mars and large asteroids) and smaller bodies such as asteroids and comets 1,2 . Four salient differences are (1) much greater global radar albedos, (2) much greater circular polarization ratios, (3) greater linear polarization ratios and (4) much weaker albedo variation with incidence angle. These differences are illustrated and quantified in Fig. 1 and Supplementary Table 1. The radar albedo is the normalized radar cross-section (NRCS). The first difference is related to the integrated areas under the curves in Fig. 1 (the global, or disk-integrated, NRCS is the integral with extra angular factors and the sum of both polarizations); the total area under the Europa curves is much greater than that of the Moon (Supplementary Table 1). The circular polarization ratio (μ C ) is the NRCS in the same circular (SC) polarization as that transmitted divided by the NRCS in the opposite circular (OC) polarization. Thus, the second difference is related to the ratio of the SC curve to that of the OC curve in Fig. 1; the icy-satellite radar backscatter is notable in that μ C > 1, whereas most observed surfaces in the solar system have μ C << 1 (Supplementary Table 1). This distinction is perhaps the most surprising and challenging for scattering models to explain. The linear polarization ratio (μ L ) is the NRCS in the opposite linear (OL) polarization as that transmitted divided by the NRCS in the same linear (SL) polarization. Both μ C and μ L are defined such that they are zero for a smooth, plane mirror. Icy satellites have global μ L ≈ 0.5, which is considerably larger than that of rocky worlds (Supplementary Table 1). The fourth difference is related to the shape of the OC curves in Fig. 1. Radar backscatter from most observed solar system surfaces is dominated by a quasi-specular component in the OC and SL polarizations that is relatively sharply peaked near 0° incidence (perpendicular to surface). A weaker diffuse component that varies much more gradually with incidence angle is observed in all polarizations. The quasi-specular component is well described by models with mirror-like reflections from nearly smooth surfaces [3][4][5][6] . The diffuse component is generally attributed to very rough surfaces and/or subsurface scattering 3-6 . Article https://doi.org/10.1038/s41550-023-01920-2 effect (CBOE, also called weak localization) from constructive interference by identical but opposite-in-direction, multiple-scattering paths, when the transmitter-target-receiver angle is zero, could fit all four radar properties, is observed in laboratories and does not require fine-tuned subsurface structures [18][19][20] . However, the upper limit of CBOE NRCS enhancement is two, and as such the utility of the CBOE model for explaining the very high global NRCS of icy satellites has been questioned (icy-satellite global NRCS can exceed that of rocky worlds by more than an order of magnitude; Supplementary Table 1), as has the consistency of assumptions in the planetary model 18 of the effect 14 . An exact solution for electromagnetic scattering by cylinders was applied to simulate ice pipes and ice lenses (approximately ice cylinders) embedded in snow in Greenland's percolation zone and shown to approximately match radar observations 10 . That model may be less applicable to icy satellites, however, because backscatter from Greenland's percolation zone has a stronger incidence-angle dependence than that of icy satellites. In addition, the formation process for terrestrial ice pipes/lenses involves solar-driven melting of snow, followed by refreezing, which should not operate on icy satellites 10 . Thus, some models may be preferred, but the specific cause(s) of the extraordinary radar properties of icy satellites, beyond simply subsurface scattering, is yet to be determined.

A continuum of icy satellites' radar properties
In this section, we demonstrate that two radar properties, global NRCS and μ C , of icy satellites have a continuum, in the sense that both properties span a range and they are correlated. Thus, the theoretical problem of identifying a scattering mechanism that explains the extraordinary icy satellite radar properties can be framed more specifically for global NRCS and μ C : rather than evaluating models on the basis of whether they can produce very high global NRCS and μ C , models should now also be evaluated on the basis of whether they can predict/fit the observed NRCS-μ C continuum. Figure 2 shows the global NRCS and μ C of icy satellites of Jupiter ( Jovian) and Saturn. The figure demonstrates that icy satellites are distributed in both global NRCS and μ C and that these two radar properties are related. We emphasize the importance of this new observational result as an improved perspective from which to understand these icy satellite radar properties: global NRCS and μ C of icy satellites are better described as variable along a continuum, which includes a very high extreme, than as an anomalous group.
A pronounced difference between the 12.6 cm disk-integrated NRCS and μ C of the leading and trailing hemispheres of Enceladus was reported 2 . Those differences are primarily due to differences in the SC polarization. No SC echo was detected from Enceladus' trailing hemisphere, whereas the leading hemisphere has an extremely high disk-integrated SC NRCS, and the OC echoes from both hemispheres were similar. The striking SC hemispherical difference is suspicious given the aforementioned similar OC echoes, lack of evidence for obvious hemispherical differences in 2.2 cm SL radar observations of Enceladus 21,22 and lack of significant hemispherical differences of radar properties of other icy satellites. However, a hemispherical difference in 2.2 cm thermal emission was observed and argued to be consistent with the 12.6 cm radar difference 23 . Nevertheless, only measurements of Enceladus' leading hemisphere from Black et al. 2 are plotted in Fig. 2; trailing-hemisphere measurements are also plotted in Supplementary Fig. 1.
In addition to the conclusion that global NRCS and μ C of icy satellites have a continuum, which is a new constraint for models, several other aspects of Fig. 2 are discussed. The figure suggests a modest difference between the NRCS-μ C relationships of the Jovian and Saturnian icy satellites. This difference is unlikely to be due to statistical noise because all four Saturnian icy satellites with similar global NRCS to a Jovian icy satellite have lower global μ C (Methods), although a In contrast to rocky surfaces, the icy satellites exhibit a gentle NRCS variation at all incidence angles in all polarizations. Icy-satellite radar backscatter is dominated by a diffuse component (similar in shape, but not magnitude, to the diffuse component from rocky surfaces) with little to no quasi-specular component.
A few other surfaces in the solar system exhibit radar properties that, to varying degrees, resemble these extraordinary radar properties of icy satellites (global-NRCS, μ C , μ L and incidence angle variation), including polar ice on Mars 7 and ice in permanently shadowed regions of Mercury's poles 8,9 . Radar scattering from the percolation zone of the Greenland ice sheet and some mountain glaciers bears some resemblance to that from icy satellites, but most terrestrial ice, including other zones of the Greenland ice sheet, does not 10,11 . We are unaware of any peer-review publications of radar measurements of surfaces beyond the Saturn system.
The markedly different radar properties of icy satellites, compared with those of rocky worlds, is very likely to be due to weaker microwave absorption by ice than rock, which permits greater radar propagation and subsurface backscatter. The specific mechanism(s) that produces such extraordinary radar properties, however, has not previously been confidently determined. At least six models have been proposed. A simulation based on large numbers of random ice facets (vacuum-ice interfaces) was claimed to account for the observed NRCS and μ C of Ganymede, although an outcome that matched both was not shown 12 . Refraction scattering, instead of reflection scattering, was shown to be capable of producing extremely large SC NRCSs; however, it may not produce the large OC NRCSs observed 13,14 . Mode decoupling during refraction scattering was shown to match the observed global NRCS, μ C and μ L ; however, it requires lenses with fine-tuned refractive indexes that generally decrease with radius, for which there is little to no other evidence and limited geologic plausibility 15 . Scattering from abundant, buried craters could match all four extraordinary radar properties of icy satellites, but similarly requires special refractive and geometric structures that may be unlikely in high abundance 16 27,36 . The substantial differences in global NRCS (related to sum of OC and SC integrals), circular polarization ratio (SC/OC, both integrated and at all incidence angles) and angular variation (shape of curve) are discernable. Although it is not perceptible on this scale, the lunar SC NRCS is similar in shape to the SC NRCS of Europa.
Article https://doi.org/10.1038/s41550-023-01920-2 systematic error is possible. Physical and/or chemical differences between the surfaces/subsurfaces of icy satellites in the two systems could also be responsible. Figure 2 also suggests that the global NRCS distribution of icy satellites is discretized rather than continuous, with Iapetus/Titan, Callisto/Dione, Ganymede/Rhea/Tethys and Europa/ Enceladus in four distinct groups. If true, this would be yet another extraordinary radar property of the icy satellites because fundamentally discretized properties are rare in remote sensing, especially for global averages. As such, we doubt the apparently distinct groups reflect physical bounds of global NRCS. We also note that, in both systems, global NRCS and μ C are anticorrelated with distance from the planet. The reversed order of Dione and Rhea is the only exception to this apparent trend. Mimas was not measured but may also depart from this trend as the 12.6 cm global NRCSs of the Saturnian icy satellites are correlated to their 2.2 cm SL radar and optical disk-integrated albedos 2,22 and both albedos for Mimas are less than those of Enceladus (and Tethys), despite Mimas' closer proximity to Saturn. Global NRCS and μ C of Saturnian icy satellites may therefore also be related to accumulated deposits from Saturn's E ring. Radar measurements of the Uranian and Neptunian satellite systems should provide further insights into the planetary system differences, global NRCS discretization and relationship to the planetary distance results of Fig. 2.
In addition to the measurements shown in Fig. 2, the Jovian icy satellites were observed with radar at a wavelength (λ) of 70 cm (ref. 24 ) and Titan was observed at λ = 3.5 cm (ref. 25 ). Those observations, however, are significantly noisier and, in some cases, resulted only in upper limits (Supplementary Information) and thus are not shown in Fig. 2. The Saturnian icy satellites were also observed with the Cassini spacecraft radar at λ = 2.2 cm in the SL polarization 21,22,26 . We are unaware of other published radar observations of icy satellites.
The two other unusual radar properties of icy satellites, large μ L and diffuse-dominated variability with incidence (generally quantified by n, where NRCS ∝ cos n (incidence)), are significantly less well constrained (including, for some radar properties, icy satellites, and/ or wavelengths, no published measurements, to our knowledge; Supplementary Table 1). Investigating similar correlations using these properties is a possible area for future research.

Coherent backscatter opposition effect as a preferred model
The continuum of icy satellites' radar properties introduced in the previous section offers a new observational test of backscatter models. All six previously published models fail to satisfactorily fit the observed continuum (Supplementary Information, including Supplementary  Fig. 2). However, one model, a scalar radiative transfer model with CBOE 18 , with a slight but important adaptation, does fit the observed continuum. Further research of the applicability of the model based on the exact solution for electromagnetic scattering by a dielectric cylinder 10 is also warranted (Supplementary Information).
We modified the backscatter model of Hapke 18 to: where NRCS OC and NRCS SC are the global NRCS in the OC and SC polarizations, NRCS S and NRCS M are the global NRCS from single and multiple scattering and E OC and E SC are the CBOE enhancement factors in the two circular polarizations. The adaptation from the model of Hapke 18 is the separation of E OC and E SC into independent parameters, instead of E OC = E SC = E, which is justified by laboratory observations that CBOE differs for each polarization 19 . NRCS S is a function of the single-scattering albedo (w) and the phase function, which includes an amplitude constant (b), whereas NRCS M is a function of w only; the equations for NRCS S and NRCS M are reproduced in the Methods. Thus, the model has four free parameters: w, b, E OC and E SC . The above model fits the icy-satellite NRCS-μ C continuum with w varying among icy satellites and constant b, E OC and E SC (Methods). The Jovian and Saturnian systems were separately fit with a modest change of b (E OC and E SC were also permitted to change but intriguingly were found to be similar for both systems). The best-fits are shown in Fig. 3 and best-fit values and uncertainties are provided in the Supplementary Information, including Supplementary Table 2.
In the context of the above CBOE model, global NRCS and μ C variations among icy satellites (and also differences between icy satellites and rocky worlds) are primarily attributed to variations of single-scattering albedo (w). Single-scattering albedo is the ratio of the scattering coefficient to the extinction (scattering plus absorption) coefficient. The large w of icy satellites (Supplementary Information) are consistent with the conclusion that icy satellites have much greater global NRCS (and μ C ) than rocky worlds primarily because of the lower radar absorptivity of ice than rock, which permits increased backscatter from multiple scattering. The variation of w among icy satellites could be due to variation of absorptivity (for example, as a result of variation of composition and/or abundance of non-ice material) and/or variation of concentration of scatterers (for example, as a result of variation of concentration of fractures and/or ice boulders). Recall that CBOE does not require fine-tuned subsurface structures; it occurs for scatterers of various shapes, sizes and separations. The curves in Fig. 3 show 0 ≤ w ≤ 1, where the limits correspond to the theoretical minimum and maximum; mathematical continuation of the curves in either direction is not physical. The best fit for Europa and Enceladus is with w of unity (Supplementary Information), which suggests that they have very little absorbing material and/or a very high concentration of scatterers. That the previous statement is true for Europa at both 3.5 cm and 12.6 cm spatial scales (wavelengths of radar observations in Fig. 2) indicates Article https://doi.org/10.1038/s41550-023-01920-2 that minimal absorption is a significant component of Europa's very high w. The dubious possibility that the icy satellites' NRCS-μ C relationship is discrete rather than continuous is again suggested by the best-fit single-scattering albedos: icy satellites in the same global NRCS groups are best fit by nearly identical w (Supplementary Information).
Furthermore, in the context of the above CBOE model, the apparent, modest difference in the radar-properties continua of Jovian and Saturnian icy satellites is primarily attributed to differences of the single-scattering phase function. The phase function describes the angular dependence of scattered power. The simple phase function of this model (Methods) has one free parameter, b, where b < 0 corresponds to greater forward than backward scattering (vice versa for b > 0) and the magnitude of the difference is proportional to |b|. The best-fit values (Supplementary Table 2) indicate that radar single scattering by observed icy-satellite surfaces is, on average, preferentially forward, and more so for Jovians than Saturnians. Phase function differences result from many causes, including differences of size, shape and composition of scatterers. It is surprising, given geologic differences among icy-satellite surfaces in each planetary system, that all observed icy satellites in each system are satisfactorily fit with the same b (and the same E OC and E SC ; Methods). That all observed icy satellites in a system can be fit with the same single-scattering phase function and CBOE enhancement factors is, within the framework of the model, the reason for the NRCS-μ C continua. The consistency in each planetary system, yet difference between the two systems, suggests a planetary-system control of the icy-satellite phase functions. Two noteworthy constraints on hypotheses for planetary-system control of phase functions are that they apply at spatial scales probed by the radar observations (including depths of ≳10-100λ (ref. 20 ) for both 3.5 cm and 12.6 cm observations), and they result in little to no leading-trailing hemispherical differences of radar properties 2,27 . Differences of impurity composition as a result of intrinsic compositional differences between the Jovian and Saturnian systems is one hypothesis for the planetary-system phase function difference that is consistent with these constraints.
The best-fit CBOE enhancement factors for the Jovian and Saturnian systems are similar (Supplementary Table 2). This similarity suggests that the values, and/or their ratio, may be universal characteristics of radar scattering from icy satellites (and possibly even of CBOE from worlds and/or possibly of CBOE from ice). Notably, the best-fit CBOE enhancement factors are considerably less than their theoretical maximum (Supplementary Table 2 and Methods). Further experimental, observational and theoretical investigations of these parameters is warranted.
The results of our icy satellite radar-properties continuum test of models ( Fig. 3 and Supplementary Information) indicate that CBOE occurs in 3.5 cm and 12.6 cm radar scattering from icy satellites. CBOE, however, can increase the NRCS only by a maximum factor of two 18 , and thus cannot be the primary reason for the approximately order of magnitude greater global NRCS of icy satellites than that of rocky worlds ( Fig. 1 and Supplementary Table 1). Recall that in equations (1) and (2), CBOE is the enhancement factor for each polarization (E OC , E SC ); the rest of the equations are one particular scalar radiative transfer model. From here onward, we consider CBOE and the radiative transfer model on which it was appended as two separate models. For example, equations (1) and (2) can fit the NRCS-μ C continuum without underlying equations for NRCS S and NRCS M (that is, free parameters NRCS S , NRCS M , E OC and E SC ) or with various underlying equations; each instance is a different model merged with CBOE. Thus, CBOE can fit the observed μ C of icy satellites and also contribute to fits of NRCS by other models. The results of our icy satellite radar-properties continuum test indicate that the CBOE model is preferred, but are less diagnostic of the Hapke 18 scalar radiative transfer model on which it was appended. On the basis of these and previous results [18][19][20] , we conclude that CBOE is the primary reason for the extraordinary circular polarization ratios (and probably linear polarization ratios) of icy satellites. The lower radar absorptivity of ice than rock, multiple scattering and CBOE all contribute to the extraordinary global NRCSs and diffuse-dominated NRCS variability with incidence angle.

Discussion
CBOE occurs whenever multiple scattering contributes to backscatter; it increases as the multiple scattering contribution increases. Thus, in the context of the CBOE model, the extraordinary radar properties of icy satellites are attributed to a very large increase of the multiple scattering contribution, compared with that from rocky worlds, rather than a fundamentally distinct scattering mechanism.
Laboratory results indicate that CBOE also differs between linear polarizations 19 . For both circularly and linearly polarized radiation, CBOE is greater for the same polarization as the incident radiation. Equations (1) and (2) are adaptable to linear polarizations; therefore, we predict that icy-satellite global linear polarization ratios are also correlated to their global radar albedos. The NRCS-μ L relation, however, is predicted to differ from the NRCS-μ C continuum because the definition of μ L is inverted relative to that of μ C (OL/SL versus SC/OC), whereas, for both cases, the CBOE enhancement factor is greater for the same polarization as that transmitted (SL and SC). If incidence angle variation likewise depends on single-scattering albedo, which it does in the Hapke 18 scalar radiative transfer model, then that variation is also predicted to be correlated to global NRCS. Thus, the conclusion that icy satellites have a NRCS-μ C continuum is predicted to extend to their extraordinary radar properties more generally. In addition to improving observational constraints of icy-satellite linear polarization ratios and incidence angle variations, another radar measurement that is likely to be distinguishing is bistatic observation (separated transmitter and receiver) of icy-satellite surfaces. The CBOE model predicts that CBOE occurs only near opposition 18,19 , a strong dependence on observation geometry. Bistatic observations of icy satellites at microwave wavelengths may be possible with the planned Europa Clipper and/or Jupiter icy moons explorer ( JUICE) missions to the Article https://doi.org/10.1038/s41550-023-01920-2 Jovian system using the spacecraft and terrestrial communications systems 28 . The success of the CBOE model implies that multiple scattering contributes considerably (dominantly in most cases) to radar backscatter from icy satellites, which indicates marked heterogeneity (that is, scatterers) of icy-satellite surfaces/subsurfaces at 3.5 cm and 12.6 cm scales (λ) to depths of ≳10-100λ (ref. 20 ). That decimetre-deep to greater than metre-deep heterogeneity on icy satellites could be chaotic, such as jumbled ice pebbles/cobbles of various sizes, or organized, such as regularly spaced faults/folds. Indeed, widespread tectonic deformation on icy satellites 29,30 may extend to the wavelength scale of the radar observations. In situ investigations of icy satellites, such as sampling beneath the surface 31 , should be designed to accommodate significant heterogeneity.
The high global radar albedos, combined with the implication of considerable multiple scattering, indicates little to no microwave absorption by the heterogeneous surfaces/subsurfaces of icy satellites, and effectively zero absorption in the extreme cases of Europa and Enceladus for the modified scalar radiative transfer model of Hapke 18 . That weak absorption suggests a high degree of ice purity. Furthermore, the weak absorption suggests that non-ice material is generally smaller than λ; if the surface/subsurface includes non-ice grains (for example, salts or silicates), they must be sparse and smaller than surrounding ice grains. We speculate that ice purity may be related to the extent of resurfacing by material from a liquid reservoir, including oceans.
We note that the radar properties of Europa and Enceladus are not the extreme maxima predicted by the modified CBOE model; they are maxima for particular phase functions and CBOE enhancement factors. Greater global-NRCS, μ C and μ L , and more exotic incidence angle variations, are theoretically possible with other phase functions, enhancement factors and underlying radiative transfer models. We speculate that worlds rich in other (non-water) ices, such as nitrogen-ice and/or methane-ice surfaces/subsurfaces at greater heliocentric distances than Saturn, may have even more extreme radar properties. Some such surfaces are known to be exceedingly bright at visible wavelengths 32 , like icy satellites 33 . Observation of such surfaces is a frontier for planetary radar.

Probability that different Jovian and Saturnian continua are noise
If the Jovian and Saturnian icy satellites had the same NRCS-μ C continuum, then a Saturnian satellite with a global NRCS similar to that of a Jovian satellite would have a 50% chance, owing to statistical noise, of having a measured global μ C of less than that of the Jovian satellite. The probability that all four Saturnian satellites with similar global NRCS to a Jovian satellite would consistently have either lesser or greater measured global μ C than those of the Jovian satellites of similar global NRCS, which would suggest separate continua, is ≈ 0.5 3 = 0.125. Although not a tiny probability, it indicates that the different continua are unlikely to be due to statistical noise. If the argument were made for the Jovians in comparison to the Saturnians, instead of the reverse as above, both the 3.5 cm and 12.6 cm measurements would be included and the probability would be conservatively estimated to be ≈ 0.5 2 × 0.5 2 = 0.0625, so 0.125 is a more conservative estimate.

Modified Hapke (1990) CBOE model
The modified Hapke 18 CBOE model is given by equations (1) and (2) with NRCS M = (2r 0 + 2r 2 where r 0 is the diffusive reflectance 34 and the other parameters are defined in the main text. Note that a factor of four has been included in the transformation from geometric albedo given in Hapke 18 to global NRCS 35 . A single-scattering phase function of p = 1 + bcos(g), where g is the phase (transmitter-target-receiver) angle, was assumed in Hapke 18 . For this phase function, b > 0 corresponds to greater backscattering and b < 0 to greater forward scattering.

Best-fit modified Hapke (1990) CBOE model
The best-fit model is defined to be the model with parameters that minimize where i denotes an individual icy satellite and the summation is of all modelled icy satellites, which for Fig. 3 is all measured icy satellites in the Jovian or Saturnian systems. NRCS i and μ C,i are the measured global radar albedo and circular polarization ratio of an individual icy satellite, and σ NRCS,i and σ μ C ,i are the uncertainties for those values, respectively. NRCS SC,i and NRCS OC,i are the modelled global NRCSs in the SC and OC polarizations, from equations (1) and (2). Equation (6) treats the global NRCS and μ C as independent, equal measurements in the sense that neither is assumed to depend on the other or is given greater weight.
Our current physical interpretation of the model bounds its mathematical parameters to 0 ≤ w ≤ 1, −1 ≤ b ≤ 1, 0 ≤ E OC ≤ 2 and 0 ≤ E SC ≤ 2. Hapke 18 specifies that 1 ≤ E ≤ 2; however, similar arguments to those for constructive interference as the cause of CBOE suggest the possibility of destructive interference, which would correspond to 0 ≤ E ≤ 1. Evidence for destructive interference away from perfect opposition (transmitter-target-receiver angle of exactly zero) has been observed in the laboratory 18 . If the conditions 1 ≤ E OC ≤ 2 and 1 ≤ E SC ≤ 2 are imposed, the best-fit models have similar equation (6) sums and appear qualitatively similar to those shown in Fig. 3, except that they extend to greater global NRCS and μ C .
To explore the effect of each of the four model parameters, w, b, E OC and E SC , on the modelled global NRCS and μ C , we plotted test examples in which each parameter varied while the other three parameters were constant. We found that varying w among icy satellites while the other three parameters are constant results in a NRCS-μ C relationship similar to that observed for icy satellites. Specifically, both global NRCS and μ C monotonically increase from zero to a maximum with increasing w (except for the limiting case of b = −1, which results in constant μ C ). The other three cases result in NRCS-μ C relationships that are significantly different from the observed icy satellite continuum. For example, in general, the other cases result in a negative correlation (inverse relationship) between global NRCS and μ C ; they also do not, in general, include the global NRCS = 0 limit. We also found that when all four parameters vary, the model exactly fits the global NRCS and μ C of every individual icy satellite. Nevertheless, it is notable that the model fits all icy-satellite global NRCS and μ C , within measurement uncertainties, of the Jovian or Saturnian systems with only one parameter, w, varying among icy satellites. Thus, we report and focus on best-fits where w varies among icy satellites and b, E OC and E SC are constant for all icy satellites in a system.
The uncertainties of the best-fit parameters are defined to be their standard deviations from 1,000 random simulations. Simulated icy-satellite global NRCS and μ C values were randomly generated, such that they follow normal distributions, with means and standard deviations given by the measurements and measurement uncertainties. The best-fit model for each simulated continuum was determined following the method described above. We verified that frequency distributions of the best-fit parameters for the simulated continua are approximately normal, with the following exceptions. The distribution Article https://doi.org/10.1038/s41550-023-01920-2 of w of Europa and Enceladus is dominated by a frequency maximum at w = 1, with a sharp decrease of frequency as w decreases (recall w ≤ 1 is imposed). For both the Jovian and Saturnian systems, the E SC distribution is bimodal with an approximately normal distribution whose maximum is near the best-fit E SC for the measured continua and another frequency maximum at E SC = 2.

Data availability
All data are from the referenced publications or provided equations.

Code availability
Code to reproduce all figures is available via email, without restriction, upon request to the corresponding author.