Leveraging the Gravity Field Spectrum for Icy Satellite Interior Structure Determination: The Case of Europa with the Europa Clipper Mission

An inverse problem methodology requires the development of a forward model. Here we describe a general structure for the interior of an icy body and detail the computation of the associated gravity potential. The modeling approach employed herein draws inspiration from the works of Koh et al. (2022) and Pauer et al. (2010). Unlike the aforementioned works, our focus lies on the gravitational potential rather than the surface acceleration produced by gravity anomalies. This choice allows for a direct comparison with the quantities typically solved for in OD solutions of deep space probes. For this reason and for the sake of self-consistency, this manuscript describes in detail our methodology to construct the interior model.

In this work we assume a five-layer model for the interior of an icy body, which is deemed sufficiently comprehensive for describing its internal differentiation. As shown in Figure 1, we account for a silicate interior and an outer hydrosphere. The silicate interior encompasses a core, a mantle, and a crust. Following Koh et al. (2022), the terms "crust" and "mantle" here refer to chemical distinctions within the silicate layer. The hydrosphere is composed of an ocean and an icy shell. Only the core is assumed spherically symmetric, while all the other layers display thickness variations, due to topography, superimposed on a spherical shape.

Zoom In Zoom Out Reset image size

The nonspherically symmetric gravitational potential is generated by the irregular interface (topography) between layers characterized by a significant density contrast. Notably, these contrasts exist between the silicate crust and the ocean, as well as between the icy shell and the vacuum of space. We note that a third contribution, associated with the degree-2 shape, is the dominant factor on the total gravitational potential. The degree-2 gravity contribution expected for a tidally deformed body is added manually, rather than being calculated from the shape, as further discussed toward the end of this section.

To generate the gravity field due to the seafloor topography, we generate a random set of spherical harmonics coefficients (C_lm , S_lm ) and rescale them to different levels of maximum expected topographic range.

We generate a random set of spherical harmonics coefficients following a prescribed power rule. In particular,

$$\begin{eqnarray}&&{P}_{l}=\displaystyle \frac{k}{{l}^{\beta }},\end{eqnarray} \tag{ 1 }$$

where P is the power spectrum of the topography field, l is the degree, and k and β are prescribed coefficients.

The power spectrum is defined as

$$\begin{eqnarray}&&{P}_{l}={\sum }_{m}^{}{C}_{{lm}}^{2}+{S}_{{lm}}^{2},\end{eqnarray} \tag{ 2 }$$

where β is a free parameter of the model and k is chosen such that the maximum (in absolute value) of the topographic height equals a prescribed value, δ h. The method used for generating the random set of spherical harmonics ensures symmetry between maximum and minimum values of the topography; thus, δ h represents half of the maximum topographic range. Renaming, for clarity, the spherical harmonics coefficients describing the topography as C_lm ^h and S_lm ^h , the corresponding gravity potential coefficients are obtained from the topography as follows (Jeffreys 1976):

$$\begin{eqnarray}&&{C}_{{lm}}^{U}=\displaystyle \frac{1}{2l+1}4\pi G{\rm{\Delta }}\rho R{\left(1-\displaystyle \frac{z}{R}\right)}^{l+2}{C}_{{lm}}^{h},\end{eqnarray} \tag{ 3a }$$

$$\begin{eqnarray}&&{S}_{{lm}}^{U}=\displaystyle \frac{1}{2l+1}4\pi G{\rm{\Delta }}\rho R{\left(1-\displaystyle \frac{z}{R}\right)}^{l+2}{S}_{{lm}}^{h},\end{eqnarray} \tag{ 3b }$$

where C_lm ^U and S_lm ^U are the gravity spherical harmonics coefficients corresponding to the topography coefficients ${C}_{{lm}}^{h},{S}_{{lm}}^{h}$. G is the universal gravitational constant, and z and R are the radial distance of the gravity potential evaluation point from the considered interface and the distance from the interface to the center of the body, respectively. In this work, we are interested in the spherical harmonics coefficients of the gravity potential at the surface of Europa. We set then z in the above equations such that the evaluation point equals the radius of Europa, and R as the radius of the layer interface.

In the following we will work with adimensional potential coefficients, obtained by normalizing Equations (3a) and (3b) with GM/R.

We account for partial elastic compensation of the silicate crust by multiplying each topographic coefficient by a factor F_l , defined as

