Enceladus’s and Dione’s floating ice shells
supported by minimum stress isostasy
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 present-day ocean within Dione.
Paper accepted for publication in Geophysical Research Letters
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 long-wavelength 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 three-fourths of the surface radius, but the crust-ocean 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 second-order 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 degree-three 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. Degree-three 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.  and McKinnon : the second-order 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 second-order tidal-rotational effects [McKinnon, 2015]. A thick shell is hardly compatible with the south polar activity, raises the issue of shell-core contact at the equator precluding isostasy, and does not match degree-three 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 long-wavelength 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.  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.  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 large-scale isostasy on Earth: thermal (Pratt) and Airy isostasy are dominant in oceanic and continental crusts, respectively, while the long-wavelength geoid is explained by mantle convection. In planetology, isostasy invariably resorts to the equal-mass 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  and Dahlen ) that the only physically meaningful isostatic prescription consists in minimizing crustal deviatoric stress in a self-consistent elastic-gravitational theory.
3.1 Hydrostatic-Isostatic Decomposition
Before explaining the three methods particular to this paper, we recall the principle of the hydrostatic-isostatic 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. , we split the harmonic components into hydrostatic and isostatic (nonhydrostatic) components:
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:
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 , we do not compute the figure of equilibrium with the method of Tricarico , which is based on nested ellipsoids, because an ellipsoidal stratification is forbidden for heterogenous hydrostatic bodies (Hamy-Pizzetti theorem, see Moritz ). 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 second-order equations for a multilayer body yields complicated expressions for the degree-two 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 :
where are the second-order corrections and is the nondimensional rotational parameter (Table S.1). Fluid Love numbers of two- and three-layer bodies are given in Zharkov  and Zharkov and Gudkova , respectively. We computed from using the second-order relations between shape and gravity (Eq. (49) of Zharkov and Gudkova ). Our formalism for the harmonic coefficients of degree two is equivalent to the one of Zharkov and Gudkova , 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 . 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 , 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 self-gravitating spherically symmetric body submitted on the one hand to a surface load and, on the other, to an internal load located at the crust-ocean 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 elastic-gravitational 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 [Greff-Lefftz et al., 2010]
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):
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  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 elastic-gravitational 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 crust-ocean boundary is given by
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 Metropolis-Hastings 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 three-layer 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 salt-rich ice is denser. For the three-layer 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 four-layer 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].
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 fluid-like (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 ice-dominated layer ( as in McKinnon ). 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. , in Van Hoolst et al. ). For a rigorous comparison, we compute the librations from the probability distribution over the parameters inferred from our gravity-shape 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.
Comparison with librations can pinpoint data biases and constrain modeling choices. First, Enceladus’s degree-three gravity favors a thinner shell than degree-two coefficients [McKinnon, 2015] (Fig. 1 and Table S.6). Degree-two 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 degree-two shape could be responsible for the disagreement: while the new ellipsoidal shape of Thomas et al.  does not affect much the results (Table S.6), the recent ellipsoidal shape of Nadezhdina et al.  predicts more degree-two compensation and thus a thinner shell. Second, we did not allow for a lot of surface porosity in our three-layer model, but we can easily do it with a four-layer 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 leading-trailing 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 [Crow-Willard 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 second-order 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. , 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 long-wavelength 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 planetary-scale isostasy.
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).
- Besserer et al.  Besserer, J., F. Nimmo, J. H. Roberts, and R. T. Pappalardo (2013), Convection-driven 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.  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.
- Crow-Willard and Pappalardo  Crow-Willard, 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  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  Dahlen, F. A. (1982), Isostatic geoid anomalies on a sphere, J. Geophys. Res., 87, 3943–3947, doi:10.1029/JB087iB05p03943.
- Dahlen and Tromp  Dahlen, F. A., and J. Tromp (1999), Theoretical Global Seismology, Princeton University Press.
- Greff-Lefftz et al.  Greff-Lefftz, 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.1365-246X.2009.04490.x.
- Gregory  Gregory, P. C. (2005), Bayesian Logical Data Analysis for the Physical Sciences, Cambridge University Press, Cambridge.
- Hammond et al.  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.  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.  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. Moragas-Klostermeyer, and R. Srama (2015), Ongoing hydrothermal activities within Enceladus, Nature, 519, 207–210, doi:10.1038/nature14262.
- Iess et al.  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.  Jaeger, J. C., N. G. W. Cook, and R. W. Zimmerman (2007), Fundamentals of Rock Mechanics, 4th edition, Blackwell.
- James et al.  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  Jeffreys, H. (1970), The Earth, Cambridge University Press, Cambridge.
- Kalousová et al.  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.1365-246X.2012.05539.x.
- Kirchoff and Schenk  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  Lambeck, K. (1980), Estimates of stress differences in the crust from isostatic considerations, J. Geophys. Res., 85, 6397–6402, doi:10.1029/JB085iB11p06397.
- Lambeck  Lambeck, K. (1988), Geophysical Geodesy, Clarendon Press, Oxford.
- McGovern et al.  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  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  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  Melosh, H. J. (2011), Planetary Surface Processes, Cambridge University Press, Cambridge.
- Moritz  Moritz, H. (1990), The Figure of the Earth, Wichmann, Karlsruhe.
- Multhaup and Spohn  Multhaup, K., and T. Spohn (2007), Stagnant lid convection in the mid-sized icy satellites of Saturn, Icarus, 186, 420–435, doi:10.1016/j.icarus.2006.09.001.
- Nadezhdina et al.  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., XLI-B4, 459–462, doi:10.5194/isprs-archives-XLI-B4-459-2016.
- Nimmo and Bills  Nimmo, F., and B. G. Bills (2010), Shell thickness variations and the long-wavelength 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.  Nimmo, F., B. G. Bills, and P. C. Thomas (2011), Geophysical implications of the long-wavelength topography of the Saturnian satellites, J. Geophys. Res., 116, E11001, doi:10.1029/2011JE003835.
- Ojakangas and Stevenson  Ojakangas, G. W., and D. J. Stevenson (1989), Thermal state of an ice shell on Europa, Icarus, 81, 220–241, doi:10.1016/0019-1035(89)90052-3.
- Park et al.  Park, R. S., A. S. Konopliv, B. Bills, N. Rambaux, J. Castillo-Rogez, 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.  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 low-degree shape of Mercury, Geophys. Res. Lett., 42, 6951–6958, doi:10.1002/2015GL065101.
- Porco et al.  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.  Postberg, F., J. Schmidt, J. Hillier, S. Kempf, and R. Srama (2011), A salt-water reservoir as the source of a compositionally stratified plume on Enceladus, Nature, 474, 620–622, doi:10.1038/nature10175.
- Sabadini and Vermeersen  Sabadini, R., and B. Vermeersen (2004), Global Dynamics of the Earth, Kluwer Academic Publishers, Dordrecht.
- Sambridge and Gallagher  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/978-90-481-8702-7_192.
- Schulson and Duval  Schulson, E. M., and P. Duval (2009), Creep and Fracture of Ice, Cambridge University Press, Cambridge.
- Tarantola  Tarantola, A. (2005), Inverse Problem Theory and Methods for Model Parameter Estimation, SIAM, Philadelphia.
- Thomas et al.  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.  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  Tricarico, P. (2014), Multi-layer hydrostatic equilibrium of planets and synchronous moons: Theory and application to Ceres and to solar system moons, Astr. J., 782, 99, doi:10.1088/0004-637X/782/2/99.
- Turcotte et al.  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.  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  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  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  Wieczorek, M. A., and M. T. Zuber (2004), Thickness of the Martian crust: Improved constraints from geoid-to-topography ratios, J. Geophys. Res., 109, E01009, doi:10.1029/2003JE002153.
- Yin and Pappalardo  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  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  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.  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/0019-1035(85)90157-5.
Enceladus’s and Dione’s floating ice shells
supported by minimum stress isostasy
List of Figures \@mkbothLIST OF FIGURESLIST OF FIGURES
List of Tables \@mkbothLIST OF TABLESLIST OF TABLES
Appendix S.1 Nonhydrostatic deviations of Enceladus and Dione
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. Second-order 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 degree-two shape and gravity ratio resulting from Eqs. (4)-(5) are
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 ).
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. ). 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]:
The second-order geoid is given by Eq. (50) of Zharkov and Gudkova  but second-order corrections are only a few meters.
The shape of the mid-size 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.  (denoted TOPA), which provides all harmonic coefficients of degree two and three. For Enceladus, we also use the ellipsoidal shape of Thomas et al.  (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.  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. ; McKinnon ) so that denote the coefficients of (error bars on are not correct in Iess et al.  but are correctly given in McKinnon ). Thomas et al.  give estimates for the semi-major axes and the mean radius of the best-fitting ellipsoid with error bars; the mean radius is defined as the radius of the sphere of equal volume; its value in Thomas et al.  should be (P. Thomas, private communication). The degree-two harmonic coefficients of the shape, to first order in the flattening, are computed from
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]:
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 second-order 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 , 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
it being understood that .
We simulate with a Metropolis-Hastings sampler [Tarantola, 2005], generating 400000 samples for each inversion and retaining every fourth sample in order to reduce the correlation within the samples.
The burn-in 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 two-parameter Bayesian confidence regions for the ice shell and ocean thicknesses, while Fig. S.2 does the same for the core radius and density.
Table S.6 shows the results of several inversions of gravity-shape 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 near-surface 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 gravity-shape 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
where and are the unnormalized associated Legendre functions. For Enceladus’s model SOL1/TOPA, the parameters in the right-hand 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 Second-order error on the isostatic model (Figure s.4)
In this paper, gravity and shape are split into hydrostatic and isostatic components. Second-order 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 second-order gravitational potential generated by the total shape (hydrostatic plus isostatic) of all density interfaces (core/ocean, ocean/crust, and surface). The second-order 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]
where the finite amplitude (FA) contribution of the shape reads
If the only nonzero shape coefficients are , the finite amplitude terms associated to these coefficients are given to second order () by
For Enceladus, the second-order error on the isostatic model is 5 to 10% of , 10 to 26% of , and 10 to 20% of (Fig. 3a). However, the degree-two isostatic components are a small fraction of the total gravity coefficient (less than 10% of and 3% of for Enceladus). Therefore the second-order error on the isostatic model is less than 1% of the degree-two total gravity coefficients and less than 20% of the degree-three 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.  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 gravity-shape 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 flexural-isostatic model used in Cadek et al.  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 degree-two 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.
Flexural-isostatic model. In our notation, the flexural-isostatic model of Cadek et al.  (their Eq. (1)) reads
where are the radius, gravitational acceleration, and nonhydrostatic topography at the crust-ocean boundary, respectively. The parameter is the degree of compensation for elastic flexure which varies between 0.75 and 0.9 in Cadek et al.  ( in the limit of Airy isostasy). With the approximation , Eq. (S.13) reduces to
where should be interpreted as the deflection of the shell due to the surface topography . Eq. (S.14) is the well-known 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]:
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]:
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. ), a kind of bottom loading:
where is defined by
The correct formula for the shell deflection due to bottom loading would be (Eq. (A10) of McGovern et al. ):
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. . 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 crust-ocean boundary: the former is measured with respect to the deflected crust-ocean boundary whereas the latter is measured with respect to the spherical surface coinciding with the crust-ocean 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. ).
Lithospheric failure. If we put aside the question of how is related to , the flexural-isostatic model of Cadek et al.  (Eq. (S.13)) can be seen as a parametric relation between surface topography and bottom topography, in which is a black-box parameter representing the effect of flexural isostasy. Following McKinnon , we now show that a thin lithosphere (less than thick) cannot support degree-two topography without breaking up, except in the limit of full isostasy with negligible deflection.
The degree-two zonal deformation of a spherical shell generates tangential membrane stresses given by the Vening-Meinesz equations (Eqs. (4.2)-(4.3) of Melosh ). We quantify the stress in the lithosphere by the maximum differential stress
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 near-complete compensation ().
The hydrostatic-isostatic 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 .
In conclusion, the nonhydrostatic degree-two 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 crust-ocean 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 (first-order) 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 higher-order approximation of thin shell theory which is not worth the effort in comparison with the minimum stress isostasy developed in our paper.
|Error bars on the differences .|
|Error bars assuming that errors on add quadratically.|
|Error bars on .|