A comprehensive revisit of select Galileo/NIMS observations of Europa

The Galileo Near Infrared Mapping Spectrometer (NIMS) collected spectra of Europa in the 0.7-5.2 $\mu$m wavelength region, which have been critical to improving our understanding of the surface composition of this moon. However, most of the work done to get constraints on abundances of species like water ice, hydrated sulfuric acid, hydrated salts and oxides have used proxy methods, such as absorption strength of spectral features or fitting a linear mixture of laboratory generated spectra. Such techniques neglect the effect of parameters degenerate with the abundances, such as the average grain-size of particles, or the porosity of the regolith. In this work we revisit three Galileo NIMS spectra, collected from observations of the trailing hemisphere of Europa, and use a Bayesian inference framework, with the Hapke reflectance model, to reassess Europa's surface composition. Our framework has several quantitative improvements relative to prior analyses: (1) simultaneous inclusion of amorphous and crystalline water ice, sulfuric-acid-octahydrate (SAO), CO$_2$, and SO$_2$; (2) physical parameters like regolith porosity and radiation-induced band-center shift; and (3) tools to quantify confidence in the presence of each species included in the model, constrain their parameters, and explore solution degeneracies. We find that SAO strongly dominates the composition in the spectra considered in this study, while both forms of water ice are detected at varying confidence levels. We find no evidence of either CO$_2$ or SO$_2$ in any of the spectra; we further show through a theoretical analysis that it is highly unlikely that these species are detectable in any 1-2.5 $\mu$m Galileo NIMS data.