$$\begin{eqnarray}&&{F}_{l}=\displaystyle \frac{1}{1+\tfrac{{\rho }_{c}-{\rho }_{w}}{{\rho }_{m}-{\rho }_{c}}{C}_{l}},\end{eqnarray} \tag{ 4 }$$

where C_l is defined as in Turcotte et al. (1981, Equation (27)) and ρ_c , ρ_m , and ρ_w are the densities of the crust, mantle, and ocean, respectively. C_l depends on the degree of the expansion, the densities of the crust and mantle, and the rigidity of the lithosphere:

$$\begin{eqnarray}&&{C}_{l}=\displaystyle \frac{1-\tfrac{3{\rho }_{m}}{(2l+1)\overline{\rho }}}{\tfrac{\sigma \left[{l}^{3}{\left(l+1\right)}^{3}-4{l}^{2}{\left(l+1\right)}^{2}\right]+\tau [l(l+1)-2]}{l(l+1)-1-{\nu }_{c}}+1-\tfrac{3{\rho }_{m}}{(2l+1)\ \overline{\rho }}},\end{eqnarray} \tag{ 5 }$$

where

$$\begin{eqnarray*}\begin{array}{ccc}\tau & = & \frac{{E}_{c}{{T}_{E}}_{c}}{2{R}_{c}^{2}g\left({\rho }_{m}-{\rho }_{c}\right)}\\ \sigma & = & \frac{\tau }{12\left(1-{\nu }_{c}^{2}\right)}{\left(\frac{{{T}_{E}}_{c}}{{R}_{c}}\right)}^{2},\end{array}\end{eqnarray*}$$

and $\overline{\rho }$ is the average density of the body and ${\nu }_{c},{R}_{c},{{T}_{E}}_{c}$ are the Poisson ratio, radius, and elastic thickness of the crust, respectively.

The very same approach can be used for the generation of the gravity field contribution from the partially compensated icy shell topography. In this case, Equations (4) and (5) can directly be used, making sure to attribute the correct value to the quantities with subscripts c, w, m, which in this case refer to the icy shell, the outer vacuum environment, and the subsurface ocean, respectively. In this work, being interested in the retrieval method, without loss of generality and for the sake of conciseness, we use a simple model for the topographic contribution to gravity, not accounting for finite-amplitude corrections and vertical attenuation of the topography at the base of the ice layer and silicate crust. The underlying interior structure model can be modified and made more complex without impacting the following of this work, provided that the basic assumptions made here continue to hold (e.g., that the topography follows a power law: Equation (1)).

Given the linearity property of spherical harmonics, the total nonspherically symmetric gravitational potential of the body (the spherically symmetric part is described by l = 0) is given by the sum of the potentials of the silicate interior and the icy shell.

Figure 2 shows the amplitude spectra of the silicate and icy shell gravity fields obtained choosing model parameters representative of Europa (see Table 1 and accompanying discussion in Section 3.2 below).

Zoom In Zoom Out Reset image size

The amplitude spectrum is defined as

$$\begin{eqnarray}&&{A}_{l}=\sqrt{\tfrac{{\sum }_{m}^{}{{C}_{{lm}}^{U}}^{2}+{{S}_{{lm}}^{U}}^{2}}{2l+1}}.\end{eqnarray} \tag{ 6 }$$

One of the most interesting features is the degree at which the icy shell becomes dominant on the total external gravitational potential. In the case shown in Figure 2, the turning point (crossing degree l_cr) is around l_cr = 30. This is an important indication, depending on the expected resolution of the gravity field recovery, of which parameters will have a substantial influence and need to be considered to explain the spacecraft measurements. Since the degree-2 gravity field is dominated by shape and rotation and this quadrupole gravity field has been measured for most of the icy bodies of interest, we limit the generation of the gravity model from degree 3 onward and set the degree-2 field to the measured values from Voyager and Galileo. Notably, A₂ in Figure 2 is not derived from the silicate and icy shell spectra but is set to the measured value (in this case from Anderson et al. 1998).

The problem of OD involves minimizing the discrepancy between actual tracking observations of a probe and the predictions generated by a comprehensive model of its dynamic environment and measurement system. A detailed discussion on the theory of OD is outside the scope of this manuscript, and interested readers are referred to seminal texts such as Tapley et al. (2004) and Milani & Gronchi (2009).

