Variations in rotation rate and polar motion of a non-hydrostatic Titan

ContentslistsavailableatScienceDirect

Icarus

journalhomepage:www.elsevier.com/locate/icarus

Variations in rotation rate and polar motion of a non-hydrostatic Titan

Alexis Coyette

^a,b,∗

, Rose-Marie Baland

^b

, Tim Van Hoolst

^b,c

a Earth and Life Institute, UCL, Place Louis Pasteur 3, Louvain-la-Neuve B-1348, Belgium

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

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

a rt i c l e i n f o

Article history:

Received 26 October 2017 Revised 29 January 2018 Accepted 2 February 2018 Available online 8 February 2018 Keywords:

Titan Libration Interior structure

a b s t r a c t

Observationoftherotationofsynchronouslyrotatingsatellitescanhelptoprobetheirinterior.Previous studiesmostlyassumethattheselargeicysatellitesareinhydrostaticequilibrium,althoughseveralmea- surementsindicatethattheydeviatefromsuchastate.Hereweinvestigatetheeffectofnon-hydrostatic equilibriumandofflowinthesubsurfaceoceanontherotationofTitan.Weconsiderthevariationsinro- tationrateandthepolarmotiondueto(1)thegravitationalforceexertedbySaturnatorbitalperiodand (2)exchangesofangularmomentumbetweentheseasonallyvaryingatmosphereandthesolidsurface.

Thedeviationofthemassdistributionfromhydrostaticitycansignificantlyincreasethediurnallibration anddecreasetheamplitudeoftheseasonallibration.Theeffectofthenon-hydrostaticmassdistribution islessimportantforpolarmotion,whichismoresensitivetoflowinthesubsurfaceocean.Byincluding alargespectrumofatmosphericperturbations,thesmallerthansynchronousrotationratemeasuredby Cassiniinthe2004–2009period (Meriggiolaetal.,2016)couldbeexplainedbytheatmosphericforc- ing.Ifourinterpretationiscorrect,wepredictalargerthansynchronousrotationrateinthe2009–2014 period.

© 2018ElsevierInc.Allrightsreserved.

1. Introduction

Changesinthespinofsolarsystembodiesprovideinsightinto their deep interior, and have, for example, recently been used to determine that the coreof Mercury isat least partiallyliquid (Margotetal., 2007)andhasa radiusofabout2000 km(Margot et al., 2012; Hauck et al., 2013; Rivoldini and Van Hoolst, 2013).

Inasimilarapproach,thelibrationsatorbitalperiodofEnceladus, detected onthebasisofCassiniopticaldata,showthat Enceladus has a global ocean below an about 20 km thick crust (Thomas etal.,2016; ˇCadeketal., 2016;VanHoolstetal.,2016),aconclu- sion that is also compatible with an examination of gravity and shape data(Beutheetal.,2016).Rotation variationscould alsobe usedtoprobetheinterioroflargeicysatellites,inparticularofTi- tanforwhichtherotation,gravityandshapehavebeenmeasured bytheCassinimission.

Theobservational determinationofrotationvariationsisbased on measurements of the shift in orientation in inertial space of Cassini radar images taken during different flybys. Cassini radar images taken between 2004 and 2009 have shown that the ro- tation period of Saturn’s moon Titan differs slightly from its or-

∗ Corresponding author at: Earth and Life Institute, UCL, Place Louis Pasteur 3, Louvain-la-Neuve B-1348, Belgium.

E-mail address: alexis.coyette@uclouvain.be (A. Coyette).

bitalperiod,althoughfirmconclusionshavebeendifficulttoobtain (Lorenz,2008;Stilesetal.,2008;2010).Thedeviationfromasyn- chronousrotationoftheiceshellis−0.024^◦± 0.018^◦/year(one-

σ

uncertainty,Meriggiolaetal.,2016).FurtheranalysisoftheCassini data, including data from the flybys performedsince 2009, may improvetheestimationsoftherotationalvariations.

Titan is assumed to be in a mean state of rotation called a Cassinistate (see e.g.Peale,1969). ItimpliesthatTitan isinsyn- chronous rotation and that the rateof precession of its rotation axisisclosetothatofthenormaltoitsorbit.Asaresult,thespin axis,thenormaltotheorbitandthenormaltotheLaplaceplane orinertialplanearenearlycoplanarandtheobliquity

η

^(theangle

betweentherotationaxisandthenormalto theorbital plane)is nearlyconstant.

Variations in the rotation rate of Titan around this mean synchronous rotation can occur for several reasons (Tokano and Neubauer, 2005; Van Hoolst et al., 2013; Richard et al., 2014).

First, Titan’s rotation changes with a period equal to Titan’s or- bitalperiodasaresultofthegravitationaltorqueexertedby Sat- urn. Current theoretical models (Van Hoolstet al., 2013; Richard et al., 2014) show that the amplitudes of the diurnal rotation variations are below the detection limit related to the posi- tion error of Cassini radar images of the order of one kilome- ter (Meriggiola et al., 2016). Second, dynamic variations in the atmosphere (and to a less extent in the hydrocarbon lakes) of https://doi.org/10.1016/j.icarus.2018.02.003

0019-1035/© 2018 Elsevier Inc. All rights reserved.

