Compositional Mapping of Europa using MCMC Modelling of Near-IR VLT/SPHERE and Galileo/NIMS Observations

We present maps of surface composition of Europa's anti-jovian hemisphere acquired using high spatial resolution IFU multi-spectral data from the SPHERE instrument on the Very Large Telescope (0.95 to 1.65$\mu$m) and the NIMS instrument on the Galileo orbiter (0.7 to 5.2$\mu$m). Spectral modelling was performed using a Markov Chain Monte Carlo method to estimate endmember abundances and to quantify their associated uncertainties. Modelling results support the leading-trailing hemisphere difference in hydrated sulphuric acid abundances caused by exogenic plasma bombardment. Water ice grains are found to be in the 100$\mu$m to 1mm range, with larger grains present on the trailing hemisphere, consistent with radiation driven sputtering destroying smaller grains. Modelling best estimates suggest a mixture of sulphate and chlorinated salts, although uncertainties derived from the MCMC modelling suggest that it is difficult to confidently detect individual salt abundances with low spectral resolution spectra from SPHERE and NIMS. The high spatial resolution offered by SPHERE allows the small scale spatial distribution (<150km) of potential species to be mapped, including ground-based detection of lineae and impact features. This could be used in combination with other higher spectral resolution observations to confirm the presence of these species.