In obtaining the OD solution, a Kaula-like regularization of the spherical harmonics coefficients is often used. The use of such regularization offers several advantages, including improved numerical stability of the solution and the ability to constrain the noise power of poorly resolved spherical harmonic coefficients. This is particularly crucial in situations where global coverage is not guaranteed, as is the case with Europa Clipper (Mazarico et al. 2023, Figure 3). However, defining a Kaula regularization requires an a priori expectation of the field power. That is, the choice of this regularization constraint affects the solution and potentially its power spectrum.

In this study, we adopt a two-step OD inversion approach. In the first step, an unconstrained estimation of the gravity field is performed. This strategy allows us to obtain a solution for the low-degree gravity field that is entirely data driven and independent of any a priori assumptions. From this initial solution, we use the well-determined lowest degrees of the gravity field to fit a power law to the spectrum. The resulting fitted power law serves as a regularization constraint in the second-step inversion of the gravity field. By employing this two-step methodology, we ensure that the higher-degree solution, which necessitates regularization, is purely data driven. It is worth noting that here we implicitly assume the power spectrum to be well approximated by a single power law in the degree range of the analysis. This is well supported by current observations of other planetary bodies (e.g., Ermakov et al. 2018, Figure 1).

Figure 3 shows a graphical example of this procedure. The gravity field derived from the first step exhibits substantial uncertainty due to the absence of regularization. In the case of Europa Clipper, we utilize degrees 3–5 to fit a power law to the estimated spectrum and infer the power of the higher degrees. As depicted in the figure, the fitted power law provides a close approximation to the true field.

Zoom In Zoom Out Reset image size

By incorporating this regularization approach into the second OD fit, we can obtain a reliable solution for higher degrees as well. We emphasize the significance of this two-step procedure, as it enables obtaining a high-quality gravity solution without relying on any prior knowledge of the body's gravity field (except the assumption of a single power-law slope).

Using the gravity field spectrum as an observable quantity in the Bayesian inversion process requires defining the associated uncertainty, in the form of a covariance matrix.

Typically, in planetary gravity studies, the uncertainty spectrum presented alongside the amplitude spectrum is the rms value (per degree) of the retrieved uncertainty on the spherical harmonics expansion coefficients:

$$\begin{eqnarray}&&{\sigma }_{{A}_{l}}=\sqrt{\tfrac{{\sum }_{m=0}^{l+1}({\sigma }_{{C}_{{lm}}}^{2}+{\sigma }_{{S}_{{lm}}}^{2})}{2l+1}}.\end{eqnarray} \tag{ 7 }$$

Although ${\sigma }_{{A}_{l}}$ provides useful information for determining the maximum resolution of the gravity field (or the degree strength; e.g., Konopliv et al. 1999), it does not represent the uncertainty of the power spectrum amplitude A_l . That is, it is representative of the typical error of a single coefficient at that degree (C_lm or S_lm ), and it is not related to the knowledge of the total power at degree l, A_l . Standard uncertainty quantification methods can be used to derive the true uncertainty ${\widetilde{\sigma }}_{{A}_{l}}$. Two equivalent methods can be used: the linear propagation method and the clone method.

The gravity OD estimation process yields the spherical harmonics coefficients of the gravity potential and the corresponding variance–covariance matrix P .

The well-known linear propagation method is defined as

$$\begin{eqnarray}&&{\sigma }_{{A}_{l}}=\sqrt{{\overrightarrow{{\rm{\nabla }}}}_{x}{\boldsymbol{P}}{\overrightarrow{{\rm{\nabla }}}}_{x}^{T}},\end{eqnarray} \tag{ 8 }$$

where ${\overrightarrow{{\rm{\nabla }}}}_{x}$ is defined as

$$\begin{eqnarray*}&&{\overrightarrow{{\rm{\nabla }}}}_{x}=\ \displaystyle \frac{\partial {A}_{l}}{\partial \overrightarrow{x}},\end{eqnarray*}$$

with x = [C₂₀, C₂₁,...,C_lm , S₂₁, S₂₂, S_lm ].

An alternative method consists in generating a large number of clones (N), where each clone represents a realization of the gravity field coefficients drawn from a normal distribution with an expected value equal to the estimated coefficients and covariance matrix P , from which the uncertainty of the spectrum can be computed. The uncertainty on the spectrum is then obtained as the standard deviation (per degree) of these N spectra.

Figure 4 compares the different uncertainty quantification approaches for a simulated gravity field recovery. See Appendix A for a similar comparison applied to the lunar gravity field measured by GRAIL (Goossens et al. 2020).

Zoom In Zoom Out Reset image size

