SUPPLEMENTARY INFORMATION
Supplementary Information

doi:10.1038/nature12400

Contents
S1 Data
S1.1 Gravity Model . . . . . . . . . . . . . . . . . . . . . . . . . . . .
S1.2 Topography Model . . . . . . . . . . . . . . . . . . . . . . . . . .

2
2
2

S2 Ice Shell Model
S2.1 Origin of long-wavelength topography . . . . . . . . . . . . . . .
S2.2 Ice shell structure . . . . . . . . . . . . . . . . . . . . . . . . . . .
S2.3 Lateral flow in the ice shell . . . . . . . . . . . . . . . . . . . . .

6
6
6
7

S3 AdmittanceS3.1 Observed Admittance . . . . . . . . . . . . . . . . . . . . . . . .
S3.1.1 Monte Carlo Analysis . . . . . . . . . . . . . . . . . . . .
S3.2 Model Admittance . . . . . . . . . . . . . . . . . . . . . . . . . .

8
8
9
9

S4 Results
15
S4.1 Degree-3 maps . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
S4.2 Uncertainty in degree-3 erosion estimates . . . . . . . . . . . . . 15
S4.3 Eﬀect on tidal Love number, k2 . . . . . . . . . . . . . . . . . . . 17
S4.4 Degree-2 admittance and fluid Love number . . . . . . . . . . . . 18
S4.5 Degree-4 predictions . . . . . . . . . . . . . . . . . . . . . . . . . 21
References

23

W W W. N A T U R E . C O M / N A T U R E | 1

RESEARCH SUPPLEMENTARY INFORMATION

S1

Data

A field F (such as gravity or topography) may be expanded on a sphere as
follows 16
F (θ, φ) =

l
∞ !
!

(Clm cos mφ + Slm sin mφ) Plm (cos θ)(S1)

l=0 m=0

where Plm (cos θ) is a fully normalised associated Legendre function, θ iscolatitude, φ is longitude and Clm and Slm are spherical harmonic coeﬃcientsof degree l and order m.

S1.1

Gravity Model

The description of Titan’s gravitational field 3 is given as non-normalised, dig′g′
, Slm
. To convert these to fully normalised
mensionless potential coeﬃcients Clm
g
g
, Slm
, we write
gravity coeﬃcients, Clmg{C, S}lm =

#− 12
"
′
(l − m)!
GM(2 − δ0m )(2l + 1)
(l + 1) 2 {C, S} glm(l + m)!
R(S2)

where the square root term does the normalisation, the (l + 1) term arises
from the diﬀerentiation associated with converting from potential to gravity, theGM
R2 term generates dimensional coeﬃcients (which we will express in terms of
mGal = 10−5 ms−2 ), and δ0m is the Kronecker delta.
Three solutions have been derived 3 for Titan’s gravity field. In the first two(SOL1a and SOL1b), results from six Cassini gravity flybys were analysed separately and then combined into multi-arc solutions. Whereas SOL1a attempts
to model the gravity field only up to degree 3, SOL1b attempts to model the field
up to degree 4, but only as a means of verifying the robustness of the degree-3
solution. It is found that the degree-3 solution diﬀers only modestly between
SOL1a and SOL1b. For the last solution (SOL2), a global model extending to
l = 3 was derived using Pioneer and Voyager data and satellite ephemerides
in addition to Cassini observations. In spite of the diﬀerent approaches, theSOL2 field closely resembles the SOL1a field. The SOL1b solution matches lessclosely, likely because it also attempts to include the degree-4 component of the
field. Although all three solutions give a consistent estimate of the periodic k2
Love number (with a central value of 0.6), the static part of the degree-4 field
is currently not well constrained due to the limited nature of the observations.
Our calculations will be based primarily on the SOL1a gravity field.

S1.2

Topography Model
′

′

h
h
The non-normalised topography coeﬃcients Clm, Slm
are derived from a combination of radar altimetry and analysis of the overlapping regions of radar
images, in a technique known as SAR topo 9,13,14 , and were used to derive shape

2 | W W W. N A T U R E . C O M / N A T U R E

SUPPLEMENTARY INFORMATION RESEARCH

solutions up to l = m = 11. Several distinct solutions were produced, depending on where the harmonic expansion was truncated. We denote these solutions
Deg4-exp, Deg5-exp, ... Deg11-exp, where the number indicates the highest
degree and order used to fit the observations.
Due to the large gaps in Cassini radar coverage (Figure S1), topographymodels with power beyond degree 6 are not adequately constrained unless an apriori restriction is applied (minimise rms deviation from best-fit sphere). Even
when a priori constraints are applied, the coeﬃcients tend to be less stable
when the model’s expansion limit exceeds degree 6 (Figure S2). Likewise, theresulting admittance estimates tend to be most stable for the topography models
with power limited to degree 6 or less (Figure S3).
90

1000

500
30
0

0

−30

Elevation (m)Latitude (degrees north)

60

−500
−60
−90

0

90

180
Longitude (degrees east)

270

360

−1000

Figure S1: Cassini radar-derived elevation data for Titan. Elevation is given
relative to the 2575-km reference sphere.
For our purposes, we prefer to use the highest resolution data available
without requiring a priori constraints in the model fits. We therefore primarily
use the Deg6-exp model in our analysis. In the main text Fig. 1b and bothparts of Fig. 3, we use the Deg6-exp topography model 9 , the coeﬃcients forwhich are given in Table S1. To convert from non-normalised coeﬃcients to
h
h
, Slm
, we write
fully normalised coeﬃcients, Clm{C, S} hlm ="− 12!
(l − m)!
h′(2 − δ0m )(2l + 1)
{C, S}lm(l + m)!

