The diurnal libration and interior structure of Enceladus Icarus

ContentslistsavailableatScienceDirect

Icarus

journalhomepage:www.elsevier.com/locate/icarus

The diurnal libration and interior structure of Enceladus

Tim Van Hoolst

^a,b,∗

, Rose-Marie Baland

^c,a

, Antony Trinh

^a

a Royal Observatory of Belgium, Ringlaan 3, Brussels B-1180, Belgium

b Instituut voor Sterrenkunde, KU Leuven, Celestijnenlaan 200D, Leuven B-3001, Belgium

c Georges Lemaître Centre for Earth and Climate Research, Earth and Life Institute, Université catholique de Louvain, Place Louis Pasteur, 3, Louvain-la-Neuve B-1348, Belgium

a rt i c l e i n f o

Article history:

Received 22 December 2015 Revised 11 May 2016 Accepted 13 May 2016 Available online 24 May 2016 Keywords:

Enceladus Interiors

Rotational dynamics

a b s t r a c t

WedetermineconstraintsontheiceshellandoceanofEnceladusfromtheobservedlibrationatorbital periodbyassessingtheeffectsofuncertaintiesinthesize,density,rigidity,andviscosityoftheinternal layersand ofthe non-hydrostaticstructure onthelibration. Theobserved librationamplitude implies thattheaveragethicknessoftheiceshellisbetween14kmand26kmandthattheoceanis21kmto 67kmthick.

© 2016ElsevierInc.Allrightsreserved.

1. Introduction

The gravitationalforce exerted by Saturn on Enceladus raises tides that are thoughtto be responsible for the observed geyser activity at the south pole of Enceladus (Porco, 2006; Spencer et al., 2006). Although the precise causeand mechanismof cry- ovolcanism isnotwell understood,thevariationoftheplumein- tensity and its correlation with the orbital position of Enceladus is strongly indicative of a tidal mechanism driving the geysers (Hedmanetal.,2013;Hurfordetal.,2007;Porcoetal.,2014).The detection of sodium-salt-rich ice grains emitted from the plume suggeststhat the plumeis connectedto asubsurface salty water reservoir incontactwithsilicaterocks (Hsu etal., 2015;Postberg etal.,2009;2011).Thesubsurfaceoceanmightexistonlybeneath the regionofactivitynearto thesouthpolarregion (Tobieetal., 2008),althoughrecentgravityandtopographydatacannotdistin- guish betweena localand a globalsubsurface ocean (Iesset al., 2014;McKinnon,2015).

The gravitational interaction with Saturn also periodically changes therotationofEnceladus,whichmightbeco-responsible fortheobservedvariationsinplumeactivity(Nimmoetal.,2014a).

The diurnallibrations ofEnceladus,representingthe variationsat orbital period in the rotation angle with respect to the steady changefora constantrotationrate, haverecentlybeenaccurately determined by a detailed analysis of images of Enceladus taken by the Cassini imaging system (Thomas et al., 2016). Although

∗ Corresponding author at: Royal Observatory of Belgium, Ringlaan 3, B-1180 Brussels, Belgium. Tel.: 32 23730668.

E-mail address: tim.vanhoolst@oma.be (T.V. Hoolst).

thelibrationamplitudeof0.120± 0.014 degree(2

σ

^{) istoosmall}

fordiurnallibration toplay an importantrole in thevariation in tidalstress neededtoexplaintheplumevariability(Nimmoetal., 2014b), it provides unique direct insight into the interior struc- turethrough thesensitivityof librationto the existenceand cer- tain characteristicsof aliquid globalinternallayer. Thomas etal.

(2016)showedthatthelibrationamplitudeprovestheexistenceof aglobaloceanbeneathanicysurface.

Hereweperformanextendedstudyofthediurnallibrationsof Enceladus inorderto determinethe best constraintson iceshell and ocean, in particular on the mean thickness of the ice shell.

We study the effect of the rigidity and visco-elastic behavior of theiceshellonthelibrations andassesstheuncertaintythey in- troduceintheinterpretationofthelibrationamplitudeintermsof thedepthto theocean.We firstusean approachin whichEnce- ladusisconsideredto beinhydrostaticequilibrium. Next,weuse theobservedtopographyandgravitationalfieldofEnceladustode- velop non-hydrostatic models of the interior structure. Two dif- ferent end-member assumptions are considered for the shape of the interfaces betweenthe core, the ocean andthe ice shell. By meansofthosemodels,wetheninvestigatethelibrationsfornon- hydrostatic models of Enceladus that are consistent withgravity and topography data (Iess et al., 2014; Nimmo etal., 2011). We show that although many uncertainties remain aboutEnceladus’

interiorstructure,themeanthicknessoftheiceshellcanbewell constrainedbythelibrationamplitude.

2. Librationwithoutglobalsubsurfaceocean

ThegravitationaltorqueofSaturn exerted onEnceladus forces therotationofEnceladustovaryslightlyduringtheorbitalmotion.

http://dx.doi.org/10.1016/j.icarus.2016.05.025 0019-1035/© 2016 Elsevier Inc. All rights reserved.

SinceEnceladusislikelyina 1:1spin-orbitresonance,thetorque willhaveperiods commensuratewiththeorbital period,whichis equaltotherotationperiodofEnceladus,andwillforcelibrations withthe sameperiods. Inaddition,Enceladus haslibrations with longer periods related to the variations in the torque associated withperiodicperturbationsoftheorbitofEnceladusduetoDione.

Thetwo mainlong-periodlibrations haveperiods ofabout4and 11years(Rambauxetal.,2010).SincetheobservationsofThomas etal.(2016)haveaperiodequaltotheorbitalperiod,weneglect thoselong-periodlibrationshere.