The clone method offers a convenient way to derive not only the uncertainty of the spectrum but also the associated variance–covariance matrix P _A . By generating a large set of clones, it becomes straightforward to compute the covariance between all the samples, resulting in the variance–covariance matrix. As elaborated in the subsequent section, both ${\widetilde{\sigma }}_{{A}_{l}}$ and P _A are necessary inputs for the definition of the Bayesian likelihood.

We can observe that the amplitude spectrum uncertainty is consistently lower than the rms uncertainty spectrum. A substantial difference between the two quantities can be expected, as they convey different information. While the latter, as stated before, relates to the resolution of the field (i.e., the capacity of measuring gravity locally and thus of resolving single coefficients), the former gauges the capability of measuring the global power of the field. Even if the single-coefficient estimate is not accurate owing to aliasing from nonuniform coverage, the spacecraft is still highly sensitive to the averaged strength (or power) of the gravity field driving its motion.

Given the well-established understanding of the Markov Chain Monte Carlo (MCMC) theory, which is extensively discussed in the literature (e.g., Robert & Casella 2004; Sharma 2017), we focus our discussion on notable specifics of our method. To mitigate scaling issues caused by the broad range of magnitudes in the parameter space, we adopt the affine-invariant ensemble sampler provided by the Python library emcee (Foreman-Mackey et al. 2013; Goodman & Weare 2010). By employing this approach, we sample the parameter space with N independent chains that can be run in parallel, reducing the computational burden. To confine the sampling to a defined parameter space region, we impose rigid bounds within the likelihood function. While MCMC processes allow for various weighting schemes across different parameter space regions, we opt for uniform priors on all parameters and reject samples that fall outside these prior bounds. In addition to the sampling strategy, the fundamental component of a Bayesian inversion is the likelihood function, which governs the acceptance or rejection of samples and influences the convergence of the method. In our study, we utilize a multivariate normal distribution:

$$\begin{eqnarray}&&{\boldsymbol{P}}=\displaystyle \frac{1}{\sqrt{{\left(2\pi \right)}^{k}}\left|{\boldsymbol{\Sigma }}\right|}\exp \left(\tfrac{1}{2}{\overrightarrow{v}}^{T}{{\boldsymbol{\Sigma }}}^{-1}\overrightarrow{v}\right),\end{eqnarray} \tag{ 9 }$$

where k is the dimension of the problem (the number of observable quantities), Σ is the variance–covariance matrix of the observations, $\left|{\boldsymbol{\Sigma }}\right|$ is the determinant of Σ, and $\overrightarrow{v}$ is the residuals vector, i.e., the vector containing the difference between the observed and computed observations. The computed observations are obtained from the forward model, while the observed ones are the result of the OD problem solution.

In the case presented here, the set of observations include the mass of the body, the moment of inertia factor (MOIF), and the gravity amplitude spectrum. We build the covariance matrix as a block matrix:

$$\begin{eqnarray}{\boldsymbol{\Sigma }}=\left[\begin{array}{ccc}{{\boldsymbol{\Sigma }}}_{s} & 0 & 0\\ 0 & {\sigma }_{M}^{2} & 0\\ 0 & 0 & {\sigma }_{\mathrm{MOIF}}^{2}\end{array}\right],\end{eqnarray} \tag{ 10 }$$

where ${\sigma }_{M}^{2}$ and ${\sigma }_{\mathrm{MOIF}}^{2}$ are the variances of the observed mass and MOIF, respectively, and Σ_s is the variance–covariance matrix of the gravity spectrum of dimensions $\left({l}_{\max }-2,\ {l}_{\max }-2\right)$.

Since the gravity amplitude spectrum spans several orders of magnitude, to avoid numerical errors that may lead to overweighting of the lowest degrees, we use log(A) as an observable instead of the spectrum A. The associated covariance matrix Σ_s , then, is computed applying the clone method (Section 2.2) to the quantity log(A).

NASA's flagship mission Europa Clipper is scheduled to launch in 2024 October and is expected to arrive in the Jupiter system in 2030. The mission is designed for conducting an extensive mapping of Europa via a series of close encounters (flybys). The gravity/radio science (G/RS) investigation is one of the 10 investigations of Europa Clipper. It will make use of the onboard telecommunication subsystem to establish a two-way Doppler link with the Earth ground stations. In this work we simulate the G/RS investigation during the nominal mission phase, consisting of 49 flybys of Europa. We use the baseline trajectory ⁵ of the nominal science phase, together with a detailed model of the solar and Jupiter system ephemerides. We integrate the trajectory of the probe and produce synthetic observables based on the spacecraft and ground station antenna scheduling. This task holds significant importance for the Europa Clipper investigation, as the flyby geometry, together with the pointing requirements on the spacecraft, constrains the onboard antenna viewing geometry.

