Compositional Mapping of Europa Using MCMC Modeling of Near-IR VLT/SPHERE and Galileo/NIMS Observations

SPHERE (Spectro-Polarimetric High-contrast Exoplanet REsearch; Beuzit et al. 2019) is an instrument on the 8 m Very Large Telescope (VLT) UT3, located at the European Southern Observatory's Paranal Observatory in Chile. It was designed primarily as an exoplanet imager, using adaptive optics to provide stable high spatial resolution observations, enabling direct imaging of exoplanets orbiting their host stars. This high spatial resolution can also be applied to other classes of objects, such as the solar system targets like Europa.

Observations of Europa were taken during SPHERE science verification in 2014 December, as summarized in Table 1 and shown in Figure 1. The observation block consists of multiple Europa observations, and then a single calibration star observation, taken immediately after the Europa observations. This ensured that the atmospheric conditions for the science and calibration observations were as similar as possible. We used SPHERE in IRDIFS_EXT mode, allowing simultaneous imaging with the Integral Field Spectrograph (IFS) and Infrared Differential Imaging Spectrometer (IRDIS) subsystems of the SPHERE instrument. This allows simultaneous spectroscopic observations covering 0.95–1.65 μm (J and H bands) and dual-band imaging at 2.1 and 2.251 μm (K band).

Zoom In Zoom Out Reset image size

The IFS is a relatively new tool in astronomy, enabling both spatial and spectral information to be recorded simultaneously on the same detector. The IFS on SPHERE (Claudi et al. 2008; Mesa et al. 2015) produces image cubes with 38 wavelength channels from 0.95 to 1.65 μm, giving a low spectral resolution of R = λ/Δλ ∼ 30. It has a high spatial resolution, with a pixel size of 7.46 mas pixel⁻¹, corresponding to ∼25 km pixel⁻¹ at Jupiter. Accounting for diffraction, this allows features ∼150 km across to be resolved. The magnitude of the diffraction is calculated by measuring the FWHM of the point spread of the calibration star, observed through VLT/SPHERE (≈6 pixels at 0.95 μm).

The IFS operates by using a lenslet array that focuses the image into a series of spaxels (spatial pixels), spots in the focal plane, each of which is from a different area of the observed field. These spaxels are used as the input into a dispersive spectrograph that transforms each spaxel into a full spectrum for its respective part of the image. These multiple spectra (one from each spaxel) are imaged by a single detector, recording both the spatial and spectral information from the observation. In the reduction process, each spectrum is fitted and transformed into the spectral dimension of the output image cube so that the final image cube has one slice for each measured wavelength (Claudi et al. 2008).

IRDIS (Dohlen et al. 2008; Vigan et al. 2010) produces simultaneous imaging through two filters, producing two images on separate parts of the same detector. Different filter combinations are available so that contrast can be maximized depending on the spectral features of individual targets. Our data use the DB_K12 filter set that has transmissions centered at 2.1 μm for the K1 filter and 2.251 μm for the K2 filter, with filter widths of 0.051 and 0.055 μm, respectively, meaning that IRDIS is sensitive to water-ice absorption at 2 μm.

The data reduction routine transforms the raw observations into cleaned and calibrated data sets that are ready for scientific use. The reductions account for detector properties, atmospheric conditions, and observation conditions to ensure that observations are calibrated and comparable.

The raw IRDIS data consist of FITS image files containing two images side by side on the same detector, one from each filter. The IRDIS reduction process is written in Python and is as follows:

• 1.  
  Perform instrument calibration by subtracting dark frame from observations and dividing by flat frame.
• 2.  
  Split the image into two images, one for each filter.
• 3.  
  Clean the images by "despiking" to remove extreme bright or dark pixels.
• 4.  
  Sum images from each observation's detector integrations into a single image for each filter in each observation.

After reduction, the IRDIS images can be corrected and mapped in a similar way to the IFS data (below).

The raw IFS data consist of FITS files containing individual spectra for each spaxel, arranged across the detector. Our reduction process consists of a series of steps to ultimately transform these individual spectra into a mapped spectral cube for the observed target. The reduction routine is written in Python and makes use of EsoRex, the ESO Recipe Execution Tool that provides standard reduction pipelines for ESO instruments. An overview of our IFS data reduction routine is as follows:

• 1.  
  Clean raw images by replacing known bad pixels with interpolated values.
• 2.  
  Run EsoRex reduction to identify spectral locations on detector, perform instrument calibrations (using detector flat and dark frames), and generate image cubes. This produces an image cube for each individual detector integration in each individual observation.
• 3.  
  Sum image cubes from each observation's detector integration into a single cube for each observation. This improves signal-to-noise ratio and improves processing efficiency, as the number of files to process is decreased significantly.
• 4.  
  Clean the image cubes by "despiking" (replacing extreme bright or dark pixels) and "destriping" to remove a prominent banding pattern (see Section 2.2.1).