Thepolargravitationaltorqueleadingtovariationsintherota- tion ratecan be expressed in terms of the external gravitational field coefficients of Enceladus (e.g. Dehant and Mathews, 2015), whichhavebeendetermined uptodegree 3byIessetal.(2014). Sixgravitational coefficients have been estimated, all five of de- greetwo (J₂, C₂₁, S₂₁, C₂₂, S₂₂) and the zonal component ofde- greethree(J₃), butonlythetorqueonC₂₂ needs tobetakeninto account.The polar torque does not depend on the zonal gravita- tionalcoefficientsbecauseofthegeometricsymmetryofthezonal harmonicswithrespecttothepolaraxis.Thecontributionsofthe tesseralcoefficientstothepolartorquecanbeneglectedsincethey areapproximatelyproportionaltotheverysmallobliquityofEnce- ladus.Theobliquityisnotknown,buttheoreticalcalculationspre- dictvaluesbelow10⁻⁵ radians (Balandetal., 2014).The S₂₂ sec- torialcoefficientisabout70timessmallerthantheC₂₂ coefficient (Iessetal.,2014),anditscontributiontothegravitationaltorqueis thereforealsoneglected.

Thetorqueforceslibrationattheorbitalperiodandatthesub- harmonic periods of it. The latter will also be neglected as the torqueatthenthsubharmonic issmallerthan that atorbitalpe- riodbyafactoreⁿ(e.g.VanHoolstetal.,2008),wheree=0.0047 istheorbital eccentricity(Thomasetal., 2016). Weintroduce the librationangle

γ

=

φ

− M,where

φ

^{istherotationangleofEnce-}

ladusandMthemeananomaly.Whenthelibrationatorbitalpe- riodof1.37daysiswrittenas

γ

=gsin(^M+

π

),thelibrationam- plitudeg_solidforanentirelysolidandelasticEnceladuscanbeex- pressedas

g_solid=6B˜− ˜A C˜

en² n²−

ω

²f

≈ 6eB˜− ˜A

C˜ , ⁽¹⁾

where

ω

fisthefrequencyofthefreelibration

ω

f=n



3

(

^B− A

)

C˜

kf− k2

k_f ⁽²⁾

and B˜=B(^kf− 5k₂/6)/kf, A˜=A(^kf− 5k₂/6)/kf and C˜= C+4k₂n²R⁵/(^9G) ^{(Van Hoolst et al., 2013). Here A} < B are theequatorialprincipal moments ofinertia,C thepolarprincipal momentof inertia, R=252.1 km the radius andG the universal gravitational constant. The mean motion of Enceladus’ orbit is denoted by n, k₂ is the classical potential Love number, and k_f the fluid Love number. The libration amplitude of the forced librationatorbitalperioddependsonafreefrequency,similarlyto aharmonicoscillator,butwe assumethatthefreelibrationmode isdampedasthedampingtimescaleisextremely shortcompared totheageoftheSolarSystem(see,e.g.Tiscarenoetal.,2009).

Forthenumericalevaluation,weusethedegree-two,order-two gravitationalcoefficientC₂₂=(^B− A)/(^4MER²)=(¹⁵⁴⁹.8± 15.6)× 10⁻⁶,with mass M_E=1.0794× 10²⁰ kg, and the estimate of the meanmomentofinertiaI=(^A+B+C)/3=(⁰.335± 0.005)^MR2(2 sigma)derivedfromRadau’sequation andthedegree-twogravity field of Enceladus (Iess et al., 2014). For an entirely solid Ence- ladus,k₂ issmallandapproximately1.5× 10⁻⁴ dependingsome- whatonthedensityandrigidityprofile (Balandetal.,2016),and k_f=(^B− A)/(^qMER²)≈ 0.989 withq=(^n2R3)/(^GME)^{the ratioof} the centrifugal acceleration to the gravitational acceleration. The

free libration period is then 5.81days (5.77 days fora homoge- neous model, Thomas et al., 2016), far enough from the forcing period of 1.37 daysfor the approximation in Eq.(1) to be valid.

Elasticity,whichhasbeenincludedhere,lengthensthefreeperiod by lessthan 0.01%since the tides of a solid Enceladus are very smallwitha radial tidaldisplacement oftheorder ofacentime- ter.Alocalsubsurfaceocean willbeable toincreasethe tides,at leastlocally,butcannotsubstantiallychangethefreelibrationpe- riodsincetheshelllibratestogetherwiththedeepersolidinterior.

The libration amplitude for a solid and elastic Enceladus is about 132 m at the equator, four times smaller than the ob- servedvalueof(528± 62)m(Thomasetal.,2016).Thomasetal.

(2016)arguedthattheincompatibilityoftheobservationwiththe librationamplitudeofasolidandrigidEnceladus isevidenceofa subsurfaceocean.Elasticitydoesnotchangethatconclusion since elasticitydecreasesthelibrationamplitude(bylessthanameter).

Any largerelasticity (smallerrigidity) thanassumedhere, forex- amplefora porouscoreassuggestedby Hsuetal.(2015),would imply an even smaller libration amplitudeand a strongerincon- sistency with the observations. Visco-elastic behavior of the ice wouldalso drainenergyfrom librationto deformation anddissi- pationandwillfurtherreducethelibrationamplitude.

Weconcludethatapartiallydecouplinglayermustexist,which allows for differential rotation between different layers of the satellite. Since measurements of the plume composition demon- strate that the water has been in contact with rocks (Postberg etal.,2011),librationstronglysuggeststhatEnceladusiscomposed ofasolidicyoutershellontopofaliquidoceanmainlycomposed ofwaterthatisincontactwithasolidcorecontainingsilicatesin its toplayers. Adetailedstudyof thelibrationof Enceladuswith anoceanispresentedinthenextsection.

3. librationwithaglobalsubsurfaceocean 3.1. Librationmodel

Inordertoassesstheeffectofvariousinteriorstructureparam- eters,such asthe thicknessofthe iceshell,on thelibration am- plitude, andto constrainsome of theseparameters basedon the observed libration, we construct a large set of interior structure modelsofEnceladuswithaglobalsubsurfaceoceanforwhichwe calculatethelibration.AllmodelssatisfythetotalmassM_Eandra- diusRandhaveameanmomentofinertiaI/(M_ER²)between0.325 and0.345, covering themoment ofinertia estimate of Iessetal.