The Doppler link must be established using different assets: the high gain antenna (HGA), the fanbeam (FBA) antennas, and the low gain antennas (LGA). FBA and LGA have a lower gain with respect to the HGA, thus reducing the amount of received power at the ground station, but the HGA is not available for tracking close to Europa owing to pointing constraints. As a result, the capability of establishing a successful Doppler link, as well as its resulting noise level, depends on a complex combination of factors involving the flyby geometry, the flyby epoch, and the spacecraft science operation planning. In this work we adhere to the latest assumptions regarding the performance of the tracking system (Mazarico et al. 2023), which assume a noise threshold, for a successful Doppler link, of 4 dB-Hz.

Using the NASA-JPL software library MONTE (Evans et al. 2018), we integrate the trajectory, considering the point-mass gravitational attraction of all main solar system and Jupiter system bodies, the nonspherical contribution of Jupiter's gravity field, and Europa's gravity field (see next section). We account also for the solar radiation pressure acting on the probe, which is modeled as a cannonball. For each Europa interior model case and thus assumed gravity field, we integrate the trajectory of the probe over the 4 hr flyby periods centered on each closest approach, simulate the tracking data, compute the partial derivatives of the observables with respect to the solved-for parameters, and run a full OD solution, iterating the process until convergence. We solve the OD problem using a multi-arc approach (see, e.g., Cascioli et al. 2021; Durante et al. 2020 for further details). The full list of parameters solved for in the OD filter includes the state of the probe (position and velocity with respect to the Jupiter system barycenter), a solar radiation pressure scale factor per flyby, Europa's body fixed frame spin axis orientation and spin rate, the libration angular amplitude (at the orbital frequency), Europa's GM, the spherical harmonics coefficients of Europa's gravity field up to degree and order 20 (in the second inversion; see Section 2.2), and the tidal Love number k₂.

One of the objectives of the G/RS investigation is the measurement of Europa's MOIF ($\tfrac{C}{{{MR}}^{2}}$, with C the polar moment of inertia and M and R the mass and radius of the body, respectively). The measurement accuracy on the MOIF depends on whether Europa satisfies hydrostatic equilibrium (Mazarico et al. 2023). While the icy shell is believed to meet the equilibrium condition because of its low rigidity, less certain predictions can be made about the silicate interior. Previous analyses of Galileo encounters of Europa implicitly assumed hydrostaticity (Anderson et al. 1998; Jacobson et al. 2000), while recent efforts have obtained solutions of Europa's gravity field relaxing the assumption (Gomez Casajus et al. 2021). The predominantly equatorial coverage of the Galileo encounters prevents a reliable estimate of both J₂ and C₂₂, resulting in substantial uncertainties, and precluding concluding whether Europa is in hydrostatic equilibrium or deviates from it. This consideration impacts the various methods employed to estimate the MOIF.

If measurements from Europa Clipper confirm that the body is in hydrostatic equilibrium, the Radau–Darwin equation, which relates J₂, C₂₂, and MOIF, can be directly used, yielding a predicted uncertainty on the MOIF of approximately 1.5 × 10⁻⁴. However, if the Radau–Darwin relationship cannot be used, alternative methods must be employed, leading to a higher uncertainty of the order of 2.3 × 10⁻³ (see Mazarico et al. 2023 and references therein for a comprehensive discussion on this topic).

As discussed in Section 2.2, the second observable quantity employed in the Bayesian inversion is the mass of the body. It serves the purpose of maintaining internal consistency within the inversion process. The uncertainty associated with the mass is derived from the estimated uncertainty of GM, multiplied by a factor of 1000. We use an uncertainty value of 8.8 × 10¹⁸ kg, which corresponds to approximately 0.01% of Europa's mass.

We establish a reference case for Europa by adopting values for key model quantities that align with those published in previous studies (Pauer et al. 2010; Koh et al. 2022). Table 1 reports the values used as a baseline. It can be argued that the chosen value of the maximum silicate layer topography (2 km) is relatively small, if compared to other differentiated silicate worlds in the solar system (e.g., Wieczorek 2015; Ermakov et al. 2018). However, it will be clear in the next section (and Appendix B) that this assumption constitutes a worst-case scenario for the methodology we propose in this work, as a larger topography would imply a stronger gravitational signal of the seafloor and thus a better overall sensitivity to it, resulting in a better determination.