TitaninducechangesinTitan’srotationwithamain periodequal to half the orbital period of Saturn. Depending on the model of thedynamicsofthe atmosphereandonthe rotationmodelused (GoldreichandMitchell, 2010;Richardetal.,2014),themaximum displacementof a givensurfacespotat theequator with respect toitsequilibriumpositionwithoutvariationsofthelength-of-day (LOD)couldbe uptoabout1km.Themaximalrotationratevari- ations associated with LOD variations predicted for a large set ofinterior models is about 0.013°/year (Van Hoolst etal., 2013), andiscompatiblewiththeobservationinthe 2−

σ

^limit.Third,

free librations with periods of the orderof a year might be ex- cited by the atmosphere of Titan, and fourth, deviations from a Keplerianorbit on different timescales introduce additional vari- abilitymostlyatlongperiods (Richard etal.,2014;Yseboodt and Van Hoolst,2014). Moreover, Titan might not exactly occupythe 1:1spin-orbit resonance so that the rotationis non-synchronous (NSR)(GreenbergandWeidenschilling,1984),eventhoughtheob- served deviation from a synchronous rotation of the ice shell is compatiblewithazero-NSRinthe2−

σ

^limit.

In addition,the gravitationaltorqueexerted bySaturn andthe atmosphericandhydrologictorquesalsoleadtofluctuationsinthe orientationofthespinaxis.Thepositionofthespinaxischanges in two ways: with respect to inertial space (precession and nu- tations,e.g. Bills and Nimmo, 2008 and Balandet al., 2011) and withrespect tothe solid surface (polar motion).Polar motion of Titan due to its atmosphere andhydrocarbon lakes has recently beenstudied byTokanoetal.(2011)andCoyetteetal.(2016).By assuming Titan(including its ocean)to be inhydrostatic equilib- rium,the atmosphere forces the spin axis to followan elliptical pathwitha typical amplitudeof about500 m inthe y-direction and200 minthe xdirectionandamain periodequaltothe or- bitalperiod ofSaturn. These valuesapply toa shell thickness of about200km.Forthinnershells,boththeamplitudeandthemain periodofthe polarmotionsensitively dependonwhetheraforc- ingperiod iscloseto theperiodofa free wobblemode ofTitan.

Forshells thinnerthan80km,the amplitudeofthepolarmotion couldreachseveraltensofkmormore.

In the existing models for the librations andpolar motion of Titanwithaninternal ocean,itisassumedthatTitanisinhydro- staticequilibrium.Withintheobservationalerrors,theratioofthe degree-twogravitationalcoefficientsagreeswiththatexpectedfor ahydrostatic Titan,suggestingthat Titan isindeed closeto a re- laxedshape (Iess etal., 2012). However, this ratiois only anec- essarybutnotasufficientconditionforasynchronoussatelliteto bein hydrostaticequilibrium. The observedshape of thesurface, whichismoreflattenedatthepolesthanexpectedforhydrostatic equilibrium(Zebkeretal.,2009), aswellasthenon-zerodegree- threegravity signal(Iessetal.,2012)clearlyindicatesome depar- ture from hydrostatic equilibrium. Here, we will therefore relax the assumption of hydrostatic equilibrium in the rotation model byconsideringthe flowinthesubsurfaceocean andtheeffectof thenon-hydrostatic surfaceofTitan on theshape ofthe internal boundaries.

The planofthepaperisasfollows.InSection2,we useAiry- like models of isostasy to calculate the shape of the interfaces between different layers of Titan based on the observations of the degree-two gravity field (Iess et al., 2012) and topography (Zebkeret al., 2009). Section 3describes theextension of the li- brationtheory developedinVan Hoolstetal.(2013) andthe po- lar motion theory developed in Coyette et al. (2016) to non- hydrostaticsatelliteswiththe inclusionof aPoincaré flow inthe subsurfaceocean.AsinCoyetteetal.(2016),theprecessionofthe rotation axis of the solid layers is assumed to be known and is thereforenot solved jointly with the polar motion.In Section 4, numericalresultsarepresentedforthediurnalandseasonalforced librationsandpolarmotion.Weconsiderthesensitivityoftheli-

brationtovarious interiorstructure parameterssuch astherigid- ity and viscosity of the ice shell in the discussion (Section 5).

WealsostudyinSection5thepossibleobservationsoflibrations, LOD andpolarmotion. Finally,concludingremarks are presented inSection6.

2. Non-hydrostaticinternalstructureofTitan 2.1. Differentiationanddensityprofile

The mean density(

ρ

=1882± 1 kgm³,from ssd.jpl.nasa.gov) and the moment of inertia (MOI) of Titan (I/MR²=0.3431± 0.0004, SOL1a of Iess et al., 2012) indicate that Titan is differ- entiated into an ice-ocean layer, a mantle (denoted by “m” in the following) and a core(c)(see alsoGrasset etal., 2000; Sohl etal., 2003;Tobieetal.,2005;Fortesetal.,2007).Fromthelarge tidal Love number k₂=0.589± 0.075 (Iess et al., 2012) and the large value oftheobliquity ofTitan (BillsandNimmo, 2008; Ba- land et al., 2011; 2014; Noyelles and Nimmo, 2014; Boué et al., 2017) it is assumed that the ice-ocean layer is divided into a shell (s) and a subsurface global ocean (o). The moment of in- ertia is derived by assuming that Titan is inhydrostatic equilib- rium. Radau’s equation then allows determining the moment of inertia fromthe degree-twocoefficients ofthegravitational field.