(2014)withinthetwosigmaprecision([0.33,0.34])andtherange [0.328,0.333]obtainedby McKinnon(2015).Therangesofvalues considered for the radii anddensities of the core, ocean andice shellare giveninTable1. Wechoosethe densityoftheiceshell to be in the interval [900, 1000] kg/m³. Thisrange is somewhat larger than the range of [900, 950] kg/m³, often considered for pure ice including some degree of porosity in studies of the ice shells of icy satellites (e.g. Nimmo et al., 2014a), to include also theeffectofthepresenceofcontaminants.Fortheocean,wecon- sideradensityrangeof[950,1100]kg/m³toincludeawiderange of salinityvalues and temperatures. We refer to Sharqawy etal.

(2010) for seawater densities on Earth and the online density ta- bleatweb.mit.edu/seawater/wherewaterdensityvaluesaregiven between943.1kg/m³and1096.2kg/m³ fortemperaturesbetween 0and120°C,atmosphericpressure,andsalinities between0and 120g/kg.

ThereferencemodelsoftheinteriorstructureofEnceladusare spherically symmetric. Tobe able to studythe libration, we cal- culatethedegree-twoshapeofthethreeinterfaces(surface,inter- face betweenthe ocean andthe iceshell, andinterface between the core and the ocean). As explained above, libration only de- pendsonthedegree-twoshapeandwethereforeconsiderthatthe

Table 1

Ranges of values of the parameters of the interior structure of Enceladus for models that satisfy the observed libration amplitude for different modeling hypotheses: rigid models in hydrostatic equilibrium (HE), elastic models in hydrostatic equilibrium with shell rigidity in the range [10 ⁹ , 5 × 10 ⁹ ] Pa, elastic non-hydrostatic models with compensation at the ocean-shell interface, and elastic non-hydrostatic models with compensation at the core-ocean interface.

The mean thickness of the shell h s and ocean h o and the mean radius of the core r c are given in km, the densities in kg m ⁻³ . The initial sampling is given in the last column. The sampling interval is 25 kg m ⁻³ for the shell and ocean density, 1 km for shell thickness and 5 km for core radius. For the non- hydrostatic case with a non-equipotential ocean-shell interface, the density of the ocean is at the high end since the density difference must be sufficiently large in order to explain the gravity data. The libration amplitudes are calculated by solving Eqs. (45) and (46) of Van Hoolst et al. (2013) .

Quantity Rigid HE Elastic HE Non-HE ocean Non-HE core Sampling

h s [16, 24] [14, 24] [14, 24] [18, 26] [5, 72]

h o [23, 65] [24, 67] [29, 66] [21, 63] [1, 82]

r c [170, 205] [170, 205] [170, 205] [170, 205] [165, 205]

ρs [90 0, 10 0 0] [90 0, 10 0 0] [90 0, 10 0 0] [90 0, 10 0 0] [90 0, 10 0 0]

ρo [950, 1100] [950, 1100] [975, 1100] [950, 1100] [950, 1100]

ρc [2158, 2829] [2158, 2829] [2166, 2829] [2158, 2829] [2156, 2946]

ρo −ρs [0, 200] [0, 200] [75, 200] [0, 200] [0, 200]

interfaces can be describedby triaxial ellipsoids. The radial dis- tancer_i(

θ

^,

φ

^{)ofapointoninterfacei with}co-latitude

θ

^andlon-

gitude

φ

^{canbeexpressedas} r_i

( θ

^,

φ )

^=ri⁰



1−2

3

α

iP₂

(

^cos

θ )

⁺¹₆

β

iP₂²

(

^cos

θ )

^cos2

φ 

, ⁽³⁾

where r⁰_i =a_i(¹−

α

i/3−

β

i/2) ^{is the mean radius of interface i} given by the spherically symmetric reference model,

α

i=[(^ai+ b_i)/2− ci]/[(^ai+b_i)/2] isthepolarflattening and

β

i=(^ai− bi)/ai

the equatorial flattening. Here, a_i, b_i, and c_i are the semi-major axes of the ellipsoid and P₂(cos

θ

^{) and P}2²(^cos

θ

) ^{are the associ-} atedLegendrefunctions. Theinterfaceofthecoreisdenotedbyc for core, the interface between the ocean andthe shell by o for ocean, and theouter surfaceby s for shell.We usedifferent as- sumptionstocalculatetheflattenings.First,inSections3.1and3.2, we considerthat Enceladus isin hydrostaticequilibriumandcal- culatethepolarandequatorialflatteningoftheinterfacesbysolv- ingClairaut’sequationforthesynchronouslyrotatingsatellite(see e.g. Hubbard et al., 1984; Moritz, 1990;Van Hoolstet al., 2008).

Second,we usethetopographyandgravitydata(Iessetal., 2014;

Nimmo etal., 2011; Thomaset al.,2016) atdegree twoto calcu- late thenon-hydrostatic flattenings ofthe interfaces fortwo dif- ferentassumptionsonthecompensationmechanism(Section3.3).

TheflatteningvaluesareshowninFig.1.

Forlargeicysatellitessuch asTitanandtheGalileanicysatel- lites, thelibration amplitudeofthe shell gs canbe well approxi- matedby

gs≈4eK₃ n²Cs

(4)

(VanHoolstetal.,2013).HereCsisthepolarmomentofinertiaof theshellandK₃istheeffectivestrengthofthetotaldiurnalforcing torqueexertedontheshell.Itessentiallyconsistsofthesumofthe gravitationaltorqueexertedby Saturnandoftheassociatedpres- suretorqueexerted bytheocean, bothonthestaticshapeofthe shellofEnceladus andthetidalbulgeoftheshell(seeVan Hoolst etal.,2013forexplicitexpressions),othertermscontributingabout 1%orless(formodelswithashellthicknessabove15km).

ForthemodelsoftheinteriorstructureofEnceladusconsidered hereandforaflatteningprofileconsistentwithhydrostaticequilib- rium,theshellhasaconsiderablysmallermomentofinertiathan the entire satellite, between6% and 63% of the total momentof inertia for shells withincreasing thickness from5 kmto 72 km.