INTRODUCTION
The Galilean moons -Io, Europa, Ganymede and Callisto -are Jupiter's four largest moons. The moons collectively form a solar system in miniature around Jupiter, with environments ranging from the volcanic and rocky Io to the icy Europa, Ganymede and Callisto.
Europa, the second Galilean moon from Jupiter, is the smallest of the four moons. It has a core of silicate rock with an outer crust of liquid water and water ice which may only be ∼20 km thick (Howell 2021). Europa is tidally locked, meaning it always presents the same face to Jupiter, and tidal heating is sufficient to maintain a liquid subsurface ocean between the surface and silicate interior Schubert et al. 2004). The presence of a liquid (i.e. current-carrying) subsurface ocean is supported by induced magnetic field measurements by Galileo (Kivelson et al. 2000). The subsurface ocean and its direct contact with Europa's silicate interior makes Europa one of the most likely candidates in the solar system to be able to support habitable conditions (Chyba 2000;Marion et al. 2003).
Europa is geologically active, with transient cryovolcanic plumes of water vapour tentatively detected in Hubble Space Telescope observations (Roth et al. 2014). It has a very smooth surface, with few impact craters, implying that the surface is geologically very young, with an average age of ∼50 Myr. This young surface age implies that recent cryovolcanic resurfacing events could leave detectable signatures from the subsurface ocean on Europa's surface, potentially allowing infrared spectroscopy to identify species that may be present in the ocean, albeit chemically processed via exposure on the surface (Pappalardo et al. 1999;Greeley et al. 2004).
Europa's surface is mainly composed of water ice, and is covered with a series of intersecting linear features or 'lineae', the largest of which are over 1000 km long. These lineae are thought to be caused by tidal stresses on Europa's surface which opens fissures, exposing the warmer ice layers beneath . 'Chaos' terrains are composed of irregular polygonal blocks of old surface material that are set in younger ice. These are also thought to be caused by endogenic processes, with some form of melting or softening of the ice crust, though the exact mechanism is not fully understood Figueredo & Greeley 2004).
Much of the current understanding of Europa's surface composition comes from infrared spectroscopic observations. Observed reflectance spectra can be compared to reference spectra measured in laboratories on Earth to identify the cause of distinctive absorption features in Europa's spectrum. The most comprehensive study of Europa comes from the Galileo orbiter mission that orbited Jupiter from 1995 to 2003 with repeated flybys of the Galilean satellites. Galileo's Near-Infrared Mapping Spectrometer, NIMS, (0.7 to 5.2 µm, R = λ/∆λ ∼ 60) (Carlson et al. 1992) confirmed a surface dominated by water ice although with significant amounts non-ice contamination. Many other spacecraft have taken observations during flybys of the Jovian system, including the New Horizons flyby in 2007 which allowed observations of Europa with the LEISA (Reuter et al. 2008) instrument (1.25 to 2.5 µm, R = 240) (Grundy et al. 2007).
One of the key areas of study is the composition of the non-ice material on Europa's surface causing distortions to the detected water-ice absorption bands. The largest component of the non-ice material is hydrated sulphuric acid, which can account for much of the contamination on Europa's trailing hemisphere. Observations by NIMS show a strong correlation between the sulphuric acid distribution and the 'bullseye' distribution jovian plasma bombardment which is centred on Europa's trailing apex (270°W, 0°N) suggesting an exogenic origin for the sulphuric acid (Carlson et al. 2005;Dalton III et al. 2013). This jovian plasma contains a high abundance of sulphur ions originating from Io's volcanos, which combines with the H 2 O present on Europa's surface to produce the detected H 2 SO 4 · nH 2 O (Carlson et al. 1999).
In addition, a number of hydrated salts have been proposed to explain the remainder of the non-ice contamination on Europa's surface, with this non-ice material likely consisting of a combination of hydrated sulphuric acid and salts (Carlson et al. 2009). McCord et al. (1998), McCord et al. (1999) and McCord et al. (2002) suggested that some mixture of magnesium and sodium sulphates and potentially sodium carbonates could explain the salts present in the NIMS spectra of Europa's trailing hemisphere. Salts were found to be concentrated in young lineae and chaos terrain, while the mixture of different salts appeared constant across the observed area, suggesting these salts may originate in the subsurface ocean. Potential resurfacing mechanisms include cryovolcanic plumes, effusive flow from cracks in the crust and lag deposits from tectonic heating (Carlson et al. 2009).
More recent advances in telescope optics have enabled studies using ground-based observatories to map compositional contrasts on Europa using at higher spectral resolution than Galileo/NIMS (R ∼ 60). Brown & Hand (2013) used Keck/OSIRIS (1.4 to 1.8 µm and 1.956 to 2.381 µm, R ∼ 2000, ∼100 km/px) to identify a small 2.07 µm absorption feature on Europa's trailing hemisphere caused by magnesium sulphate salts. The apparent correlation of this feature with exogenic radiation products suggested that magnesium sulphate is itself a radiation product, rather than a completely endogenic constituent of Europa's sub-surface ocean. Observations with VLT/SINFONI (1.452 to 2.447 µm, R = 1500, ∼70 km/px) (Ligier et al. 2016) however suggested that magnesium chlorinated salts provide better fits to the overall spectrum than sulphate salts. These salts were found to be correlated with geological units, suggesting they have an endogenous origin from the sub-surface ocean. Trumbo et al. (2019) found that the 450 nm sodium chloride absorption is strongly correlated with chaos terrain on Europa's leading hemisphere, again suggesting an endogenous origin with chlorinated salts in the sub-surface ocean.
The Juno mission has enabled spectroscopic observations with the JIRAM spectrometer (2 to 5 µm) from within the Jupiter system. Studies into >2 µm water ice absorption bands suggest a mix of amorphous and crystalline ice with grain sizes ranging from tens to hundreds of microns (Filacchione et al. 2019;Mishra et al. 2021) and the temperature dependent reflectance peak around 3.6 µm suggests a maximum ice temperature of 132 K (Filacchione et al. 2019).
None of the previous ground-based studies and very few of the NIMS observations cover the J band (<1.4 µm) spectral range that is covered by SPHERE. This wavelength range covers spectral features such as water ice absorptions around 1.1 µm and 1.3 µm and various hydrated salt absorption bands around 1.2 µm (see section 4.3). For ground-based observations, there is also a trade-off between spectral and spatial resolution, with the previous studies using instruments which offer higher spectral resolution at the expense of spatial resolution. SPHERE (R ∼ 30, ∼25 km/px), on the other hand, offers a much higher spatial resolution, complementing these previous studies.
The icy Galilean moons, (Europa, Ganymede and Callisto), are due to be studied in the coming decade by ESA's Jupiter Icy Moons Explorer (JUICE) and NASA's Europa Clipper. JUICE will carry out flybys of Callisto and Europa, and will then orbit Ganymede, providing detailed global mapping data of the entire moon (Grasset et al. 2013). Europa Clipper will also carry out flybys of Europa, Ganymede and Callisto, and will have a specific focus on Europa, performing 45 flybys of the moon to produce effectively global data coverage. The modelling of SPHERE and NIMS observations in this paper is a precursor to the analysis that will be possible with the MISE and MAJIS spectrometers on Europa Clipper and JUICE respectively.
In this paper, we will discuss our reduction and analysis of near-infrared reflected sunlight observations of Europa using the SPHERE instrument on the groundbased Very Large Telescope, VLT/SPHERE (Beuzit et al. 2019). These observations will be used to analyse the composition of the Europa's surface, in relation to the physical and chemical processes shaping these worlds. We will demonstrate that ground-based observations are capable of reproducing and extending those from visiting spacecraft, thereby providing key support for future missions.  (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 highspatial-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 December 2014, as summarised 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) sub-systems of the SPHERE instrument. This allows simultaneous spectroscopic observations covering 0.95 to 1.65 µm (J and H bands) and dual band imaging at 2.1 µm and 2.251 µm (K band).
The Integral Field Spectrograph (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/px, corresponding to ∼25 km/px at Jupiter. Accounting for diffraction, this allows features ∼150 km across to be resolved. The magnitude of the diffraction is calculated by measuring the full width at half maximum of the point spread of the calibration star, observed through VLT/SPHERE (≈6 px at 0.95 µm).
The IFS operates by using a lenslet array that focusses 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 maximised depending on the spectral features of individual targets. Our data uses the DB_K12 filter set that have transmissions centred at 2.1 µm for the K1 filter and 2.251 µm for the K2 filter, with filter widths of 0.051 µm and 0.055 µm respectively, meaning IRDIS is sensitive to water ice absorption at 2 µm.

Data reduction
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 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:  Table 2). The shaded orange region shows the area of the 'E6' observation where NIMS has full spectral coverage of the SPHERE wavelength range, which is used for full spectral modelling. Regions with extreme emission and incidence angles (> 75°) have large residual photometric errors, so are not used in this study or shown on this map. The background visible light reference image of Europa is from Becker (2013). 0°W is the sub-jovian longitude and 180°W is the anti-jovian longitude. 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 consists 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: 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.