Since the measured gravity and topography show that Titan de- viates from hydrostatic equilibrium, we consider a range ofMOI values.GaoandStevenson(2013)showedthatdeviationsfromthe hydrostaticequilibriumvalueof10%orevenmorearepossible.As inBalandetal.(2014),weconsiderthatthetruemomentofinertia ofTitanliesbetween0.30(Ganymede-like, e.g. Sohletal., 2003) and0.36(Callisto-like,e.g. Fortes,2012).Theupperlimitalsocor- respondstothevalueproposedbyHemingwayetal.(2013)based onananalysisofthegravitytotopographyadmittance.Thiswide range has the advantage of not excluding possible although less likelydensityprofiles andatthesame timeallows better explor- ingconsequencesofanon-hydrostaticmomentofinertia.

We assume that all layers are homogenous in composition.

Pressure induced density variations within a layer are expected to be less than a few % for Titan’s core and below 2% for Ti- tan’s iceshell.Since thesevariations aresmaller than theuncer- tainty on the mean density of the layers and density variations haveasmalleffectontheamplitudesoftheforcedlibrations(see, e.g. Dumberry etal., 2013forthe case oflibrations ofMercury), weassumeeachlayerktobeofuniformdensity

ρ

k.Aspherically symmetricreferencemodelofTitanisthendeterminedbyspecify- inginadditiontheouterradii(R_k)ofthelayers.

Sinceareferenceinteriorstructuremodelisuniquelyspecified by 8parameters (the densities

ρ

k andradii R_k ofall fourlayers) andonly two quantities(the surfaceradius R andthetotal mass M_T) are known, we consider physically plausible ranges of den- sitiesand interfaces radii asin Balandet al.(2014). The method used here to construct the interior models differs from Baland etal.(2014)aswewillchoosethedensitiesandradiirandomlyin- side therangesinstead ofpickingequally-spacedvalues.Thishas theadvantageofgivingafinerexplorationoftheparametersspace.

We consider the density of the ice shell to be between 920kg/m³,thedensityofpureiceIhatambientpressureandtem- perature(e.g.Sotinetal.,1998)and1065kg/m³,correspondingto contaminatediceand/ordense clathrates(e.g.Fortesetal.,2007).

Depending on the composition (ammonia-water or salted water) andpressure,theoceandensitymaytypicallyvaryfrom950kg/m³ to 1350 kg/m³ (Croft et al., 1988; Vance and Brown, 2013). The mantle can be made of high-pressure water ices possibly con- taminatedwithrocky materialswithdensitybetween1300kg/m³ (e.g.Grassetetal., 2000) andabout2000kg/m³.Thecorecanbe made of hydrated silicates or rocks mixed with ice and/or iron

withadensitybetween2500and4500kg/m³ (Sohletal., 2003).

Wealsoimposethedensitytodecreasewithincreasingradiusand consideraniceshellthicknesshsbetween10kmand200kmand aminimalthicknessof5kmforallthelayers.Themeanradiusof Titan’ssurfaceis2574.73± 0.09km(Zebkeretal.,2009).

The eight interior parameters (density and radius of the four layers)are thenrandomly chosen insidethe rangesdescribe here above.Ifaninteriormodeliscompatiblewiththemeandensity

ρ

andthechosen rangeofmomentofinertia (0.30<I/(MR²)<0.36), wecomputethetidalLovenumberofthismodelfromtheanalyt- icalsolutionforhomogeneousincompressiblelayers(Sabadiniand Vermeersen,2004).Forthestandardmodels,weusethefollowing rigidities:3.3GPafortheshell,4.6GPaforthemantleand100GPa forthecore.InSection5,weinvestigatethedependenceonrigid- ityandviscosity.IfthetidalLovenumberoftheinteriormodellies inthepossiblerangededucedfromtheobservationofthegravity field ofTitanby Cassini(k₂=0.589± 0.075,Iessetal.,2012),we keepthemodelasapossibleinteriorstructureofTitan.Ifnot,the model is rejected. In total, we retained 50,000 interior structure modelsofTitan.

2.2. Departurefromhydrostaticequilibrium

The observed shape of the surface of Titan does not corre- spond to the hydrostaticshape derived fromthe observed grav- ity field. In hydrostatic equilibrium, the surface of Titan derived fromtheobservedgravityfieldisexpectedtobeatriaxialellipsoid with radii a=2574.97 km, b=2574.66 km and c=2574.56 km (Iess et al., 2010). However, the observed surface shape corre- spondsto an ellipsoidwithradiigivenby a=2575.15± 0.02km, b=2574.78± 0.06 km and c=2574.47± 0.06 km (Zebker et al., 2009). Theobservedshape ofTitanisthereforemoreflattened at the poles than expected foran hydrostatic shape. In addition to theobserved shapeofTitan,thepresence ofanon-zerodegree-3 inthegravity signal(Iessetal.,2010;2012)alsoclearlyindicates theexistenceofadeparturefromthehydrostaticequilibrium.The non-hydrostaticdeparturemustbepartiallycompensatedinorder to be consistent with the gravity field. Similarly as Nimmo and Bills (2010),Hemingwayetal.(2013)andBalandetal.(2014),we use here a compensation acting through variations in the thick- ness of the ice shell (Airy-type model). This leads to a thinner shellatthepoles, inagreementwiththemaximumtidaldissipa- tionthere(Tobieetal.,2005;Beuthe,2013).Asecondanddifferent waytoresolvethenon-hydrostaticdeparture wouldhavebeento consider a compensation actingthrough variationsin densityin- steadofvariationsinthickness(Pratt-typemodel).Choukrounand Sotin(2012)havedevelopedsuchPratt-typemodelbyconsidering a 8%increase iniceshelldensitynearthe polesasaresultofan accumulationofdenseethane-richclathratesatthatplacesdueto the larger precipitationof ethane-dominated liquids. As obliquity resultsforPratt-typemodelsdo nosignificantlydifferfromthose obtainedwithAiry-typemodels(Balandetal.,2014),wewillhere onlyconsideranAiry-typemodel.