INTRODUCTION
Europa's young surface and geological features like chaos terrains suggest an active exchange of materials between the surface and subsurface ocean (Carr et al. 1998;Moore et al. 2009;Kattenhorn & Prockter 2014;Trumbo et al. 2019a). Accessing the ocean directly to determine its composition will require a sustained program of scientific and technological developments (Hand & German 2018), and thus our best near-term prospect for determining the ocean composition and constraining Europa's habitability is to understand its surface composition with available spectra. Europa's icy surface is, at present, our primary window into the chemistry of the putative subsurface ocean.
Most of our knowledge of Europa's surface composition comes from reflectance spectroscopy, using ground-and space-based telescopes, as well as instruments onboard spacecraft. The near-infrared reflectance data collected by the Galileo mission, specifically from the Near-Infrared Mapping Spectrometer (Carlson et al. 1992), has revolutionized our understanding of Europa's surface (e.g. McCord et al. 1998;Carlson et al. 1999b;McCord et al. 1999;Hansen & McCord 2004;Carlson et al. 2005;Hansen & McCord 2008;Dalton et al. 2012;Shirley et al. 2016). Spectroscopy in the near-infrared has helped reveal the majority of species identified on Europa's surface (Carlson et al. 2009), including hydrated sulfuric acid (Carlson et al. 1999a(Carlson et al. , 2005, hydrated sulfates (McCord et al. 1998;Dalton 2007), chlorinates (Strazzulla 2011;Ligier et al. 2016) and oxidants (Hansen & McCord 2008;. Although serendipitous observations of Europa by Juno (Filacchione et al. 2019;Mishra et al. 2021) have provided a new window into this world, the Galileo NIMS dataset still holds great potential to further improve our understanding of Europa's surface composition, by the application of new, comprehensive methods, as elucidated below.
A standard approach to finding a species in reflectance data is through a spectral fitting analysis (e.g. Carlson et al. 2005;Clark et al. 2012). Using a radiative-transfer model, such an analysis also constrains a subset of physical properties like abundances and average grain-sizes. For the NIMS data of Europa, much work has been done with models that fit the data with a linear combination of laboratory spectra of various species (e.g. Dalton 2007;Dalton et al. 2012). However, this is a rough approximation to the highly non-linear radiative-transfer process through a planetary regolith (e.g. Shirley et al. 2016). As a result, linear-mixing type analyses do not allow us to disentangle key properties of the regolith such as average grain-size, porosity, observational geometry parameters, and species abundances, which are highly degenerate with each other.
A more physically motivated alternative is a thorough radiative-transfer model that uses fundamental optical properties of each end-member, alongside parameterizations for physical properties of the regolith (e.g. Hapke 2012a). Such a model takes into account various inherent degeneracies in the radiative-transfer process and thus yields more conservative estimates of end-member properties, such as abundances and grain-sizes. Physically motivated radiative transfer schemes, such as the Hapke model (Hapke 2012a), require that optical constants (e.g. refractive indices and extinction coefficients) be measured for several species of interest in the relevant temperature (∼ 100 K) and wavelength regime (∼ 1−5 µm for most NIMS data). Such laboratory measurements can be challenging, which limits the potential Europan species we can explore using this technique. Thankfully, such optical constants exist for certain key species: amorphous and crystalline water ice at 120 K Mastrapa et al. 2009), CO 2 at 179 K (Quirico & Schmitt 1997;Quirico et al. 1999;Schmitt et al. 2018), SO 2 at 125 K (Schmitt et al. 1994(Schmitt et al. , 2018 and sulfuric-acid-octahydrate (hereafter SAO) at 77 K (Carlson et al. 2005).
SAO is a major product of the radiolytic sulfur cycle on the trailing hemisphere of Europa (Carlson et al. 2009). Carlson et al. (2005) used NIMS data and a simple two-component reflectance model of crystalline ice and SAO to map the distribution of SAO on the leading and trailing hemispheres of Europa, finding concentrations as high as 90% near the trailing hemisphere apex. However, an important form of water ice in the upper-most layers of Europa's surface may be amorphous ice (e.g Hansen & McCord 2008), which was not included in the work by Carlson et al. (2005) due to the unavailability of its optical constants at that time. Since the study of Carlson et al. (2005), a comprehensive database of cryogenic water-ice optical constants, in both amorphous and crystalline forms, was published by Mastrapa et al. (2009).
The trace oxides CO 2 and SO 2 are astrobiologically significant oxidants on Europa's surface and might play an important role in creating an environment of chemical disequilibrium in its subsurface ocean (e.g Chyba & Phillips 2001;Hand et al. 2007). Moreover, CO 2 has also been linked to young geological regions on Europa like the chaos terrains (Trumbo et al. 2019b), which indicates a possibly endogenic origin of the carbon that fuels a carbon-cycle on Europa's surface (Carlson et al. 2009;Hand & Carlson 2012). Published estimates of CO 2 abundances exist for the leading side of Europa, where the 4.25 µm feature in NIMS spectra (Hand et al. 2007) was used to constrain the abundance to be ∼ 360 ppm. For SO 2 , a rough abundance estimate of ∼ 0.2% comes from disk-integrated UV observations of Europa (Hendrix et al. 2008;Carlson et al. 2009) -although a 4.0 µm feature has been identified in NIMS spectra (Hansen & McCord 2008). Theoretical studies of CO 2 and SO 2 clathrate formation indicate that the bulk ice-shell of Europa could have an oxidant concentration of up to ∼ 7% (Hand et al. 2006). It should be noted that while sulfur and its allotropes are also present in significant amounts on the trailing side of Europa -due to the Io-genic bombardement (Carlson et al. 2009) -they lack strong features in the NIR wavelength region and hence can be assumed to not affect the NIMS spectra.
To date, only water-ice and SAO have been simultaneously considered in a radiative-transfer based analysis of NIMS data. However, Europa-specific optical constants are now available for amorphous ice, crystalline ice, SAO, CO 2 , and SO 2 . Hence, there is an opportunity to model all five species simultaneously in a revised analysis of Galileo NIMS observations of Europa.
In this work we re-analyze a subset of NIMS spectra using the Hapke reflectance model (Hapke 2012a) in a Bayesian inference framework. Bayesian inference has gained popularity in planetary surface spectroscopic inversion in recent years (e.g. Fernando et al. 2013;Schmidt & Fernando 2015;Fernando et al. 2016;Lapotre et al. 2017;Rampe et al. 2018;Belgacem et al. 2020). The goal of this work is to demonstrate a spectroscopic analysis approach rooted in Bayesian inference (Mishra et al. 2021) that permits several advances over previous analyses of Europan data, specifically: (1) include five species (amorphous ice, crystalline ice, SAO, CO 2 , and SO 2 ) in a radiative-transfer analysis of Europan spectra; (2) include parameters that account for physical effects, such as radiation-induced band shifts, porosity of the regolith, and calibration uncertainty of the instrument; (3) yield statistical constraints on parameters, including confidence intervals; and (4) provide detection significances for each individual species with a statistical metric via Bayesian model comparisons. The fitting analysis in our approach enables identification of trace species, taking into account their effect on the spectrum over the entire wavelength range, instead of relying on a few sharp features.
We apply this framework to three NIMS spectra of the trailing hemisphere of Europa, spanning 1.0-2.5 µm. These data were available in their reduced form, along with all the observation geometry parameters needed for our analysis, through previous work by Carlson et al. (2005). The data and their properties are described in section 2. Section 3 lays out the details of our analysis framework, including the Hapke reflectance model (section 3.1) and the Bayesian inference methodology (section 3.2). We discuss the results from the application of this framework to the three spectra in section 4, followed by a theoretical study (section 5) to gauge the detectability of CO 2 and SO 2 as a function of their abundance in synthetic NIMS spectra over the 1.0-2.5 µm range. We conclude with a discussion (section 6) of our results, placing them in the context of existing knowledge of Europa's surface composition, before summarizing (section 7) the key takeaways from our work.

DATA
During the Galileo spacecraft's first orbit of Jupiter, the NIMS instrument acquired several spectral maps of Europa, the first being a map of the northern hemisphere in the trailing-side anti-Jupiter region and termed the G1ENNHILAT01A observation. Approximately 300 spectra were obtained in this observation and projected as a point-perspective spectral radiance-factor map. Three types of terrain (highly hydrated, icy, and intermediate) were chosen, for illustration of the diversity in composition, and a single pixel spectrum for each region was extracted and fit using experimentallyderived indices of refraction (Carlson et al. 2005; note the interchange of radiance factor and coefficient therein). This work focuses on those three spectra, which are shown in Figure 1, and will henceforth be referred to as S1, S2 and S3. The locations on Europa's surface corresponding to these three spectra are shown in Figure 2, and details about the spacecraft observation geometry is shown in Table 1. Carlson et al. (2005) found that S1, S2 and S3 represent a range of compositions of Europa's trailing hemisphere, from highly hydrated to minimally hydrated. Hence, they also allow us to explore trends in the physical and chemical parameters included in our model across regions of differing water-ice content. Table 1. NIMS spectra from the G1ENNHILAT01 observation of Europa being analyzed in this work and referred to as S1, S2 and S3. These spectra are shown in Figure 1.  (Carlson et al. 2005). Latitude 0 • and Longitude 0 • correspond to the sub-Jovian point on Europa's surface, with longitude increasing positively westward towards the leading hemisphere. S1, S2 and S3 are all from the trailing hemisphere of Europa.
a The original data for this orbit was written to a CD and was submitted to the PDS as Volume 'go 1104' and file 'g1e002ci.qub'. The PDS filename is now 'g1e003' followed by 'cr' or 'ci', with 'cr' and 'ci' denoting radiance and radiance factor data respectively

S1
S2 S3 Figure 1. The three Galileo NIMS observations, referred to as S1, S2 and S3, being used in this work. The spectra come from observation G1ENNHILAT01 by Galileo that mapped the trailing hemisphere of Europa. Details for the three spectra are provided in Table 1 and the their locations on Europa's surface is shown in Figure 2. S1, S2 and S3 also correspond to the bottom, middle and top spectra, respectively, shown in Figure 3 of Carlson et al. (2005).
The compositional discriminative power of a spectral analysis method like Bayesian inference depends sensitively on the uncertainties in the data (see section 3 of Mishra et al. (2021)). For most NIMS spectra, including S1, S2 and S3, accurate uncertainties could not be calculated because Galileo's satellite observations were limited in temporal duration by the relative spacecraft-satellite velocity and by the solar illumination incidence angles. In addition, data transmission to the Earth was severely limited, consequently repeated measurements of a particular region on Europa's surface were precluded. It is the intrinsic variation among the repeated measurements that informs the uncertainties on the data. Estimates of the uncertainties on data of the kind considered here are usually obtained from repeated measurements of a particular region of interest. Because the spacecraft's flybys when the data was collected were very short in duration, such repeated measurements could not be taken.
While rough estimates on the SNR (signal-to-noise ratio) of Galileo data in the literature put the number between 5 and 50 ), the noise characteristics vary with wavelength and depend on instrumental and external factors that change from observation to observation (Carlson et al. 1992). The noise in the NIMS instrument includes thermal noise (which comes from the blackbody IR emission of the interior instrument housing), detector current noise that is approximately proportional to the inverse of the solar spectral irradiance, a small amount of digitization noise, and radiation noise that produces pulses of detector current with widely ranging amplitudes. The obser-  Figure 2. The approximate locations on Europa's surface, marked as red, blue and green stars, corresponding to spectra S1 (7.5 • N, 261.4 • W), S2 (19.9 • N, 225.6 • W) and S3 (65.5 • N, 254.6 • W) respectively (see Table 1 for details). This map is centered on 180 • W, where the 0 • W corresponds to the sub-Jovian longitude. The global base map of Europa was downloaded from USGS 'Map a Planet' a , derived from Galileo and Voyager images, and is shown in cylindrical projection. a https://astrocloud.wr.usgs.gov/index.php?view=map2 vations we are using did not suffer excessive radiation noise, but a small but non-quantifiable amount.
We're assuming we're working in the photon-noise limited regime, which allows us to measure all the sources of "error" via the scatter they induce on the signal. We obtain estimates of this inherent scatter in each of S1, S2 and S3, caused by a combination of all the factors listed above, by adopting the following process: 1. Assume a SNR of 50 across all wavelength channels.
2. Perform a fitting analysis with a full model (all five species included) and evaluate the best-fit or the maximum-likelihood solution. The details of the analysis framework are discussed in the next section (section 3).
3. Calculate the residuals of this fit to the data.
4. Calculate the standard deviation of these residuals or the SDR.
5. Repeat the steps for SNR of 5.
The larger of the two SDR values, from the SNR=5 and SNR=50 cases, is chosen to be a conservative estimate of the noise in a given spectrum and is assigned as the uncertainty across all channels. This analysis results in data uncertainties of 0.0081, 0.0063, and 0.0047 (in radiance factor units), which correspond to an average SNR (averaged across the wavelength channels) of 39.51, 44.43 and 28.35 Wavelength in m top-right and bottom-right panels show the k spectra convolved to the NIMS resolution. CO 2 and SO 2 k spectra have been individually shown in the bottom panels for clarity. Amorphous and crystalline ice optical constants, at 120 K, were published and included in Mastrapa et al. (2009); SAO optical constants at 77 K were measured by Carlson et al. (2005); CO 2 optical constants, at 179 K, are from Quirico & Schmitt (1997); Quirico et al. (1999); SO 2 optical constants, at 125 K, are from Schmitt et al. (1994Schmitt et al. ( , 1998 for S1, S2 and S3 respectively. These values fall in the range of 5 to 50 that is generally assumed to be the average SNR of NIMS spectra of Europa ). Finally, an important dataset needed for our model are the optical constants for the species considered, which are shown in Figure 3 (sources are mentioned in the figure caption). We use amorphous and crystalline water-ice optical constants at 120 K, similar to Carlson et al. (2005). For SAO, CO 2 and SO 2 , we have selected available optical constants that were derived at temperatures close to the Europan temperature range of 80-130 K (Spencer et al. 1999). We note that optical constants for the species considered here are temperature dependent, thus both the measured and theoretical reflectance spectra employed in this study will be sensitive to temperature as well (Dalton 2010). Future work that expands on either the species or wavelengths consider here would be best served by additional laboratory-derived optical constants obtained at Europan temperatures.

METHODS
We use a Bayesian inference framework to analyze S1, S2 and S3, which has been described in detail in Mishra et al. (2021). Briefly, the two major components of this framework are the Forward model and the Bayesian posterior sampling, described in the next two sections. The Bayesian posterior sampling efficiently samples the parameter space we wish to explore, generating millions of instances of Forward model spectra, and selects the samples for which the corresponding model instances fit the data well. From this collection of samples, one can get marginalized posterior distributions of individual parameters, pair-wise 2D distributions of parameters that highlight correlations, and a marginal probability of the model itself, also known as the Bayesian evidence, that is used for model comparison (e.g. Trotta 2017).

Forward model
To a good approximation, the bidirectional reflectance of a planetary regolith (Hapke 2012a) can be written as: Here • r F is the radiance factor, which is the ratio of bidirectional reflectance of a surface to that of a perfectly diffuse surface (Lambertian) illuminated and observed at an incidence angle, i = 0, relative to the surface normal.
• µ is the cosine of the emergence angle e.
• µ 0 is the cosine of the incidence angle i.
• g is the phase angle.
• K is the porosity coefficient.
• ω is the single scattering albedo.
• P is the particle phase function.
• H is the Ambartsumian-Chandrasekhar function that accounts for multiple scattered component of the reflection (Chandrasekhar 1960).
Here, we have ignored parameters related to backscattering and other opposition effects as they become important below small phase angles for the coherent backstatter effect and typically less than a few tens of degrees for the shadow-hiding opposition effect (Hapke 2012b(Hapke , 2021, whereas the phase angles of all three NIMS spectra used here are much higher (∼ 30 • , see Table 1). In addition, we ignore the photometric contribution of macroscopic roughness, which Hapke characterizes with a mean topographic slope angle, θ. Our available range of photometric geometries for S1, S2, and S3 is too limited to constrain θ, which, even at a single phase angle, requires sufficient coverage of incidence and emission angles to characterize the limb-darkening behavior of the surface. Moreover, global average values of θ derived from Voyager and Galileo imaging data vary little from 16 • to 22 • among different works (Verbiscer et al. 2013). At the photometric geometry of the NIMS spectra for S1, S2, and S3, the effects should be small enough to be ignored. However, we note that the detectability of photometric roughness is albedo-dependent such that the measurable values of θ diminish as the surface albedo approaches unity because of multiply-reflected light among topographic surface facets. Since the particle spectral albedos vary with wavelength, it is reasonable to expect that there may be a wavelength-dependent variation in the effects of macroscopic roughness on our spectra. None of our reflectance values approach close enough to unity for this to be a significant concern.
The equations for the single scattering albedo, ω, come from the equivalent-slab approximation (Hapke 2012c). For the phase function P , we use a two-parameter Henyey-Greenstein function (Henyey & Greenstein 1941), which has been shown to be representative of a wide variety of planetary regoliths, including the icy particles on Europa (Hapke 2012c). A detailed description of all the parameters in eq. 1 can be found in the Appendix of Mishra et al. (2021).
We are also assuming that the regolith in the heavily bombarded and sputtered trailing hemisphere of Europa, where S1, S2 and S3 come from, is well-mixed or 'churned' (Carlson et al. 2009;Shirley et al. 2016). Hence, particles of different species are mixed homogenously together in a intimate mixture. The averaging process in this intimate mixing model is over the individual particle, and ω and P in Eq. 1 become volume averages of different materials in the mixture: where for a species of type j, N j is the number of particles per unit volume, σ j (= πD 2 /4) is the cross-sectional area, Q E,j is the volume average extinction efficiency, ω j is the the single scattering albedo, P j is the phase function, and f j is the number density fraction, equivalent to N j / N j . For more details about the mixing equation parameters, please refer to section 2.2.2 of Mishra et al. (2021).

Bayesian posterior sampling
The parameters we are fitting for, which define the parameter space are: • log 10 (X i ), where X i is the abundance (or number density fraction) of each species • log 10 (D i ), where D i is the grain-size of each species • filling factor φ, which is the volume fraction occupied by particles, equivalent to 1 -(porosity of the regolith). φ is related to the porosity coefficient K in eq. 1 through the relation • a multiplicative factor c for the model spectrum. In the 1-2.5 µm region, the absolute calibration of the NIMS instrument, or the accuracy of its measurements, is uncertain to around ±10% (Carlson et al. 1992(Carlson et al. , 2005. • a wavelength-shift parameter ∆ wav , to shift the absorption band centers of the model spectrum by ∆ wav , mimicking radiation induced shifts. Wavelength in m Figure 4. Best-fit or maximum likelihood solutions of a SAO and amorphous ice model fit to S1 (blue points). The left panel shows the solution (orange) for the case where a wavelength-shift parameter is not included, while the right panel shows the solution (orange) for the case where that parameter is included (this turns out to be the preferred model for S1 as shown in section 4). In the left panel, the model is slightly shifted to longer wavelengths as compared to the data, around the band-edges near 1.4 and 1.9 µm. The reduced χ 2 values of the fits are also noted, which illustrate that the fit in the right panel is better.
We have included the ∆ wav parameter because Europa's trailing hemisphere spectra, like S1, S2 and S3, have water absorption-band centers that are shifted to slightly shorter wavelengths, by about tens of nanometers. This is a known effect of particle irradiation on fundamental absorption bands of hydrates, sulfates and silicates (e.g. Dybwad 1971;Nash & Fanale 1977). For Europa specifically, it has been emphasized that comparisons of lab and numerically modelled spectra to Europa spectra must allow for the radiation-induced shifts of hydrate bands, so as to not bias our inference of the chemical composition from the spectra (Carlson et al. 2009). One way to circumvent the problem would be to discard the data points near band-edges (e.g Carlson et al. 2005). However, we find that our wavelength-shift parameter ∆ wav works sufficiently well and results in excellent fit to the data. An example is shown in Figure 4, where S1 is fit with a model with and without the ∆ wav parameter, with a much superior fit in the latter case. Moreover, measuring this wavelength shift for different spectra can also inform us about trends in the level of irradiation experienced by the regoliths corresponding to the spectra. Hence, we include this parameter in all subsequent analysis.
To sample this multi-dimensional parameter space, we use the nested sampling algorithm (Skilling 2006). Specifically, we employ the Python implementation called dynesty 1 , which implements dynamic nested sampling (Higson et al. 2019), a more computationally accurate version of the standard nested sampling algorithm. Schematically, the iterative parameter exploration by nested sampling proceeds as follows: 1. A random instance of a set of parameters or a parameter vector is drawn from the parameter space, defined by the prior distribution function. For abundances, we use a Dirichlet prior (e.g. Benneke & Seager 2012;Lapotre et al. 2017), which ensures that the abundance parameters being sampled satisfy the constraint that they should sum up to 1. For the rest of the parameters, which do not have any such constraints, we use a uniform distribution, which is relatively uninformative (as compared to a Gaussian prior, for example), leading the data to drive the solution. The prior function also serves the purpose of defining the boundary of the parameter space to explore. Table 2 lists all the free parameters, their prior function type and bounds 2 .
2. Using these parameters, with the optical constants of the component species and the observation geometry parameters as the main inputs, the forward model (section 3.1) generates a model spectrum and then convolves it with the instrument response function (or bins it to the instrument's resolution) to produce simulated data points. Each NIMS channel has a triangular response function (Carlson et al. 1992(Carlson et al. , 2005, spanning the channel-width of ∼ 0.026 µm.
3. These predicted data points are compared with observed data points to compute the posterior probability of this particular instance of parameters, defined by the Bayes' theorem as follows: where L(d|θ, M i ) is the likelihood probability function, π is the prior, Z is the Bayesian evidence and p(θ|d, M ) is the posterior probability function. Since we assume the errors on our data to be Gaussian and independent, the likelihood function is also a Gaussian equal to where N obs is the number of observed data points (i.e., number of wavelength channels/data points in the observed spectrum), χ 2 is the familiar goodness-of-fit metric, and σ is the standard deviation.
4. The posterior probability decides whether this particular parameter set is saved or selected and, in turn, informs the choice of the next parameter set as well. Details about this process/algorithm can be found in Higson et al. (2019). Note-φ's lower limit has been set to 0.01 because K, which is the model parameter dependent on φ (see eq. 1), plateaus at a value of 1 for φ 0.01. The upper limit of 0.52 comes from Hapke (2008) where it is described as a critical value above which the medium is too tightly packed, such that coherent effects become important and diffraction can't be ignored.
After this process converges, the final set of samples can be marginalized to get individual or pairwise parameter distributions. Another useful quantity that the nested sampling process returns is the Bayesian evidence of the model or Z (eq. 4), which is discussed in the next section. It quantifies the probability of a model given the data, integrated over the entire parameter space (eq. 4). Hence, two models can be readily compared using Z, in a process called Bayesian model comparison. We adapt this into a nested Bayesian model comparison process, where we compare the evidence for a model with all species (the full model) with models where each species is removed one at a time. Assuming no a priori preference of either the full model or the model without species X, if the Bayesian evidence is higher for the latter, we conclude that X is not needed to explain the data and hence there is no evidence for it. This process helps us select the simplest model, i.e., with the fewest number of species, compatible with the observations. The comparisons are also used to quantify the confidence in the presence of each species remaining in the simplest model, through the σ-significance metric (see section 2.2.3 of Mishra et al. (2021) for a full description of statistical aspects of the Bayesian inference methodology). As per empirically calibrated evidence thresholds known as Jeffrey's scale (Jeffreys 1939), a detection significance exceeding around 2.1σ, 2.7σ and 3.6σ levels are referred to as 'weak', 'moderate' and 'strong' evidence respectively. We use these metrics to guide our interpretation of our model results.
Over the course of analyzing S1, S2 and S3, we performed around 25 Bayesian inference analyses with models containing different combinations of the five species -amorphous ice, crystalline ice, SAO, CO 2 and SO 2 . A full-model analysis (with all the five species) amounts to 12 free parameters and computes ∼ 4 x 10 5 model spectra. Using this methodology, for each of the three spectra, we now proceed to describe the detection significance of all five species, the simplest model that best explains the spectrum, and the constraints on the parameters of that model.

RESULTS
The results for the nested model comparisons performed for S1, S2 and S3 are shown in Table 3. For each spectrum, we first perform a Bayesian inference analysis with a full model, i.e., with all five species -amorphous ice, crystalline ice, SAO, CO 2 and SO 2 . We then remove each species one at a time, and rerun the analysis with these alternative models. Comparing each of their results against the reference model with all five species, we can quantify the evidence for the presence of each species. The retrieved parameter constraints for the models with highest evidence (i.e., the simplest models that explain the data) are shown in Table 4. Note-For each of the three spectra, S1, S2 and S3, the Bayesian evidence ln(Z i ) and the minimum reduced chi-squared χ 2 r,min , i.e, chi-squared of the maximum-likelihood solution are listed for all models. The 'All species' model is the reference model to compare to and contains amorphous water-ice, crystalline waterice, SAO, CO 2 and SO 2 . The other four alternative models are each without one of the five species as specified. In brackets next to the ln(Z i ) values, an n-σ detection indicates the degree of preference for the reference model over the alternative model, which can also be interpreted as the significance of detection of the species excluded from the alternative model. 'N.A.' indicates that no evidence for the particular species has been found.
For S1, we find that there is a very strong evidence for SAO (at > 30σ confidence level). This is expected as the octahydrate, which is taken to generically represent hydrates in our model, is the dominant component in the trailing hemisphere of Europa, as established in previous studies (e.g. Carlson et al. 2005;Ligier et al. 2016). Weak evidence for amorphous ice at a 2.4σ confidence-level is also found. We don't find any evidence for crystalline ice, CO 2 or SO 2 , which means that the simplest model that explains S1 consists of only SAO and amorphous ice. Indeed, a nested model comparison with the SAO and amorphous ice mixture as the reference model (not shown in Table 3) provides very strong evidence for both: > 30σ for SAO and 8.6σ for amorphous ice. Detailed results from the analysis with the SAO and amorphous ice mixture model are presented in Figure 5. The median retrieved spectrum provides an excellent fit to S1, with a reduced χ 2 value of 1.02. Also shown in the figure are the posterior probability distributions of the model parameters. SAO as 1 − φ). The wavelength shift parameter ∆ wav has also been well constrained to -0.0198 +0.0015 −0.0016 , slightly less than the NIMS channel width of 0.026 µm. This means that the water band centers at around 1.5 and 2.0µm in the model spectrum needed to be shifted Table 4. Derived parameter values, along with the 1-σ upper and lower bounds (68% confidence intervals), for models with highest evidence for S1, S2 and S3. The parameters are described in towards shorter wavelengths by around 0.02 µm or 20 nm to account for the radiation-induced shift seen in the data. Next, for S2, we once again find very strong evidence (> 30σ) for SAO. While no evidence for amorphous ice is found, S2 shows a strong (5.9σ) evidence for crystalline ice. No evidence for either CO 2 or SO 2 is found, which means that the simplest model that fits the data comprises of a mixture of SAO and crystalline ice. Nested model comparisons with SAO and crystalline ice mixture model as the reference provides very high evidence for both SAO (> 30σ) and crystalline ice (29.7σ). The results of the analysis with this preferred model are shown in Figure 6. The median retrieved spectrum provides an excellent fit to the data, with a reduced χ 2 value of 1.04. All the model parameters have well constrained Gaussian distributions. Crystalline ice is present in a very minor amount (0.24 +0.0012 −0.0007 %), but with very large grains of average size 382.44 +78.19 −70.66 microns. SAO heavily dominates the composition with an abundance of 99.76 +0.0007 −0.0012 %, but smaller grains of average size of 44.98 +3.23 −3.12 microns. The filling factor φ is low (0.18 +0.06 −0.06 ), indicating a porous regolith. Finally, a slight band-center shift is once again detected and the wavelength-shift parameter ∆ wav is constrained to -0.0098 +0.001 −0.001 , which is about half the channel-width of NIMS. Finally, for S3, we find very strong evidence for SAO (> 30σ) and crystalline ice (4.76σ), whereas a moderate to strong evidence for amorphous ice (3.36σ). Once again, no evidence is found for either CO 2 or SO 2 . The preferred model is thus a SAO, amorphous ice and crystalline ice mixture, whose results are shown in Figure 7. Nested comparisons for this model provides strong detectionlevels for all three components: > 30σ for SAO, 3.8σ for amorphous ice, and 7.2 σ for crystalline ice. The median retrieved spectrum, as shown in Figure 7, provides a good fit to the data, with a reduced χ 2 value of 0.79. The retrieved parameter distributions tell us that SAO is once again the dominant component (71.15 +13.82 −16.38 %), albeit with an abundance slightly lower as compared to the  Table 2. The dashed purple lines are the median values along with the 1-σ upper and lower bounds (68% confidence intervals). All the values are listed in Table 4. Note that the distribution for f am has not been shown here as it is dependent on f sao and not one of the free parameters. The derived values for f am are shown in Table 4.
values for S1 and S2. Amorphous ice is also a major species, with an abundance of 28 microns respectively. The regolith is constrained to be porous once again, with a significantly low value for φ of (0.03 +0.03 −0.01 ). Finally, the spectral band-centers are negligibly shifted, as the ∆ wav parameter is retrieved to be -0.0043 +0.0015 −0.0015 µm. Along with the individual parameter distributions that are shown in Figures 5-7, our analysis also yields pair-wise distributions of the various parameters (Figures 10-12), which highlight the correlations between them. Some correlations are expected, such as the correlation between grainsize and abundances. In all three cases, we see that the grain-size of any species (SAO, amorphous ice or crystalline ice) is anti-correlated with its abundance. This makes sense as the 'weight' of a species in the mixing equation in Hapke's model is proportional to the product of its abundance and its grain-size squared. Hence, decreasing the abundance of a species should have the same overall effect as increasing its grain-size.  On the other hand, Figures 10-12 also illuminate complex parameter correlations, e.g., the correlation between the filling factor φ and the calibration parameter c is puzzling -we see that φ and c show no correlation for S1, positive correlation for S2 and negative correlation for S3. The filling factor φ enters the Hapke radiative transfer equation (eq. 1) as the porosity coefficient K, via the relation K = −ln(1 − 1.209φ 2/3 )/(1.209φ 2/3 ). K affects eq. 1 as a multiplicative factor and via the multiple scattering function H. When φ and c are allowed to vary in isolation (the abundance and grain-size parameters are kept fixed), they are highly anti-correlated, since both essentially act like multiplicative factors. However, in our exercise where abundance and grain-size parameters are free as well, φ is not the only parameter affecting the multiple scattering function H, and hence the correlations between φ and c become complex.

ESTABLISHING DETECTABILITY LIMITS FOR CO 2 AND SO 2
Since neither CO 2 nor SO 2 were detected in any of the three spectra, it was of interest to perform a theoretical study that establishes the minimum amount of CO 2 and SO 2 that can be detected with NIMS data, in the 1-2.5 µm range. Figure 3 shows that in the 1-2.5 µm range, both CO 2 and SO 2 possess sharp features in their k (imaginary part of refractive index) spectra. CO 2 especially has numerous features, with a few at around 2.0 microns matching the strength of water-ice features. Since k is a proxy for the absorption coefficient (α = 4πk/λ), the sharp/broad nature and position of peaks in the k spectrum also determine the absorption features in the reflectance spectrum.
We investigate the effect of CO 2 and SO 2 on the net reflectance in a mixture with the dominant species SAO, amorphous ice and crystalline ice also present. An example is shown in Figure 8, where theoretical spectra of a mixture of SAO, amorphous ice, crystalline ice and CO 2 /SO 2 were generated. We took the solution of the SAO, amorphous ice and crystalline ice mixture from the S3 analysis (Figure 7) and injected a fourth species -either CO 2 or SO 2 -in increasing amounts (SAO's abundance was being reduced to ensure the abundance fractions add up to 1). For CO 2 , only at an unrealistically high abundance of around 60% does its absorption features near 2.0 µm start to appear in the high-resolution spectrum (resolution ∼ 0.001µm). These features do not appear in the coarser NIMS resolution spectrum (resolution ∼ 0.026µm). For SO 2 , no sharp features appear at native resolution even at abundance values approaching 60%. Figure 8 also illustrates that both CO 2 and SO 2 affect the overall continuum level of the spectra as well. However, this effect is degenerate with parameters like the filling factor φ, which also changes the overall amplitude of the spectrum.
A Bayesian model comparison analysis shows that significant evidence (> 3σ) for CO 2 and SO 2 only starts emerging at abundances of 60% or higher. Hence, given the low expected abundance of around a few percent for CO 2 and SO 2 on Europa's surface (Hand et al. 2007;Carlson et al. 2009), and the relatively coarse spectral resolution of NIMS, detection of CO 2 , SO 2 and possibly other oxides like H 2 O 2 in the 1-2.5 µm wavelength range Galileo/NIMS spectra is unlikely.

SUMMARY AND DISCUSSION
We applied a novel Bayesian inference framework to three Galileo NIMS observations of Europa (Figure 1), referred to as S1, S2 and S3. Our analysis framework successfully evaluated evidence for and provided constraints on abundances and average grain-sizes of important trailing hemisphere species -amorphous water ice, crystalline ice, sulfuric-acid-octahydrate (SAO), CO 2 , and SO 2 . These parameter constraints account for the degeneracies with physical parameters of the regolith, like porosity and radiation-induced spectral band-center shift. We also employ Bayesian model compari- son to find the simplest model that explains each spectrum, for which we also evaluate the statistical significance of detection of each species. The model comparisons (Table 3) show that S1 is best explained by a mixture of SAO and amorphous ice, S2 by a mixture of SAO and crystalline ice, and S3 with a mixture of SAO, amorphous ice, and crystalline ice. The retrieved species abundances and grain-sizes for these preferred models are shown in Table 4. No evidence for either CO 2 or SO 2 is found for any of the three spectra, but given the limitations of the Galileo NIMS data we cannot set useful limits on the abundances of these two species. Oxides, in trace amounts, can also be trapped in the matrix of the major species. Hence, the use of optical constants for CO 2 or SO 2 trapped in water-ice, rather than those of pure CO 2 or SO 2 ice, may be more appropriate. However, because of the minimal impact of CO 2 and SO 2 in the current dataset (see section 5 and Figure 8), the choice of optical constants for the two oxides should have minimal impact on our main conclusions. SAO dominates the composition (> 30σ detection) for all three spectra, with S3 having the lowest amount of SAO among the three. This is consistent with S3's location at a high latitude of 67.9 • , where the abundance of SAO is expected to be lower than nearer the equator (Carlson et al. 2005;Ligier et al. 2016;. Carlson et al. (2005) analyzed the same spectra and derived SAO's abundance to be 89% for S1, 54% for S2 and 30% for S3. While the abundance we retrieve for S1, 87.11 +7.27 −16.98 %, is consistent with their result to within 1-σ, our numbers for S2 (99.76 +0.0007 −0.0012 %) and S3 (71.15 +13.82 −16.38 %) do not agree with their estimates. Our average-grain-size estimates for SAO are in better agreement with their results. Carlson et al. (2005) used log-normal distributions for the icy and acidic particle sizes and reported the radii, r m , corresponding to the mean values of ln(r). The effective scattering radius, r ef f , is the cross-section weighted radius (Hansen & Travis 1974) and for the log-normal case is r ef f = r m exp[(5/2)σ 2 ] ≈ 3.0r m , using the variance σ 2 = 0.44 found for Antarctic ice grains (Aitchison & Brown 1957). The corresponding SAO diameters they obtain for the S1, S2 and S3 cases are 37, 42, and 57 microns, respectively, and generally agree with the values we obtain (S1 : 42.21 +5.96 −4.61 microns, S2 : 44.98 +3.23 −3.12 microns, and S3 : 97.83 +11.87 −10.75 microns). There can be many reasons for the discrepancy between our abundance estimates for SAO and those of Carlson et al. (2005). Firstly, Carlson et al. (2005) used a model with only two components: sulfuric-acid-octahydrate and crystalline water ice. Amorphous water ice was not included due to the unavailability of its optical constants when they performed the work. Secondly, the reflectance model used in that work assumed spherical particles, while the radiative-transfer model we are using (see section 3.1) is better suited for irregular particles (Hapke 2012c) and is more likely to apply to airless planetary surfaces like Europa's. Our best direct evidence for irregular particles on airless bodies comes from microscopic studies of lunar regolith particles (e.g Slyuta 2014). As for Europa, detailed surface photometry and Hapke modeling show that the average single particle phase function of Europan regolith particles are irregular in shape (Domingue & Hapke 1992;Domingue & Verbiscer 1997;Hapke 2012c;Verbiscer et al. 2013). Finally, to navigate around the problem of radiationinduced shift of the absorption band centers, which could have biased their results, Carlson et al. (2005) discarded the band-edge data near 1.4 and 1.9µm. They also discarded data below 1.3 µm for all three spectra, to ensure consistent comparisons with other datasets used in their work that did not go below 1.3 µm. We deal with the band-center shift problem by introducing a parameter that can perform that shift in wavelength, and hence are able to use all the data points in each spectrum.
The strong detection of amorphous water ice in S3 (3.36σ level) fits with our expectation. The prevalence of amorphous ice on the surface of Europa is well established, especially in regions of high particle irradiation (e.g. Strazzulla et al. 1992) and/or colder temperatures (e.g. Hansen & McCord 2004). S3 comes from the colder plain-terrains region in the high latitudes of the trailing hemisphere of Europa (see Figure 2). While particle irradiation is weaker at higher latitudes as compared to the equatorial latitudes (Nordheim et al. 2018), dominance of water at higher altitudes (Fanale et al. 1999;Hansen & McCord 2004;Carlson et al. 2005) along with colder temperatures might enhance amorphous ice abundance. S1 shows weak evidence of amorphous ice (2.4σ level) and comes from a region quite close to the equator. Perhaps the higher average temperatures in these regions -that cause the annealing of amorphous ice to crystalline ice -are dominating over the opposing effect of irradiation induced amorphization, resulting in only a minor presence of amorphous ice. S2 on the other hand, which also comes from equatorial latitudes, doesn't show any evidence for amorphous ice. This might be explained by the fact that S2 belongs more to the anti-Jovian hemisphere, where particle irradiation is weak beyond ∼ 20 • latitudes, as compared to near the trailing hemisphere apex (Nordheim et al. 2018). Moreover, S2 also comes from a region containing numerous linea (image not shown here), implying high tidal stresses and perhaps heating from differential motions (Doggett et al. 2009), which would favour crystallization and work against irradiation induced amorphization. This potential form of enhanced heating could also explain why we see a very high SAO abundance of ∼ 99% for S2, as heating of the surface on Europa has also been linked to the formation of sublimational lag deposits (Fagents et al. 2000;Carlson et al. 2005). Water has a higher vapor pressure as compared to sulphur allotropes and sulfuric acid, and even modest heating can saturate the surface with sulfurous material.
For crystalline ice, we obtain strong evidence for its presence in S2 (5.9σ) and S3 (4.76σ), consistent with the 1.65 µm absorption feature corresponding to crystalline ice being clearly visible in both the spectra (Figure 1). Evidently, at the locations of S2 and S3, water-ice has not been amorphized to the depths being probed by reflectance spectroscopy in the 1-2.5 µm wavelength region. S1 on the other hand comes from the visibly darkest regions of the trailing hemisphere (Figure 2), quite close to the apex, where amorphization probably extends to the full sub-mm depths being probed by NIMS. This also relates well to the fact that only S3 shows strong evidence (> 3σ) for both forms of water ice, probably because of its location in the high-latitude regions where the water-ice abundance is significantly higher than the equatorial latitudes in the trailing hemisphere. Hence, we are seeing both amorphous ice, formed due to irradiation and/or colder high-latitude temperatures, in the top-most layers of the regolith, and crystalline ice that is likely present just below it. Interestingly, although we retrieve very low abundances of crystalline ice for both S2 (0.24 +0.0012 −0.0007 %) and S3 (0.36 +0.0009 −0.0006 %), its effect on the spectra mainly comes from the very large average grain-size values of 382.44 +78.19 −70.66 microns for S2 and 919.14 +58.26 −96.88 microns for S3. These large grain-sizes are consistent with the range of 10s-100s of microns that have been found in previous spectroscopic and polarimetric studies (e.g. Dalton et al. 2012;Ligier et al. 2016;Poch et al. 2018). It should be noted that it is more appropriate to think of these large grain-sizes as a collection of crystallites or grains, rather than one big crystal. While the initial production of hydrate would be uniform, sublimation of the more volatile H 2 O would produce localized regions of water-rich and hydrate-rich material. The scale of these concentrations depends on photon path-lengths inside and outside the regolith, which is not easy to estimate. However, one can assume that the transport inside this regolith dominates, and hence the sizes of these regions can be parameterized as a pseudo-particle-size.
The wavelength-shift parameter, ∆ wav , shows a clear decreasing trend from S1 to S3. The shift is almost a channel-width: -0.0198 +0.0015 −0.0016 µm for S1, about half of that at -0.0098 +0.001 −0.001 µm for S2, and negligible at -0.0043 +0.0015 −0.0015 for S3 (the negative sign indicates a shift towards shorter wavelengths). As mentioned earlier, irradiation-induced wavelength shifts, especially of water absorption bands, are a well-studied phenomenon in the laboratory (e.g. Nash & Fanale 1977;Poston et al. 2017;Brunetto et al. 2020). They are primarily caused by changes in the cubic structure of water ice (Baratta et al. 1991). If the shifts we infer are indeed correlated with the amount of irradiation experienced, the location corresponding to S1 has perhaps experienced the most irradiation among the three locations considered here. It should be noted that our simplified ∆ wav parameter shifts the whole spectrum, rather than just the regions near the water absorption bands, the latter being much more complicated to implement. However, as Figure 4 shows, the shift mostly affects the fit near the water bands. The 1-2.5 µm wavelength range covered by our data has very little continuum and, although the fit near the continuum gets a little worse when the model is shifted, the overall impact is minimal. For data with a broader wavelength coverage, a more complicated shift treatment may be warranted.
Finally, for the filling factor parameter φ, we retrieve a value of 0.15 +0.07 −0.06 for S1, 0.18 +0.06 −0.06 for S2 and 0.03 +0.03 −0.01 for S3. These values indicate regoliths of very high porosity, (1 -φ), specifically 0.85 +0.06 −0.07 for S1, 0.81 +0.06 −0.06 for S2 and 0.97 +0.016 −0.032 for S3, respectively. For comparison to familiar terrestrial snow examples, fresh wind-blown snow has measured porosity values in the range of 83% to 85% (Clifton et al. 2007), and the topmost layer of frigid, dry snow deposits can reach almost 91% porosity (Fu et al. 2018). Despite the complex correlation between φ and other parameters discussed earlier, particularly with the calibration parameter c, the retrieved distributions for φ are well constrained in all three cases (Figures 5-6). Reported values for Europa's porosity in the literature are limited and have never employed spectroscopic analysis. Our relatively high values of porosity are generally consistent with measurements of the low thermal-interia of airless icy bodies as well as Europa's surface in particular (Ferrari & Lucas 2016;Rathbun & Spencer 2020;Rathbun et al. 2010). Nelson et al. (2018) show that laboratory polarization studies of airless solar system body surface analogues indicate that the surface of Europa can have a porosity as high as 90%. On the other hand, Poch et al. (2018) performed polarimetric studies of ice particles, with their analysis showing that Europa most likely has sintered grains, leading to a more compact or less porous regolith. However, they caution that, unlike their pure water-ice lab samples, Europa's surface contains a significant amount of dark non-icy materials, which can significantly affect the degree of polarization. Moreover, they did not directly study the effect of porosity on polarimetric phase curves. Fits of the most recent Hapke (2012a) model to combined Earth-based telescopic and spacecraft observations of Europa's visual photometric phase curve simultaneously constrain the porosity from the angular-width of the shadow-hiding opposition effect as well as from the way increased compaction states can amplify the reflectance overall. The global-average fit reported in Verbiscer et al. (2013) constrains Europa's porosity at around 0.61±0.08, 3 somewhat less porous than what we have obtained in our individual regions. However, our own results suggest that the porosity may vary significantly over different terrain-types and regions of Europa's surface.
Since reflectance spectroscopy at the Galileo NIMS wavelengths is probing just the upper 100s of microns-millimeter of the regolith, it is reasonable to think that the surface we are seeing is highly porous or 'fluffy' in a thin layer at the top, but with compactness varying with depth and perhaps increasing. Another explanation of the high porosities we retrieve might be poorly calibrated data. As Figure 9 shows, our model can only match the low reflectance level of S3 with very low values of φ (and extremely high values of the crystalline ice grain-size). Finally, deficiencies in the Hapke model itself might be forcing the retrieved porosity to become extreme. Helfenstein & Shepard (2011) tested Hapke's photometric model's ability to constrain porosity of lab samples through their phase curve data, finding preliminary evidence that the effects of porosity may be difficult to detect in high-albedo surfaces. However, laboratory work is required to asses how well the K parameter constrains the porosity of an icy regolith using spectroscopy.

CONCLUSIONS
The highly non-linear radiative-transfer process through a planetary regolith calls for a radiativetransfer model, such as the Hapke model, to lie at the heart of spectroscopic analyses of planetary reflectance data. Here, we took the most up-to-date version of the Hapke model and make several quantitative improvements, relative to prior implementations, for Europan data analysis: (i) the simultaneous inclusion of amorphous and crystalline water ice, sulfuric-acid-octahydrate (SAO), CO 2 , and SO 2 ; (ii) the use of physical parameters like regolith porosity; and (iii) consideration of a radiation-induced band-center shift and NIMS calibration-uncertainty. To thoroughly explore the complex multidimensional parameter space of the Hapke model, we employed a Bayesian inference framework. Most importantly, the Bayesian approach folds in both data uncertainties and model degeneracies into all parameter constraints, providing conservative estimates on key parameters of interested, such as species abundances. The key findings from the application of our framework to the three Galileo NIMS spectra S1, S2 and S3 (Figures 1 & 2) are as follows: • Sulfuric-acid-octahydrate or SAO is strongly detected (at > 30σ confidence level) and dominates the composition in all three spectra (Table 3). This matches our expectations, as all three spectra come from the trailing hemisphere of Europa where sulfur-based hydrates are the dominant product of the surface radiolysis reactions. Among the three spectra, S3 has the lowest amount of SAO (Table 4), which agrees with its physical location in the water-ice rich northern plain-terrains of the trailing hemisphere.
• The presence of amorphous vs. crystalline water ice in the three spectra can possibly be understood through a balance of the opposing forces of heating (higher temperatures favor crystalline ice) and irradiation (higher irradiation levels favor amorphous ice). S1 shows weak evidence of amorphous ice (2.4σ) and no evidence of crystalline ice, which might indicate that the irradiation induced amorphization dominates over the higher than average temperatures near the trailing-hemishpere's apex. S2 shows no evidence for amorphous ice but a strong evidence for crystalline ice (5.9σ), consistent with the irradiation being weaker at S2's location in the anti-Jovian hemisphere, despite the equatorial latitudes. Moreover, S2's location has numerous linea, implying heating from differential motions, which favors crystallization. Finally, S3 shows strong evidence for both forms of ice, perhaps because it belongs to the cold northern plain-terrains where water ice is more abundant. The strong evidence for crystalline ice in S2 and S3 is also consistent with the presence of the 1.65 µm absorption feature seen clearly in the two spectra ( Figure 1). We note that this picture of the dominance of amorphous over crystalline ice and vice-versa, as a balance of particle irradiation and temperature, is may be a simplified picture of a complex set of factors that are yet to be well-understood.
• The grain-sizes retrieved for crystalline ice in S2 (382.44 +78.19 −70.65 microns) and S3 (919.14 +58.26 −96.88 microns) are rather large. We interpret these large grain-sizes as being representative of the scale of water ice concentrations formed due to the active sublimation of water-ice. Also, since our model retrieves these large grain-sizes, the corresponding abundances of crystalline ice for S2 (0.24 +0.0012 −0.0007 %) and S3 (0.36 +0.0009 −0.0006 %) are fairly low. These two parameters are usually highly anti-correlated, as evidenced in the 2D parameter distribution plot for S2 (Figure 11).
• The wavelength-shift parameter, ∆ wav , shows a clear decreasing trend from S1 to S3. The shift, towards shorter wavelengths, is almost a channel-width (0.0198 +0.0015 −0.0016 µm for S1, about half of that at 0.0098 +0.001 −0.001 µm for S2, and negligible at 0.0043 +0.0015 −0.0015 for S3). This fits well with the general expectation that S1 comes from the most irradiated, visibly dark, equatorial regions of Europa's trailing hemisphere, while S3 is from the northern latitudes where irradiation is much weaker.
• We retrieve high porosities for all three spectra. In particular, S3 requires a high porosity for the model to match the low reflectance levels of the data (Figure 9). While our porosities are consistent with previous studies that use polarimetric and photometric phase curves to constrain the porosity of Europa, more laboratory work is needed to test the ability of the K parameter (porosity coefficient) in Hapke's model to constrain the porosity of an icy regolith using spectroscopic data.
Our work illustrates that the application of physically motivated reflectance models, in state-ofthe-art statistical frameworks, can yield new insights from the archival Galileo NIMS data of Europa. The technique of quantifying the evidence of a given species from the data, using Bayesian model comparison, holds promise for uncovering trace species abundances like organics on Europa, which is one of the primary goals for next generation missions like Europa Clipper (Blaney et al. 2017) and JUICE (Grasset et al. 2013). Crucially, like many studies before (e.g. Dalton 2010;Dalton et al. 2012;Ligier et al. 2016;Shirley et al. 2016;Filacchione et al. 2019;Mishra et al. 2021), our study highlights the need for laboratory measurements of cryogenic optical constants relevant to Europan conditions. The optical constants of irradiated samples are also important to accurately model the conditions on Europa. For example, we now know that irradiated sodium chloride exists on Europa due to laboratory work that demonstrated the development of color-centers in the salt when irradiated (Trumbo et al. 2019a). As more optical constants of species relevant to Europa become available, we plan to incorporate them into our framework and revisit the Galileo NIMS dataset, as well as NIR Europan data collected by other missions like Cassini (Brown et al. 2003) and New Horizons (Grundy et al. 2007), where plenty of new knowledge about the complex Europan surface composition might be awaiting.  log 10 (D cr ) Figure 9. Sensitivity of our best model for S3 to φ and D cr . φ and D cr are varied in the left and right panels respectively, while all the other parameters are kept fixed to the retrieved median solution for S3 (Table 4). The black points in both panels are the S3 data. Only for very small values of φ and very large values of D cr , is the model able to match the low reflectance levels of the data. . Each diagonal plot also shows the corresponding median values and the associated 1σ lower and upper limits for the distribution. The remaining plots are pair-wise 2D distributions, that illustrate the correlations between parameters. The contours in these 2D distributions correspond to 0.5,1,1.5 and 2.0σ intervals. Parameter symbols are described in Table 2.