Enceladus’s and Dione’s floating ice shells
supported by minimum stress isostasy
Abstract
Enceladus’s gravity and shape have been explained in terms of a thick isostatic ice shell floating on a global ocean, in contradiction of the thin shell implied by librations. Here we propose a new isostatic model minimizing crustal deviatoric stress, and demonstrate that gravity and shape data predict a thick ocean beneath a thick shell agreeing with – but independent of – libration data. Isostatic and tidal stresses are comparable in magnitude. South polar crust is only thick, facilitating the opening of water conduits and enhancing tidal dissipation through stress concentration. Enceladus’s resonant companion, Dione, is in a similar state of minimum stress isostasy. Its gravity and shape can be explained in terms of a thick isostatic shell overlying a thick global ocean, thus providing the first clear evidence for a presentday ocean within Dione.
Paper accepted for publication in Geophysical Research Letters
http://agupubs.onlinelibrary.wiley.com/agu/journal/10.1002/(ISSN)19448007/
1 Introduction
Saturn’s moon Enceladus is celebrated for its huge south polar fractures venting jets of water vapor and ice particles [Porco et al., 2006] while its neighbor Dione is more discreet, though essential in maintaining Enceladus’s eccentricity through a 2:1 orbital resonance. The composition of Enceladus’s plume shows that it originates in a subsurface ocean in contact with a silicate core [Postberg et al., 2011; Hsu et al., 2015], but a more detailed picture of the interior must be based on geodesy. Inferences about the internal structure of Enceladus and Dione were first based on their longwavelength shape [Thomas et al., 2007; Nimmo et al., 2011] and gravity field [Iess et al., 2014; Hemingway et al., 2016]. These data reveal the presence of a hydrated silicate core with a radius of about threefourths of the surface radius, but the crustocean partition of the outer layer has remained controversial. The recent measurement of Enceladus’s librations [Thomas et al., 2016; Nadezhdina et al., 2016] demonstrates that the crust is a thin ice shell floating on a global ocean.
For both satellites, the task of building interior models compatible with the gravity and shape data is complicated by strong deviations from hydrostatic equilibrium. For synchronously rotating satellites, nonhydrostatic deviations can be estimated either from the gravity ratio or from the analog ratio for the shape (Text S.1). If the satellite is hydrostatic, these shape and gravity ratios are both equal to 10/3 to first order in the flattening [Zharkov et al., 1985]. Using the shape and gravity data of Cassini and taking into account secondorder corrections, one concludes that the observed shape and gravity ratios of Enceladus deviate from hydrostaticity by thirty and ten percent, respectively, whereas the corresponding deviations for Dione are twice as large (Text S.1). Such differences between the nonhydrostatic components of shape and gravity are indicative of isostasy, in which surface topography is mechanically supported and gravitationally compensated by a subsurface mass anomaly in such a way that below a certain depth, called the compensation depth, pressure is everywhere hydrostatic [Lambeck, 1980, 1988]. Enceladus’s degreethree zonal gravity harmonic () also points to isostasy: it is fully nonhydrostatic and nonzero within its error interval, but only a third of what is expected from the corresponding shape harmonic. Degreethree compensation is attributed to isostatic support of the south polar depression [Iess et al., 2014].
In this paper, we reexamine the two main aspects of the method of Iess et al. [2014] and McKinnon [2015]: the secondorder figure of equilibrium and the isostatic compensation model. Our computation of the figure of equilibrium is more rigorous but does not make a significant numerical difference. By contrast, our new theory of isostasy has important implications for the structure of the layer. We put our conclusions on a firm footing by doing a Bayesian inversion of the gravity data taking into account the uncertainties on the shape and gravity.
2 The Problem with Classical Isostasy
Isostasy can occur either through crustal density variations (Pratt) or through variations in crustal thickness (Airy). For Pratt isostasy, the only plausible scenario involves porosity variations close to the surface, but compensation is too high [McKinnon, 2015]. With Airy isostasy, Enceladus’s ice shell was initially estimated to be on average 30 to thick [Iess et al., 2014], but this value was revised to due to secondorder tidalrotational effects [McKinnon, 2015]. A thick shell is hardly compatible with the south polar activity, raises the issue of shellcore contact at the equator precluding isostasy, and does not match degreethree compensation [McKinnon, 2015]. Furthermore, the result conflicts with libration models predicting that the shell is half as thick [Thomas et al., 2016; Van Hoolst et al., 2016]. Compared to Enceladus, Dione’s case is even more serious because the thick ice shell implied by Airy isostasy does not fit into the thick layer [Hemingway et al., 2016] (see also Fig. 1). Although a solution is possible within uncertainties, it entails that Dione’s shape deviates even more from hydrostaticity.
The easy way out is to suppose that an elastic lithosphere partly supports the load (flexural isostasy) so that the subsurface mass anomaly can be smaller and located closer to the surface. Flexural support of longwavelength loads, however, generates stresses that are not only larger than the tensional strength of intact ice (1 to [Schulson and Duval, 2009]) but also much above the Coulomb failure criterion for a pervasively cracked thin lithosphere [McKinnon, 2013]. Cadek et al. [2016] nevertheless argue for partial support of Enceladus’s topography by a thin lithosphere, without addressing the question of lithospheric failure. Besides, they overestimate the elastic thickness required for top loading by a factor of ten [Turcotte et al., 1981; McGovern et al., 2002; Kalousová et al., 2012]. We give a detailed rebuttal of Cadek et al. [2016] in the Supporting Information (Text S.7).
Classical Airy (or Pratt) isostasy is plagued by ambiguities about the isostatic prescription, i.e. the constraint relating the subsurface mass anomaly to the surface load [Vening Meinesz, 1946]. At the longest wavelengths, various isostatic prescriptions lead to geoid anomalies differing by more than a factor of two [Dahlen, 1982]. This question has never been settled because of the complexity of largescale isostasy on Earth: thermal (Pratt) and Airy isostasy are dominant in oceanic and continental crusts, respectively, while the longwavelength geoid is explained by mantle convection. In planetology, isostasy invariably resorts to the equalmass prescription applied to conical columns [Lambeck, 1988]. This prescription however neither takes into account horizontal stresses nor geoid perturbations due to the loads themselves. Here we consider instead (following Jeffreys [1970] and Dahlen [1982]) that the only physically meaningful isostatic prescription consists in minimizing crustal deviatoric stress in a selfconsistent elasticgravitational theory.
3 Methods
3.1 HydrostaticIsostatic Decomposition
Before explaining the three methods particular to this paper, we recall the principle of the hydrostaticisostatic decomposition. Let and denote the cosine real harmonic coefficients of degree and order of the shape and nondimensional gravity potential, respectively. Following Iess et al. [2014], we split the harmonic components into hydrostatic and isostatic (nonhydrostatic) components:
(1)  
(2) 
The hydrostatic components with and are determined by computing the figure of equilibrium (the other hydrostatic components vanish). At each harmonic degree , the isostatic components are related by the isotropic admittance (or equivalently by the nondimensional compensation factor ) characterizing the isostatic model:
(3) 
where is the ice shell density, is the bulk density and is the surface radius.
3.2 Figure of Equilibrium
Enceladus’s rapid rotation leads to significant deviations from the first order figure of equilibrium. Contrary to McKinnon [2015], we do not compute the figure of equilibrium with the method of Tricarico [2014], which is based on nested ellipsoids, because an ellipsoidal stratification is forbidden for heterogenous hydrostatic bodies (HamyPizzetti theorem, see Moritz [1990]). We work instead with the classical theory of equilibrium figures extended from first order (Clairaut’s equation) to second order in the flattening [Zharkov, 2004].
Solving the secondorder equations for a multilayer body yields complicated expressions for the degreetwo shape and gravity coefficients (the latter normalized with respect to the mean radius). These complete solutions to second order are closely approximated by compact formulas depending on the fluid Love numbers and :
(4)  
(5) 
where are the secondorder corrections and is the nondimensional rotational parameter (Table S.1). Fluid Love numbers of two and threelayer bodies are given in Zharkov [2004] and Zharkov and Gudkova [2010], respectively. We computed from using the secondorder relations between shape and gravity (Eq. (49) of Zharkov and Gudkova [2010]). Our formalism for the harmonic coefficients of degree two is equivalent to the one of Zharkov and Gudkova [2010], except that their Eqs. (51)(52) and (57)(58) are wrong and must be replaced by the two equations above.
The error on Eqs. (4)(5) is similar to the error on the method of Tricarico [2014]. It is of third order in the flattening if the satellite is homogeneous. If not, the error is formally of second order, but numerically of third order. Compared to the method of Tricarico [2014], Eqs. (4)(5) are simpler to use and more easily applicable to multilayer bodies. For the Bayesian inversion, we do not use Eqs. (4)(5) but the complete solutions of the classical approach to second order (i.e. the error is strictly of third order).
3.3 Minimum Stress Isostasy
The principle of minimum stress isostasy is based on the idea that, over time, the crust has reached the state of minimum deviatoric stress compatible with the observed topography through internal deformation, lithospheric failure and viscoelastic relaxation [Dahlen, 1981]. The last mechanism depends on the possible viscoelastic processes and the loading timescale. These questions cannot be answered without knowing more about the crust (rheology and thermal state) and about the origin of the putative shell thickness variations. If desired, one can relax the assumption of minimum stress so as to include less bottom loading and more flexural support (flexural isostasy), but the model is not predictive if we do not know the elastic thickness of the lithosphere. Our results will show that such additional parameter is not needed, and we tend to favor the simplest model that can successfully explain the observations. We thus assume here that the crust has reached the state of minimum stress isostasy at the largest scale.
In order to minimize crustal stresses, we need to know the deformations, gravity field perturbations, and stresses due to loads acting on the top and bottom of the crust. In terms of standard geophysical techniques, this means computing the elastic Love numbers of a selfgravitating spherically symmetric body submitted on the one hand to a surface load and, on the other, to an internal load located at the crustocean boundary. The two solutions are linearly combined with an arbitrary loading ratio, which is then fixed by minimizing the second invariant of the deviatoric stress tensor. The result takes the form of Eq. (3).
We solve the elasticgravitational problem with the incompressible propagator matrix method [Sabadini and Vermeersen, 2004]. The surface and internal loads are modeled as thin layers with densities (per unit surface) and , respectively. The surface load Love numbers and the internal load Love numbers are obtained as boundary values of the full solution. By definition, the perturbation of the gravity potential (dimensional here: ) and the radial displacement of the surface read [GreffLefftz et al., 2010]
(6) 
where and are the loading potentials ( is the surface gravity and where is the shell thickness). The former potential depends on the topography through where , while the latter can initially be written as , where the loading ratio is a negative number of order unity to be fixed by the isostatic prescription. The compensation factor follows from Eq. (3):
(7) 
The final step consists in computing the loading ratio with the chosen isostatic prescription. The magnitude of deviatoric stresses can be measured with the second invariant of the deviatoric stress tensor [Dahlen, 1982] which is proportional to the shear (or distortional) strain energy density [Jaeger et al., 2007]. For each harmonic degree , we minimize the total shear energy of the crust , which depends quadratically on crustal deformations (Eq. (8.128) of Dahlen and Tromp [1999] restricted to the crustal volume). Since crustal deformations are linear combinations of the deformations due to the surface and internal loads, the total shear energy of the crust can be decomposed as where the terms do not depend on . Crustal stresses are minimum for the loading ratio , which is a computable quantity once the elasticgravitational problem has been solved for surface and internal loads.
In minimum stress isostasy, crustal deformation is very small: with . If deflection is neglected, the topography of the crustocean boundary is given by
(8) 
This extension to finite amplitude topography introduces negligible errors in the gravitational potential (Text S.6). The local shell thickness is determined from the average shell thickness and the coefficients (Eq. (S.9) in Text S.5).
3.4 Bayesian Inversion
The figure of equilibrium and the isostatic model define the forward problem: given the interior structure and the isostatic load (), we can predict the shape and gravity coefficients. What interests us more is the inverse problem which is nonlinear, under or overconstrained depending on the parameter, and based on uncertain gravity and shape data. It is thus well suited to a Bayesian inference method [Sambridge and Gallagher, 2011]. The result of the Bayesian inversion is the posterior probability density function, i.e. the conditional probability for the parameters given the data [Tarantola, 2005; Gregory, 2005]. We simulate the posterior probability density function with a MetropolisHastings sampler. From the generated samples, we compute for each parameter the probability density function, the mean value and the Bayesian confidence intervals (Text S.4).
Enceladus and Dione are modeled as threelayer incompressible bodies made of an elastic core, an ocean, and an elastic shell. The model parameters are the densities (, ) and thicknesses (, ) of the ice shell and ocean (the radius and density of the core are derived parameters). The prior information is described by uniform probability density functions subjected to the constraint that the shell thickness at the south pole is not negative (Text S.4 and Fig. S.1). For Enceladus, the uncorrelated prior ranges are for and for . For Dione, the prior ranges are for and for . Enceladus’s ocean is probably similar in salinity to the ice grains in the plume ( salt by mass [Postberg et al., 2011]) while the detection of silica nanoparticles sets a 4% upper bound on the salinity [Hsu et al., 2015]. The prior range for is ( salt by mass) and we assume the same for Dione. Given the low pressure within the crust, the density of pure ice varies between and depending on the temperature, but porous ice is lighter whereas saltrich ice is denser. For the threelayer body, the prior range for extends from (10% porosity in the upper half of the crust) to (no porosity, similar salt content as the ocean maximum). For Enceladus, we model a very porous crust with a fourlayer model in which the crust is made of a bottom layer with fixed density () and a thick upper layer with a prior density range of [Besserer et al., 2013].
4 Results
Given the compensation factor, minimum stress isostasy requires a shell nearly half as thick as in classical isostasy (Fig. 1). The compensation factor is most sensitive to the shell thickness which plays the role of compensation depth. The shear moduli of the shell and core (radius ) are set to (pure ice) and (hydrated silicates), respectively. The compensation factor is nearly independent of , , and as long as the shell and core are not fluidlike (i.e. for Enceladus). As a simple example, we solve Eqs. (1)–(5) for the two parameters and the eight hydrostatic/isostatic components of degree two, assuming zero data uncertainty and an icedominated layer ( as in McKinnon [2015]). For Enceladus, and , while and for Dione. Fig. 1 shows that classical Airy isostasy predicts a very thin ocean for Enceladus and no ocean at all for Dione.
For Enceladus, our reference data are the coefficients of the gravity solution SOL1 (Table S.2), combined with the shape TOPA (Table S.4). For this model, the inversion yields , , , and at (see Fig. S.1 for posterior distributions). Shell and ocean densities are not constrained. The shell and ocean thicknesses are inversely correlated (Fig. 2) because the core radius is well determined (Fig. S.2). Our estimates for Enceladus’s shell thickness overlap with those of librations ( in Thomas et al. [2016], in Van Hoolst et al. [2016]). For a rigorous comparison, we compute the librations from the probability distribution over the parameters inferred from our gravityshape inversion assuming a rigid shell with nonhydrostatic boundaries [Van Hoolst et al., 2016]. The predicted distribution ( at ) is wider than the distribution of observed librations ( at ) (Fig. 3). Thus librations put tighter bounds on the average shell thickness of Enceladus, although they do not constrain the other interior parameters.
For Dione, the inversion of gravity and shape data (Tables S.3 and S.5) yields , , , and (at ). Errors are much larger than for Enceladus (Fig. 2) because of the large relative error on the shape (Tables S.5 and S.7). The gravity and shape data can thus be explained if there is a global ocean deep under the surface, whose past existence was already suggested by observation of ridge flexure [Hammond et al., 2013] and thermal history models [Multhaup and Spohn, 2007]. We predict that Dione undergoes librations of amplitude at (Fig. 3), one order of magnitude below Enceladus and thus not detectable in Cassini images.
5 Discussion
Comparison with librations can pinpoint data biases and constrain modeling choices. First, Enceladus’s degreethree gravity favors a thinner shell than degreetwo coefficients [McKinnon, 2015] (Fig. 1 and Table S.6). Degreetwo gravity however predicts a thinner shell (in agreement with librations) if is higher than its SOL1 central value, as suggested by the alternative gravity solution SOL2 (Table S.2). Alternatively, the degreetwo shape could be responsible for the disagreement: while the new ellipsoidal shape of Thomas et al. [2016] does not affect much the results (Table S.6), the recent ellipsoidal shape of Nadezhdina et al. [2016] predicts more degreetwo compensation and thus a thinner shell. Second, we did not allow for a lot of surface porosity in our threelayer model, but we can easily do it with a fourlayer model. The estimated shell thickness increases with porosity because surface topography contributes less to the gravity signal and must be less compensated (Table S.6 and Fig. S.3). Consistency with librations suggests however that porosity is not an important factor.
Enceladus’s shell thickness varies mainly in latitude from at the equator (zonal average) to and at the north and south poles, respectively (Fig. 4 and Text S.5). Longitudinal variations are either subdominant, with shell thickening along the tidal axis, or could be absent altogether as suggested by the latest estimates of the ellipsoidal shape (Table S.4). The very thin south polar crust facilitates the passage of water from the ocean to the surface and increases the concentration of tidal heating in the area. The variation in shell thickness could be due to crustal tidal heating which is indeed highest at the poles and lowest along the tidal axis [Ojakangas and Stevenson, 1989]. For Dione, shell thickness varies by less than five percent with a minimum at the poles and a maximum along the leadingtrailing axis. Similarly to Enceladus, the zonal variation in shell thickness could be due to tidal heating, but we have no ready explanation for the longitudinal variation anticorrelated with tidal heating. Both satellites are more tectonized or resurfaced in their leading and trailing hemispheres than close to their tidal axis [CrowWillard and Pappalardo, 2015; Kirchoff and Schenk, 2015].
The average crustal stress is about for Enceladus, which is half of the average topographic stress, as expected in isostasy [Melosh, 2011]. It is comparable in magnitude to tidal stresses [Nimmo et al., 2007a] and could trigger the formation of the south polar terrain margins by gravitational spreading [Yin and Pappalardo, 2015]. We also evaluated isostatic stresses in the core in order to check that they are always smaller than isostatic stresses in the crust. Therefore, the core can be modeled as an elastic layer regarding isostatic loading. A final contentious point is that isostasy is only valid to first order in the flattening, contrary to the figure of equilibrium. The secondorder error on the isostatic model, however, changes the total gravity potential [Wieczorek and Phillips, 1998] by less than the data uncertainty (Text S.6).
Beyond Enceladus and Dione, our new take on isostasy is applicable to large icy satellites with global oceans, such as Europa [Nimmo et al., 2007b], Titan [Nimmo and Bills, 2010], and particularly Ganymede whose gravity and shape will be measured by the JUICE mission. According to Park et al. [2016], gravity and shape data from the Dawn mission suggest isostasy on Ceres, but the case is far from clear because compensation does not occur for all gravity components. Finally, isostasy plays a crucial role in understanding the longwavelength gravity and shape as well as estimating the crust thickness of the planets Mars [Wieczorek and Zuber, 2004], Venus [James et al., 2013], and Mercury [Perry et al., 2015]. Thanks to the simultaneous availability of gravity/shape and libration data, Enceladus’s case constitutes the first validation of planetaryscale isostasy.
Acknowledgments
All data used here are freely available in the literature. M.B. is supported by the Brain Pioneer contract BR/314/PI/LOTIDE. A.R. is supported by the Belgian PRODEX program managed by the European Space Agency in collaboration with the Belgian Federal Science Policy Office. A.T. received support from the ‘Supplementary Researcher’ programme managed by the Belgian Federal Science Policy Office, and from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 670874).
References
 Besserer et al. [2013] Besserer, J., F. Nimmo, J. H. Roberts, and R. T. Pappalardo (2013), Convectiondriven compaction as a possible origin of Enceladus’s long wavelength topography, J. Geophys. Res., 118, 908–915, doi:10.1002/jgre.20079.
 Cadek et al. [2016] Cadek, O., G. Tobie, T. Van Hoolst, M. Massé, G. Choblet, A. Lefèvre, G. Mitri, R.M. Baland, M. Behounkova, O. Bourgeois, and A. Trinh (2016), Enceladus’s internal ocean and ice shell constrained from Cassini gravity, shape and libration data, Geophys. Res. Lett., 43, 5653–5660, doi:10.1002/2016GL068634.
 CrowWillard and Pappalardo [2015] CrowWillard, E. N., and R. T. Pappalardo (2015), Structural mapping of Enceladus and implications for formation of tectonized regions, J. Geophys. Res., 120, 928–950, doi:10.1002/2015JE004818.
 Dahlen [1981] Dahlen, F. A. (1981), Isostasy and the ambient state of stress in the oceanic lithosphere, J. Geophys. Res., 86, 7801–7807, doi:10.1029/JB086iB09p07801.
 Dahlen [1982] Dahlen, F. A. (1982), Isostatic geoid anomalies on a sphere, J. Geophys. Res., 87, 3943–3947, doi:10.1029/JB087iB05p03943.
 Dahlen and Tromp [1999] Dahlen, F. A., and J. Tromp (1999), Theoretical Global Seismology, Princeton University Press.
 GreffLefftz et al. [2010] GreffLefftz, M., L. Métivier, and J. Besse (2010), Dynamic mantle density heterogeneities and global geodetic observables, Geophys. J. Int., 180, 1080–1094, doi:10.1111/j.1365246X.2009.04490.x.
 Gregory [2005] Gregory, P. C. (2005), Bayesian Logical Data Analysis for the Physical Sciences, Cambridge University Press, Cambridge.
 Hammond et al. [2013] Hammond, N. P., C. B. Phillips, F. Nimmo, and S. A. Kattenhorn (2013), Flexure on Dione: Investigating subsurface structure and thermal history, Icarus, 223, 418–422, doi:10.1016/j.icarus.2012.12.021.
 Hemingway et al. [2016] Hemingway, D. J., M. Zannoni, P. Tortora, F. Nimmo, and S. W. Asmar (2016), Dione’s internal structure inferred from Cassini gravity and topography, Lunar Planet. Sci. Conf. 47, abstract 1314.
 Hsu et al. [2015] Hsu, H.W., F. Postberg, Y. Sekine, T. Shibuya, S. Kempf, M. Horányi, A. Juhász, N. Altobelli, K. Suzuki, Y. Masaki, T. Kuwatani, S. Tachibana, S.I. Sirono, G. MoragasKlostermeyer, and R. Srama (2015), Ongoing hydrothermal activities within Enceladus, Nature, 519, 207–210, doi:10.1038/nature14262.
 Iess et al. [2014] Iess, L., D. J. Stevenson, M. Parisi, D. Hemingway, R. A. Jacobson, J. I. Lunine, F. Nimmo, J. W. Armstrong, S. W. Asmar, M. Ducci, and P. Tortora (2014), The gravity field and interior structure of Enceladus, Science, 344, 78–80, doi:10.1126/science.1250551.
 Jaeger et al. [2007] Jaeger, J. C., N. G. W. Cook, and R. W. Zimmerman (2007), Fundamentals of Rock Mechanics, 4th edition, Blackwell.
 James et al. [2013] James, P. B., M. T. Zuber, and R. J. Phillips (2013), Crustal thickness and support of topography on Venus, J. Geophys. Res., 118, 859–875, doi:10.1029/2012JE004237.
 Jeffreys [1970] Jeffreys, H. (1970), The Earth, Cambridge University Press, Cambridge.
 Kalousová et al. [2012] Kalousová, K., O. Souček, and O. Čadek (2012), Deformation of an elastic shell with variable thickness: a comparison of different methods, Geophys. J. Int., 190, 726–744, doi:10.1111/j.1365246X.2012.05539.x.
 Kirchoff and Schenk [2015] Kirchoff, M. R., and P. Schenk (2015), Dione’s resurfacing history as determined from a global impact crater database, Icarus, 256, 78–89, doi:10.1016/j.icarus.2015.04.010.
 Lambeck [1980] Lambeck, K. (1980), Estimates of stress differences in the crust from isostatic considerations, J. Geophys. Res., 85, 6397–6402, doi:10.1029/JB085iB11p06397.
 Lambeck [1988] Lambeck, K. (1988), Geophysical Geodesy, Clarendon Press, Oxford.
 McGovern et al. [2002] McGovern, P. J., S. C. Solomon, D. E. Smith, M. T. Zuber, M. Simons, M. A. Wieczorek, R. J. Phillips, G. A. Neumann, O. Aharonson, and J. W. Head (2002), Localized gravity/topography admittance and correlation spectra on Mars: Implications for regional and global evolution, J. Geophys. Res., 107, 5136, doi:10.1029/2002JE001854.
 McKinnon [2013] McKinnon, W. B. (2013), The shape of Enceladus as explained by an irregular core: Implications for gravity, libration, and survival of its subsurface ocean, J. Geophys. Res., 118, 1775–1788, doi:10.1002/jgre.20122.
 McKinnon [2015] McKinnon, W. B. (2015), Effect of Enceladus’s rapid synchronous spin on interpretation of Cassini gravity, Geophys. Res. Lett., 42, 2137–2143, doi:10.1002/2015GL063384.
 Melosh [2011] Melosh, H. J. (2011), Planetary Surface Processes, Cambridge University Press, Cambridge.
 Moritz [1990] Moritz, H. (1990), The Figure of the Earth, Wichmann, Karlsruhe.
 Multhaup and Spohn [2007] Multhaup, K., and T. Spohn (2007), Stagnant lid convection in the midsized icy satellites of Saturn, Icarus, 186, 420–435, doi:10.1016/j.icarus.2006.09.001.
 Nadezhdina et al. [2016] Nadezhdina, I. E., A. E. Zubarev, E. S. Brusnikin, and J. Oberst (2016), A libration model for Enceladus based on geodetic control point network analysis, Int. Arch. Photogramm. Remote Sens. Spatial Inf. Sci., XLIB4, 459–462, doi:10.5194/isprsarchivesXLIB44592016.
 Nimmo and Bills [2010] Nimmo, F., and B. G. Bills (2010), Shell thickness variations and the longwavelength topography of Titan, Icarus, 208, 896–904, doi:10.1016/j.icarus.2010.02.020.
 Nimmo et al. [2007a] Nimmo, F., J. R. Spencer, R. T. Pappalardo, and M. E. Mullen (2007a), Shear heating as the origin of the plumes and heat flux on Enceladus, Nature, 447, 289–291, doi:10.1038/nature05783.
 Nimmo et al. [2007b] Nimmo, F., P. C. Thomas, R. T. Pappalardo, and W. B. Moore (2007b), The global shape of Europa: Constraints on lateral shell thickness variations, Icarus, 191, 183–192, doi:10.1016/j.icarus.2007.04.021.
 Nimmo et al. [2011] Nimmo, F., B. G. Bills, and P. C. Thomas (2011), Geophysical implications of the longwavelength topography of the Saturnian satellites, J. Geophys. Res., 116, E11001, doi:10.1029/2011JE003835.
 Ojakangas and Stevenson [1989] Ojakangas, G. W., and D. J. Stevenson (1989), Thermal state of an ice shell on Europa, Icarus, 81, 220–241, doi:10.1016/00191035(89)900523.
 Park et al. [2016] Park, R. S., A. S. Konopliv, B. Bills, N. Rambaux, J. CastilloRogez, C. A. Raymond, A. T. Vaughan, A. Ermakov, M. T. Zuber, R. R. Fu, M. J. Toplis, C. T. Russell, A. Nathues, and F. Preusker (2016), A partially differentiated interior for (1) Ceres deduced from its gravity field and shape, Nature, 537, 515–517, doi:10.1038/nature18955.
 Perry et al. [2015] Perry, M. E., G. A. Neumann, R. J. Phillips, O. S. Barnouin, C. M. Ernst, D. S. Kahan, S. C. Solomon, M. T. Zuber, D. E. Smith, S. A. Hauck, S. J. Peale, J.L. Margot, E. Mazarico, C. L. Johnson, R. W. Gaskell, J. H. Roberts, R. L. McNutt, and J. Oberst (2015), The lowdegree shape of Mercury, Geophys. Res. Lett., 42, 6951–6958, doi:10.1002/2015GL065101.
 Porco et al. [2006] Porco, C. C., P. Helfenstein, P. C. Thomas, A. P. Ingersoll, J. Wisdom, R. West, G. Neukum, T. Denk, R. Wagner, T. Roatsch, S. Kieffer, E. Turtle, A. McEwen, T. V. Johnson, J. Rathbun, J. Veverka, D. Wilson, J. Perry, J. Spitale, A. Brahic, J. A. Burns, A. D. Del Genio, L. Dones, C. D. Murray, and S. Squyres (2006), Cassini observes the active south pole of Enceladus, Science, 311, 1393–1401, doi:10.1126/science.1123013.
 Postberg et al. [2011] Postberg, F., J. Schmidt, J. Hillier, S. Kempf, and R. Srama (2011), A saltwater reservoir as the source of a compositionally stratified plume on Enceladus, Nature, 474, 620–622, doi:10.1038/nature10175.
 Sabadini and Vermeersen [2004] Sabadini, R., and B. Vermeersen (2004), Global Dynamics of the Earth, Kluwer Academic Publishers, Dordrecht.
 Sambridge and Gallagher [2011] Sambridge, M., and K. Gallagher (2011), Inverse Theory, Monte Carlo Method, in Encyclopedia of Solid Earth Geophysics, edited by H. K. Gupta, pp. 639–644, Springer Netherlands, Dordrecht, doi:10.1007/9789048187027_192.
 Schulson and Duval [2009] Schulson, E. M., and P. Duval (2009), Creep and Fracture of Ice, Cambridge University Press, Cambridge.
 Tarantola [2005] Tarantola, A. (2005), Inverse Problem Theory and Methods for Model Parameter Estimation, SIAM, Philadelphia.
 Thomas et al. [2007] Thomas, P. C., J. A. Burns, P. Helfenstein, S. Squyres, J. Veverka, C. Porco, E. P. Turtle, A. McEwen, T. Denk, B. Giese, T. Roatsch, T. V. Johnson, and R. A. Jacobson (2007), Shapes of the saturnian icy satellites and their significance, Icarus, 190, 573–584, doi:10.1016/j.icarus.2007.03.012.
 Thomas et al. [2016] Thomas, P. C., R. Tajeddine, M. S. Tiscareno, J. A. Burns, J. Joseph, T. J. Loredo, P. Helfenstein, and C. Porco (2016), Enceladus’s measured physical libration requires a global subsurface ocean, Icarus, 264, 37–47, doi:10.1016/j.icarus.2015.08.037.
 Tricarico [2014] Tricarico, P. (2014), Multilayer hydrostatic equilibrium of planets and synchronous moons: Theory and application to Ceres and to solar system moons, Astr. J., 782, 99, doi:10.1088/0004637X/782/2/99.
 Turcotte et al. [1981] Turcotte, D. L., R. J. Willemann, W. F. Haxby, and J. Norberry (1981), Role of membrane stresses in the support of planetary topography, J. Geophys. Res., 86, 3951–3959, doi:10.1029/JB086iB05p03951.
 Van Hoolst et al. [2016] Van Hoolst, T., R.M. Baland, and A. Trinh (2016), The diurnal libration and interior structure of Enceladus, Icarus, 277, 311–318, doi:10.1016/j.icarus.2016.05.025.
 Vening Meinesz [1946] Vening Meinesz, F. A. (1946), The indirect isostatic or Bowie reduction and the equilibrium figure of the Earth, Bull. Géod., 1, 33–107, doi:10.1007/BF02519025.
 Wieczorek and Phillips [1998] Wieczorek, M. A., and R. J. Phillips (1998), Potential anomalies on a sphere  Applications to the thickness of the lunar crust, J. Geophys. Res., 103, 1715–1724, doi:10.1029/97JE03136.
 Wieczorek and Zuber [2004] Wieczorek, M. A., and M. T. Zuber (2004), Thickness of the Martian crust: Improved constraints from geoidtotopography ratios, J. Geophys. Res., 109, E01009, doi:10.1029/2003JE002153.
 Yin and Pappalardo [2015] Yin, A., and R. T. Pappalardo (2015), Gravitational spreading, bookshelf faulting, and tectonic evolution of the South Polar Terrain of Saturn’s moon Enceladus, Icarus, 260, 409–439, doi:10.1016/j.icarus.2015.07.017.
 Zharkov [2004] Zharkov, V. N. (2004), A theory of the equilibrium figure and gravitational field of the Galilean satellite Io: The second approximation, Astr. Lett., 30, 496–507, doi:10.1134/1.1774402.
 Zharkov and Gudkova [2010] Zharkov, V. N., and T. V. Gudkova (2010), Models, figures and gravitational moments of Jupiter’s satellite Io: Effects of the second order approximation, Planet. Space Sci., 58, 1381–1390, doi:10.1016/j.pss.2010.06.004.
 Zharkov et al. [1985] Zharkov, V. N., V. V. Leontjev, and V. A. Kozenko (1985), Models, figures, and gravitational moments of the Galilean satellites of Jupiter and icy satellites of Saturn, Icarus, 61, 92–100, doi:10.1016/00191035(85)901575.