Image cleaning
A strong regular banding pattern (see Figure 2a) 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 cross-talk 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.
To remove this pattern, its specific frequencies were identified by analysing 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, then converting back to image space using where F is the Fourier transform (image domain to frequency domain) and 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.
The destriping routine removes the pattern for wavelengths where it occurs, and 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 wavelengthdependent systematic errors.
Images are 'despiked' by calculating the mean and variance of all pixels within a 2 px 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 whilst 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.

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 (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 centred on the calibration star that have radii calculated from the size of the Airy disc. The star aperture is cho- Figure 2. Image cleaning process for IFS data. The top images show observations of Europa at 1.37 µm, the wavelength where the striping pattern is the worst. The bottom images show the observations' corresponding Fourier Transform magnitudes on a logarithmic scale. The frequencies that cause the striping pattern in the raw image (the bright dots away from the centre of the Fourier transform) are set to zero to remove the artificial pattern. The image is then despiked to remove any extreme bright or dark pixels, which are replaced by interpolated values. sen so that the first airy diffraction maximum is contained within the aperture while still minimising the background flux within the aperture. The background flux is summed from a 10 px 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).

Mapping
In order to map the observations, the exact location, size and orientation of the moon's disc in the observed image must be known. The size of the moon (i.e. its pixel radius) and its orientation (i.e. the angle between its north pole and the image vertical) can be calculated from ephemeris data and the known plate scale and orientation of the detector.
The pointing information for the telescope (RA and Dec) are not accurate enough to use for determining the location of the moon's disc in the image for mapping purposes, as accurate mapping and photometric correction requires the disc location to be known to 1 px for the best results. Therefore, the disc location is calculated from the image itself. This is done by applying a threshold to the image (such that all dark pixels are 0 and bright pixels are 1), and then calculating the centroid of this image (its centre of brightness). The threshold step ensures that any variations in brightness across the moon do not significantly bias the results, as all pixels within the disc of the moon will have the same value. This method identifies the centre of the disc within a few pixels.
The calculated disc location is then corrected to reduce errors caused by the moon's phase. This is done by creating a photometric model of the target and calculating its apparent position through the same routine. The difference between the actual centre of the photometric model and its calculated centre is then used to correct the calculated centre of the observation. This step ensures the dark limb of a target observed at nonzero phase angles does not affect the calculation of the its centre. For our Europa data, observed at a phase angle of ∼10°, this step corrects the disc location by ∼2 px.
Once the position of the disc is known, the latitude, longitude, incidence angle and emission angle are calculated for each pixel using ephemeris data for the target and the rotation information for the telescope. The incidence and emission angles are used to photometrically correct the image (section 2.2.4). Finally, the image is transformed to an equirectangular map projection using the calculated coordinates for each pixel. Each wavelength of each individual observation is mapped separately to ensure there were no errors introduced by Europa's rotation between observations or any wavelength dependent drifts.

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 towards the edge of the disc (see Figure 4a). The photometric effect must be accounted for to correct the reflectance for areas of an observed disc away from the centre. The correction takes the form Corrected image = Image Photometric model (3) where the photometric model gives the ratio of the brightness of the disc 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 was both the sub-solar point and the sub-observer point (i.e. where the incidence and emission angles are zero).
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 towards the edge of the disc, overcompensating for the limb darkening and producing unphysical reflectance values (> 1) (see Figure 4b).
Therefore, we used the Oren-Nayar reflectance model (Oren & Nayar 1994), which modifies on 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 the surface is composed of a series of V-shaped cavities, with facets that are themselves Lambertian scatterers. The cavities are modelled to have a Gaussian distribution of slopes, with a mean slope of µ = 0 and a standard deviation of σ that parametrises the surface roughness. The Oren-Nayar model reduces to standard Lambertian scattering for σ = 0.
For a surface roughness σ and albedo ρ, the modelled Oren-Nayar surface brightness is where (i, φ i ) and (e, φ e ) are the polar and azimuth angles of the incident and emitted rays respectively. L 1 is the contribution from single scattering and L 2 is the contribution from multiple scattering with The single scattering coefficients C n are dependent on the surface roughness σ and the viewing geometry, C 3 = 0.125 σ 2 σ 2 + 0.09 4αβ π 2 2 (12)  Incorporating roughness increases the modelled brightness at the edge of the disc 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 towards centre of Europa's disk. Large values of σ produce undercorrected observations where the limb is too dark (e.g. Figure 4a) and small values of σ produce overcorrected observations where the limb is too bright (e.g. Figure 4b). σ = 33° (Figure 4c) was found to produce the minimum brightness variation 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°.
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 modelled 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 disc model before being used to correct the image. The radius of the airy disc is proportional to the wavelength of the specific image being corrected, with the radius set to r airy = 10 px for 1.65 µm, as measured from observations of the calibration star (see Figure 4d for 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 disc. This is higher than previous studies using Lambertian models that typically extend to 50°-60° (Brown & Hand 2013;Grundy et al. 2007;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 disc 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.

Galileo/NIMS
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 Near-Infrared Mapping Spectrometer (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/px 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/px).
The NIMS instrument used different detectors to cover different wavelength ranges, one of which (detector 3, covering 0.99 to 1.26 µm) failed early in the Galileo mission, during Galileo's 4th orbit of Jupiter, so there are no Europa observations with full spectral coverage from 1997 onwards. 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 datasets given in Table 2 were used in this study. These datasets 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 datasets 17EN-GLOBAL01A and 17ENGLOBAL02A were taken after the failure of NIMS detector 3, so lack the 0.99 to 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 datasets using the same mapping and photometric correction code as the SPHERE dataset.
For our analysis, we selected NIMS observations taken at a low phase angle with a relatively high spatial resolution ( 100 km/px) and similar spatial coverage to our SPHERE observation (spatial resolution ∼150 km), allowing direct comparison between the two datasets.

WATER ICE ABSORPTION BANDS
The most prominent features in Europa's near-IR spectra are a series of absorption bands caused by water ice ). These bands can be used to provide a simple qualitative indication of the water ice distribution on the surface of icy bodies like Europa, with stronger water ice absorption corresponding to higher water ice abundance (Hapke 1993). As Europa's crust is mainly composed of water ice , the inverse of the water ice distribution provides an indication of the spatial distribution of non-ice contaminants on Europa's surface. Figure 6 shows the strength of the mapped 1.5 µm band in the SPHERE IFS and Galileo/NIMS datasets using the 1.50/1.36 µm spectral ratio. The use of the ratio (rather than simply measuring the reflectance at 1.50 µm) means that variations in this relative band depth are mainly caused by variations in water ice abundance rather than varying grain size or ice crystallinity (see graph in Figure 6). Figure 7 shows the strength of the 2 µm water ice band, as directly imaged by the two filters of the IRDIS instrument.
Both bands and instruments show consistent distributions, with the only major differences due to the higher spatial resolution of the NIMS dataset. Water ice appears more abundant at high latitudes and there is significant contamination of Europa's trailing hemisphere, particularly in Dyfed Regio. This contamination is centred on the trailing apex and is consistent with the bullseye distribution of exogenous plasma bombardment increasing the non-ice fraction of the trailing hemisphere at low latitudes (Carlson et al. 2005).
The shape of the contaminated regions follow the structures of geological units, such as the clear outlines of Dyfed Regio and Powys Regio. The effect of Manannán and Pwyll impact craters are visible due to the impacts exposing less contaminated ice leading to an increased water ice signature in and around the craters. Signatures of some of Europa's large lineae are also visible in both the SPHERE and NIMS data around Cadmus Linea, as an arc of lower ice content material from the north of Dyfed Regio. This is the first time that the lineae have been directly observed using a ground-based telescope, demonstrating the significant improvement in spatial resolution that SPHERE enables. Figure 6. Normalised 1.50/1.36 µm spectral ratio for SPHERE and NIMS observations. The NIMS data are a combination of the datasets in Table 2; in regions where the observations overlap, the median reflectance value is used. Blue areas in the map indicate stronger absorption at 1.50 µm, generally implying a higher water ice abundance and red areas indicate reduced water ice abundance. The graph on the left shows reflectance spectra for a range of water ice grain sizes (darker lines indicate smaller ice grains, see Figure 8) normalised to unity at 1.36 µm. The choice of ratio wavelengths measures the full relative depth of the 1.5 µm absorption and produces a map ratio which is mainly affected by ice abundance (rather than grain size). The Oren-Nayar photometric correction allows our mapping to extend to higher emission angles than typical for ground-based observations. Here, we assume data below an emission angle of e = 75°is useful, as above this angle the spatial resolution of the mapped data degrades significantly. This is an improvement on previous studies, which typically discard data above e = 50°to 60°, allowing us to usefully observe and map higher latitudes with ground-based observations.

Linear unmixing
The mapped spectral cubes are analysed by fitting to laboratory spectra from reference cryogenic libraries (described in section 4.3). Our fitting routine uses linear spectral modelling to fit the observed spectra, where the modelled spectrum, M , is 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 modelled spectrum, and are subject to the constraint which ensures the fractional abundance of the different endmembers sums to 100% and 0 ≤ w i ≤ 1 which ensures 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 non-linear scattering that occurs, and may therefore be responsible for some small residual errors after fitting (Ligier et al. 2016;Shirley et al. 2016).

Markov Chain Monte Carlo modelling
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 theorem (Bayes 1763) states that the posterior probability of a set of parameter values w given an observation o can be calculated as where P ( o| w) gives the likelihood of measuring the observation o given the parameters w. P ( w) 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 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 modelling 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 towards and then remains in a region of parameter space which 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 modelling, MCMC enables the simulation of reflectance spectra, and the calculation of the posterior probability distribution of endmembers abundances being 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 modelling 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 confi-dence 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 modelling 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 modelling 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 which would imply unphysical abundance values (i.e. less than 0% or greater than 100%) are rejected by the MCMC process.
We assume all possible combinations of abundance values are equally likely so 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 one used in Lapotre et al. (2017). The sum-to-one 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 then setting the final endmember's abundance to be w N = 1 − N −1 i w i . This is mathematically equivalent to rejecting all solutions where the abundances do not sum to one, but 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 optimised fitting routine to perform accurate fits across the whole observed disc within a reasonable timescale: 1. Perform a simple linear fit to calculate optimized parameter values, w LF i .
2. Generate initial positions, pos0, for MCMC chains as a random Gaussian ball centred on w LF i with a standard deviation of initialisation_sd. These initial chain positions are corrected to ensure 0 ≤ w i ≤ 1 and i w i = 1. Initialising around the linear fit values significantly decreases the required chain length, as the chain is initialised 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 initialisation. This 'burns in' the chain to ensure the chain position is relatively independent of the initialisation.
4. Continue running emcee in a series of runs with run_steps steps. Each run is initialised 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 it remains accurate and the choice of different parameters do not influence the final calculated values. For example, initialisations with completely random values for pos0 ultimately converge to identical values as initialisations using linear fit values, but 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 Table 4. 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 which 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%.

Spectral library
The full list of endmembers in our spectral library are given in Table 3 and selected spectra are shown in Figure 8. These cryogenic reference spectra include hydrated sulphuric acid (Carlson et al. 1999) and a variety of hydrated salts (Hanley et al. 2014;Dalton & Pitman 2012;Dalton III 2007). The sulphuric acid spectra (Carlson et al. 1999) do not cover wavelengths below 1 µm, so this limits our modelling to cover the 1 to 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).
Our water ice reflectance spectra are calculated from measured refractive indices (Grundy & Schmitt 1998) using the Hapke bidirectional reflectance model (Hapke 1993). Modelling 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 extreme 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 modelled 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 modelling 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 modelling (Hapke 1993).

REGIONS OF INTEREST
Initial modelling was performed on regions of interest in Powys Regio (170°W, 5°N) and Dyfed Regio (240°W, 30°N). These locations both have high non ice content, as indicated by the weakness of the 1.5 µm water ice absorption in Figure 6. These also sample both the leading (Powys Regio) and trailing (Dyfed Regio) hemispheres, so provide a useful comparison between the non-ice material in the two hemispheres and their associated plasma and radiation environments. In addition to comparing the composition of the two locations, we can compare the SPHERE and NIMS fits for the same location to understand the limitations of the datasets.
The MCMC fit results for both locations are summarised in Figure 9, where the violin plots on the right hand side show the fitted abundance distributions for each endmember and family of endmembers. Both the SPHERE and NIMS datasets show consistent results, with a roughly even mixture of ices, acids and salts. As expected, Dyfed Regio on the trailing hemisphere has a higher acid abundance, and there are differences in the salt mixtures for each location (discussed in more detail in section 6.3).
It is notable that the uncertainty on the fitted abundances is larger for some endmembers in the SPHERE dataset. For example SPHERE is less able to discriminate between different water ice grain sizes in Powys Regio, and is less able to rule out the detections of individual salts. This demonstrates the utility of the additional spectral range covered by NIMS in helping to lift degeneracies between different endmembers and increase the confidence of their detections. However, the best estimate values are consistent for both datasets, suggesting that the missing spectral range for SPHERE is not leading to any spurious results.
The MCMC results can also be summarised using corner plots that show the relationship between posterior distributions of pairs of endmembers from the modelling of an individual spectrum. Figure 10 shows the corner plot result for the SPHERE Powys Regio fit. Generally, most of the 2D-histograms in Figure 10 are roughly circular, implying relatively little relationship between the individual endmember abundances distributions. However some of the 2D-histograms, particularly between water ice grains (top left), are clearly skewed, implying correlation between the endmember abundance distributions. For the case of small and large water ice grains, this negative covariance implies that much of the uncertainty in individual endmember abundance is due to the slight degeneracy between grain sizes. This can be seen in the violin plots in Figure 9 where there is relatively large uncertainty on some individual ice grain size endmembers, but these cancel out, leaving a smaller uncertainty on the total water ice abundance.
6. COMPOSITIONAL MAPS MCMC fitting was performed for each observed spectrum for both the SPHERE and NIMS datasets. This produces posterior abundance distributions for each location that were sampled to create maps of best estimate abundances for different species, showing their spatial distributions.

Hydrated sulphuric acid
The hydrated sulphuric acid distribution shown in Figure 11 is mainly concentrated towards the trailing Table 3. Spectral library hemisphere, with lower abundances towards 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 bullseye distribution centred 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 sulphuric acid hydrate present on Europa's surface.
The peak abundance is 68% for SPHERE (1-σ range 60% to 75%) and 65% for NIMS (1-σ range 60% to 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 datasets shown in Figure 11 only cover some of the eastern part of the trailing hemisphere, so the values are not directly comparable.

Water ice
The modelled 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 towards the trailing apex. The water ice distribution is highly anticorrelated with the sulphuric acid distribution, with the combined hydrated sulphuric acid and ice abundance ∼ 80% when averaged over the observed area. This anti-correlation 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.
As found in previous studies (Dalton III et al. 2012;Ligier et al. 2016), the majority of the surface is dominated by 100 µm to 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 localised 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.

Salts
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 III 2007;Shirley et al. 2010), sulphates appear to be the main salts present on Europa's surface (see Figure 13). Magnesium bearing salts (MgSO 4 · nH 2 O, MgSO 4 brine and Na 2 Mg(SO 4 ) 2 · 4 H 2 O) appear correlated with dark ter- rain, particularly Powys Regio. As shown in the lower panels of Figure 13, the lower bounds on the abundances of magnesium sulphates are all very low and close to zero. This suggests that the uncertainties and degeneracies between the salt spectra mean it is not possible to positively identify any individual magnesium sulphate salts with the SPHERE and NIMS spectral resolutions.
Mirabilite (Na 2 SO 4 · 10 H 2 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 modelled salts, and the corresponding non-zero lower bounds suggest that it is likely to be present on the leading hemisphere. The modelled 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 sulphate salts and all have lower bounds close to zero, implying it is difficult to positively identify individual salts, particularly with SPHERE. The spatial distribution of the chlorinated salts appears more globally uniform than the sulphate salts with slightly higher abundances around geological units.
Hydrated magnesium chlorate (Mg(ClO 3 ) 2 · 6 H 2 O), perchlorate (Mg(ClO 4 ) 2 · 6 H 2 O) and sodium perchlorate (NaClO 4 ) 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  Table 5. The shaded 2D-histograms give the posterior abundance distribution with brighter colours indicating a higher density. In low density regions towards the edge of the distribution, individual points of the distribution are plotted. For clarity, some endmembers have been combined. Circular distributions imply no strong correlation between the abundance distributions of the two endmembers whereas skewed distributions imply a correlation between the abundances of those endmembers. Figure 11. Hydrated sulphuric acid spatial distribution for the SPHERE and NIMS datasets. The values shown here are the best estimates that are calculated as the median of the posterior abundance distribution for each location. The acid abundance is highest towards the trailing apex (270°W, 0°N) and lower towards the leading hemisphere and high latitudes. The slight increase in acid abundance towards the edge of the SPHERE dataset is likely to be caused by residual errors from the photometric correction. The typical uncertainty on the abundance values is ±10 percentage points for SPHERE and ±7 percentage points for NIMS. 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, so is consistent with our low abundances.
Hydrated magnesium chloride (MgCl 2 · nH 2 O) appears the most abundant chlorinated salt, with best estimate abundances up to 10% in the NIMS dataset. 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.

DISCUSSION
Compositional contrasts on Europa's surface appear to be primarily driven by exogenic plasma bombardment and the resulting longitudinal gradients in rates of radiation driven processes. As with previous studies (Carlson et al. 2005;Dalton III et al. 2012;Ligier et al. 2016), the leading-trailing contrast in acid abundance strongly suggests an exogenic origin for the sulphuric acid as the acid abundances appear correlated with sulphur ion bombardment. Figure 12. Water ice spatial distributions for SPHERE (top) and NIMS (bottom). The total water ice abundance (left) is lowest towards the trailing apex and highest at high latitudes where abundances are > 50%. The grain size column (right) shows the ice size endmember with the highest individual abundance. The ice grain size varies with longitude, with larger (300 µm to 1 mm) grains more abundant on the trailing hemisphere and smaller grains (100 µm to 300 µm) more abundant on the leading hemisphere.
Similarly, the leading-trailing hemisphere dichotomy of ice grain sizes suggests that exogenic processes are likely to be responsible for the contrast in grain sizes. As described in Clark et al. (1983), more intense radiation on the trailing hemisphere will lead to a faster sputtering rate for water ice grains compared to the leading hemisphere. This sputtering erodes grains at a rate proportional to the surface area of the grain, preferentially destroying smaller grains due to their higher surface area to volume ratio. Therefore, the leading hemisphere's lower sputtering rate will increase the lifetime of smaller grains, thus reducing the average grain size on the leading hemisphere .
In addition to the longitudinal variation in water ice grain size, the SPHERE data also show some latitudinal contrast, with larger grains extending further eastwards in the northern hemisphere. The cause of this latitudinal variation is unclear, especially given that it is not repeated in the southern hemisphere, suggesting that thermal effects are unlikely to be the cause as these would be expected to be symmetric around the equator. Therefore, this region of larger ice grains may be caused by some endogenic process, potentially related to the age of the surface.
Cratered areas and ridged plains, particularly at high latitudes, show high water ice abundance. Chaos areas and bands on the other hand appear more contaminated, especially Dyfed Regio which is located near the trailing apex and therefore receives the most intense exogenic plasma bombardment. Best estimate salt abundances are highest in the chaos areas in Dyfed and Powys Regio.
The uncertainties on the modelled salt abundances demonstrate how the low spectral resolution SPHERE and NIMS data are consistent with a variety of different  salt mixtures. The lack of unique distinguishing spectral features at these low spectral resolutions means that the endmember spectra are partially degenerate, meaning that the observed spectra generally cannot be used to positively identify individual salt species (see Figure 13 and Figure 14). This makes detailed analysis of the salt distributions and their potential origins difficult due to the wide variety of potential mixtures. However, the high spatial resolution offered by SPHERE and NIMS can be used to infer the likely spatial distribution the salts would have if they are present. These spatial distributions can then be used in comparison with high spectral resolution observations that have more confidently identified specific species in a general area.
The clearest trend in our data was the higher mirabilite (Na 2 SO 4 · 10 H 2 O) abundance in the leading hemisphere (highest abundance at ∼ 130°W) contrasting with the higher magnesium sulphate abundances towards and into the trailing hemisphere (highest abundance at ∼ 160°W). This is very similar to the longitudinal contrasts in hexahydrite (MgSO 4 · 6 H 2 O) and mirabilite abundances found in Shirley et al. (2010). Mc-Cord et al. (2001) found that magnesium sulphates are more stable than sodium sulphates against thermal dehydration and radiolysis, meaning that higher radiation on the trailing hemisphere will destroy mirabilite at a faster rate. Therefore, as with acids and ices, the exogenic radiation environment may explain these observed longitudinal compositional contrasts.
Magnesium ions have not been detected in the Io plasma torus (Carlson et al. 2009), so the most plausible source for Mg ions is Europa's sub-surface ocean. McCord et al. (2002) suggested that magnesium sulphates may originate directly from the sub-surface ocean whereas Brown & Hand (2013) suggested that magnesium sulphate is a radiation product of exogenic sulphur ions and magnesium salts already present on Europa's surface. Our observed magnesium sulphate distributions ( Figure 13) appear broadly consistent with both hypotheses, with magnesium sulphates appearing mainly in geological units in regions where Brown & Hand (2013) detected magnesium sulphate. SPHERE observations with increased longitude coverage would allow more detailed investigation of the radiation product hypothesis, although the detection of a distinctive magnesium sulphate absorption feature mainly on the trailing hemisphere in high spectral resolution observations by Brown & Hand (2013) appears compelling. This 2.07 µm absorption feature is very narrow, so is not detectable with NIMS' spectral resolution.
Chlorine bearing salts, particularly magnesium chloride, appear correlated with geological units such as lin-eae and darker areas ( Figure 14). The magnesium chloride abundances are low ( 10%), but appear consistent with the abundances and spatial distributions found in Ligier et al. (2016). The strong correlation with geological units suggests an endogenic origin for these salts. If chlorine ions are present in Europa's sub-surface ocean, magnesium chloride would form directly in the ocean. Therefore, this observed magnesium chloride may have originated directly from the ocean (Ligier et al. 2016) and would provide the magnesium source for any radiation produced magnesium sulphate (Brown & Hand 2013).
Magnesium perchlorate appears correlated with both the magnesium chloride distribution and the exogenic plasma bombardment centred on the trailing apex (Figure 14). This is consistent with the endogenic magnesium chloride being converted to magnesium perchlorate under the high radiation ion bombardment in the trailing hemisphere around Dyfed Regio (Carrier & Kounaves 2015;Ligier et al. 2016).
As demonstrated by the uncertainties in salt abundances produced by the MCMC fitting routine, broader wavelength coverage and higher spectral resolution data will be necessary to positively identify or rule out the presence of many individual salts on Europa's surface. Future Juno observations of Europa will be able to use JIRAM's high spectral resolution and the very high spatial resolution offered by Juno's position in the Jupiter system to study and constrain Europa's non-ice surface composition at small scales (Mishra et al. 2021). The next generation of instruments such as ELT/HARMONI and JWST/NIRSPEC will enable spectroscopy to identify narrow features, such as the 2.07 µm magnesium sulphate absorption in Brown & Hand (2013), whilst also acquiring high spatial resolution mapping to identify small-scale compositional contrasts.

CONCLUSION
The near-infrared spectral observations with SPHERE and NIMS have enabled compositional mapping of Europa's anti-jovian hemisphere. The two instruments' datasets appear generally consistent, demonstrating the utility of SPHERE's low spectral resolution and high spatial resolution ground-based spectroscopy for solar system targets.
Our reduction pipeline has enabled the use of SPHERE to observe extended targets by removing the banding pattern likely caused by crosstalk between different lenslets. The use of the Oren-Nayar reflectance model (Oren & Nayar 1994), combined with a diffraction model, allows photometric correction of observations to emission angles > 70°. This allows mapping to extend to higher latitudes than previous studies, observing > 90% of the observed disc. The limiting factor is generally now the low spatial resolution caused by the extreme viewing angles at the edge of the disc.
The use of MCMC modelling for spectral fits enables detailed uncertainty estimation on final fitted abundances. This is particularly useful for the relatively low-abundance hydrated salts, where the MCMC uncertainties quantitatively show that confident detection of individual salt endmembers is generally not possible with these datasets due to the significant degeneracies between hydrated salt spectra. High spatial resolution observations, such as with SPHERE, can be used to complement other observations by spatially mapping species which have been more confidently detected with higher spectral resolution, lower spatial resolution observations. Future studies may use further statistical techniques such as the Beysian information criterion (Schwarz 1978) to help to select between different potential models of Europa's surface composition.
The identified spatial distributions of species appears consistent with previous studies of Europa, with high acid abundance and larger ice grains on the trailing hemisphere and a variety of potential hydrated salts. Future studies will help to further constrain Europa's sur-face composition by expanding the spatial and spectral coverage and resolution of near-infrared observations. Higher spectral resolution observations (e.g. JWST and ELT/HARMONI) and laboratory reference spectra measurements will help to constrain salt abundances by identifying narrow characteristic features in the spectra, while high spatial resolution spacecraft observations will enable accurate geolocation of compositional features. Table 5. Abundances values for case studies shown in Figure 9.