• 5.  
  Calibrate the image cubes to remove telluric contamination and calculate reflectance spectra, using the calibration star observation (see Section 2.2.2).
• 6.  
  Map the observations by performing a photometric correction and transforming the images to an equirectangular map projection (see Section 2.2.3).
• 7.  
  Combine the maps from individual observations into a single mapped cube for each target. This gives a complete spectrum for each observed point on the map.

2.2.1. Image Cleaning

A strong regular banding pattern (see Figure 2(a)) remains after the EsoRex reduction for extended sources such as Europa, so we developed a method using a Fourier transform filter to identify and remove this artificial pattern. This pattern is likely to be produced from crosstalk between different lenslets not being fully accounted for in the EsoRex step in the reduction pipeline. The pattern is regular and fixed in position for different wavelengths and observations, so it can be systematically removed using a Fourier transform filter to remove the specific frequencies that generate the pattern. The strength of the pattern is wavelength dependent and is strongest at ∼1.35 μm.

Zoom In Zoom Out Reset image size

To remove this pattern, its specific frequencies were identified by analyzing the Fourier transform of IFS images. These frequencies are removed from the image by converting the image to Fourier space, setting the frequencies to zero, and then converting back to image space using

$$\begin{eqnarray}&&\mathrm{Destriped}={{ \mathcal F }}^{-1}\left\{\mathrm{Mask}\times { \mathcal F }\left\{\mathrm{Image}\right\}\right\},\end{eqnarray} \tag{ 1 }$$

where ${ \mathcal F }$ is the Fourier transform (image domain to frequency domain) and ${{ \mathcal F }}^{-1}$ is the inverse Fourier transform (frequency domain to image domain). The mask used to remove the unwanted frequencies is shown in Figure 3.

Zoom In Zoom Out Reset image size

The destriping routine removes the pattern for wavelengths where it occurs, and it has negligible effect for wavelengths where the pattern does not occur. Therefore, it can be safely applied to every wavelength of an image cube to avoid introducing any wavelength-dependent systematic errors.

Images are "despiked" by calculating the mean and variance of all pixels within a 2-pixel radius of a given pixel. If the pixel value is greater than two standard deviations from the mean of the surrounding pixels, it is assumed to be a "bad" pixel, and its value is replaced with the mean value. This removes extreme spurious values while preventing loss of actual data, as diffraction means that any variation in the data would be seen across multiple pixels. The despiking routine is run repeatedly until there are no bad pixels to replace. It is based on the sigma_filter routine in the IDL Astronomy User's Library (Landsman 2016).

See Figure 2 for a summary of this image-cleaning process.

2.2.2. Spectral Calibration

The spectral calibration calculates the reflectance spectrum of the observed moon from the measured flux. The calibration star observation is used to radiometrically calibrate the spectrum, removing any telluric contamination in the process. This process, based on the telluric correction step described in Ligier et al. (2016), uses the known spectrum of the calibration star and the known solar spectrum to calculate the reflectance

$$\begin{eqnarray}\begin{array}{rcl}R({\lambda }) & = & \displaystyle \frac{{F}_{{\rm{m}}{\rm{o}}{\rm{o}}{\rm{n}}}({\lambda })-B({\lambda })}{{F}_{{\rm{s}}{\rm{t}}{\rm{a}}{\rm{r}}}({\lambda })}\\ & & \times \displaystyle \frac{{T}_{{\rm{s}}{\rm{t}}{\rm{a}}{\rm{r}}}}{{T}_{{\rm{m}}{\rm{o}}{\rm{o}}{\rm{n}}}}\times \displaystyle \frac{1}{{\rm{\Omega }}\times {10}^{\tfrac{{m}_{{\rm{s}}{\rm{t}}{\rm{a}}{\rm{r}}}({\lambda })-{m}_{{\rm{S}}{\rm{u}}{\rm{n}}}({\lambda })}{2.5}}},\end{array}\end{eqnarray} \tag{ 2 }$$

where R(λ) is the reflectance spectrum of Europa. F_moon and F_star are the measured fluxes from Europa and the calibration star, respectively, B(λ) is the background flux, T is the integration time of an observation, Ω is the solid angle subtended by each detector pixel, m_star(λ) is the apparent magnitude of the calibration star at Earth (Cutri et al. 2003), and m_Sun(λ) is the apparent magnitude of the Sun at Europa (Blanton & Roweis 2007).

The stellar and background fluxes are calculated by summing the flux in two circular apertures, both centered on the calibration star, that have radii calculated from the size of the Airy disk. The star aperture is chosen so that the first airy diffraction maximum is contained within the aperture while still minimizing the background flux within the aperture. The background flux is summed from a 10-pixel-wide annular aperture around the same calibration star, which reduces the effect on the calibration of any slight variation in sensitivity across the detector. The average spectral slope of the SPHERE data is finally adjusted to be equal to the spectral slope from the NIMS data to remove any slight variations caused by differing spectral slopes from the calibration star spectrum and the solar spectrum. This removes the need to allow the continuum spectral slope to vary as part of the spectral fitting process as in Ligier et al. (2016).