SUPPLEMENTARY INFORMATION
Enceladus’s and Dione’s floating ice shells
supported by minimum stress isostasy
Contents\@mkbothCONTENTSCONTENTS
\@afterheading\@starttoc
toc
List of Figures \@mkbothLIST OF FIGURESLIST OF FIGURES
\@afterheading\@starttoc
lof
List of Tables \@mkbothLIST OF TABLESLIST OF TABLES
\@afterheading\@starttoc
lot
Appendix S.1 Nonhydrostatic deviations of Enceladus and Dione
The degreetwo shape and gravity ratios are defined by [Zharkov et al., 1985; Nimmo et al., 2011]
(S.1) 
where and denote the cosine real harmonic coefficients of degree and order of the shape and nondimensional gravity potential, respectively. If the body is in hydrostatic equilibrium, to first order in the flattening, whatever the density stratification. Secondorder corrections change the shape and gravity ratios of a hydrostatic body by up to a few percent with respect to 10/3, with a slight dependence on the internal structure. The degreetwo shape and gravity ratio resulting from Eqs. (4)(5) are
(S.2) 
If the body is homogeneous, for Enceladus and for Dione. Cassini observations yield for Enceladus (Tables S.2 and S.4) and for Dione (Tables S.3 and S.5). To second order in the flattening, the nonhydrostatic shape deviations (with respect to the homogeneous body) are thus for Enceladus and for Dione, while the nonhydrostatic gravity deviations are for Enceladus and for Dione. The density stratification decreases these values by less than one percent ( and for Enceladus if ).
Appendix S.2 Gravity data (Tables s.2 & s.3)
Two gravity solutions are available for Enceladus [Iess et al., 2014]. SOL1 includes the five coefficients of degree two and the zonal coefficient of degree three, i.e. the least set of parameters able to fit the data at the noise level and compatible with the topography. SOL2 includes the five coefficients of degree two and the seven coefficients of degree three. Given the uncertainties, SOL2 does not contain more physical information than SOL1. On the other hand, the errors on SOL1 cannot reflect all the biases due to the sparse data (only three flybys) or resulting from the assumptions involved in spacecraft trajectory modeling (discussed in the Supplementary Material of Iess et al. [2014]). In particular, the central value of changes from SOL1 to SOL2 by more than while changes by more than and becomes compatible with a zero value within (the shift in is negligible). The error bars being much larger on SOL2, the shift of may be meaningless. Nevertheless we retain both gravity solutions in order to assess the impact of gravity field uncertainties on interior modeling. Table S.2 gives the coefficients differing from zero at (at least for SOL1): , , and . Table S.3 gives the coefficients of the preliminary gravity solution for Dione [Hemingway et al., 2016].
The geoid, to first order in the figure of equilibrium, is given by , except for the two components affected by the rotation and by the permanent static tide [Zharkov, 2004; Zharkov and Gudkova, 2010]:
(S.3) 
The secondorder geoid is given by Eq. (50) of Zharkov and Gudkova [2010] but secondorder corrections are only a few meters.
Appendix S.3 Shape data (Tables s.4 & s.5)
The shape of the midsize Saturnian satellites has been determined by fitting an ellipsoid [Thomas et al., 2007, 2016; Nadezhdina et al., 2016] and by directly estimating the spherical harmonic components of the topography [Nimmo et al., 2011]. For Enceladus and Dione, we use the shape of Nimmo et al. [2011] (denoted TOPA), which provides all harmonic coefficients of degree two and three. For Enceladus, we also use the ellipsoidal shape of Thomas et al. [2016] (denoted TOPB), which is based on more recent data, from which the harmonic coefficients of degree two can be derived. TOPB is less flattened than TOPA and its coefficient follows the geoid (see Table S.4). The new ellipsoidal shape of Nadezhdina et al. [2016] has a polar flattening () close to the one of TOPA but its equatorial flattening () is nearly hydrostatic (similarly to TOPB).
We use here spherical harmonic coefficients in unnormalized form (as in Iess et al. [2014]; McKinnon [2015]) so that denote the coefficients of (error bars on are not correct in Iess et al. [2014] but are correctly given in McKinnon [2015]). Thomas et al. [2016] give estimates for the semimajor axes and the mean radius of the bestfitting ellipsoid with error bars; the mean radius is defined as the radius of the sphere of equal volume; its value in Thomas et al. [2016] should be (P. Thomas, private communication). The degreetwo harmonic coefficients of the shape, to first order in the flattening, are computed from
(S.4) 
Appendix S.4 Basic equations of Bayesian inversion
The result of the Bayesian inversion is the posterior probability density function, i.e. the conditional probability for the parameters given the data [Tarantola, 2005]:
(S.5) 
represents the prior information on while is the likelihood function associating a probability to the data given parameter values. is related to by Eqs. (1) to (5) and (7), though we replace in practice Eqs. (4)(5) by the complete secondorder solutions for the figure of equilibrium. We assume that the gravity data () are normally distributed and uncorrelated (Tables S.2 and S.3). The model parameters are the densities (, ) and thicknesses (, ) of the ice shell and ocean (the radius and density of the core are derived parameters). The prior information on is described by uniform probability density functions (ranges given in Section 3.4) subjected to the constraint that the shell thickness at the south pole is not negative (Fig. S.1 shows the resulting prior distributions).
If the data includes the shape coefficients (), the isostatic shape coefficients must be treated as additional model parameters. We choose instead to include the shape coefficients with their uncertainties in the model itself which becomes probabilistic. This is similar to the problem of fitting a straight line to data with errors in both coordinates. According to Section 4.8.2 of Gregory [2005], the likelihood function is the same as if the shape was exactly known except that the standard deviation of the gravity data is replaced by
(S.6) 
where the and are the uncertainties on the gravity (Tables S.2 and S.3) and shape (Tables S.4 and S.5) coefficients, respectively. The likelihood function reads
(S.7) 
where
(S.8) 
it being understood that .
We simulate with a MetropolisHastings sampler [Tarantola, 2005], generating 400000 samples for each inversion and retaining every fourth sample in order to reduce the correlation within the samples.
The burnin length is equal to 1000.
From the selected samples, we compute for each parameter the probability density function (Fig. S.1), the mean value and the Bayesian confidence intervals (Tables S.6 and S.7).
The correlation between parameters is illustrated by Fig. 2 showing the twoparameter Bayesian confidence regions for the ice shell and ocean thicknesses, while Fig. S.2 does the same for the core radius and density.
Appendix S.5 Bayesian inversion results (Tables s.6 & s.7)
Table S.6 shows the results of several inversions of gravityshape data for Enceladus, corresponding either to different datasets or to different interior models. SOL1 and SOL2 refer to the gravity data of Table S.2 while TOPA and TOPB refer to the shape data of Table S.4. The SOL1/TOPA column shows the results of the reference case discussed in the text (Fig. S.1 shows the prior and posterior distributions). The SOL2/TOPA column shows the results with the alternative gravity solution SOL2. In that case, the mean value of the shell thickness () is slightly higher than with SOL1 () because the larger data uncertainties lead to a distribution extending to much larger values of the shell thickness, though the mode is (Fig. S.1).
The next two columns of Table S.6 show the results with different shape data. The shape given by TOPB is only of degree two, thus we compare inversions for SOL1/TOPA and SOL1/TOPB without the third degree. The last column of Table S.6 shows the results for an alternative interior model (POROUS) in which the nearsurface porosity is a free parameter (see Section 3.4). As expected, the estimated shell thickness is significantly higher because the surface topography contributes less to the gravity coefficients. Fig. S.3 shows the dependence of the estimated shell thickness on crustal porosity (for that figure, inversions were made for given values of the porosity).
Table S.7 shows the results of the inversion of gravityshape data for Dione. SOL1 and TOPA refer to the gravity and shape data of Tables S.3 and S.5. The important effect of the shape uncertainties is illustrated by an inversion in which the shape is precisely known (TOPB).
The distribution of the local shell thickness at the colatitude and longitude is obtained by computing for each sample
(S.9) 
where and are the unnormalized associated Legendre functions. For Enceladus’s model SOL1/TOPA, the parameters in the righthand side have mean values and uncertainties given by , , , and . These parameters are correlated so that the error on can only be estimated from the complete distribution.
Appendix S.6 Secondorder error on the isostatic model (Figure s.4)
In this paper, gravity and shape are split into hydrostatic and isostatic components. Secondorder corrections are included in the figure of equilibrium but not in the isostatic model. In order to estimate the error due to this approximation, we compute the secondorder gravitational potential generated by the total shape (hydrostatic plus isostatic) of all density interfaces (core/ocean, ocean/crust, and surface). The secondorder error on the isostatic model is estimated by the difference .
Consider first a spherically symmetric body of radius stratified into layers of uniform density. Next, add topographic relief to the interface of radius between two layers which have a positive density contrast . Powers of the relief can be expanded in real spherical harmonics: , where the first subscript denotes cosine () or sine () components. The gravitational potential caused by this relief at the surface is given by [Wieczorek and Phillips, 1998]
(S.10) 
where the finite amplitude (FA) contribution of the shape reads
(S.11) 
If the only nonzero shape coefficients are , the finite amplitude terms associated to these coefficients are given to second order () by
(S.12) 
For Enceladus, the secondorder error on the isostatic model is 5 to 10% of , 10 to 26% of , and 10 to 20% of (Fig. 3a). However, the degreetwo isostatic components are a small fraction of the total gravity coefficient (less than 10% of and 3% of for Enceladus). Therefore the secondorder error on the isostatic model is less than 1% of the degreetwo total gravity coefficients and less than 20% of the degreethree gravity coefficient, which is comparable to the uncertainty on the data. Finite amplitude corrections are smaller for Dione.
Appendix S.7 Comment on a paper by Cadek et al. (2016)
Outline. A recent paper by Cadek et al. [2016] claims that Enceladus’s lithosphere provides flexural elastic support so that the shell can be thinner than predicted by classical isostasy. We show here that their claim is wrong.
In that paper, the classical Airy relation between the surface and compensating topography is replaced by one depending on an additional parameter: the elastic thickness of the lithosphere (see Eq. (S.13) below). Because of this new free parameter, the gravityshape problem becomes underdetermined (see their Fig. 2) and the average shell thickness must be constrained by independent information provided by librations. Cadek et al. tune so that the shell thickness matches the one required in the libration model (). Thus it is not surprising that they draw the same conclusions about a thin south polar crust as do the libration models with nonhydrostatic crustal boundaries [Van Hoolst et al., 2016].
We will make two points:

The flexuralisostatic model used in Cadek et al. [2016] is wrong: it is based on an equation for top loading, but the degree of compensation is computed with a formula valid for bottom loading. If flexural isostasy results from top loading, should be about instead of as stated in the paper.

Supporting more than of topography with a thick lithosphere stretches belief. More generally, whatever the exact thickness of the lithosphere (as long as it is thin), the degreetwo deflection required to provide elastic support induces stresses that lead to lithospheric failure.
We round off the section by discussing the problems raised by isostasy in thin shell theory.
Flexuralisostatic model. In our notation, the flexuralisostatic model of Cadek et al. [2016] (their Eq. (1)) reads
(S.13) 
where are the radius, gravitational acceleration, and nonhydrostatic topography at the crustocean boundary, respectively. The parameter is the degree of compensation for elastic flexure which varies between 0.75 and 0.9 in Cadek et al. [2016] ( in the limit of Airy isostasy). With the approximation , Eq. (S.13) reduces to
(S.14) 
where should be interpreted as the deflection of the shell due to the surface topography . Eq. (S.14) is the wellknown compensation equation for top loading [Turcotte et al., 1981] (without geoid effects), and we expect a choice for that is consistent with top loading [Turcotte et al., 1981]:
(S.15) 
where and ( is Poisson’s ratio). The bending resistance is neglected here so that there is only one parameter measuring the rigidity of the shell [Turcotte et al., 1981]:
(S.16) 
where is Young’s modulus and is the elastic thickness of the lithosphere.
Instead, Cadek et al. define by a formula which was proposed for the modeling of dynamic topography (Eq. (25) of Kalousová et al. [2012]), a kind of bottom loading:
(S.17) 
where is defined by
(S.18) 
The correct formula for the shell deflection due to bottom loading would be (Eq. (A10) of McGovern et al. [2002]):
(S.19) 
where should be interpreted as the thickness of the bottom load and as the shell deflection. Clearly, there is a contradiction in using Eq. (S.14) with defined by Eq. (S.17). The correct way to proceed would be to use Eq. (S.14) with the degree of compensation for top loading, Eq. (S.15).
If and , and is much larger than (Fig. S.5). In particular, corresponds to , which is the elastic thickness favored in Cadek et al. [2016]. By contrast, using the correct formula for the degree of compensation () yields . It is hardly credible that more than of topography can stand on such a thin lithosphere without breaking it (see below).
Flexural isostasy with top loading is not the only possibility.
Alternatively, Eqs. (S.17) and (S.19) can be used if the deflection of the shell (or crust) results from a subsurface density anomaly.
Beware that the thickness of the bottom load is not identical to the topography of the crustocean boundary: the former is measured with respect to the deflected crustocean boundary whereas the latter is measured with respect to the spherical surface coinciding with the crustocean boundary before deflection.
If top and bottom loads are both present, the flexure problem cannot be solved without specifying another free parameter, which is the ratio of subsurface to surface loading (Eq. (B5) of McGovern et al. [2002]).
Lithospheric failure. If we put aside the question of how is related to , the flexuralisostatic model of Cadek et al. [2016] (Eq. (S.13)) can be seen as a parametric relation between surface topography and bottom topography, in which is a blackbox parameter representing the effect of flexural isostasy. Following McKinnon [2013], we now show that a thin lithosphere (less than thick) cannot support degreetwo topography without breaking up, except in the limit of full isostasy with negligible deflection.
The degreetwo zonal deformation of a spherical shell generates tangential membrane stresses given by the VeningMeinesz equations (Eqs. (4.2)(4.3) of Melosh [2011]). We quantify the stress in the lithosphere by the maximum differential stress
(S.20) 
where is the crustal shear modulus and is the shell flattening. The latter is related to the harmonic component of the shell deflection () by . We neglect here the equatorial flattening because Enceladus’s nonhydrostatic deformation is mostly zonal.
If compensation results from top loading, is equal to the deflection of the shell given by in Eq. (S.13), which is much larger in magnitude than the topographic load itself because of the nearcomplete compensation ().
The hydrostaticisostatic decomposition yields so that the maximum differential stress is .
If compensation results from bottom loading, is equal to the shell deflection ( in Eq. (S.19)), so that the maximum differential stress is .
Both estimates are much larger than the tensional strength of intact cold ice (1 to ) [Schulson and Duval, 2009].
The lithosphere actually fails at lower stress levels because of its finite brittle strength, as discussed in detail in McKinnon [2013].
In conclusion, the nonhydrostatic degreetwo topography on Enceladus cannot be supported by lithospheric flexure.
The only way out consists in distributing the loads on the top and bottom of the crust so as to minimize the shell deflection and the resulting membrane stresses.
But this state of full isostasy does not depend on the elastic thickness of the lithosphere.
Isostasy in thin shell theory. Applying thin shell theory to isostasy is problematic for three reasons. First, thin shell theory is not valid if the shell is thicker than ten percent of the shell radius. Thus one must suppose a priori that the ice shells of Enceladus and Dione are less than and thick, respectively. Enceladus’s ice shell seems to satisfy, but just barely, the thin shell threshold whereas Dione’s ice shell does not. Thin shell theory could be valid for Dione (or Enceladus) if the lithosphere is much thinner than the whole ice shell. It is not obvious, however, how the subsurface loading applied on the crustocean boundary is transmitted to the lithosphere itself (see third reason below). Second, the geoid cannot be included in a consistent way: the topographic loads are measured with respect to the geoid which is affected, in turn, by the deflection of the shell as well as by the surface and subsurface loads themselves. This effect is crucial at long wavelengths. Third, the standard (firstorder) formalism of thin shell theory does not take into account the difference in area between the inner and outer shell surfaces. It is thus not possible to introduce ad hoc factors and as was done by Cadek et al. (compare Eqs. (S.13) and (S.14)), since factors of similar magnitude are neglected elsewhere in thin shell theory. This problem must be solved by working with a higherorder approximation of thin shell theory which is not worth the effort in comparison with the minimum stress isostasy developed in our paper.
Parameter  Symbol  Enceladus  Dione  Unit 

Surface radius  252.1  561.4  km  
Gravitational parameter  
Bulk density  1610  1478  
Surface gravity  0.113  0.232  
Period  1.37  days  
Rotational parameter   
Coeff.  SOL1  SOL2  

Coeff.  SOL1  

TOPA  TOPB  GEO1  GEO2  Unit  
km  
km  
km  
km  
5  3  m  
2  4  m  
1      m  
5  5        
Error bars on the differences .  
Error bars assuming that errors on add quadratically. 
TOPA  GEO1  Unit  
km  
km  
km  
km  
8  m  
9  m  
12      
Error bars on . 
Parameter  Prior range  SOL1/TOPA  SOL2/TOPA  SOL1/TOPA  SOL1/TOPB  POROUS 

(km)  
(km)  
(km)  
()  
()  
() 
Parameter  Prior range  SOL1/TOPA  SOL1/TOPB 

(km)  
(km)  
(km)  
()  
()  
() 