Weuseheretheobservednon-hydrostaticshapeoftheiceshell and consider that the mantle and core have hydrostatic shapes caused by the centrifugal, tidal and gravitational potentials. The flattenings of the interior and the mantle can be obtained from the classical theoremof Brunsthat statesthat the topographyof an equipotential surface is relatedto the ratio of the perturbing potential andthegravitationalacceleration(e.g.Moritz,1990,see alsoAppendixA.inBalandetal.,2014).Wethereforehave,forthe core(k=c)andmantle(k=m)

−2

3

α

kRk=



²⁰k

(

^Rk

)

g

(

^Rk

)

⁽¹⁾

1

6

β

kR_k=



²²_k

(

^Rk

)

g

(

^Rk

)

^{, (2)}

where

α

k=[(^ak+b_k)/2− ck]



[(^ak+b_k)/2] isthe polarflattening and

β

k=(^ak− bk)/ak istheequatorialflatteningoftheoutersur- face of layer k with principal axes a_k, b_k, c_k and ²⁰_k (^Rk) ^and

²²_k (^Rk) ^{are the degree-2 order-0 and degree-2 order-2 parts of} the total perturbing potential at the outer surface of thelayer k (withradiusR_k)(seeAppendixAforthemathematicalexpressions of²⁰_k (^Rk)^and²²_k (^Rk)^).

Tocompletethesystem,weusethesurfaceflatteningsderived fromtheobservedshapeofTitanfromZebkeretal.(2009)

α

^{s=19.2236× 10−5 (3)}

β

^{s=14.36381× 10−5, (4)}

and the observed C₂₀ and C₂₂ gravity field coefficients (C₂₀=

−33.599× 10⁻⁶ and C22=10.121× 10⁻⁶, from Iess et al., 2012), whichcanbeexpressedasafunctionoftheflatteningsofthedif- ferentlayers(seee.g.VanHoolstetal.,2013)as

C20=− 1 MTR²

8

π

15

 ρ

^s



R⁵_s

α

^{s− R5}o

α

^o



+

ρ

^o



R⁵_o

α

^{o− R5}m

α

^m



+

ρ

^m



R⁵_m

α

^{m− R5}c

α

^c



+

ρ

^cR5c

α

^c



(5)

C22= 1 4MTR²

8

π

15

 ρ

s



R⁵_s

β

s− R⁵o

β

o



+

ρ

o



R⁵_o

β

o− R⁵m

β

m



+

ρ

^m



R⁵_m

β

^{m− R5}c

β

^c



+

ρ

^cR5c

β

^c



, (6)

whereM_TisthemassofTitan.

From theseeight equations(Eqs. (1)and(2) forthe coreand the mantle and Eqs. ((3)–(6))), we obtain the flattenings of the core, mantle, ocean and shell (see Fig. 1). They can be com- pared with hydrostatic flattenings obtained through integration of Clairaut’s equation (e.g. Van Hoolst et al., 2008). The non- hydrostaticflatteningsofthecoreandmantledonotdeviatemuch fromthehydrostaticones.Incontrast,thenon-hydrostaticflatten- ingsoftheoceanandshellstronglydifferfromthehydrostaticflat- tenings. The non-hydrostaticpolar flattening

α

^{o isnegativesince}

itcompensatesforthelarge deviationofthesurfacefromhydro- staticity. Mostnon-hydrostaticinteriormodels are alsocharacter- izedby anegativeequatorialflattening

β

^{o.Ourflatteningsofthe}

oceandifferfromthoseofBalandetal.(2014)duetothefactthat we choose to constrain the tidal Love number k₂ in the 1−

σ

range.Usinga2−

σ

^{range,wewouldhaveobtainedthesamebe-}

haviorasinBalandetal.(2014) witha broaderrangefor

α

^oand

β

o in the presence of thick iceshells. By construction, the non- hydrostaticflattenings ofthe shell are the same forall theinte- riormodels.Forthehydrostaticmodels,theflatteningsoftheshell presentasmallvariabilityrelatedtotherangeinmeanmomentof inertia ofthe interiormodels. The strongdifference between the non-hydrostaticandhydrostaticflatteningsof theoceancan have alargeimpactontheMOI difference(^Bs− As),especiallyforthin shells(seeFig.2)andthereforeontherotationoftheshell.

3. Librationandpolarmotionforanon-hydrostaticTitan 3.1. Angularmomentumequations

We study the rotation of Titan from the angular momentum equation, or Euler–Liouville equation, that expresses the basic physicalprinciple thatthetime-derivativeoftheangularmomen- tumH oftherotatingbodyisequaltothetotalappliedtorque ^.