The total torque onthe shell isbetween61% and79% of theto- taltorqueontheentiresatellite.Therelativereductioninthetotal torqueontheshellwithrespecttothetorqueontheentiresatel- liteisthereforemuchlessthantherelativereductioninthepolar moment ofinertia of theshell withrespect to that ofthe entire satellite, in particularfor thinshells. As a result, the comparison of Eqs. (1)and (4) showsthat the libration of the shell above a

globaloceanis largerthan the librationifEnceladus were not to haveaglobaldecouplingocean.Wenote thatinsharpcontrastto thelargeicysatellitesforwhichtidaldeformationcanreduce the planet torque on the shell by up to one order ofmagnitude for thinshells, the torques related tothe tidal deformation of Ence- ladusdecreasethetotaltorqueontheshellwithrespecttoarigid satelliteby11%atmostforthemodelswiththethinnestshellasa resultofthemuchsmallertidaldeformation.Forashellthickness of20km,thereductionisatmost3%.

The libration amplitude is further increased by a resonance witha freemode, inparticular forthinshells. Twofree libration modesexistfora satellitewitha globalsubsurfaceocean: (1)an in-phaselibrationof thesolid coreandshell, (2)anout-of-phase libration(Van Hoolstetal.,2013).Theperiodofthefirst modeis closetothat foran entirelysolid satellitesince theocean isthin (Fig.2),justifyingneglectingthesquareofthisfreefrequencywith respectto n²,asdone inEq.(4).The periodofthe second mode decreasesfrom about4.7daysto about1.5days,closetothe or- bitalforcingperiod,forshellthicknessesdecreasingfrom72kmto 5km(Fig.2).Wethereforeextendapproximation(4)forthelibra- tionamplitudeto

gs≈ 4eK3

Cs

(

ⁿ²− K1/Cs

)

^{. (5)}

Here,thesquare ofthesecond freefrequencyisapproximatedby K₁/C_s,asfollowsfromEq.(64)ofVanHoolstetal.(2013),seealso their Eq. (47) forthe definition of the couplingstrength K₁. The ratioofK₁/Cston² increasesfrom0.08to0.78forshellthickness decreasing from 72km to 5 km, showing that the amplitude of librationcanincreasebyafactorofupto4forthethinnestshells asaresultoftheresonancewiththesecondfreemode.

InFig.2,thelibration amplitudeis representedforallinterior structure models of Enceladus. The libration amplitude strongly dependson the thickness of the ice shell. Fora shell rigidity of 3.3× 10⁹ Pa (Shojiet al., 2013), only models witha mean shell thickness in the range [16,23] km are consistent with the ob- served libration amplitudeat the equatorialsurface ofEnceladus of[466,590]m(2

σ

uncertaintyofThomasetal.,2016),seeFig.2. Thomasetal.(2016) concludedthat theiceshellmust beon av- erage between 21 km and 26 km thick based on the libration modelofVanHoolstetal.(2009)forrigidsolidlayers,about4km thickerthanourresult.Thedifferenceisnotduetoelasticity,since therange forrigid solid layers is extended by only 1km to [16, 24] km,butmainly dueto the differentsampling ofthe interior structure models. If the ocean would havea highsalinity (

ρ

o >

1100kg/m⁻³ andupto

ρ

^o=1200kg/m⁻³)dueto exchangewith the underlying core or if the shell would have a high porosity (850kg/m⁻³<

ρ

s<900kg/m⁻³), thethicknessoftheshellcould beupto26km,avaluecorrespondingtotheupperlimitgivenby Thomas etal.(2016),who consider very low densities ofthe ice

Fig. 1. Equatorial flattenings of the surface (top), ocean-shell interface (middle) and core (bottom) for the interior structure models of Enceladus as a function of the shell thickness. Three cases are considered: hydrostatic equilibrium (in grey), a non-equipotential ocean-shell interface (left), and a non-equipotential core-ocean interface (right).

The scale of the vertical axes is not always the same. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

shell.Thelower limitweobtainforthemeanthicknessoftheice shellisalsoduetoourwidersamplinginwhichwealsoconsider smalldensitydifferencesbetweentheoceanandtheshell.Wenote thatthesurfaceflatteningsinThomas etal.(2016)correspondto those of the observed topography, and are therefore not hydro- static.Sincetheflatteningsoftheinternalinterfacesdependonthe non-hydrostatic surface flattenings in the calculation of Thomas etal.(2016)(see theirEqs.(15) and(16)),theinternal flattenings cannotbe trulyhydrostaticflattenings.Byconvention, hydrostatic flatteningsdescribethedegree-twoshapeofequipotentialsurfaces relatedtoself-gravity,thecentrifugalpotentialandthestatictidal potential (Hubbard etal., 1984; Moritz, 1990). As such they can- not depend on non-hydrostatic surface flattenings. Since the ob- servedequatorial surfaceflattening ofThomas etal.(2016) (

β

s= 0.0187± 0.0007) almostoverlapstherange[0.0178,0.0191]ofhy- drostatic flattenings for our set of interior structure models, the effectontheiceshellthicknessisneverthelessbelow1km.

Because of the differentflatteningsof thesurface andthe in- terfacebetweenthe ocean andiceshell, thethickness ofthe ice shellisvariable.Fig.3showstheiceshellthicknessasafunction

ofco-latitudeandlongitudeforamodelwithameanshellthick- nessof20kmthatsatisfiestheobservedlibrationamplitude.Since wehereconsiderhydrostaticequilibriumandtheshellisthin,the flatteningsofthe surfaceandoftheinterface betweenthe ocean andthe shell arevery similar (see Fig. 1). As aresult, the varia- tions intheshellthicknessaresmallandlessthan1km(Fig.3).

This mayseem surprisingsince both theactivity andthe gravity andtopographydataindicatethattheshellismuchthinneratthe southpolethanelsewhere.WewillshowinSection3.3.thatmuch largervariationsintheshellthicknessandamuchthinnershellat the south poleis possible when we relax the assumption of hy- drostaticequilibriumandtakeintoaccountthegravityandtopog- raphydata.

3.2. Effectoficerigidityandviscosity

InFig. 2,the rigidityof theiceshell is fixedto 3.3× 10⁹ Pa, closetothevalueofHelgerud(2009)atzeropressureandzerode- greeCelsius,andtherigidityofthecoreistakentobe50× 10⁹Pa.