Some authors, moreover, have predicted the possible presence of a thick (3–10 km) sedimentary layer of gypsum on the seafloor (Melwani Daswani et al. 2021). Two possible end-member models, as proposed by Koh et al. (2022), would consist of a flattening of the seafloor due to the gypsum deposits filling the topography, or alternatively a uniform blanket that would uniformly increase the seafloor topography. While the second end-member does not have any significant effect on the gravity signal of the seafloor (as shown by Koh et al. 2022), the first end-member can substantially change its amplitude. Although not explicitly introduced here, we argue that this end-member case can be equivalently modeled by either a modification of the density contrast between the silicate crust and the ocean or a reduction of the topographic profile (as a consequence of Equations (3a) and (3b)). Among the end-member models we tested in this work (table in Figure 5), we have considered ample variations in the topographic profiles, thus implicitly accounting for the possible presence of a sedimentary deposit on the base of the ocean.

The gravity model based on parameters from Table 1 is shown in Figure 6. The difference between the gravity fields of the silicate ocean floor and the icy shell is particularly evident in the scale of the associated gravity anomalies.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The icy shell's low rigidity leads to isostatic compensation at the longest wavelengths, resulting in the suppression of large-scale anomalies and the prominence of smaller scales. Conversely, the higher rigidity of the seafloor enables the persistence of significant power in longer spatial scales. This fundamental mechanism allows us to effectively distinguish and separate the contributions of the two.

The morphology and amplitude of the resulting gravity spectrum are highly dependent on the choice of the physical parameters defining the body's interior. Figure 7 shows the effect of varying the ocean depth. Variations in the ocean depth result in strong variations in the amplitude and shape of the gravity spectrum and variations in the MOIF larger than the measurement uncertainty.

Zoom In Zoom Out Reset image size

To ensure a comprehensive exploration of the sensitivity of the Bayesian inversion process, we recognize the need to consider multiple cases rather than relying on a single interior structure scenario. Figure 7 serves as an illustrative example, demonstrating that the amplitude of the gravity spectrum depends on the interior structure configuration. This variation implies differing degrees of sensitivity for the spacecraft measurements. Specifically, a deep ocean scenario results in a reduced amplitude and lower sensitivity to seafloor variations, while a shallow ocean scenario gives greater sensitivity, potentially leading to a more robust determination of the ocean depth and seafloor parameters. Consequently, we conducted an extensive examination of the parameter space, encompassing a wide range of cases that represent the extreme ends of the spectrum. We explored the parameter space, varying the main parameters of the model (see Appendix B for the full set of figures). We were thus able to determine the subset of parameters whose plausible variations have a limited impact on the spectrum and moment of inertia, allowing us, then, to fix them and not solve for them. From the other, more important parameters, we defined six different interior structure scenarios. The table in Figure 5 details these six cases we will explore in the following.

As the figure shows, we chose to fix d_crust, ρ_ocean, ${\rho }_{\mathrm{core}}$, and ρ_crust because of the reduced sensitivity of the observables to these parameters. By reducing the number of solved-for parameters, we mitigate correlations among them and minimize the dimensionality of the parameter space, thereby improving computational efficiency. However, to test the solution's robustness against potential biases resulting from fixed parameter assumptions, we define subcases where these parameters are still held fixed during inversion but perturbed off the actual, correct value (Table 2). The amount of perturbation is determined by assessing the parameter space (Appendix B) and identifying the variations corresponding to a 3σ offset from the computed average MOIF value (with 1σ representing the predicted uncertainty in the hydrostatic case). This approach ensures a thorough evaluation of the solution's sensitivity to variations in these parameters.

Table 2. Definition of Subcases

	Model Value	3σ Distance	Subcase a	Subcase b	Subcase c
dcrust (km)	30	20	+	⋯	⋯
ρocean (g cc–1)	1.01	0.03	+	+	+
ρcrust (g cc–1)	2.9	0.12	+	⋯	+
${\rho }_{\mathrm{core}}$ (g cc–1)	5.61	0.3	⋯	+	⋯

Note. Plus and minus signs indicate that the 3σ distance is added or subtracted to the model value. The 3σ distance is defined as the parameter value corresponding to a distance equivalent to three times the predicted formal uncertainty on the MOIF in the hydrostatic case.

Download table as:  ASCIITypeset image