Fig. 1. Polar ( αon the left) and equatorial ( βon the right) flattenings of the core, mantle, ocean and shell (from top to bottom) of a large set of interior models of Titan described in Section 2.1 , as a function of the ice shell thickness. The hydrostatic models are represented in lighter blue and the non-hydrostatic ones in darker blue. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 2. MOI difference (B s − A s) (and normalized MOI difference (B s − A s) / (M T R ²) at the right) as a function of the equatorial flattening of the ocean βo for the hy- drostatic (lighter blue) and non-hydrostatic (darker blue) models of Titan. (For in- terpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

AsweconsiderinteriormodelsofTitancomposedoftwosolidre- gions separatedby the ocean,we havethree angularmomentum equations,onefortheshell(s),onefortheocean(o)andone for the interior (denoted by i and comprising the core andmantle).

Theangularmomentumequationsfortheshell,theoceanandthe interiorare(seeMathewsetal.,1991)

dHs

dt +

ω

^{s× Hs=}

^s (7)

dHo

dt +

ω

^{s× Ho=}

^{o (8)}

dHi

dt +

ω

^{i× Hi=}

i. (9)

Wechoosetoexpresstheequationsfortheshellandfortheocean intheBodyFrame(BF)oftheshellandtheequation fortheinte- riorintheBFoftheinterior.TheBFofasolidlayeristheframere- latedtotheprincipalmomentsofinertiaA_k<B_k<C_kofthatlayer.

HereH_kand

ω

k=(

ω

^xk,

ω

^yk,

ω

^zk)^{aretheangularmomentumandthe} rotation vector ofthe layerk,respectively. (x, y,z) isthe generic namegiventothecoordinatesoftheBFcorrespondingtothethree principal axesofasolid layer(shell orinterior).Theangularmo- mentaaredescribedinSection3.1.1.

IncomparisonwithrotationstudiesofMathewsetal.(1991)for the Earth or Dumberry and Wieczorek (2016) for the Moon, we herestudyasystemofonlythreeequations.Wedonotstudythe relativemotionoftheBFoftheinteriorwithrespecttotheBF of the shell,asinMathewsetal.(1991) (see their fourthequation).

Wealsodonotstudytheprecessionandnutationoftheshellro- tationaxisinspace, whichwouldcorrespondtothefifthequation introduced in Dumberry and Wieczorek (2016). As Euler’s kine- maticequations(e.g.Lambeck,1980)show,polarmotionisalways associatedwithnutation,andvice-versa.However,theseequations alsoshow that long-periodpolarmotionisassociatedwithmuch smaller nutation andthat long-periodnutation in space leadsto muchsmallerpolarmotionamplitude.Long-periodpolarmotionis thereforealmostdisconnectedfromnutation,makingitpossibleto studypolarmotionindependentlyfromnutation.Incontrast,diur- nalpolarmotionintheBFiscoupledwithnutationandisnotcon- sideredhere.Wetestedthatnutation onlysignificantlyinfluences polarmotionforperiodsshorterthanabouttenTitandays.These short-periodswillthereforenotbe consideredintheatmospheric forcingofpolarmotion,andalsointheexternalgravitationalforc- ingby Saturn.Wewillseethatthisamountssimplytoneglecting

theobliquityofthesolidlayers(seeSection3.1.1).Inasimilarap- proximation,Cassinistatemodelscanbedevelopedindependently byneglectingthepolarmotion(e.g.Balandetal.,2011,foran hy- drostaticsatellite).Librationsandlength-of-dayvariationsarenot affectedbyourchoiceofsettingobliquitiestozero.

Due to the synchronous rotation, the rotation rateof a solid layercan be writtenas

ω

^z_k=n+

γ

˙k wheren isthe meanmotion ofTitanand

γ

k(t) thelibration,orsmalldeviationfromtherota- tion angle of the layer k from the rotation angle determined by theconstant rotation raten.The two smallcomponents (

ω

^x_k,

ω

_k^y) describethepolarmotionofthat layer.Thefirsttwo components ofeachangularmomentumequation(Eqs.(7)–(9))willbeusedto studythepolarmotionofthesolidlayersofTitanwhilethelibra- tionsorlength-of-day(LOD)variationswillbeextractedfromtheir thirdcomponent.

The total torques k on the right hand side of Eqs. ((7)–(9)) are the sum of the gravitational torques ^e_k ^{and }_k^{g exerted on} layer k by Saturn and by the other internal layers of Titan, re- spectively,onlayerk. kalsoincludesthepressuretorque _k^pex- ertedbytheliquidoceanontheinterfaceswithlayerk.Thispres- suretorqueisdividedintoahydrostaticpressuretorque _k^p,Hand a non-hydrostatic pressure torque _k^{p,NH. In addition, the atmo-} sphere(andlakes) ofTitanalsoexertsatorque ^{Atm ontheshell.}

InSection3.2,wedetermineexpressionsforallthesetorques.

3.1.1. Angularmomenta

Theangularmomentumofanysolidlayerk canbewritten,in theBFofthatlayer,as

H_k=I_k

ω

k, (10)

where

I_k=

⎛

⎜ ⎝

A_k+c¹¹_k c¹²_k c_k¹³ c¹²_k B_k+c²²_k c²³_k c_k¹³ c²³_k C_k+c³³_k

⎞

⎟ ⎠

⁽¹¹⁾

istheinertiatensoroflayer kwithprincipalmoments A_k,B_kand C_k andincrementsc^{i j}_k (seeAppendix B). As all large icysatellites in synchronous rotation,Titan has a main ellipsoidal shape with momentsof inertiaA<B<C duetorotation andstatictides.The small inertia increments c^{i j}_k represent deviations from this ellip- soidalshapeoflayerkthatariseduetorotationalandtidaldefor- mationsofTitan.

Tothefirstorderinthesmallquantities

ω

_k^x,

ω

^y_k,

γ

˙kandc^{i j}_k,we have

Hk=

⎛

⎜ ⎝

Ak

ω

^x_k+nc¹³_k B_k

ω

^y_k+nc_k²³ C_k

(

ⁿ⁺

γ

^˙k

)

^+nc33k

⎞

⎟ ⎠

^{. (12)}

We divide the ocean into two parts: a bottom part (that we designatebythesubscriptob)thatisalignedwiththeinteriorand atoppart(subscriptot)alignedwiththeshell.Sincethetopand bottompartsoftheoceanarealignedwiththeshellandtheinte- rior,respectively, theprincipal momentsof thetop oceanAot,Bot

andCotaredefinedwithrespecttothesameprincipalaxesasthe shell while the principal moments of the bottom ocean A_ob,B_ob andC_obaredefinedwithrespecttothesameprincipalaxesasthe interior.ByusingtherotationmatrixR_BFfromtheBFoftheinte- rior tothe BF of the shell(see Eq. (C.12)), the total angularmo- mentumoftheoceanisthengiven,intheBFoftheshell,by Ho

BFs=

(

^Iot+RBFIobR⁻¹_BF

) ω

^o, (13) where

ω

^{ois therotation vector ofthe oceanwith respectto the} BF ofthe shell.Bytakinginto account thefact that Titan isin a

Cassini state (see Eqs. (12)–(15)of Coyette et al., 2016), we then find,correctuptothefirstorderinsmallquantities

ω

_k^x/yandobliq-

uity

η

k

Ho

BFs=Hot+H_ob+C_ob

⎛

⎜ ⎝

 ω

^x− n

 η

^sin

(

^M+

ω )  ω

^{y− n}

 η

^cos

(

^M+

ω )

0

⎞

⎟ ⎠

^{, (14)}

where we use the notation 

η

=

η

^s−

η

i, 

ω

^x=(

ω

^xs−

ω

^xi) ^and 

ω

^y=(

ω

^ys−

ω

^y_i)^.HereHot andH_obaretheangularmomentumof thetopandbottomparts oftheocean expressedintheBF ofthe adjacentsolidlayer:

Hot=

⎛

⎜ ⎝

Aot

ω

^xo+nc¹³_ot Bot

ω

^yo+nc²³_ot Cot

(

ⁿ+

γ

˙^o

)

+nc³³_ot

⎞

⎟ ⎠

⁽¹⁵⁾

H_ob=

⎛

⎜ ⎝

Aob

( ω

o^x−

 ω

^x+n

 η

^sin

(

^M+

ω ))

+nc_ob¹³ B_ob

( ω

^yo−

 ω

^y+n

 η

^cos

(

^M+

ω ))

+nc²³_ob

C_ob

(

ⁿ⁺

γ

^˙o

)

^+nc33ob

⎞

⎟ ⎠

^{. (16)}

The angularmomentum of the ocean dependson the polar mo- tion and obliquity difference betweenthe shell andthe interior.

The dependenceon the polar motioncomponents is essential to ourmodel while the obliquity terms can be neglected aswe do notstudyherethediurnalpolarmotion.Thediurnalpolarmotion hasaninfluence smallerthan 3%on thepolarmotionamplitude, aninfluencethatissmallerthantheuncertaintiesduetothelayer’

densities and radii. Therefore, neglecting the diurnal terms does notsignificantly influenceour resultsforthepolarmotionofthe shell.

3.2.Torques

3.2.1. Pressuretorque

Thepressureofthesubsurfaceoceanactingonitstopandbot- tomsurfacesexertsatorqueontheshellandontheinteriorgiven by



s^p=



So

(

^ro× ˆno

)

^P

(

^ro,

ϕ

,

λ )

^{dS (17)}



i^p=−



Si

(

^ri× ˆni

)

^P

(

^ri,

ϕ

,

λ )

^dS. (18)

Thepressuretorqueexertedontheoceanisgivenby



o^p=−

s^p− 

i^p. (19)

InEqs. (17) and(18), S_o and S_i are the outer and inner surfaces boundingthe ocean andnˆo and nˆ_i are the outward unit normal onthese surfaces.The top andbottom surfaces of theocean are oblatespheroidsdeformedbytides androtation. Thelocalradius

r_kofalayerkisgivenbyEq.(C.3).Thepressuregradientinsidethe subsurfaceoceancanbeobtainedfromtheNavier–Stokesequation expressedas

∇

^P=−

ρ

^o



∇

^W+

∇ 

int+

ω

s/i×

( ω

^s/i× r

)

+

ω

˙s/i× r+2

ω

s/i× 

v

r+d

v

r

dt



. (20)

Here thefirst andsecond termsare related to thepressure field inducedbyWandint,theexternalandinternalgravitationalpo- tential,respectively.ThethirdtermofEq.(20)representsthecen- trifugaltermduetotherotationofthereferenceframecorotating

with the adjacentsolid layer witha rotation speed

ω

s/i and the fourthtermarisesduetothetimevariationofthisrotationspeed.

ThefifthtermrepresentstheCorioliseffectduetherelativeveloc- ityofthefluid

v

r.Finally,thelasttermrepresentsthetotalderiva- tiveofthefluidvelocity.

If thesubsurface ocean is assumed to be in hydrostaticequi- librium, the right-hand side of Eq. (20) reduces to the first two components of thepressure gradient which are already included in Coyette etal. (2016). In that case, the pressure torque can be rewritten,usingGauss’theorem,as

s^p,H=−



Vot

r×

ρ

^o

( ∇ 

int

(

^r

)

+

∇

^W

(

^r

))

^{dV (21)}

i^p,H=



Vi

r×

ρ

^o

( ∇ 

int

(

^r

)

+

∇

^W

(

^r

))

^dV, (22)

whereV_otandV_iarethevolumeofthetop partoftheoceanand the of the interior, respectively. The part of the pressure torque that dependson the external potential W(^r) ^{will be referred to} astheexternalhydrostaticpressuretorque s^pe,Handthe partde- pending on the internal potential int(^r) ^{will be referred to as} the internal hydrostatic pressure torque  s^pg,H. From a compari- sonofexpressions(21)and(22)forthepressuretorquesandEqs.

(27)and(34)forthegravitationaltorque,thehydrostaticpressure torque due to the internal gravitational potential exerted by the differentlayers onthe shell (or interior)and to theexternal po- tentialexertedby Saturnon theshell(orinterior) isequaltothe gravitationaltorque exerted onthe adjacenttop (orbottom)part oftheocean (see e.g.Buffett, 1996).Thesetorqueswilltherefore be incorporated in the gravitational torque in Sections 3.2.2 and 3.2.3.

We herealso includethe flow inthe subsurface oceanin the formof a Poincaré flow,a uniform-vorticity flow witha residual flowrequiredinordertohaveatangentialflowatthetopandbot- tomoceanboundaries(Poincaré,1910,seealsoDehantandMath- ews,2015).Thisflowisthesimplestflowthatcanbeincorporated inananalytical approachandisthe sameastheflow inthe core of the Earth used in studies of the rotation of the Earth. In the presence of such a flow, Eq.(20) can be used to obtain the fol- lowing expression for the total torque acting on the ocean (see Mathewsetal.,1991andVanHoolstandDehant,2002forthetri- axialcase)

o=

ω

^{o× Ho. (23)}

Since the total gravitational torque (including the hydrostatic part of the pressure torque) acting on the ocean vanishes (see Van Hoolst et al., 2009 and Sections 3.2.2 and 3.2.3), the total torque is equal to the non-hydrostatic pressure torque acting on theocean o^p,NH.UsingEq.(14)fortheangularmomentumofthe ocean,wethenfind

o^p,NH=n

⎛

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎝

(

^{Cot− Bot}

) ω

^yo+

(

^Cob− Bob

)( ω

^yo−

 ω

^y

+n

 η

^cos

(

^M+

ω ))

− n

(

^c23ot +c²³_ob

)

−

(

^Cot− Aot

) ω

^xo−

(

^Cob− Aob

)( ω

^xo−

 ω

^x

+n

 η

^sin

(

^M+

ω )

+n

(

^c13ot +c_ob¹³

)

0

⎞

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎠

. (24)

This equation is similar to Eq. (34) of Van Hoolst (2007) (or Eq. (35) of Van Hoolst (2015) with a typo in the sign of the

^2c13^f term of the y-component) where the fluid layer was as- sumedto extendtothecenter. ThepartofEq.(24)that depends on the moment of inertia of the top of the ocean must be due to the shell while the part that depends on the moment of in- ertia of the bottom originates at interface with the interior. As

o^p,NH=− s^p,NH−  _i^p,NH,wethenhave

s^p,NH=−n

⎛

⎜ ⎝

(

^Cot− Bot

) ω

^yo− nc²³_ot

−

(

^Cot− Aot

) ω

^xo+nc¹³_ot 0

⎞

⎟ ⎠

⁽²⁵⁾

_i^p,NH=−n

⎛

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎝

(

^Cob− Bob

)( ω

o^y−

 ω

^y

+n

( η

^s−

η

i

)

^cos

(

^M+

ω ))

^{− nc23}ob

−

(

^Cob− Aob

)( ω

^xo−

 ω

^x

+n

( η

s−

η

i

)

^sin

(

^M+

ω ))

+nc¹³_ob 0

⎞

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎠

. (26)

The z-component of Eqs. (24)–(26) is zero, indicating that the Poincaré flow in the subsurface ocean only influences the polar motionandnotthevariationsinrotationrate.

3.2.2. Externalgravitationaltorques

ThegravitationaltorqueexertedbySaturnonthelayerkisde- finedby

^ek=−



Vk

r×

ρ (

^r

) ∇

^W

(

^r

)

dV, (27)

whereW(^r) ^{is the} gravitationalpotential ofSaturn atpositionr from the mass center of Titan andV_k is the volume of layer k. Thisexpressioneasilycompareswiththeexpressionfor _k^pe,H,the hydrostaticpressuretorqueduetotheexternalpotential(seeEqs.

(21)–(22)),andshowsthatthetorqueontheoceanvanisheswhen onlythehydrostaticpartofthepressureistakenintoaccount.The externaltorqueonlayerkcanalsobeexpressedastheoppositeof thetorqueoflayerkonSaturn

^ek=Mprp×

∇

^Wk

(

^rp

)

, (28)

whereM_p isthe massofthe centralplanetandW_k(^rp) ^thegrav- itationalpotential oflayerk atpositionrp oftheplanetfromthe mass centerofTitan. Inthe BFof layerk,andrestrictingW_k ex- pressed intermsofthecomponentsofthe inertiatensoroflayer k (MacCullagh’s theorem) to degree two, we then have(see also Balandetal.,2016)

^ek=3n²



_a

d



3

⎛

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎜ ⎝

(

^Ck− Bk+c³³_k − c²²_k

)

^r_k^yrk^z+c²³_k

(

^ry_k²− r_k^z2

)

+

(

^c13kr_k^y− c¹²_k r^z_k

)

^rk^x

(

^Ak− Ck+c¹¹_k − c_k³³

)

^r_k^xr_k^z− c¹³_k

(

^r_k^x2− r^z2_k

)

+

(

^c12_k ^rz_k− c²³_k r^x_k

)

^r_k^y

(

^Bk− Ak+c²²_k − c¹¹_k

)

^rxkr_k^y+c¹²_k

(

^rxk²− r^y_k²

)

+

(

^ck²³r^x_k− c¹³_k r^y_k

)

^rzk

⎞

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎟ ⎠

,

(29) whererˆ_k=(^rx_k,r^y_k,r^z_k)^{istheunitvectorpointingfromthecenterof} massofTitantoSaturn,expressedintheBFofthelayerk.disthe distancefromthe centerofmassofSaturn tothecenterofmass ofTitanandathesemi-majoraxisoftheorbitofTitan.

In an orbital framewhose x-axisischosen inthedirection of Saturn,thisunitvectorrˆ_kisexpressedbyconstructionas(1,0,0).

Byperformingtransformationsbetweenthisorbitalframeandthe BF ofthe layerkandbytakingintoaccounttheCassini state,we can expressthe unitvector rˆ_kasafunction oftherotationvaria- tions(librationsandpolarmotion),obliquityandorbiteccentricity

as(seeFig.3andCoyetteetal.,2016)

⎛

⎜ ⎝

r^x_k r^y_k r^z_k

⎞

⎟ ⎠

⁼

⎛

⎜ ⎝

1 2esinM−

γ

k

η

ksin

(

^M+

ω )

−^ω_n^xk

⎞

⎟ ⎠

^{, (30)}

where2esinMistheequation ofthecentercorrectuptothefirst orderinorbiteccentricitye.Here

ω

^{istheargumentofperiapsisof}

TitanandMitsmeananomaly.Eq.(30)showsthattheunitvector pointsapproximatelytowardthex-directionoftheBFasexpected in the Cassini state. The y-component of unit vector rˆ_k depends ontheorbiteccentricityandonthelibration angleasthisvector isnearly in theorbital androtational equators.The z-component dependsontheobliquityangleandpolarmotion.

Theexternalgravitationaltorque(Eq.(29))isthengiven,correct uptothefirstorderine,

η

k,

γ

k,

ω

^x_k^andtheincrementalinertiac_k^{i j} by

k^e=3n²

⎛

⎜ ⎜

⎝

0

(

^Ak− Ck

) 

η

ksin

(

^M+

ω )

−^ω_n^xk



− c¹³_k

(

^Bk− Ak

)(

^2esinM−

γ

k

)

+ c¹²_k

⎞

⎟ ⎟

⎠

^{. (31)}

Thex-componentoftheexternal torqueis zerocorrectupto the firstorderasthex-axisoftheBFapproximatelypointstowardSat- urn.We are then left withonlytorque componentsin they and zdirectionsthatbothtendtorealignthelong-axisofTitaninthe directionofSaturn.

Includingthe hydrostaticpressuretorque relatedto theexter- nalgravitationalpotentialofSaturn(seeSection3.2.1),wefindthe followingexpression

s^e+

s^pe,H=3n²

⎛

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎝

0

(

^As− Cs+Aot− Cot

) 

η

ssin

(

^M+

ω )

−^ω_n^xs



−c¹³s − c¹³ot

(

^Bs− As+Bot− Aot

)(

^2esinM−

γ

^s

)

+c¹²_s +c¹²_ot

⎞

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎠

,

(32) fortheexternalgravitationaltorqueactingontheshelland

i^e+

i^pe,H=3n²

⎛

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎜ ⎝

0

(

^Ai− Ci+A_ob− Cob

) 

η

isin

(

^M+

ω )

−^ω_n^xi



−c¹³_i − c¹³_ob

(

^Bi− Ai+B_ob− Aob

)(

^2esinM−

γ

ⁱ

)

+c¹²_i +c¹²_ob

⎞

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎟ ⎠

,

(33) for the external gravitational torque acting on the interior. Here

 s^pe,H and  _i^{pe,H represent the part of the} hydrostatic pressure torquethatisduetotheexternalpotentialW(seeEqs.(21)–(22)).

The external gravitational torque (including hydrostatic external pressuretorque)actingontheoceanis o^e+ o^pe,H=0.

3.2.3. Internalgravitationaltorques

The presenceof a subsurfaceocean implies that the principal axesofthesolidinteriorcanbemisalignedwiththoseoftheshell.

Eachlayerthereforeexertsagravitationaltorqueontheothermis- alignedlayers.Thegravitationaltorqueactingonalayerkisgiven by

k^g=−



Vk

r×

ρ

k

∇ (

^r

)

^{dV (34)}