Therigidityofthecorehasanegligibleeffectonthelibrationdue

Fig. 2. Amplitude of libration of the ice shell and periods of the two free libration modes as a function of shell thickness. Shell rigidity is 3.3 × 10 ⁹ Pa and the rigidity of the core is 50 × 10 ⁹ Pa. We here assume that the shape of the interfaces be- tween the different layers is determined from the reaction of a hydrostatic satellite to the centrifugal and static tidal potentials. The amplitude calculated with the ap- proximated solution (4) differs from the amplitude represented here and obtained by solving the full set of equations governing libration as determined in Van Hoolst et al. (2013) by less than 0.9%. The period of the in-phase mode has a nearly con- stant value close to that for an entirely solid Enceladus (5.81 days). For shell thick- nesses decreasing from 72 km to 5 km, the period of the out-of-phase mode period decreases from about 4.7 days to about 1.5 days, close to the orbital forcing period of 1.37 days (red line), allowing a resonant amplification of the libration amplitude which increases with decreasing shell thickness. For thick ice shells, the libration amplitude is close to the libration of a solid Enceladus (132 m). Only models with a shell thickness in the range [16, 23] km are consistent with the observed libration amplitude at the equatorial surface of Enceladus of [466, 590] m (blue region). The yellow marker ( • ) denotes the model used in Figs. 3 and 4 , which has a mean shell thickness of 20 km. Models for non-hydrostatic shapes and other shell rigidity and viscosity values show a qualitatively similar dependence on the mean thickness of the ice shell. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

tothemuchsmallertidaldisplacementofthecorebelowthesub- surfaceoceanthanoftheiceshell.Hereweconsidernominalshell icerigidity values between10⁹ Pa and5× 10⁹ Pa.Lowervalues fortherigidityarenotimpossiblebutrequiresubstantialporosity oftheice(e.g.Wahretal.,2006).Thelibrationamplitudeincreases byabout10%(50m)overthisrange(Fig.4),fromlowtohighice rigidity.Forrigidityvalueswithinthenominalrangethethickness estimate ofthe ice shell onlymarginally changes to [14, 24] km (Table1).Shellrigiditiesbelow10⁸Paarenotconsistentwiththe observedlibrationamplitudesincethepredictedamplitudeforall models isthen belowthelowest limitof theobservations(shells thinner than 5 kmwould be required to resonantly increase the amplitudetotheobservedrange).

Viscosity of the ice reduces the libration amplitude (Jara- Orué and Vermeersen, 2014). We extend the libration model of VanHoolstetal.(2013)toincludevisco-elasticdeformationbycal- culatingthetidaldeformationandLovenumbersforavisco-elastic rheology. We describe the complex,frequency-dependent rigidity

μ

^{softheshellforaMaxwellmodelby}

μ

^s

(

ⁿ

)

=i n

μη

μ

⁺ⁱⁿ

η

^{, (6)}

where i is the imaginary number,

μ

^{the real rigidity of the ice}

shell, and

η

^{the iceviscosity(e.g. Tobieetal., 2005; Vermeersen}

andSabadini,2004).Sinceviscositystronglydependsontempera- ture,wetake alargerviscosity(largerthan10¹⁵Pas)forthecold top layer than forthe hotter bottom layer of the ice shell (total temperaturedifference ofabout180K).We assumethat thebot- tomhalfoftheiceisconvectingandhasanearlyconstanttemper- ature,whereas in thetop halfheat is transportedby conduction.

Thebottomiceisconsideredtohaveaviscosityclosetothemelt-

Fig. 3. Map of the local thickness of the ice shell for hydrostatic equilibrium (top), a non-equipotential ocean-shell interface (middle), and a non-equipotential core- ocean interface (bottom) for an interior structure model of Enceladus with shell thickness of 20 km, shell density of 925 kg m ⁻³ , ocean density of 1050 kg m ⁻³ , and core radius of 190 km, shown in Fig. 2 by a yellow marker ( • ). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

ingpointviscositybetween10¹³Pasand10¹⁵Pas(valuedepend- ing onthe grain size, Sotinet al., 2009). The librationamplitude doesnot depend on theexact value of thetop viscosity because thetopicebehaves elasticallyfortheMaxwellmodel considered.

TheMaxwelltimescale,separatingelasticfromfluidbehavior,cor- respondstothelibrationperiodfora viscosityof3.9× 10¹⁴Pas forthe nominalrigidity value. The libration amplitudeis almost constant forbottom viscositiesbetween 10¹⁴ Pas and 10¹⁵ Pa s, decreasesduetomorefluid-likebehaviorwithdecreasingbottom viscositybetween10¹⁴ Pasand10¹³ Pas andisalmostconstant forviscositiessmallerthan10¹³Pas(Fig.4).Thevariationinlibra- tionamplitudeislessthan20moverthenominalviscosityrange considered,wellbelowtheobservationalprecisionofThomasetal.

(2016).Weconcludethat theuncertaintyintheiceviscositydoes

8.0 8.5 9.0 9.5 10.0 10.5 11.0 250

300 350 400 450 500 550

Log(ice shell rigidity/Pa)

librationamplitude[m]

8.0 8.5 9.0 9.5 10.0 10.5 11.0

Log(ice shell rigidity/Pa)

12.0 12.5 13.0 13.5 14.0 14.5 15.0

Log(ice shell viscosity/(Pa s))

12.0 12.5 13.0 13.5 14.0 14.5 15.0

Log(ice shell viscosity/(Pa s))

Fig. 4. Amplitude of libration as a function of shell rigidity (blue) and bottom shell viscosity (red) for an interior structure model of Enceladus with shell thickness of 20 km, shell density of 925 kg m ⁻³ , ocean density of 1050 kg m ⁻³ , and core ra- dius of 190 km, shown in Fig. 2 by a yellow marker ( • ) Viscosity is neglected for the study of the dependence on shell rigidity and rigidity is set to 3.3 GPa for the study of the dependence on shell viscosity. Only bottom shell viscosity is varied, top viscosity is set to 10 ¹⁵ Pa s. Thick lines indicate results for nominal values of rigidity and bottom ice viscosity. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