2.2.3. Mapping

2.2.4. Photometric Correction

A sphere of uniform albedo illuminated by a point-like source (i.e., the Sun) appears to have varying brightness, with brightness dropping off toward the edge of the disk (see Figure 4(a)). The photometric effect must be accounted for to correct the reflectance for areas of an observed disk away from the center. The correction takes the form

$$\begin{eqnarray}&&\mathrm{Corrected}\ \mathrm{image}=\displaystyle \frac{\mathrm{Image}}{\mathrm{Photometric}\ \mathrm{model}},\end{eqnarray} \tag{ 3 }$$

where the photometric model gives the ratio of the brightness of the disk to the brightness of a surface viewed and illuminated at a normal incidence angle. Therefore, the values of the corrected image give the reflectance each location would have if it were both the subsolar point and the subobserver point (i.e., where the incidence and emission angles are zero).

Zoom In Zoom Out Reset image size

The simplest photometric correction assumes a Lambertian surface—a surface that is a perfect diffuse scatterer and therefore appears equally bright from all directions (Hapke 1993; Lambert 1760). Therefore, the brightness is simply proportional to the projected area of the surface to the illumination source, $\cos (i)$, where i is the angle of the source from the surface normal. The Lambertian correction improves on the basic images but accumulates significant errors toward the edge of the disk, overcompensating for the limb darkening and producing unphysical reflectance values (>1) (see Figure 4(b)).

Therefore, we used the Oren−Nayar reflectance model (Oren & Nayar 1994), which modifies the Lambertian model by accounting for surface roughness, a significant improvement on the implicit assumption of a perfectly smooth surface in the Lambertian model. The Oren−Nayar model assumes that the surface is composed of a series of V-shaped cavities, with facets that are themselves Lambertian scatterers. The cavities are modeled to have a Gaussian distribution of slopes, with a mean slope of μ = 0 and a standard deviation of σ that parameterizes the surface roughness. The Oren−Nayar model reduces to standard Lambertian scattering for σ = 0.

For a surface roughness σ and albedo ρ, the modeled Oren−Nayar surface brightness is

$$\begin{eqnarray}&&{f}_{r}=\cos (i)\left({L}_{1}+{L}_{2}\right)\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}\begin{array}{rcl}{L}_{1} & = & \displaystyle \frac{\rho }{\pi }({C}_{1}+\cos (\phi )\tan (\beta ){C}_{2}\\ & & +(1-| \cos (\phi )| )\tan \left(\displaystyle \frac{\alpha +\beta }{2}\right){C}_{3})\end{array}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{L}_{2}=0.17\displaystyle \frac{{\rho }^{2}}{\pi }\left(\displaystyle \frac{{\sigma }^{2}}{{\sigma }^{2}+0.13}\right)\left(1-\cos (\phi ){\left(\displaystyle \frac{2\beta }{\pi }\right)}^{2}\right),\end{eqnarray} \tag{ 6 }$$

where (i, ϕ_i ) and (e, ϕ_e ) are the polar and azimuth angles of the incident and emitted rays, respectively. L₁ is the contribution from single scattering, and L₂ is the contribution from multiple scattering, with

$$\begin{eqnarray}&&\phi ={\phi }_{e}-{\phi }_{i}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\alpha =\mathrm{Max}(e,i)\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\beta =\mathrm{Min}(e,i).\end{eqnarray} \tag{ 9 }$$

The single scattering coefficients C_n are dependent on the surface roughness σ and the viewing geometry,