(S3)

W W W. N A T U R E . C O M / N A T U R E | 3

RESEARCH SUPPLEMENTARY INFORMATION

40

Coefficient value (m)

20

0

−20

−40
C30
C31
S31
C32
S32
C33
S33

−60

−80

3

4

5

6

7

8

9

10

11

12

Model expansion limit

Figure S2: Degree-3 normalised topography model 9 coeﬃcients (with 1-σ errorbars) as a function of the maximum spherical harmonic degree allowed when
fitting the data.

50Admittance estimate (mGal/km)

−5

−10

−15

−20

−25

−30

−35

−40

−45

3

4

5

6

7

8

9

10

11

12

Model expansion limit

Figure S3: Degree-3 admittance estimates corresponding to the topographymodel coeﬃcients shown in Figure S2 and the SOL1a gravity field 3 . Admittanceestimates are obtained as described in section S3.1, using a Monte Carlo analysis(error bars illustrate the standard deviations of each Monte Carlo distribution).

4 | W W W. N A T U R E . C O M / N A T U R E

SUPPLEMENTARY INFORMATION RESEARCH

Table S1:
metres).

Fully normalised Deg6-exp topography model 9 coeﬃcients (inTerm

Estimate

1−sigma error

C00
C10
C11
S11
C20
C21
S21
C22
S22
C30
C31
S31
C32
S32
C33
S33
C40
C41
S41
C42
S42
C43
S43
C44
S44
C50
C51
S51
C52
S52
C53
S53
C54
S54
C55
S55
C60
C61
S61
C62
S62
C63
S63
C64
S64
C65
S65
C66
S66

2574750.0
0.0
23.1
18.5
−169.5
−17.8
31.8
120.8
20.1
−4.9
2.8
−21.3
−8.8
−2.9
0.0
14.3
−39.3
22.1
75.9
4.5
−26.8
16.7
−50.2
0.0
−47.3
28.6
−26.9
−14.0
37.1
−12.4
0.0
0.0
0.0
0.0
0.0
0.0
−2.5
8.9
−61.0
−24.1
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0

6.0
5.8
5.2
5.8
4.5
5.4
7.0
4.6
4.6
3.8
3.7
5.6
5.9
5.9
7.2
7.2
3.7
3.2
4.2
4.5
4.5
0.0
0.0
0.0
0.0
3.3
2.3
4.7
6.2
6.2
0.0
0.0
0.0
0.0
0.0
0.0
3.1
2.5
3.8
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0

W W W. N A T U R E . C O M / N A T U R E | 5

RESEARCH SUPPLEMENTARY INFORMATIONS2

Ice Shell Model

S2.1

Origin of long-wavelength topography

Titan’s shape is primarily determined by tides and rotation; the former resulting
in elongation along the tidal axis, and the latter in an equatorial bulge (flattening
at the poles). The degree-2 and degree-4 topography will also be aﬀected byvariations in ice shell thickness that arise due to tidal heating 7 . As discussed
below, non-Newtonian flow within the lower part of the ice shell could cause
degree-3 shell thickness variations to develop from a pattern that is initiallyconfined to degrees 2 and 4. This may, in part, explain the source of the observed
degree-3 topography. Since ongoing tidal heating will support the maintenanceof shell thickness variations at degrees 2 and 4, those variations could persist
even as lower crustal flow (discussed below) continues to generate shell thickness
variations at degree 3.
Heterogeneities in the ice shell could also contribute to a departure from thepurely degree-2 and -4 pattern predicted from tidal heating and could thus beresponsible for part of the shell thickness variations, and therefore topography,
at degree 3.

S2.2

Ice shell structure

We follow the common strategy of modelling the ice shell as an elastic layer
overlying an inviscid layer. Roughly speaking, ice will undergo a transition
from elastic to viscous behaviour at temperatures in the range 160 − 180 K,
depending on the exact strain rate and grain size assumed 30 . For a conductive
ice shell with basal temperature Tb = 270 K, the elastic thickness (T ) will thenbe 38-50% of the total shell thickness (d), while if Tb = 210 K, then T will be 5875% of d. While the transition from elastic to viscous behaviour will occur oversome finite region, that region will be thin because of the very strong variationin viscosity with depth. Hence, our two-layer model is a good approximation.
Assuming a heat flux of F ≈ 4 mW/m2 through Titan’s (conductive) ice
shell 7 , and allowing thermal conductivity to vary with temperature 31 , we canestimate the shell’s elastic thickness (T ) according to T = 567 ln (Tz /Ts ) /F ,where Tz is the temperature at which the shell transitions from elastic to viscous
behaviour. The resulting estimated elastic thickness is T ≈ 82 − 98 km. Despite
the highly approximate nature of this analysis, it yields an elastic thickness thatis consistent with our estimates (see main text Fig. 3b).
Throughout our analysis, we assume the ice shell to be in an equilibrium state
where the various forces (flexure within the elastic part of the shell, weight ofthe overlying topography, and buoyancy of the root) are in balance. This isreasonable because the vertical response time of the shell (analogous to postglacial rebound on Earth) should be fast compared with the loading timescale.

6 | W W W. N A T U R E . C O M / N A T U R E

SUPPLEMENTARY INFORMATION RESEARCH

S2.3

Lateral flow in the ice shell

Shell thickness variations lead to flow in the lowermost, low-viscosity part ofthe shell, which will tend to smooth out any such variations. For a Newtonian
fluid, the timescale (τ ) for removal of variations is given by 26τ=ηbg∆ρδ 3 k 2(S4)

where ηb is the viscosity at the base of the shell, g is the acceleration due togravity, ∆ρ is the density contrast between the shell and the fluid underneath, δis the eﬀective channel thickness in which flow occurs and k is the wavenumber(k = l/R, where l is spherical harmonic degree and R is the planetary radius).
Assuming a linear temperature gradient and a thermal conductivity 31 whichgoes as c/T , where c = 567 W/m, the eﬀective channel thickness (δ) is given
by 26
δ=

R g Tb d
Q ln (Tb /Ts )(S5)

where Rg is the gas constant, d is the shell thickness, Q is the activationenergy and Tb and Ts are the basal and surface temperatures, respectively.
Finally, the viscosity ηb is given by
"
#$
!
1
1
Q
ηb = ηref exp
(S6)
−
Rg Tb
Tref

where the viscosity of ice is ηref at a reference temperature Tref . Table S2
gives τ as a function of Tb for spherical harmonic degree 3. Here, we haveassumed ηref = 1014 Pa s at Tref = 273 K, Q = 60 kJ/mol, Ts = 90 K, d =100 km, g = 1.35 m/s2 , ∆ρ = 80 kg/m3 and R = 2575 km. For the range ofTb values explored, δ = 3.4 km. Rheological parameters 32 are subject to someuncertainty; nonetheless, the results of Table S2 serve to illustrate the main
conclusion, which is that flow is slow if the ocean is suﬃciently cold (Tb ! 220 K).
A temperature of 220 K corresponds to 25 wt% ammonia in a simple NH3 −H2 O
system 33 .
Table S2: Viscosity and timescale for removal of degree-3 ice shell thickness
variations as a function of temperature.
Tb (K) ηb (Pa s) τ (Myr)273
1.0 × 1014
0.55
250
1.1 × 1015
6.34
16
230
1.4 × 1077.6
210
2.8 × 1017
1491
In practice, the rheology of ice may be non-Newtonian, in which case ourflow timescales will be underestimates 26 , permitting larger values of Tb . An
important consequence of non-Newtonian flow is that mode-coupling occurs: an

W W W. N A T U R E . C O M / N A T U R E | 7

RESEARCH SUPPLEMENTARY INFORMATION

initially degree-2 or degree-4 pattern (e.g., due to tidal heating) will develop adegree-3 component as flow proceeds, thus potentially explaining the observed
degree-3 signal.
Finally, we note that if tidal heating is indeed occurring, a balance may
develop wherein shell thickness variations are being generated by tidal heating
just as quickly as lateral flow is removing those variations. Such an equilibrium
situation could be stable even if the relaxation timescales are relatively short,again permitting higher values of Tb .

S3

Admittance

The admittance, Zl , may be thought of as the ratio between gravity and topography at a particular wavelength 34 , and is typically measured in mGal/km.
An equivalent method makes use of the coherence between the Bouguer gravity and topography 35 . However, as discussed below, the likelihood of gravity
disturbances originating in the silicate core leads us to prefer the admittancetechnique in this case.

S3.1

Observed Admittance

We may define the cross-power spectrum Dij (l) between two fields i and j as
Dij (l) =

l
!

j
j
i
i
Clm
Clm
+ Slm
Slm

(S7)

m=0

With this definition, the estimated admittance, Z(l), and the correlation,
γ(l), between the gravity and surface topography, represented by subscripts g
and h, are as follows 15
Z(l) =

Dhg (l)Dhh (l)Dhg (l)γ(l) = "Dhh (l)Dgg (l)

(S8)
(S9)

If some fraction of the gravity signal is not correlated with the surface topography, then the coherence (γ 2 ) will be less than one. However, the crucial
advantage of equation (S8) is that any such gravity noise does not aﬀect theestimated admittance, Z. For the case of Titan, contributions to gravity fromdeeper interfaces (such as the silicate interior) are likely to be important, while
contributions to the surface topography from these processes are likely negligible. An approach like that embodied in equation (S8), which is unaﬀected bynoisy gravity, is essential for interpreting the limited observations available atTitan.

8 | W W W. N A T U R E . C O M / N A T U R E

SUPPLEMENTARY INFORMATION RESEARCH

S3.1.1Monte Carlo Analysis

The admittance estimates illustrated in the main text Fig. 2 were obtained
through a Monte Carlo analysis. For each of the three gravity and three topography models (nine combinations in total), the correlations and estimatedadmittances were obtained from a distribution of 100,000 distinct admittanceand correlation estimates, each of which was based on gravity and topography
coeﬃcients that were generated randomly according to the 1-σ errors in themodel coeﬃcients. Individual correlations and admittance estimates were computed according to equations (S8) and (S9).
The admittance estimated based on the Monte Carlo analysis will have aslightly smaller magnitude than the admittance estimated directly from thecoeﬃcients (i.e., when uncertainties are ignored). This is because, as long as
there is uncertainty in the topography coeﬃcients, the mean of the distributionof Dhh (see equation S7) will always be greater than the value of Dhh obtained
directly from the estimated topography coeﬃcients (because if x is normally2
distributed, then E(x2 ) > [E(x)] ). For example, if uncertainties are ignored,
the admittance computed directly from the SOL1a gravity 3 coeﬃcients andthe Deg6-exp topography 9 coeﬃcients is −39 mGal/km, whereas when uncertainties are accounted for using a Monte Carlo analysis, the mean estimatedadmittance is −32 mGal/km. We adopt the latter value because it is moreconservative—more negative admittances would require higher magnitudes of
erosion and/or larger elastic thicknesses.

S3.2

Model AdmittanceIn order to interpret the observed admittance, we require a model for how the
topography is supported. Here, we will assume that the topography is supported
by some combination of shell thickness variations and flexure 7 .
Gravity Anomaly
The gravity anomaly at degree l due to a thin surface layer of amplitude hl anddensity ρc is given by 36∆glt =(l + 1)4πGhl ρc
(2l + 1)(S10)

and similarly, the gravity anomaly due to a thin layer (a "root") of thickness
rl and density contrast ∆ρ = ρm − ρc at the base of the shell, is given by(l + 1)4πGrl ∆ρ∆glb = −(2l + 1)

"l+2
!
d
1−
R

(S11)

where the mean thickness of the shell is d, the radius of the body is R,
and ρm is the density of the material underlying the shell (i.e., the subsurface
ocean). In the short-wavelength limit (l ≫ 1), equation (S10) reduces to the

W W W. N A T U R E . C O M / N A T U R E | 9

RESEARCH SUPPLEMENTARY INFORMATION

usual flat-plate formula, as required. When the net gravity anomaly and thesurface topography (hl ) are known, the theoretical admittance is given byZl =∆glt + ∆glb
hl

(S12)

In the remainder of this development, we drop the subscripts from both hand r and take it as understood that these parameters correspond to a specificwavelength.
In practice, it will be diﬃcult to observe r and therefore to compute ∆glb
according to equation (S11). Instead, we would like to find an expression for r
in terms of h, which can be more readily observed. This is generally possible
because, for a finite elastic thickness, there will be a pressure balance between
the overburden of positive surface topography (ρc gh), the buoyancy of the root(∆ρgr), and the restoring forces due to flexure.
Lithospheric Deflection, Cartesian case
If the deflection of the initial ice shell is w, in a Cartesian system, and if geoid
undulations are neglected for the moment, this pressure balance can be writtenasD∇4 w = ∆ρgr − ρc gh(S13)Here, D represents flexural rigidity and is given byD=

ET 3
12(1 − ν 2 )(S14)

where T is the eﬀective elastic thickness, E is Young’s modulus and ν is
Poisson’s ratio. We treat w as positive upward bending, h as positive upward
relief above the reference ice shell surface, and r as positive downward relief
from the base of the ice shell. The relationship between r and h depends on theelastic properties of the shell and the thickness of loads applied at the top (dt )and bottom (db ) of the shell (Figure S4); dt is the thickness of material added at
the surface (a negative value would indicate erosion), and db is the thickness ofmaterial added at the base of the ice shell (a positive value would indicate basal
freezing, a negative value, basal melting). Our model represents the equilibrium
state achieved after the lithosphere has finished deflecting in response to theapplied load(s). The model also assumes that the ice shell properties do notchange over time.
Our formulation is similar to previous work 20,21 , however, our sign convention diﬀers slightly and, for simplicity, we assume that material added at the
top or bottom of the shell is also of density ρc . Our formulation also diﬀers
in that we handle top and bottom loads simultaneously with w being the total
deflection resulting from the combined eﬀects of top and bottom loading.

1 0 | W W W. N A T U R E . C O M / N A T U R E

SUPPLEMENTARY INFORMATION RESEARCH

h

w

dt

d

w

r

db

Figure S4: Illustration of the influence of top loading (pale pink material, dt )and bottom loading (pale blue material, db ) on ice shell flexure (w), surface
relief (h) and root thickness (r).
From Figure S4, we haveh = w + dt

(S15)

r = db − w(S16)Then, assuming the loads are periodic and in-phase, we can solve (S13) forw, obtainingw=∆ρdb − ρc dt
ρm + µ(S17)

Here, we have introduced a parameter, µ, which will serve as a shorthandfor the flexural rigidity at a particular wavelength and gravity. In a Cartesiansystem,µ(k) =ET 3 k 4
12(1 − ν 2 )g

(S18)

1
1 + µρ

(S19)

where k is a wavenumber. The advantage of using this shorthand will become
clear when we move from a Cartesian to a spherical system.
It is useful to define a compensation function, C(ρ), that expresses the degreeof compensation under flexural support compared with the case of pure isostasy.
This can be defined as the ratio of the deflection, w, according to (S17) tothe zero-rigidity deflection, w0 , obtained from (S17) when µ = 0. That is,
C = w/w0 , orC(ρ) =

When the elastic thickness, T , is zero, C = 1 (fully compensated). The
parameter ρ is the density contrast that is resisting the flexure (i.e., related to

W W W. N A T U R E . C O M / N A T U R E | 1 1

RESEARCH SUPPLEMENTARY INFORMATION

buoyancy, overburden pressure, or both). In the isostatic limit, ∆ρr = ρc h,
whereas in the top loading case, ∆ρr = ρc hC(∆ρ), while in the bottom loading
case, C(ρc )∆ρr = ρc h, as we will see. The theoretical value of C = 0 corresponds
to the zero compensation case which occurs when the ice shell is infinitely rigid(i.e., as µ → ∞). In this case, deflection (w) becomes zero (equation S17) andso, from (S15) and (S16), h = dt and r = db . In this scenario, h and r are
independent of one another and so both dt and db must be specified in orderto predict admittance. However, as long as C > 0, there will be some finite
deflection and it will be possible to obtain r as a function of h.
If C > 0 and both h and dt are specified, then from (S15), (S16), (S17) and(S19), it can be shown that!
"ρc h 1 − dht
dt
+
(S20)
r=
∆ρ C(ρc )h
We have factored out ρc h/∆ρ in order to facilitate direct comparison with
the isostatic case and because it will be convenient to do so when calculating
admittance using (S10), (S11) and (S12). In the case where no material has been
added to the surface (i.e., dt = 0, so that loading is purely from the bottom),this expression reduces to C(ρc )∆ρr = ρc h. Similarly, it can be shown that if
loading occurs purely from the top, ∆ρr = ρc hC(∆ρ).
Lithospheric Deflection, spherical case
The foregoing gives correct values for r in the Cartesian case, which is appropriate for short wavelengths. However, in order to interpret admittance at very
long wavelengths, we must consider the spherical case. Assuming the icy crust
behaves as a thin elastic shell of radius R, equation (S13) becomes$#D∇6 w+4D∇4 w+ET R2 ∇2 w+2ET R2 w = R4 ∇2 + 1 − ν (∆ρgr − ρc gh + ρm ghg )(S21)
where D is as in (S14). Here, we have adopted a modified version of previous
approaches 20,21 , which themselves follow an earlier derivation 19 . The final termin (S21) accounts for the elevation or depression of the geoid (hg , which we treat
as positive upward) that occurs with loading of the ice shell. Here, we will adopt
$
#
d l+2
≈ 1, andan approximation 20 to obtain hg , namely, we assume that 1 − Rthat the geoid and gravitational acceleration do not change with depth in theshell (as has been pointed out 37 , this was an implicit assumption of the previous
work 20 ).
Having obtained hg , and using (S15) and (S16), we rewrite (S21) as
⎡

'

⎣ 1−

3ρm
(2l + 1)ρ1 2 | W W W. N A T U R E . C O M / N A T U R E(−1⎛ T 2 ∇6 +4∇4
⎞
⎤
(
)
2
+
∇
+
2
ET ⎝ R2 12(1−ν 2 )⎠ + ρm ⎦ w = ∆ρdb − ρc dt
R2 g∇2 + 1 − νSUPPLEMENTARY INFORMATION RESEARCH

where ρ is the mean density of the body. If w is expressed in sphericalharmonics, we can replace ∇2 with −l(l + 1) and solve for w, recovering equation (S17),w=∆ρdb − ρc dt
ρm + µbut with the flexural rigidity parameter now being!
µ(l) = 1 −3ρm(2l + 1)ρ"−1

⎛ T 2 l3 (l+1)3 −4l2 (l+1)2
⎞
[
]
+
l(l
+
1)
−
2
2
2
ET ⎝
R 12(1−ν )⎠ (S22)R2 g
l (l + 1) − (1 − ν)For a spherical system, it is also necessary to account for the ratio of surface
areas at the top and bottom of the shell (since the buoyancy of the root depends
on its volume, not its thickness). This eﬀect complicates the derivations but it
can be shown that, if dt and h are specified, then the surface area correction
leads to
'
(!
"−2
dt
d
ρc h 1 − dht
+
1−
(S23)
r=
∆ρ C(ρc )hR
As required, this expression reduces to the Cartesian equivalent as R →∞. This correction also partially relaxes the simplifying assumption we made*)d l
≈ 1 rather thanto obtain the geoid height so that we now assume 1 − R
*
)
l+2
d
≈ 1.
1− R
Having obtained an expression for r as a function of h (which is possible aslong as C > 0), we can now substitute (S23) into (S11) and combine with (S10)
and (S12) to get an expression for admittance that depends on h and dt , but
not r
,!
'
+
"l (
1 − dhtdt
d(l + 1)4πGρc 1 −
+
1−
(S24)
Z(l) =(2l + 1)
C(ρc )
hR
Equation (S24) implicitly accounts for the root thickness (r), bottom load
thickness (db ), and mantle density (ρm ), such that these terms do not appear
in the final expression. Model admittance is sensitive to mantle density only
insofar as the ratio ρm /ρ influences the geoid, the eﬀect of which is captured inC(ρc ) via (S22) and (S19).
When C = 1 (pure isostasy), (S24) reduces to
'
"l (
!
d
(l + 1)4πGρc 1 − 1 −(S25)Z(l) =(2l + 1)Rland admittance will always be positive since (1 − d/R) must always be lessthan 1.

W W W. N A T U R E . C O M / N A T U R E | 1 3

RESEARCH SUPPLEMENTARY INFORMATION

100

80

Admittance (mGal/km)

60

40

20

0

−20

−40

−60

−80

−100
−1500

−1000

−500

0

500

1000

1500

Top load, dt (m)Figure S5: Admittance as a function of top load (dt , where negative dt indicateserosion) when h = 66 m, T = d = 200 km and assuming the properties listed inTable S3.
In general (i.e., when 0 < C < 1), admittance depends on topography(h), the amount of top and/or bottom loading (dt and/or db ), as well as theelastic thickness, T , and the mean shell thickness, d. For a given wavelength,mean shell thickness (d), elastic thickness (T ), and a fixed, positive h, admittance!is a positive linear function of dt (FigureS5) and crosses zero when"#−1 "$
#&
%
1
d −l
1− R
− 1 . When both top and bottom
dt = h 1 − C(ρc ) − 1

loading have taken place, and if h is known, admittance may be either positiveor negative, and is uniquely defined if either dt or db is specified.
It can be shown that, when loading occurs purely from the top (i.e., db = 0),
admittance is independent of surface relief, h, and is necessarily positive. Basedon (S24), when loading occurs purely from the bottom (i.e., dt = 0), admittanceis, again, independent of surface relief, h, and is positive as long as C(ρc ) >l(1 − d/R) . For degree 3, this is always true for the parameters given in Table S3.
Hence, degree-3 admittance may be negative only if dt /h is negative (i.e., whenerosion has occurred at topographic highs) or if substantially diﬀerent parametervalues are adopted. To obtain a negative degree-3 admittance with pure bottom
loading would require an increase of ∼ 30% in the ratio of Young’s modulus to
the ice shell density (E/ρc ). It is also possible to obtain a negative admittancewithout erosion for suﬃciently large elastic thicknesses and suﬃciently shortwavelengths (e.g., for l = 6 when T > 350 km, or for l = 9 when T > 200 km).

1 4 | W W W. N A T U R E . C O M / N A T U R E

SUPPLEMENTARY INFORMATION RESEARCH

Table S3: Parameter values assumed for admittance calculations.
Parameter
Symbol Assumed Value
Poisson’s ratio for iceν
0.25
Young’s modulus for ice
E
9 GPa
Crustal (ice shell) density
ρc
920 kg/m3
Mantle (subsurface ocean) density
ρm1000 kg/m3
Titan’s mean densityρ1880 kg/m3
Titan’s radius
R
2575 km
Acceleration due to gravity at the surfaceg1.35 m/s2

S4

Results

S4.1

Degree-3 maps

Our results suggest that the negative admittance we observe at degree-3 isthe result of negative gravity anomalies from large roots dominating over the
positive gravity anomalies from the associated topography. We tested this scenario by computing, everywhere over the surface of Titan, the gravity anomaly implied by the observed topography and then comparing the result withthe observations (Figure S6c). The gravity anomaly is obtained by multiplying the observed degree-3 topography (Figure S6a) by equation (S24), assuming
T = d = 200 km and a degree-3 erosion amplitude of 293 m (i.e, 293 m of erosion
at the topographic peaks and 293 m of deposition in the valleys). This is theamount of erosion required to produce −39 mGal/km, the admittance obtained
directly (i.e., neglecting uncertainties) from the SOL1a gravity 3 and the Deg6exp topography 9 (see section S3.1.1). Figure S6b shows the resulting gravity
anomaly, computed everywhere over the surface. The gravity field predicted
through this procedure resembles the observed field (compare panels (b) and(c) in Figure S6). Assuming the mantle density given in Table S3, the impliedroot thickness amplitude is ∼ 1.4 km.

S4.2

Uncertainty in degree-3 erosion estimates

As illustrated in Figure S5, admittance is approximately a direct linear functionof top load. Conversely, the top load required to produce a given admittancecan be obtained by solving equation (S24) for dt :%&#−1 $"#"#−ldZ(l) (2l + 1)
1
−1
1−
1−
dt = h 1 −−1
C(ρc )
4πGρc (l + 1)
R
(S26)
Based on (S26), Figure S7 shows how the implied erosion (negative dt ) varies with the estimated admittance given various combinations of shell thickness(d) and elastic thickness (T ). The dashed black line corresponds to Z(3) =!

"W W W. N A T U R E . C O M / N A T U R E | 1 5

RESEARCH SUPPLEMENTARY INFORMATION

(a)

(b)

(c)
Figure S6: Degree-3 topography and gravity maps: (a) Deg6-exp topography 9 ;(b) gravity computed as described in section S4.1 (assuming T = d = 200 km
and 293 m of erosion); (c) SOL1a gravity 3 .

1 6 | W W W. N A T U R E . C O M / N A T U R E

SUPPLEMENTARY INFORMATION RESEARCH−32 mGal/km, the admittance estimate obtained from the Monte Carlo analysis (i.e., accounting for uncertainties) based on the SOL1a gravity 3 field andthe Deg6-exp topography 9 solution. This is also the admittance assumed ingenerating Fig. 3b in the main text. The ±11 mGal/km uncertainty in that
admittance estimate (see main text Fig. 2) translates to ±81 m uncertainty inthe erosion estimate when T = d = 200 km. Diﬀerent combinations of T and dlead to slightly diﬀerent uncertainties, but roughly ±30% is typical.
Although, as we argued in section S2.2, the ice shell is not likely to be
entirely elastic, adopting T = d leads to more conservative estimates of themagnitude of erosion. For example, the magnitude of erosion required to giverise to Z(3) = −32 mGal/km is ∼ 241 m when we assume T = d = 200 km and∼ 577 m when we assume T = 100 km, d = 200 km.
500

0T=200km

Top load, dt (m)

−500

−1000

d=200km

0km
d=10
0km
km
T=10
0
0
2
d=
0km
T=10

0km

10

−1500

0km

d=

5
T=

−2000

d=

m

0k

20

m
0k

5

T=

−2500

−3000
−70

−60

−50

−40

−30

−20

−10

0

10

Admittance (mGal/km)

Figure S7: Magnitude of implied erosion as a function of admittance given
several combinations of mean shell thickness (d) and elastic thickness (T ). Negative top loads correspond to surface erosion. The dashed black line indicatesthe mean admittance estimate obtained from the Monte Carlo analysis basedon the SOL1a gravity 3 and the Deg6-exp topography 9 (−32 mGal/km).

S4.3

Eﬀect on tidal Love number, k2

The tidal Love number, k2 , of Titan has been measured 3 , with a 2-σ lower
bound on k2 of 0.413-0.439. It is therefore important to check that the kindof rigid elastic lid that we are proposing does not contradict the observations.
To do so, we constructed a highly simplified model for the interior of Titan(Table S4).
This is not meant to be realistic, but suﬃces to demonstrate our results.
The density structure approximately satisfies the nominal moment of inertia

W W W. N A T U R E . C O M / N A T U R E | 1 7

RESEARCH SUPPLEMENTARY INFORMATION

Table S4: Simple model of Titan’s interior, used to determine the eﬀect of arigid shell on the tidal Love number, k2 .
Layer
Outer
Rigidity
Viscosity
DensityRadius(GPa)(Pa s)(kg/m3 )(km)
Solid Core
2110
3
1021
2600
21
High-Pressure Ice
2275
3
10
1000
Ocean
2575 − d001000
Rigid Outer Shell
2575
3
1021
920
constraint, while the low rigidity in the inner layers is designed to reproduce theobserved k2 . We followed previous work 38 in calculating the model k2 values.
For rigid shell thicknesses of d =5, 100 and 200 km we obtained k2 values of
0.568, 0.519 and 0.413, respectively. Hence, adding a rigid shell of thickness 100
km or 200 km reduces k2 by 9% or 27%, respectively—not enough to contradict
the observed k2 values 3 .

S4.4

Degree-2 admittance and fluid Love number

Degree-2 admittance analysis of Titan’s ice shell is complicated by the fact thatthe body is tidally and rotationally distorted. Tidal/rotational distortion dominates the degree-2 gravity signal and also makes a large contribution to thedegree-2 topography. If we assume that basal freezing, uplift and erosion processes act similarly at degrees 2 and 3, then ice shell thickness variations should
be responsible for a portion of the observed degree-2 gravity and topography
signals.
Separating the degree-2 gravity and topography signals into their hydrostatic (i.e., tidal/rotational) and non-hydrostatic (i.e., due to ice shell thickness
variations) parts is necessarily an iterative process. We begin by estimating
Titan’s fluid Love number, h2f , from the observed degree-2 gravity signal (J2 )according to
6 g
J2
(S27)
5 Rω 2where R is Titan’s mean radius, g is the mean surface gravity, and ω is theangular frequency of rotation. Based on the observed gravity field 3 , we obtainh2f ≈ 2.0. This allows us to predict the expected hydrostatic topography, hT
according to
!
$
#"
# 5
R2 ω 2 1 "3 cos2 φ + 1 1 − cos2 θ −(S28)
hT = h2fg26
h2f = 1 +which we then subtract from the observed topography to get the non-hydrostatic
topography, hshell (i.e., that which is due to variations in ice shell thickness). We
then multiply hshell by equation (S24) to estimate the gravity signal due to an ice

1 8 | W W W. N A T U R E . C O M / N A T U R E

SUPPLEMENTARY INFORMATION RESEARCH

shell with anomalously deep roots (the large root size will be forced implicitly
by our choices of C and dt , both of which will be a function of mean shell
thickness, d and elastic thickness, T ). Having obtained this gravity anomaly
due to ice shell thickness variations (Δgshell ), we conclude that the portion of
the gravity signal that is due to tidal distortion is ∆gtidal = ∆gtotal − ∆gshell .
Finally, we use the newly obtained tidal gravity field to get an updated estimate
for h2f , again using (S27). After 3-4 iterations, our estimate for h2f converges
to the fourth decimal place, allowing us to separate, in a self-consistent way,
the tidal/rotational and ice shell thickness contributions to the degree-2 gravity
signal. The final estimate for h2f depends on the assumed mean shell thickness(d) and elastic thickness (T ), as illustrated in Figure S8. If we assume thatT = 100 km and d = 200 km, then we obtain h2f ≈ 2.15 (moment of inertia
factor ∼ 0.36).
2.2

2

3

2.1

2.

2.1

2.16
2.14

400

2.12
14

2.1

2.

300

5

1
2.

17

2.

16

8
.1

2

9

2.1

2.08

2.

200

2.06
654321987654321987654321987654321987654321987654321987
221211.1.1.1100.1109.8079.68057.4063.5024.1302.0918.796.857.46352413654321908721654321987654321987654321987654321987

Fluid Love number, h2f

1

09
2.

11

2.

500

Shell Thickness, d (km)

2.18

08
2.

2.0

7

600

2.04

100

2.02

50

100

150

200

250

300

350

400

2

Elastic Thickness, T (km)Figure S8: Estimate for Titan’s fluid Love number, h2f , as a function of meanshell thickness (d) and elastic thickness (T ).
Using this fluid Love number, we can estimate the non-hydrostatic portionsof the degree-2 topography and gravity signals. Figure S9 shows how the observed degree-2 gravity field (a) compares with the predicted field (b), where
the predicted field is the sum of the estimated hydrostatic gravity (c), basedon h2f = 2.15, and the gravity anomaly expected from the estimated ice shell
thickness variations (d), assuming an erosion amplitude of 577 m (obtained fromthe main text Fig. 3b assuming T = 100 km and d = 200 km). The amplitude ofthe estimated non-hydrostatic gravity is ∼ 2 mGal while the estimated hydrostatic gravity amplitude is ∼ 21 mGal. This is one measure of Titan’s departurefrom hydrostatic equilibrium.

W W W. N A T U R E . C O M / N A T U R E | 1 9

RESEARCH SUPPLEMENTARY INFORMATION

(a)

(b)

(c)

(d)
Figure S9: Degree-2 gravity maps centered on the sub-Saturnian point (180°
longitude): (a) SOL1a gravity 3 ; (b) total predicted gravity signal; (c) gravity
signal caused by tidal/rotational distortion assuming h2f = 2.15; (d) gravity
signal caused by ice shell thickness variations assuming T = 100 km and d =200 km, and therefore 577 m of erosion.

2 0 | W W W. N A T U R E . C O M / N A T U R E

SUPPLEMENTARY INFORMATION RESEARCH

S4.5

Degree-4 predictions

Assuming once again that the observed topography is the result of uplift due
to basal freezing and that the surface has experienced a similar magnitude oferosion at degrees 3 and 4, we can predict the admittance at degree 4. We
first obtain the magnitude of degree-3 erosion (dt ) over a range of values for T
and d from the main text Fig. 3b. Using this same value of dt for the degree-4
erosion amplitude, we then use equation (S24) to estimate the admittance overthe same range of values for T and d. Figure S10 illustrates that the degree-4
admittance should be negative if the elastic thickness, T , accounts for most of
the total shell thickness, d.
If, for example, T = d = 200 km, and we assume the same magnitude oferosion at degrees 3 and 4 (in this case, 241 m), then based on the observed
topography 9 (Figure S11a), we obtain a degree-4 admittance of −5.1 mGal/km
and we can compute the implied degree-4 gravity anomaly everywhere overthe surface (Figure S11b). Although the amplitudes are similar, our result is
spatially unlike the reported degree-4 gravity field 3 (Figure S11c). However, as
noted in section S1.1 as well as in the main text, the degree-4 gravity field is not
currently regarded as reliable; future gravity flybys are expected to improve thedetermination of the degree-4 field by a factor of two, providing a better test of
our prediction. Note that if we assume instead that T = 100 km and d = 200 km(the corresponding erosion amplitude being 577 m), the greater compensation
leads to a muted gravity signal and an admittance that approaches zero. Hencea weak degree-4 gravity signal, or a degree-4 admittance that is positive or
only weakly negative, may be an indication that the elastic layer accounts for asmaller portion of the total shell thickness (Figure S10).
20

600
16
14

15
12

10
2

8

400

10

0

6
4

−6

−2

5

−4

300

−−420

0

200

Degree−4 Admittance (mGal/km)

Shell Thickness, d (km)

500

−5

100
−6

50

100

150

200

250

300

350

400

−10

Elastic Thickness, T (km)Figure S10: Degree-4 admittance predicted for a range of elastic thicknesses (T )and total shell thicknesses (d).

W W W. N A T U R E . C O M / N A T U R E | 2 1

RESEARCH SUPPLEMENTARY INFORMATION

(a)

(b)

(c)
Figure S11: Degree-4 topography and gravity maps centered on the subSaturnian point (180° longitude): (a) Topography 9 ; (b) Predicted gravity signal
caused by ice shell thickness variations assuming T = d = 200 km and therefore 241 m of erosion; (c) SOL1b gravity 3 (not currently considered reliable atdegree 4).

2 2 | W W W. N A T U R E . C O M / N A T U R E

SUPPLEMENTARY INFORMATION RESEARCH

References
30.

Nimmo, F., Pappalardo, R. T. & Giese, B. Effective elastic thickness andheat flux estimates on Ganymede. Geophysical Research Letters 29, 1158
(2002).

31.

Klinger, J. Influence of a Phase Transition of Ice on the Heat and Mass
Balance of Comets. Science 209, 271–272 (1980).

32.

Goldsby, D. L. & Kohlstedt, D. L. Superplastic deformation of ice:
Experimental observations. Journal of Geophysical Research 106, 11017–
11030 (2001).

33.

Kargel, J. S. Physical Chemistry of Ices in the Outer Solar System. Solar
System Ices 3–32 (Kluwer: 1998).

34.

Lewis, B. T. R. & Dorman, L. M. Experimental Isostasy 2. An Isostatic
Model for the USA Derived from Gravity and Topography Data. Journalof Geophysical Research 75, 3367–3386 (1970).

35.

Forsyth, D. W. Subsurface Loading and Estimates of the Flexural
Rigidity of Continental Lithosphere. Journal of Geophysical Research 90,
12623–12632 (1985).

36.

Jeffreys, H. The Earth: its origin, history, and physical constitution. 535–
566 (Cambridge University Press: New York, 1976).

37.

Belleguic, V., Lognonné, P. & Wieczorek, M. Constraints on the Martianlithosphere from gravity and topography data. Journal of GeophysicalResearch 110, E11005 (2005).

38.

Moore, W. B. & Schubert, G. The tidal response of Ganymede andCallisto with and without liquid water oceans. Icarus 166, 223–226(2003).

W W W. N A T U R E . C O M / N A T U R E | 2 3