notaffecttheinterpretationofthelibrationamplitudeintermsof theshellthicknessatthelevelofprecisionoftheobservations.

3.3.Effectofnon-hydrostaticstructure

MeasurementsofthegravityfieldandtopographyofEnceladus bytheCassini spacecraftduringcloseflybys showthat Enceladus deviates from hydrostatic equilibrium. The second-degree grav- itycoefficients(C₂₂ and J₂=(⁵⁴³⁵.2± 34.9)× 10⁻⁶) yield a ratio J₂/C22=3.51± 0.05 (Iess et al., 2014) that is significantly larger thanthevalueofabout3.25expectedforahydrostaticEnceladus rotatingat its currentrotation speed (McKinnon,2015). Also the ellipsoidalshape ofEnceladusasderivedfromlimbprofiles isin- compatiblewithhydrostaticequilibrium(Nimmoetal.,2011).

A non-hydrostaticshape ofthe outer surfaceaffectsthe libra- tionofEnceladusbychangingtheprincipalmomentsofinertiaof theiceshellandthereforealsothetorquesexertedontheshell.In additiontoa non-hydrostaticshapeof thesurface,also theother interfaces(core-oceanandocean-shell)canbeofanon-hydrostatic shape, further modifying the libration. We assume herethat the densitywithin each layer is homogeneous. The observed gravity coefficients J₂ and C₂₂ can then only be explained by assuming that atleast one of the internal interfaces has a non-hydrostatic shape.Wehereconsidertwoextremecases:either(1)weexplain the degree-two gravity coefficients by assuming that the bottom ofthe iceshell isnot an equipotential surfaceor (2)we assume thattheinterfacebetweenthecoreandtheoceandoesnotcorre- spondtoanequipotentialsurface.Theothersurfaceisassumedto bean equipotential surfaceof theself-gravitational potential,the centrifugalpotential and the statictidal potential. We determine the equatorial and polar flattenings of both surfaces by solving fourequations.Twoequationsexpressthatoneinterfaceisequipo- tentialinitsdegree-two,order-zeroanditsdegree-two,order-two components,theothertwoequationsexpressthatthedegree-two, order-zeroand thedegree-two, order-two gravity coefficientsare equalto their observed values.In the first caseof compensation

throughshellthicknessvariations,wehave 2

3

α

^crc⁰gc= 5

6n²

(

^r0c

)

²+8

π

^G

15

 α

^s

ρ

^s+

α

^o

( ρ

^o−

ρ

^s

)

+

α

^c

( ρ

^c−

ρ

^o

) 

(

^r0c

)

^{2, (7)}

1

6

β

cr_c⁰gc= 1

4n²

(

^r0c

)

²+2

π

^G

15 [

β

s

ρ

s+

β

o

( ρ

o−

ρ

s

)

+

β

c

( ρ

c−

ρ

o

)

^]

(

^rc⁰

)

², ⁽⁸⁾ wheregcisthelocalgravityatthecoreinterface,and

C20 =− 1 MER²

8

π

15

 ρ

s

(

^R5

α

s−

(

^ro⁰

)

⁵

α

o

)

+

ρ

^o

((

^ro⁰

)

⁵

α

^o−

(

^r0c

)

⁵

α

^c

)

+

ρ

^c

(

^r0c

)

⁵

α

^c



, (9)

C22 = 1 4MER²

8

π

15

 ρ

s

(

^R5

β

s−

(

^r0o

)

⁵

β

o

)

+

ρ

^o

((

^ro⁰

)

⁵

β

^o−

(

^r0c

)

⁵

β

^c

)

+

ρ

^c

(

^r0c

)

⁵

β

^c



, (10)

(seealsoBalandetal.,2014).Forthesecondcaseinwhichweas- sumethatthecore-oceaninterfaceisnotanequipotentialsurface, wecalculatetheflatteningsbysolving Eqs.(11) and(12) together with

2

3

α

or⁰_ogo= 5

6n²

(

^ro⁰

)

²+8

π

^G

15



[

α

s

ρ

s+

α

o

( ρ

o−

ρ

s

)

^]

(

^r0c

)

²

+

α

^c

( ρ

^c−

ρ

^o

) (

^rc⁰

)

⁵

(

^ro⁰

)

³



, (11)

1

6

β

^or0ogo= 1

4n²

(

^ro⁰

)

²⁺²

π

^G

15



[

β

^s

ρ

^s+

β

^o

( ρ

^o−

ρ

^s

)

^]

(

^rc⁰

)

²

+

β

^c

( ρ

^c−

ρ

^o

) (

^rc⁰

)

⁵ r_o⁰

)

³



, (12)

where g_o is the local gravity at the ocean-shell interface. The degree-twoshapeofthesurfaceistakenfromNimmoetal.(2011). These models are therefore compatible with both the observed long wavelengthtopographyandthe degree-twogravity. We also consideredthedegree-twoshapeofThomasetal.(2016),butonly quote here the results obtained for the values of Nimmo et al.

(2011) since those differmore from hydrostaticequilibrium than thevaluesofThomasetal.(2016).

We adapt the libration model ofVan Hoolstet al.(2013) de- velopedforhydrostaticequilibriumtonon-hydrostaticinteriorsby replacingtheflatteningsby theirnon-hydrostaticvaluesinallthe equations.Weonlyretainmodelsforwhichtheshell(ocean)thick- ness in any direction is at least2.5 km (1 km). In case(1), the non-hydrostaticstructure almostdoesnot changetheestimate of thethicknessoftheiceshell(Table1)although thelibrationam- plitudeofindividualmodelsthatsatisfytheobservedlibrationam- plitudecanbeupto160msmallerwithrespecttothecorrespond- ingmodelsinhydrostaticequilibriumfornominalicerigidityval- ues. The assumption of a non-equipotential core boundary shifts the interval ofpossible iceshellthicknesses compatible withthe observed librationamplitudeto larger valuesof[18, 26] km.The librationamplitudeoftheindividualmodelsisbetween49mand 129mlargerwithrespecttothehydrostaticmodels.