$$\begin{eqnarray}&&{C}_{1}=1-0.5\displaystyle \frac{{\sigma }^{2}}{{\sigma }^{2}+0.33}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}{C}_{2}=\left\{\begin{array}{ll}0.45\displaystyle \frac{{\sigma }^{2}}{{\sigma }^{2}+0.09}\sin (\alpha ) & \mathrm{if}\ \cos (\phi )\geqslant 0\\ 0.45\displaystyle \frac{{\sigma }^{2}}{{\sigma }^{2}+0.09}\left(\sin (\alpha )-{\left(\displaystyle \frac{2\beta }{\pi }\right)}^{3}\right) & \mathrm{otherwise}\end{array}\right.\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{C}_{3}=0.125\left(\displaystyle \frac{{\sigma }^{2}}{{\sigma }^{2}+0.09}\right){\left(\displaystyle \frac{4\alpha \beta }{{\pi }^{2}}\right)}^{2}.\end{eqnarray} \tag{ 12 }$$

Incorporating roughness increases the modeled brightness at the edge of the disk when compared to the Lambertian model, as shown in Figure 5. The surface roughness parameter σ was selected by photometrically correcting the observation with a range of σ values and visually comparing the corrected brightness at the limb and toward the center of Europa's disk. Large values of σ produce undercorrected observations where the limb is too dark (e.g., Figure 4(a)), and small values of σ produce overcorrected observations where the limb is too bright (e.g., Figure 4(b)). A value of σ = 33° (Figure 4(c)) was found to produce the minimum brightness variation and so was selected as our surface roughness parameter. Some small residual photometric errors are still present at the very edge of the disk with σ = 33°, so we conservatively constrain the region of useful data to e < 75°.

Zoom In Zoom Out Reset image size

The only other free parameter required by the model is the surface's albedo. Variation in the albedo has a very negligible effect on the modeled photometry, so precise fitting of this value was not required and a value of ρ = 0.5 was used.

To account for the diffraction of the telescope optics, the photometric model is convolved with an airy disk model before being used to correct the image. The radius of the airy disk is proportional to the wavelength of the specific image being corrected, with the radius set to r_airy = 10 pixels for 1.65 μm, as measured from observations of the calibration star (see Figure 4(d) for the final photometric model and corrected observation).

This photometric correction allows accurate correction of the images to emission angles e > 70°, corresponding to ∼90% of the observed disk. This is higher than previous studies using Lambertian models that typically extend to 50°–60° (Grundy et al. 2007; Brown & Hand 2013; Ligier et al. 2016). Improvements beyond ∼75° are unlikely to provide much more useful data, as the extreme viewing angles at the edge of the disk produce very low spatial resolutions once these areas are mapped. Therefore, we conservatively use e = 75° as the limit of the reliable mapped region of the data.

The Galileo orbiter was launched in 1989 and orbited Jupiter from 1995 to 2003. Galileo's orbit included a series of flybys of Europa, enabling detailed imaging and mapping of Europa's surface.

The NIMS (Carlson et al. 1992) instrument on Galileo performed spectroscopy from 0.7 to 5.2 μm with a spectral resolution of R = λ/Δλ ∼ 60. During flybys, a series of scans by the NIMS instrument would map the infrared spectra of Europa's surface beneath the spacecraft. The spatial resolution of the mapping was dependent on the distance between Galileo and Europa, with spatial resolutions <1 km pixel⁻¹ for the closest flybys. However, there are also large regions of Europa's surface with no NIMS data or only very low spatial resolution data (>500 km pixel⁻¹).

The NIMS instrument used different detectors to cover different wavelength ranges, one of which (detector 3, covering 0.99–1.26 μm) failed early in the Galileo mission, during Galileo's fourth orbit of Jupiter, so there are no Europa observations with full spectral coverage from 1997 onward. Therefore, there are large regions of Europa with no observations or only very low spatial resolution observations at these wavelengths. The SPHERE IFS wavelength range includes this spectral range, so SPHERE can fill this missing gap in the near-infrared spectra of the Galilean moons.

The NIMS data sets given in Table 2 were used in this study. These data sets were taken at a low solar phase angle and collectively cover Europa's anti-Jovian hemisphere, allowing direct comparison to our SPHERE observations of the same region. The data sets 17ENGLOBAL01A and 17ENGLOBAL02A were taken after the failure of NIMS detector 3 and so lack the 0.99–1.26 μm spectral range.

All NIMS observations of Europa were downloaded from the NASA Planetary Data System and reduced using a pipeline similar to that for our ground-based observations. The I/F NIMS data cubes were processed using ISIS3 (U.S. Geological Survey 2020) to calculate the viewing angles (solar incidence angle, emission angle, and phase angle) and coordinates (latitude and longitude) for each observed pixel in the NIMS images. These angles and coordinates were used to map and photometrically correct the data sets using the same mapping and photometric correction code as the SPHERE data set.

For our analysis, we selected NIMS observations taken at a low phase angle with a relatively high spatial resolution (≲100 km pixel⁻¹) and similar spatial coverage to our SPHERE observation (spatial resolution ∼150 km), allowing direct comparison between the two data sets.

The mapped spectral cubes are analyzed by fitting to laboratory spectra from reference cryogenic libraries (described in Section 4.3). Our fitting routine uses linear spectral modeling to fit the observed spectra, where the modeled spectrum, M, is

$$\begin{eqnarray}&&M(\lambda )=\sum _{i}{w}_{i}{E}_{i}(\lambda )\times S,\end{eqnarray} \tag{ 13 }$$

where w_i are the weights of the different endmembers E_i . The parameter S accounts for any continuum spectral slope that may remain in the data. The weights w_i give the fractional abundance of each endmember in the modeled spectrum and are subject to the constraint

$$\begin{eqnarray}&&\sum _{i}{w}_{i}=1,\end{eqnarray} \tag{ 14 }$$

which ensures that the fractional abundance of the different endmembers sums to 100%, and

$$\begin{eqnarray}&&0\leqslant {w}_{i}\leqslant 1,\end{eqnarray} \tag{ 15 }$$

which ensures that the individual abundances are all physically realistic values (0% to 100%).

The use of linear unmixing is generally valid for ground-based observations, as the relatively low spatial resolutions (∼100 km) mean that the observed spectrum for each pixel is naturally a linear combination of the spectra of different geological units within that pixel (Ligier et al. 2016). This method will not, however, account for any nonlinear scattering that occurs, and it may therefore be responsible for some small residual errors after fitting (Ligier et al. 2016; Shirley et al.2016).

Markov Chain Monte Carlo (MCMC) simulates systems by sequentially randomly sampling a multidimensional parameter space. MCMC uses a series of "walkers" that sequentially model the system in a chain of different parameter values, ultimately building a posterior distribution of parameter values consistent with the observed data.

Bayes's theorem (Bayes 1763) states that the posterior probability of a set of parameter values w given an observation o can be calculated as

$$\begin{eqnarray}&&P\left({\boldsymbol{w}}| {\boldsymbol{o}}\right)\propto P\left({\boldsymbol{o}}| {\boldsymbol{w}}\right)P\left({\boldsymbol{w}}\right),\end{eqnarray} \tag{ 16 }$$

where $P\left({\boldsymbol{o}}| {\boldsymbol{w}}\right)$ gives the likelihood of measuring the observation o given the parameters w . $P\left({\boldsymbol{w}}\right)$ gives the prior probability of a set of parameter values and is used to include any information known before taking any observations (e.g., the constraint that abundances must sum to 100%).

At a given iteration of the MCMC chain, a walker will have a specific set of free parameter values, which can be described as the vector w _n . These parameter values can be related to the observed data using a cost function

$$\begin{eqnarray}&&C({\boldsymbol{w}})\propto -\mathrm{log}\left(P\left({\boldsymbol{w}}| {\boldsymbol{o}}\right)\right),\end{eqnarray} \tag{ 17 }$$

where parameters that fit the observed data will have a low cost and parameters that have a poor fit have a high cost. In the spectral modeling used in this study, the free parameters are the abundances of the endmembers (w_i in Equation (13)).

At each step in the MCMC chain, the walker selects a new set of parameter values to test t _n and moves to the new parameters with a probability proportional to the ratio of the cost functions C( w _n )/C( t _n ). This means that if t _n+1 is a better fit to the data, the walker will likely move to those parameters ( w _n+1 = t _n ); otherwise, it will likely remain in the same place ( w _n+1 = w _n ). Therefore, each walker generally moves toward and then remains in a region of parameter space that has a good fit to the data.

After an initial "burn-in" period, the probability of the walker occupying a given set of parameter values is proportional to the posterior probability of the set of parameter values being consistent with the observed data. Therefore, the positions of a large ensemble of walkers can be sampled to produce a posterior probability distribution of the parameter values.

When applied to spectral modeling, MCMC enables the simulation of reflectance spectra, and the calculation of the posterior probability distribution of endmember abundances is consistent with an observed spectrum. The posterior distribution for each endmember can be sampled to calculate the most likely abundance value and its associated uncertainty.

The use of MCMC techniques allows more robust modeling of reflectance spectra than simple linear optimization by quantifying the uncertainties of and correlations between endmember abundances (Lapotre et al. 2017). This ultimately allows us to measure the confidence of different detections, which is especially important for observations with relatively low spectral resolutions and potentially degenerate fit results where multiple different compositional mixes are possible (Brown & Hand 2013).

Our MCMC modeling uses the "emcee" library to directly calculate the walker Markov Chains. See Foreman-Mackey et al. (2013) for a more detailed discussion of MCMC and the specific emcee implementation.

The only free parameters in our MCMC modeling are the abundances of each endmember (w_i in Equation (13)). All abundances are assumed to be equally likely, so a constant prior between 0 and 1 is used. This means that any potential solutions that would imply unphysical abundance values (i.e., less than 0% or greater than 100%) are rejected by the MCMC process.

We assume that all possible combinations of abundance values are equally likely so that we have a constant prior for all physically possible mixtures. Therefore, we have no prior bias between mixtures consisting of a single endmember and mixtures containing many endmembers. This is equivalent to the Dirichlet prior with all concentration parameters equal to 1 used in Lapotre et al. (2017). The sum-to-1 constraint, i.e., that the total abundance of all N endmembers should be 100% (Equation (14)), is directly enforced by allowing the abundances of N − 1 endmembers to vary and then setting the final endmember's abundance to be ${w}_{N}=1-{\sum }_{i}^{N-1}{w}_{i}$. This is mathematically equivalent to rejecting all solutions where the abundances do not sum to 1, but it is computationally much more efficient.

Due to its underlying randomness, MCMC requires simulations using large numbers of walkers, each with long chains to produce a stable output. This makes MCMC computationally intensive, typically requiring several minutes on a single processor core to run an MCMC fit for a single spectrum. Therefore, we used an optimized fitting routine to perform accurate fits across the whole observed disk within a reasonable timescale:

• 1.  
  Perform a simple linear fit to calculate optimized parameter values, ${w}_{i}^{\mathrm{LF}}$.
• 2.  
  Generate initial positions, pos0, for MCMC chains as a random Gaussian ball centered on ${w}_{i}^{\mathrm{LF}}$ with a standard deviation of initialisation_sd. These initial chain positions are corrected to ensure 0 ≤ w_i ≤ 1 and ∑_i w_i = 1. Initializing around the linear fit values significantly decreases the required chain length, as the chain is initialized in a realistic region of the parameter space that is likely to be close to the region it will converge to.
• 3.  
  Run emcee for burn_in_steps steps using pos0 as the initialization. This "burns in" the chain to ensure that the chain position is relatively independent of the initialization.
• 4.  
  Continue running emcee in a series of runs with run_steps steps. Each run is initialized with the final position of the previous run (i.e., it directly continues the chain). The abundance values are recorded for each successive run by calculating the median value over all walkers in the run.
• 5.  
  If the median values of every parameter vary by less than convergence_difference over convergence_runs successive runs, the MCMC chain is assumed to have converged and the simulation is stopped. Otherwise, the simulation is automatically stopped at max_steps steps to prevent indefinite calculations (though in practice this limit is never reached).
• 6.  
  The final abundance values and associated uncertainties are calculated from the posterior distribution of values produced by the final run. We take the median of the posterior distribution as the best estimate of the abundance and use the 16th and 84th percentiles of the distribution to calculate the 1σ uncertainty on this best estimate.

This fitting routine has been tested over a wide range of inputs to ensure that it remains accurate and the choice of different parameters does not influence the final calculated values. For example, initializations with completely random values for pos0 ultimately converge to values identical to initializations using linear fit values, but they just take significantly longer to do so. Therefore, using the linear fit values is an effective shortcut to reduce calculation time without sacrificing accuracy. The values of the parameters used for our MCMC fitting routine are given in the Appendix.

The abundances of a single class of endmembers (e.g., all ice endmembers) can be combined into a single abundance distribution by summing the abundances of the individual endmembers at each time step for each walker. This produces a posterior distribution for the whole class of endmembers that often has a much smaller uncertainty than the individual endmember abundances that contribute to it. For example, if two specific endmembers have abundances of 10% and 40%, respectively, at one time step and then 30% and 20% at another time step, their individual uncertainties appear to be relatively large. However, the uncertainty on the combination of these two endmembers is much smaller, as at both time steps their summed abundance is 50%.

The full list of endmembers in our spectral library is given in Table 3, and selected spectra are shown in Figure 8. These cryogenic reference spectra include hydrated sulfuric acid (Carlson et al. 1999) and a variety of hydrated salts (Dalton 2007; Dalton & Pitman 2012; Hanley et al. 2014). The sulfuric acid spectra (Carlson et al. 1999) do not cover wavelengths below 1 μm, so this limits our modeling to cover the 1–2.5 μm spectral range. Our spectral library includes all available laboratory spectra covering the SPHERE wavelength range measured in Europa-like conditions, and particular care was taken to ensure the inclusion of the various non-ice species detected in prior works (Carlson et al. 1999; McCord et al. 2002; Brown & Hand 2013; Ligier et al. 2016; Trumbo et al. 2019).

Zoom In Zoom Out Reset image size

Our water-ice reflectance spectra are calculated from measured refractive indices (Grundy & Schmitt 1998) using the Hapke bidirectional reflectance model (Hapke 1993). Modeling was performed using a variety of simulated ice grain sizes ranging from 1 μm to 1 cm to fully explore the grain size parameter space. Initial models did not identify any extremely small (<30 μm) or large (>3 mm) grains, so our final model presented here covers grains ranging from 30 μm to 3 mm.

The final modeled water-ice spectra use four size bins ranging from 30 μm to 3 mm, where the spectrum is a blend of grain sizes within each bin. Each blended spectrum is calculated by modeling 250 discrete grain sizes linearly spaced within the size bin and then calculating the mean of these 250 spectra. This provides a more physically realistic simulation of the grain sizes, as we would not expect there to only be a single discrete ice grain size present on Europa's surface. This averaging also removes any residual Mie oscillations in the calculated spectrum that can occur when a single exact grain size is used for modeling (Hapke 1993).

The hydrated sulfuric acid distribution shown in Figure 11 is mainly concentrated toward the trailing hemisphere, with lower abundances toward the leading hemisphere and at higher latitudes. There is no strong correlation with geological units, and the distribution for both SPHERE and NIMS follows the expected bull's-eye distribution centered on the trailing apex, with a clear distinction between the leading and trailing hemispheres. This spatial distribution is consistent with exogenic plasma bombardment being the dominant source of the sulfuric acid hydrate present on Europa's surface.

Zoom In Zoom Out Reset image size

The peak abundance is 68% for SPHERE (1σ range 60%–75%) and 65% for NIMS (1σ range 60%–70%), both of which occur around the trailing apex. These values are similar to the ∼65% maximum abundance in Ligier et al. (2016) though lower than the ∼90% reported in Carlson et al. (2005). However, it is important to note that these previous studies have full spatial coverage of the trailing apex, whereas the SPHERE and especially NIMS data sets shown in Figure 11 only cover some of the eastern part of the trailing hemisphere, so the values are not directly comparable.

The modeled water-ice distribution in Figure 12 shows the same spatial pattern as the 1.5 μm (Figure 6) and 2 μm (Figure 7) absorption bands, with high abundance at high latitudes and lower abundance toward the trailing apex. The water-ice distribution is highly anticorrelated with the sulfuric acid distribution, with the combined hydrated sulfuric acid and ice abundance ∼80% when averaged over the observed area. This anticorrelation is a real surface feature (and not a spectral degeneracy), where the spatial distribution of water-ice abundance is driven by the absence of contaminants, mainly acids in the trailing hemisphere and to a lesser extent salts in the leading hemisphere. The highest water-ice abundances occur at high latitudes, where the exogenic plasma bombardment intensity is significantly lower.

Zoom In Zoom Out Reset image size

As found in previous studies (Dalton et al. 2012; Ligier et al. 2016), the majority of the surface is dominated by 100 μm–1 mm grains, with larger grains (>300 μm) on the trailing hemisphere and northern latitudes and smaller grains on the leading hemisphere. There also appears to potentially be a small localized area of >1 mm grains in Dyfed Regio; however, the overall ice abundance in this region is very low, so the confidence of this detection of larger grains is lower.

A variety of hydrated salts were identified in contaminated regions, especially around Dyfed Regio and Powys Regio. As with previous analyses of NIMS observations (McCord et al. 2002; Dalton 2007; Shirley et al. 2010), sulfates appear to be the main salts present on Europa's surface (see Figure 13). Magnesium-bearing salts (MgSO₄ · nH₂O, MgSO₄ brine, and Na₂Mg(SO₄)₂ . 4H₂O) appear correlated with dark terrain, particularly Powys Regio. As shown in the bottom panels of Figure 13, the lower bounds on the abundances of magnesium sulfates are all very low and close to zero. This suggests that the uncertainties and degeneracies between the salt spectra mean that it is not possible to positively identify any individual magnesium sulfate salts with the SPHERE and NIMS spectral resolutions.

Zoom In Zoom Out Reset image size

Mirabilite (Na₂SO₄ · 10H₂O) appears more abundant in eastern areas of the observations in Figure 13, with a relatively clear distinction between low abundances in the trailing hemisphere and higher abundances into the leading hemisphere. Mirabilite has the highest abundances (∼20%) of any of the modeled salts, and the corresponding nonzero lower bounds suggest that it is likely to be present on the leading hemisphere. The modeled abundances suggest an upper bound of ∼30% for mirabilite abundance on the anti-Jovian hemisphere.

Chlorinated salts (Figure 14) generally have lower best-estimate abundances than the sulfate salts, and all have lower bounds close to zero, implying that it is difficult to positively identify individual salts, particularly with SPHERE. The spatial distribution of the chlorinated salts appears more globally uniform than the sulfate salts, with slightly higher abundances around geological units.

Zoom In Zoom Out Reset image size

Hydrated magnesium chlorate (Mg(ClO₃)₂ · 6H₂O), perchlorate (Mg(ClO₄)₂ · 6H₂O), and sodium perchlorate (NaClO₄) all have very low abundances and appear slightly correlated with geological units. Magnesium perchlorate has a single region of relatively high abundance in the NIMS observation of Dyfed Regio, which is consistent with the location of the highest abundance found in Ligier et al. (2016). Sodium chloride (NaCl) and sodium perchlorate also have a low abundance in our observations, with Powys Regio and Dyfed Regio the only areas with tentative detections. UV absorptions attributed to sodium chloride have been previously detected around Europa's leading apex using HST (Trumbo et al. 2019); however, this region is outside of our observed area and so is consistent with our low abundances.

Hydrated magnesium chloride (MgCl₂ · nH₂O) appears the most abundant chlorinated salt, with best-estimate abundances up to 10% in the NIMS data set. The distribution of magnesium chloride appears correlated with disrupted geological units, particularly Dyfed Regio and the northern lineae in the SPHERE observation and Argadnel Regio in the NIMS observation. This distribution appears consistent with Ligier et al. (2016), who found abundances of ∼10% around Dyfed Regio.

The data sets used in this study can be found at 10.5281/zenodo.6034735.

The list of parameters and their values used for our MCMC fitting routine are given in Table 4.

The abundance values for the case study sites are given in Table 5.

Table 5. Abundances Values for Case Studies Shown in Figure 9

	Powys Regio		Dyfed Regio
Endmember	SPHERE	NIMS	SPHERE	NIMS
H2O 30–100 μm	${5.6}_{-4.1}^{+7.9} \%$	${4.1}_{-2.9}^{+4.8} \%$	${1.3}_{-1.0}^{+2.0} \%$	${1.2}_{-0.8}^{+1.4} \%$
H2O 100–300 μm	${15.7}_{-8.8}^{+8.1} \%$	${21.7}_{-5.7}^{+5.7} \%$	${2.9}_{-2.2}^{+3.7} \%$	${2.6}_{-2.0}^{+3.3} \%$
H2O 300 μm–1 mm	${5.8}_{-4.1}^{+6.4} \%$	${2.2}_{-1.6}^{+2.9} \%$	${12.5}_{-7.3}^{+6.3} \%$	${7.6}_{-4.8}^{+5.5} \%$
H2O 1–3 mm	${2.7}_{-2.1}^{+3.8} \%$	${0.7}_{-0.5}^{+1.4} \%$	${3.6}_{-2.7}^{+4.0} \%$	${3.9}_{-2.9}^{+4.5} \%$
H2SO4 · 6.5H2O	${5.9}_{-4.3}^{+6.4} \%$	${2.2}_{-1.6}^{+3.2} \%$	${1.9}_{-1.5}^{+3.0} \%$	${1.5}_{-1.1}^{+2.2} \%$
H2SO4 · 8H2O 5 μm	${1.7}_{-1.2}^{+2.5} \%$	${0.7}_{-0.5}^{+1.1} \%$	${0.8}_{-0.6}^{+1.3} \%$	${0.6}_{-0.4}^{+1.0} \%$
H2SO4 · 8H2O 50 μm	${17.6}_{-9.0}^{+8.6} \%$	${23.3}_{-5.7}^{+5.7} \%$	${45.1}_{-6.6}^{+6.4} \%$	${43.2}_{-4.8}^{+4.7} \%$
Mg(ClO3)2 · 6H2O	${1.3}_{-1.0}^{+2.1} \%$	${0.5}_{-0.4}^{+1.0} \%$	${1.7}_{-1.3}^{+2.5} \%$	${2.6}_{-1.9}^{+3.9} \%$
Mg(ClO4)2 · 6H2O	${1.5}_{-1.1}^{+2.2} \%$	${0.6}_{-0.4}^{+1.2} \%$	${1.8}_{-1.3}^{+3.0} \%$	${4.2}_{-2.9}^{+4.7} \%$
MgCl2 · nH2O	${5.6}_{-2.5}^{+3.5} \%$	${2.6}_{-1.2}^{+1.5} \%$	${6.5}_{-3.0}^{+4.8} \%$	${10.2}_{-4.6}^{+7.4} \%$
MgSO4	${2.1}_{-1.6}^{+3.5} \%$	${1.1}_{-0.9}^{+1.8} \%$	${1.1}_{-0.8}^{+1.5} \%$	${0.9}_{-0.7}^{+1.4} \%$
MgSO4 . nH2O	${4.6}_{-2.4}^{+3.5} \%$	${2.7}_{-1.5}^{+2.1} \%$	${2.3}_{-1.3}^{+2.1} \%$	${2.3}_{-1.3}^{+2.0} \%$
Na2Mg(SO4)2 · 4H2O	${10.2}_{-4.5}^{+6.0} \%$	${13.3}_{-5.1}^{+5.1} \%$	${4.7}_{-2.2}^{+3.1} \%$	${5.8}_{-2.3}^{+3.4} \%$
Na2SO4 · 10H2O	${5.7}_{-4.0}^{+6.3} \%$	${18.4}_{-4.0}^{+3.9} \%$	${2.6}_{-2.0}^{+3.9} \%$	${2.3}_{-1.7}^{+2.8} \%$
NaClO4	${0.8}_{-0.6}^{+1.5} \%$	${0.3}_{-0.3}^{+0.5} \%$	${0.7}_{-0.5}^{+1.1} \%$	${0.5}_{-0.4}^{+0.8} \%$
NaClO4 · 2H2O	${1.3}_{-0.9}^{+1.9} \%$	${0.5}_{-0.4}^{+0.8} \%$	${1.5}_{-1.2}^{+2.4} \%$	${1.9}_{-1.4}^{+3.0} \%$
NaCl	${1.0}_{-0.7}^{+1.6} \%$	${0.4}_{-0.3}^{+0.6} \%$	${0.7}_{-0.6}^{+1.3} \%$	${0.6}_{-0.4}^{+0.9} \%$
All ices	${33.0}_{-5.3}^{+5.8} \%$	${30.4}_{-3.7}^{+3.6} \%$	${22.4}_{-5.5}^{+5.2} \%$	${17.5}_{-4.0}^{+4.2} \%$
All acids	${27.1}_{-7.5}^{+6.6} \%$	${27.1}_{-5.1}^{+5.1} \%$	${49.1}_{-6.0}^{+5.1} \%$	${46.1}_{-4.6}^{+3.9} \%$
All salts	${39.7}_{-6.7}^{+7.7} \%$	${42.2}_{-5.2}^{+5.5} \%$	${29.0}_{-5.0}^{+5.2} \%$	${36.6}_{-4.1}^{+4.2} \%$

Download table as:  ASCIITypeset image