Fig.3showsthelargeeffectofthe non-hydrostaticmodelson thelocaliceshellthickness.Theresultsareillustratedforonein- teriorstructuremodelwithameanshellthicknessof20kmthat satisfiestheobservedlibrationamplitude.Thesouthpolarregionis onlyafewkilometerthickcomparedtothicknessesofover30km in the equatorial regions when it is assumed that the interface between the ocean and the shell is a non-equipotential surface.

Within this compensation model, the small densitycontrast be- tween the shelland theocean leadsto pronounced variationsin

thepositionoftheocean-shellinterfacewithaninterfacecloserto thesurfaceattopographic lowssuch asthesouthpole region,as can beseenfromthenegativesignoftheequatorialflatteningof theocean-shellinterface inFig.1(seealsoFig.5ofthethickness of theshell andocean inMcKinnon,2015). The localthicknesses showninFig.3arebasedonthedegree-twoshapesofthesurface andtheocean-shellinterface.Thisimpliesthatthenorthpolarre- gion has an equally thin shell. Higher-degreecomponents of the shapewillaffecttheseestimates,buthavenoinfluenceontheli- bration amplitude (see Section 3.1). When the interface between the core and the ocean is considered to be a non-equipotential surface, the variationsin thethickness ofthe iceshell are much smallerbecauseofthelargerdensitycontrastofthecorewiththe H₂Olayer,andoftheorderof3km.Onlyforanon-equipotential ocean-shell interfacecan theshell bebelow10 kmthick.Forthe nominalrangeiniceshellrigidity([1,5]GPa),thethicknessofthe shellatthesouthpoleisthenbetween2.5kmand17km.

4. Discussionandconclusions

We confirm and strengthen the results of Thomas et al.

(2016)thattheobservedlibrationamplitudeatorbitalperiodisev- idenceofaglobalsubsurfaceocean.Theexistenceofaglobalocean has also been suggested on the basis of gravity and topography data(McKinnon,2015)andisalsoconsistentwiththeobservedor- bitalmodulationoftheplumebrightnessduetotides(Behounková etal., 2015). Thosestudies predictan iceshellthicknessofabout 50kmormoreandathinoceanlessthan20kmthick. Thelibra- tion amplitudesuggests theopposite: athin shellwith a thicker ocean. Althoughthe gravityandtopography dataare very useful, andarealsousedinouranalysisofthelibration,itmustbekept inmind that theirinterpretation intermsof thethicknessofthe shellisbasedonseveraluncertainassumptions.Intheanalysesof Iesset al.(2014) andMcKinnon(2015), Airyisostasy is assumed.

Insuchamodel,themasswithinanycolumnofmatterabovethe compensationdepth isconsidered tobe thesame andduetolo- calthickeningorthinningoftheiceshell.Compensation,however, canalsobeduetoelasticanddynamicsupport(Wieczorek,2015).

Dependingonthephysicaloriginofthemassanomalies,theshell thicknessmaybe stronglydifferentfromwhat isassumedwithin theassumptionofAiryisostasy.Notealsothat thelowest-degrees twoandthreeareused,forlackofhigherdegreegravitationaldata, andthosecanbe morestronglyinfluencedbythedeeper interior thanhigherdegrees.

Table 1 summarizes the mean thicknesses of the shell and ocean for different modeling hypotheses and model parameter ranges.Theresultsshowthat thelibrationamplitudeimpliesthat the ice shell (ocean) is between 14 km and 24 km(24 km and 67 km) thick for both hydrostatic models and models that sat- isfy the observed degree-two gravity and topography by a non- equipotentialocean-shellinterface.Ifinstead,topographyandgrav- itywouldbe dueto anon-equipotential coreboundary, theshell (ocean) thickness is between 18 km and 26 km (21 km and 63 km). Other parameters of the interior structure are not con- strained by the libration. The thicknessof the shell atthe south poleisbetween2.5kmand17km(16.4kmand24.3km)forthe non-hydrostatic modelswithnon-equipotential ocean-shell (core- ocean)interface.Thesethicknessvariationsarebasedonlyonthe degree-two components of the shape of the surface and of the ocean-shellinterface.

Thickershellscouldbepossibleforlargerdensitycontrastsbe- tween the ice andthe ocean than considered here. For example, foroceandensitiesupto1200kg/m³,theupperlimitforthemean thicknessoftheiceshellincreasesby 3kmfortheinteriorstruc- turemodelsinhydrostaticequilibriumwiththereferencevalueof therigidityoftheiceshell.Alower-densityiceshellalsoincreases

themean shellthickness corresponding tothe observed libration amplitudesince itincreasesthelibration amplitude.Forexample, for an extreme case with a shell porosity of about 50% (densi- tiesconsidered between450and500kg/m³),themeanthickness of the shell can be up to 34 kmfor hydrostatic models. Such a highporositycan decreasetherigidity byaboutan orderofmag- nitude(Hessingeretal.,1996).Sinceadecreaseintheshellrigidity alsodecreasesthe librationamplitude(see Fig.4)the meanshell thicknesscorresponding to the observed libration amplitude will besmallerthan34kmwhenwe takeintoaccount thisreduction inicerigidity.Whenwetakearigidityvalueof0.33GPa,themean thicknessoftheiceshellcanbeupto31km.

Acknowledgments

RMBisfundedby aFSRgrant fromUCL,ATreceived financial supportoftheBelgianPRODEXprogrammanagedbytheEuropean SpaceAgencyincollaborationwiththeBelgianFederalSciencePol- icyOffice.

References

Baland, R.M. , Tobie, G. , Lefèvre, A. , et al. , 2014. Titan’s internal structure inferred from its gravity field, shape, and rotation state. Icarus 237, 29–41 .

Baland, R.-M. , Yseboodt, M. , Van Hoolst, T. , 2016. The obliquity of enceladus. Icarus 268, 12–31 .

Behounková, M. , Tobie, G. , Cadek, O. , et al. , 2015. Timing of water plume eruptions on Enceladus explained by interior viscosity structure. Nat. Geosci. 8, 601–604 . Dehant, V. , Mathews, P.M. , 2015. Precession, Nutation, and Wobble of the Earth.

Cambridge University Press, UK .

Hedman, M.M. , Gosmeyer, C.M. , Nicholson, P.D. , et al. ,2013. An observed correlation between plume activity and tidal stresses on enceladus. Nature 500, 182 . Helgerud, M.B. , Waite, W.F. , Kirby, S.H. , et al. , 2009. Elastic wave speeds and moduli

in polycrystalline ice lh, sl methane hydrate, and sll methane–ethane hydrate. J.

Geophys. Res. 114, B02212 .

Hessinger, J. , White, B.E. , Pohl, R.O. , 1996. Elastic properties of amorphous and crys- talline ice films. Planet. Space Sci. 44, 937–944 .

Hsu, H.W. , Postberg, F. , Sekine, Y. , et al. , 2015. Ongoing hydrothermal activities within enceladus. Nature 519, 207210 .

Hubbard, W.B. , Waite, W.F. , Kirby, S.H. ,1984. Planetary Interiors. Van Nostrand, New York .

Hurford, T.A. , Helfenstein, P. , Hoppa, G.V. , et al. , 2007. Eruptions arising from tidally controlled periodic openings of rifts on enceladus. Nature 447, 292 .

Iess, L. , Stevenson, D.J. , Parisi, M. , et al. , 2014. The gravity field and interior structure of enceladus. Science 344, 78–80 .

Jara-Orué, H.M. , Vermeersen, B.L.A. , 2014. The forced libration of Europa’s de- formable shell and its dependence on interior parameters. Icarus 229, 31–44 . McKinnon, W.B. , 2015. Effect of Enceladus’s rapid synchronous spin on interpreta-

tion of Cassini gravity. Geophys. Res. Lett. 42 (7), 2137–2143 .

Moritz, H. , 1990. The Figure of the Earth. Herbert Wichmann Verlag GmbH, Karl- sruhe .

Nimmo, F. , Bills, B.G. , Thomas, P.C. , 2011. Geophysical implications of the long-wave- length topography of the Saturnian satellites. J. Geophys. Res. 116, E11 . Nimmo, F. , Porco, C. , Mitchell, C. , 2014a. Tidally modulated eruptions on Enceladus:

CASSINI ISS observations and models. Astrophys. J. 148, 46 .

Nimmo, F. , Porco, C. , 2014b, Encyclopedia of the solar system. In: Spohn, T., Breuer, D., Johnson, T.V. (Eds.). Elsevier, pp. 851–860 .

Porco, C. , DiNino, D. , Nimmo, F. , 2014. How the geysers, tidal stresses, and thermal emission across the south polar terrain of enceladus are related. Astrophys. J.

148, 45 .

Porco, C.C. , Helfenstein, P. , Thomas, P.C. , et al. , 2006. Cassini observes the active south pole of enceladus. Science 311, 1393–1401 .

Postberg, F. , Kempf, S. , Schmidt, J. , et al. , 2009. Sodium salts in e ring ice grains from an ocean below the surface of enceladus. Nature 459, 1098–1101 . Postberg, F. , Schmidt, J. , Hillier, J. , et al. , 2011. A salt-water reservoir as the source of

a compositionally stratified plume on enceladus. Nature 474, 620–622 . Rambaux, N. , Castillo-Rogez, J.C. , Williams, J.G. , et al. , 2010. Librational response of

enceladus. Geophys. Res. Lett. 37, L04202 .

Sharqawy, M.H. , Lienhard, J.H. , Zubair, S.M. , 2010. Thermophysical properties of sea- water: A review of existing correlations and data. Desalination Water Treat. 16, 354?380 .

Shoji, D. , Hussmann, H. , Kurita, K. , et al. , 2013. Ice rheology and tidal heating of enceladus. Icarus 226, 10–19 .

Sotin, C. , Tobie, G. , Wahr, J. , et al. , 2009. Tides and tidal heating on Europa. In: Pap- palardo, R.T., McKinnon, W.B., Khurana, K.K. (Eds.), Europa. University of Arizona Press, Tucson, pp. 85–117 .

Spencer, J.R. , Pearl, J.C. , Segura, M. , et al. , 2006. Cassini encounters enceladus: Back- ground and the discovery of a south polar hot spot. Science 311, 1401–1405 . Thomas, P.C. , Tajeddine, R. , Tiscareno, M.S. , et al. , 2016. Enceladuss measured physi-

cal libration requires a global subsurface ocean. Icarus 264, 37–47 .

Tiscareno, M.S., Thomas, P.C., Burns, J.A., 2009. The rotation of janus and epimetheus. Icarus 204, 254–261. doi: 10.1016/j.icarus.2009.06.023 .

Tobie, G. , Cadek, O. , Sotin, C. , 2008. Solid tidal friction above a liquid water reservoir as the origin of the south pole hotspot on enceladus. Icarus 196, 642–652 . Tobie, G. , Mocquet, A. , Sotin, C. , 2005. Tidal dissipation within large icy satellites:

Applications to Europa and titan. Icarus 177, 534?549 .

Van Hoolst, T. , Baland, R.M. , Trinh, A. , 2013. The effect of tides on the longitudinal librations of large synchronously rotating icy satellites. Icarus 226, 299–315 . Van Hoolst, T. , Rambaux, N. , Karatekin, O. , et al. , 2009. The effect of gravitational

and pressure torques on titan’s length-of-day variations. Icarus 200, 256–264 .

Van Hoolst, T. , Rambaux, N. , Karatekin, O. , et al. , 2008. The librations, shape, and icy shell of Europa. Icarus 195/1, 386–399 .

Vermeersen, B. , Sabadini, R. , 2004. Global Dynamics of the Earth. Kluwer Academic Publishers, Dordrecht .

Wahr, J.M., Zuber, M.T., Smith, D.E., et al., 2006. Tides on Europa, and the thickness of Europa’s icy shell. J. Geophys. Res. 111, E12005. doi: 10.1029/2006JE002729 . Wieczorek, M. , 2015. Treatise on geophysics. 2nd edition. In: Gerald, S. (Ed.), Planets

and Moons, vol. 10. Oxford: Elsevier . Section 10.05.