LETTERS
https://doi.org/10.1038/s41550-019-0958-x

Cascading parallel fractures on Enceladus
Douglas J. Hemingway! !1,2*, Maxwell L. Rudolph! !3 and Michael Manga! !2
Active eruptions from the south polar region of Saturn’s~500-km-diameter moon Enceladus are concentrated alonga series of lineaments known as the ‘tiger stripes’1,2, thought
to be partially open fissures that connect to the liquid water
ocean beneath the ice shell3,4. To date, no study simultaneously explains why the tiger stripes should be located onlyat the south pole, why there are multiple approximately parallel and regularly spaced fractures, what accounts for their
spacing of about 35!km, and why similarly active fissures havenot been observed on other icy bodies. Here we propose thatsecular cooling, which leads to a thickening of the ice shelland building of global tensile stresses5,6, causes the first fracture to form at one of the poles, where the ice shell is thinnest owing to tidal heating7. The tensile stresses are thereby
relieved, preventing a similar failure at the opposite pole.
The steadily erupting water ice loads the flanks of the open
fissure, causing bending in the surrounding elastic plate andfurther tensile failure in bands parallel to the first fracture—aprocess that may be unique to Enceladus, where the gravity is
too weak for compressive stresses to prevent fracture propagation through the thin ice shell. The sequence of fissuresthen cascades outwards until the loading becomes too weakor the background shell thickness becomes too great to permit
through-going fractures.
To explain the polar location of the tiger stripes, we connected
the shell structure predictions from a tidal heating model7 with athermal evolution scenario and associated tensile-failure model5,6.
To explain the regular spacing of the tiger stripes, we then modelled
the elastic bending stresses resulting from loading at the edge of thefirst fracture and the consequent initiation of parallel fractures.
Given that the erupted ice grains appear to be sampling the internal liquid-water ocean8–10, the tiger stripes may be taken as open
fissures that fully penetrate the ice shell11. Accordingly, although
the south polar region bears evidence of a complex geologic history12–15, most studies presume an extensional origin for the tiger
stripes. Whereas Yin et al.14,15 suggest a strike-slip origin for the tiger
stripes, attributing their regular spacing to a stress shadowing effect,
we regard an extensional origin as more likely because it can create
open conduits directly to the underlying liquid ocean, allowing theeruptions to be initiated and sustained by controlled boiling11,16.
Over time, changes in the semi-major axis and eccentricity ofthe orbit of Enceladus will alter the efficacy of tidal heating—theprimary source of internal heat production for Enceladus17. Even if
Enceladus is close to a steady-state configuration, in which viscous
relaxation at the base of the ice shell is counterbalanced by ongoingfreezing and melting18,19, any period of gradual secular cooling will
result in some net freezing of the ocean. As long as there exists alayer within the ice shell that is free of active faults and capable ofaccumulating tensile stresses, the volume increase associated withthe phase change from water to ice will result in ocean pressurization5 (Methods). As freezing progresses, tangential stresses build
until the tensile failure limit is exceeded somewhere in the ice shell(Fig. 1 and Supplementary Fig. 1). Since tidal heating should cause
the shell to be thinnest at the poles7,20, tangential stresses are maximized at the poles such that the initial failure should occur at oneof the poles, with either being equally likely. In contrast to previous
work21–23, this mechanism thus accounts for the polar concentration of the eruptive activity without requiring true polar wander.
Although periods of secular heating and net melting of the ice shellhave also probably occurred, especially if Enceladus began cold,
this would have generated compressive stresses and thrust faulting.
Because such faults do not readily produce open fissures, however,
we do not favour compressive failure as the mechanism primarily
responsible for the formation of the tiger stripes.
Whereas the cold upper part of the ice shell behaves elastically,
the warmer ice toward the base of the shell behaves viscously on
long timescales (Methods). Provided that the ductile portion of theice shell is not too thick, however, the fracture can penetrate theentire ice shell, establishing an open pathway directly to the underlying ocean6 (Supplementary Fig. 2). Crucially, once the first fissure
forms, the ocean overpressure is relieved, removing the mechanism
for generating large, global tangential stresses. That is, it is no longer
possible for a similar fracture to develop at the opposite pole, or
anywhere else. In contrast with previous models involving diapirs
or impact events21–23, our model thus accounts for the concentrationof activity at a single pole by virtue of the fact that the critical pressure- and stress-relieving failure can occur only once.
We therefore suggest that Baghdad Sulcus, which cuts directly
through the geographic south pole, was the first fracture to form
and that the remaining fractures formed through a distinct, thoughrelated, process (see below). Baghdad’s orientation of ~30° from thetidal axis approximately maximizes normal tensile stresses arising
due to diurnal tidal deformation24. Although these tidal stresses
are weaker (~14–85 kPa), when combined with the isotropic background tensile stress field resulting from ocean pressurization, thetotal could be sufficient to cause failure. This may explain the orientation of the tiger stripes, provided that the ice shell has not experienced non-synchronous rotation since their formation.
Once the first open fissure has formed, liquid water rises to thelevel of neutral buoyancy (the lower 90% becomes filled with water)and begins to boil off at the top, where it is exposed to the vacuum. Whereas ocean overpressure was important for forming thefirst fracture, it is relieved once that fissure is open, and does nottherefore contribute to driving the eruptions themselves (Methods),
which are instead driven by controlled boiling16. Turbulent dissipation associated with tidally driven cyclic flushing and refilling of
the liquid-filled fissure can prevent it from freezing shut4. Provided
these (or similar) mechanisms can maintain the stability of such

Department of Terrestrial Magnetism, Carnegie Institution for Science, Washington, DC, USA. 2Department of Earth & Planetary Science,
University of California Berkeley, Berkeley, CA, USA. 3Department of Earth & Planetary Sciences, University of California Davis, Davis, CA, USA.
*e-mail: dhemingway@carnegiescience.edu
1

NATURE ASTRONOMY | www.nature.com/natureastronomy

LETTERS

NATURE ASTRONOMY

ocean-to-surface pathways, the eruptions may persist for extended
periods. Some 91% of the erupted solids fall ballistically back to thesurface of Enceladus25, accumulating preferentially on the flanksof the open fissures26, possibly explaining the origin of the ridges
observed along the tiger stripes27.
If the cold upper part of the ice shell behaves elastically, both thematerial accumulated on the flanks of the open fissure and the lossof buoyancy associated with any localized ice-shell thinning act likedownward forces applied near the edge of the broken elastic plate,causing bending stresses to develop in the surrounding ice shell28,29(Supplementary Fig. 3). Given a load acting at the edge of the broken plate, it can be shown that the maximum bending stresses occur
at a distance from the first fracture
xm = απ4ab3ρw g (1 − ν 2 )6π
α e −π∕ 4sin
2
4
Te

Bending
stresses

Eruptions(2)

where E and ν are the Young’s modulus and Poisson’s ratio for theicy shell, g is the acceleration due to gravity, ρw is the density of theocean, and Te is the effective elastic-layer thickness. If the typicaltiger-stripe spacing of ~35 km is governed by the position where
maximum bending stresses occur, then, assuming E = 9 GPa,
ν = 0.25 and ρw = 1,020 kg m− 3, and with g = 0.113 m s− 2, we obtainTe≈ 5.2 km (Supplementary Fig. 4). Given the expected temperature
structure of the ice shell, this corresponds to a total local shell thickness of roughly 8–9 km (Methods), consistent with the south-polarshell thickness inferred from gravity, topography and librations7.
Approximating the load as a line load, V0, acting at the edge of thebroken plate, it can be shown (Methods) that the resulting bending
causes maximum tensile stresses at xm to reachσmax = V0

xm

Incipient secondary
fracturec

α4 =Ridge loading(1)

where α is the characteristic length scale for elastic flexure, given by
ETe3

First fracture(3)

Bending stresses can thus initiate secondary fractures in bands
parallel to the first fracture once the load acting at the edge of theplate is sufficient to cause σmax to exceed the tensile failure limit
for ice, σcrit. Although localized ice-shell thinning could, in principle, contribute to the bending of the plate, the process may be
self-limiting because the melting ice is opposed by viscous closureand the freezing that results from reduced local turbulent dissipation, making this effect too small to produce the necessary loading4(Methods and Supplementary Fig. 5). A more effective source ofedge loading may be the accumulation of erupted material onto
the ridges flanking the open fissure. The eruption rate and the fraction of erupted material falling back onto the ridges determines
the time required before the bending stresses cause tensile failure
and the initiation of a secondary fracture (Supplementary Fig. 6).
If, for example, the eruption rate of the solids from a single fissureis 20 kg s−1 and if 91% of this material goes into the formation of
the flanking ridges, then, taking the fissure length to be 100 km,
the loading on each ridge reaches the critical value after ~1 Myr
assuming σcrit = 1 MPa (Methods).
As the secondary fracture begins to propagate, the broken portionof the elastic plate can no longer contribute to supporting the bending moment and the effective elastic thickness is reduced. Assumingthat the load is not changing on the timescale of fracture propagation, the stress profile in the elastic plate becomes increasingly steep,Fig. 1 | Proposed sequence for tiger-stripe formation. a, The global
tensile stress field (orange shading) arising from secular cooling and theconsequent ocean pressurization leads to tensile failure at one of the poles,
where the ice shell is thinnest. b, Following the first fracture’s formation,
the erupting solids then accumulate preferentially in ridges flanking theopen fissure, loading the edge of the broken ice shell and producing bending
stresses in the surrounding elastic plate (orange and blue shades indicate
tensile and compressive stresses, respectively, in the elastic layer; theductile layer is not shown). The bending stresses eventually become large
enough to initiate a set of secondary fractures parallel to the first and at a
distance, xm, determined by the ice shell’s elastic properties, according toequation (1). c, Once open, the secondary fractures then develop in a fashion
similar to the first, resulting in a cascading sequence of parallel fissures.

with the crack-tip stress being always tensile and increasing steadily
as the fracture propagates. This effect is countered by the compressive stresses from overburden pressure, which increase with depth(Supplementary Fig. 7). For Enceladus, however, gravity is so weakthat such compressive stresses are not important and the crack-tip
stresses are always increasing as the crack propagates, meaning thatthe crack should penetrate rapidly through the entire elastic layer.
For larger icy-ocean worlds, the gravity is generally too strong toallow the crack-tip stress to build in this way, potentially explaining why similarly active fissures are not observed on other bodies(Fig. 2 and Methods).
Although roughly the lower 40% of the ice shell is ductile onlong timescales (Methods), it may behave elastically on the timescale of fracture propagation, such that the secondary fracture canpenetrate all the way through to the ocean—though future work is
required to model the full viscoelastic nature of this problem. Oncea through-going fracture is established in this way, it evolves in a
fashion similar to the first fracture, allowing the sequence to cascade outward from the original fracture in symmetric pairs (Fig. 1):
after Baghdad, Cairo and Damascus, then Alexandria and the feature informally named ‘E’ by Yin and Pappalardo14. Although thedriving processes and geometries differ, loading-induced stresses
NATURE ASTRONOMY | www.nature.com/natureastronomy

LETTERS

NATURE ASTRONOMY
1,500

Enceladus

1,000

10

500

0

1
Europa–500

Bending-induced cracks
can penetrate ice shell
0.1

0

0.5

1.0

–1,000

1.5

Crack-tip stress gradient (Pa m–1)

Elastic-layer thickness (km)Crack propagation
hindered by gravity

Titan
Ganymede

Callisto

100

Pex =

–1,500Acceleration due to gravity (m s )Fig. 2 | Crack-tip stress gradient. Change in tip stress as the crack
propagates, as a function of surface gravity and effective elastic-layer
thickness, given by equation (31). Positive values (green) correspondto conditions that permit the crack-tip stress to become increasingly
tensile as the crack propagates. Negative values (magenta) correspondto conditions in which the compressive stresses due to overburdenpressure build too rapidly to permit crack propagation immediately
following the initial bending-stress-induced failure (Supplementary Fig. 7).
The dashed black contour represents the transition between the two
regimes. The vertical black lines represent the estimated ranges ofelastic-layer thicknesses for icy-ocean worlds.

have similarly been invoked to explain the spacing of hotspot volcanoes on Earth30,31.
There are several ways in which the cascade may be arrested. If
the eruption rate becomes too slow, bending stresses may not reachthe critical value on a timescale for which the ice behaves elastically.
If the eruption rate through each fissure decreases as new fissures
are formed, those forming later may not be loaded quickly enoughfor the resulting bending stresses to generate additional fractures.
Fractures forming farther from the pole, where the background iceshell thickness is greater, may also have a more difficult time propagating through the thicker ductile part of the ice shell.
Finally, although our model may be able to account for the
initial formation of the tiger stripes, other processes, such as
tidal and gravitational stresses, may continue to rework thefissures, forming strike-slip features14,15 or the funiscular plains13,
for example. We emphasize that, in spite of the youthful appearance of the south-polar terrain, and our conclusion that the total
south-polar shell thickness at the time of tiger-stripe formation isclose to the present-day value of ~9 km, the continuing eruptive
activity, along with ongoing tidal reworking, could keep the surface young even if the tiger-stripe fractures did not form recently.
Likewise, our results are not incompatible with stratigraphicallyolder features having formed through similar or distinct processes
that may have operated prior to the formation of the currently
active fissures as part of the south polar region’s evidently complex
geologic history12–15.

( )h 1− ρiw
3⎛
1 ⎛R⎞ ⎞⎜⎜
1 + ⎜⎜ ⎟⎟⎟ ⎟33
β (ri − rc )2 ⎜⎝ ξ ⎟⎠ ⎟⎟⎟ξ⎜
+ E ⎜1 + 2ν
⎟⎟
3
3ri2
⎜⎜
⎛ R ⎞⎟⎜ ⎟ −1 ⎟⎟⎟⎠⎜⎜⎝⎜⎝ ξ ⎟⎟⎠(5)

where ξ is the radius corresponding to the base of the elastic layer (ξ = R − Te) and
β = 4 × 10−10 Pa−1 is the compressibility of water. This excess ocean pressure drives
tangential stresses in the ice shell that, when evaluated at r = R, are given by
σt =

3
2

Pex

( ) −1Rξ

3(6)

Supplementary Fig. 1 shows that, starting with an unbroken ice shell that is capableof supporting global-scale tangential stresses, a few hundred metres of ocean
freezing is sufficient to cause tensile failure in the ice shell.
Manga and Wang5 assumed a uniform ice shell such that there was no preferred
location for the tensile failure. For the purposes of the above calculation, we havemade the same assumption. However, whereas a non-uniform ice-shell thickness
should not meaningfully affect the above result34, it does affect the preferred
location for the tensile failure, as stresses will concentrate in the thinnest parts ofthe shell. Even before the initiation of eruptive activity, the ice shell is expected
to have been thinnest at the poles due to the non-uniform distribution of tidal
heating7. Hence, tangential stresses resulting from secular cooling should have been
concentrated at both poles, with tensile failure being equally likely to occur at either
pole. Although they are an order of magnitude smaller, tidal stresses and stresses
associated with topographic anomalies will also affect the stress field and may
contribute to determining the precise location and orientation of the first failure.
Propagation of initial fracture. Once initiated at the surface (Fig. 1a), tensile
cracks will propagate downward on a timescale that is rapid compared with theMaxwell time. We develop a model based on linear elastic fracture mechanics to
calculate the maximum depth of penetration of fractures. As the subsurface ocean
is pressurized, tensile stresses are generated globally in the elastic layer (we assume
that stresses do not accumulate in the lower, warmer part of the ice), and if thesestresses exceed the tensile strength of ice, a fracture will initiate at the surface
and propagate downward. The model from the previous section always predicts
maximum tensile stresses at the surface, rather than at intermediate depths within
the ice shell. If the tip of the fracture extends below the stressed elastic layer,
the upper region of the crack is under tension but the tip of the crack is undercompression due to overburden pressure. We model the downward propagation of
the fracture under the assumption that the entire ice shell behaves as a linear elastic
solid on the timescale of crack propagation but that only the elastic layer supports
tensile stresses that encourage fracture penetration. We use a boundary-element
code based on the displacement discontinuity method35, extended to include acrack-tip element6 to more accurately resolve the stress field near the crack tip. In
each calculation, we initiate a short crack near the surface, extending to one-halfthe depth (dc) where lithostatic compression balances the applied tensile stress (σt)
Tedc

∫0 σt dz = ∫0 ρigz dz(7)where z is depth from the surface. We note that this formula is similar to the
maximum depth of fracture penetration from Qin et al.36, except that we accountfor the presence of a viscous ice layer. We solve for dc, obtaining

Methods

Temperature structure. The total ice-shell thickness (d) and the effective
elastic thickness (Te) are related by the temperature structure of the ice shell.
For a conductive ice shell, accounting for the temperature-dependent thermal
conductivity, which goes as T −1, the temperature structure is given by32

dc =

z ∕d

NATURE ASTRONOMY | www.nature.com/natureastronomy

Tensile-stress accumulation due to secular cooling. Beginning with equations(2)–(5) by Manga and Wang5, we obtain a set of analytical expressions that capture
the effects of ocean overpressure and ice-shell tensile-stress accumulation resulting
from secular cooling and the corresponding ice-shell thickening. Following their5
notation, we assume an ice shell of outer radius R and inner radius ri above a global
liquid-water ocean over a rocky core of radius rc. Taking the ice and water densities
to be ρi and ρw, respectively, it can be shown that the excess pressure in the ocean
resulting from freezing a thickness of water h is given by
ρ

–2

⎛T ⎞
T (z ) = Ts ⎜⎜ b ⎟⎟⎟⎜⎝ Ts ⎟⎠

where z is the depth below the surface and where Ts and Tb are the surface and basal
temperatures, respectively. On long timescales, the warmer parts of the ice shell
behave viscously while the coldest parts of the ice remain elastic33. If we assume thatthe ductile layer corresponds to the ice that is warmer than 160–180 K, then from (4),and assuming Ts = 75 K and Tb = 273 K, we obtain Te∕d ≈ 0.59–0.68.

(4)

2Teσt
ρi g

We calculate the mode-I stress intensity factor KI , describing the magnitudeof stresses near the crack tip. Next, we incrementally extend the crack length,computing a new solution for displacements and KI for each successive crack(8)

LETTERS

NATURE ASTRONOMY

length. The crack is arrested if KI < 0, which is equivalent to assuming zero fracture
toughness. Our numerical and analytical results are in good agreement with one
another (Supplementary Fig. 2) and are consistent with earlier work6.
Bending stresses. The load acting at the edge of the broken plate causes bending
stresses to develop in the elastic part of the ice shell28,29,37. We assume that the ice
shell is completely broken and that there are no remaining horizontal loads—weneglect membrane stresses since they cannot be transmitted across the open
fissure. Approximating the problem in a Cartesian geometry, and assuming
cylindrical bending, the bending moment at any given point is given by
M (x ) = D

d2w
dx 2

(9)

where w is the deflection, x is the horizontal distance increasing away from thefracture and D is the flexural rigidity, given by
D=

ETe3

(10)

12(1 − ν 2)

We assume that the flexural rigidity does not vary with x . The load is related to thedeflection byD

d4 w+ ρw gw = q (x )
dx 4

(11)

(we assume only vacuum exists above the ice shell).
For a line load acting at the edge of the broken plate, it can be shown28 that the
deflection is given byw=

V0α 3 −x ∕α
x
e
cos
2D
α

(12)

where V0 has units of N m−1 and where α is given by
α4 =

4D
ρw g

(13)

It can be shown that the bending moment per unit length along the fracture (unitsof N) is then given by
M (x ) = V0α e−x ∕α sin

x
α

(14)

The maximum bending moment occurs where dM ∕dx = 0, at
xm = απ
4(15)The fibre stresses within the plate can be written⎞
12M (x ) ⎛⎜ Te⎟
σxx (x , z ) =⎜ − z ⎟⎟⎟⎠
Te3 ⎜⎝ 2(16)

where z is the vertical position measured downward from the top of the plate(Supplementary Fig. 3). The maximum fibre stresses occur at x = x m, where thestresses are tensile in the upper half ( z < Te∕2) and compressive in the lower half
( z > Te∕2) of the deflected elastic plate.
This distance, x m, varies as a function of the elastic properties of the shelland is a 3/4-power function of Te and a 1/4-power function of E (Supplementary
Fig. 4). Assuming uniform properties across the ice shell, tensile failure will
occur at this distance and parallel to the first fracture. Given E = 9 GPa, we findthat the maximum stresses occur at a distance of 35 km from the fracture whenTe≈ 5.2 km (Supplementary Figs. 3 and 4). Assuming a predominantly conductive
ice shell, and taking the elastic layer to correspond to the uppermost part of the ice
shell where the temperature is <160 K, this elastic thickness implies a total shell
thickness (elastic plus ductile layers) of ~8.9 km, consistent with the polar iceshell thickness inferred from shape, gravity and libration observations7. Whereas~5.2 km may have been the relevant effective elastic thickness at the time of tigerstripe formation, subsequent heating or reworking and faulting of the ice shell
may have reduced the modern effective elastic thickness to perhaps <2 km, inline with the low end of estimates based on studies of local flexural and tectonic
features38,39—though we emphasize that these studies provide only a lower boundover the history of Enceladus.
The elastic properties similarly determine how the magnitude of the fibre
stresses is related to the magnitude of the load. Given that the maximum tensilestress occurs at x = x m and z = 0, from equation (16), its magnitude isσmax =

6M (x m)Te2

(17)which can also be written in terms of V0 as
σmax = V06π
α e−π∕4sin
4
Te2

(18)

Tensile failure occurs when the load is sufficient to make the bending moment at
x = x m equal to the critical value ofσ T2
Mcrit = crit e
6(19)

where σcrit is the tensile failure limit for cold ice, which we take to be 1 MPa for
intact ice40 or 100 kPa for previously weakened ice.
Ice-shell thinning. Here, we assess the possible effect of loss of buoyancy due
to thinning ice in the vicinity of an open fissure. We calculate the steady-state
temperature distribution around a fissure using a radial-basis-function finitedifference approach. We use an isothermal surface boundary condition of Ts = 75 K,
and an isothermal boundary condition of Tm = 273 K on the crack wall and along
the ice–ocean interface (Supplementary Fig. 5a). For the portion of the crackthat extends above the level of neutral buoyancy, the temperature varies linearly
between the surface temperature and Tm at the upper surface of the water. We
applied a far-field insulating boundary condition at a distance of 30 km fromthe crack. We determine the temperature-dependent thermal conductivity usingk = c ∕T , where c = 651 W m−1 is an experimentally derived constant41. For a verticalcrack in an ice shell with thickness 12 km, the steady-state conductive heat flow is~750 W m−1 per side.
Once open, the tendency for the fissure to narrow due to freezing is opposedby dissipation as water is cyclically flushed in and out of the fissure4. Since thetemperature gradient, and therefore the rate of conductive heat loss away from thefissure, is greatest near the surface and decreases towards the base of the fissure, thefreezing will be most rapid at the top of the fissure and slower toward its base. We
idealize the crack as a wedge shape and estimate the dissipation within the crack as
well as the heat conducted away from the crack.
We estimate the steady-state opening angle, θ, for the melt-back wedge
illustrated in Supplementary Fig. 5a. Assuming that the crack opens and closes
periodically42, we can write the width of the crack, b, as a function of depth, z, andtime, t , asb (z , t ) = b0 + 2z tanθ + Asin(ωt )(20)

where b0 is the mean width of the crack at the surface, A is the amplitude of theoscillations and ω is the angular frequency. Conservation of mass requires that thevelocity of the liquid in the crack, u, be related to the variations in crack width∂ (ub) ∂b
=
∂z
∂t

(21)

The energy dissipated by oscillatory flow in a fracture is found by relating theresolved shear stress on the wall, τ, to the mean flow velocity, u , calculated fromequation (21), using the Darcy–Weisbach equation1
τ=f ρ u 22(22)

where f is a friction factor and ρ is the fluid density. Because the eruptiondriven net outward flow of water is many orders of magnitude slower than thepeak velocities in the oscillatory flow, it does not contribute meaningfully to thedissipation and is therefore neglected here. While highly turbulent, the average
properties of the flow are in quasi-steady state, and the net acceleration of the wateris negligible. The wall shear stress is thus balanced by a dynamic pressure loss perunit depth as2τ = b

dP
dz

(23)

In turn, the rate of energy dissipation per unit volume (Qv) is related to the rate ofdecrease in dynamic pressure
Qv = u

dP
dz

(24)

We multiply Qv by the crack width to obtain the energy dissipated per unit area ofcrack wall
qd = fρ u 3(25)

Dissipation decreases very rapidly with increasing crack width (increasing opening
angle). In Supplementary Fig. 5b, we show predicted values of the dissipation
NATURE ASTRONOMY | www.nature.com/natureastronomy

LETTERS

NATURE ASTRONOMY
integrated along the wall of the crack for a crack with a surface width of 1 m and an
oscillation amplitude of 0.9 m, assuming a friction factor f = 0.01 (representativeof fully developed turbulent flow in a smooth channel). Owing to the large
uncertainties in the crack width, amplitude of oscillation, and friction factor, thedissipation is also very uncertain, but our result nevertheless demonstrates that,
if the crack begins to freeze shut, dissipation in the crack becomes much larger
than the conducted heat flow, causing melting. On the other hand, dissipation
alone cannot produce enough melt-back to achieve opening angles larger than afew hundredths of a degree. Thus, the dissipation mechanism, while capable ofpreventing the fissure from freezing out4, does not produce a melt-back wedge wide
enough to contribute significantly to the bending stresses discussed above.
Ridge accumulation. Erupted material can accumulate in ridges at the surface,loading the edge of the broken elastic plate from the top. Such a load, per unit
length along the fracture, is
mg
V0 =
L

mġ
t
L

(27)

V0 = Aρr g(28)

where ρr is the assumed density of the ridge. Assuming ρr = 900 kg m− 3, for example,
we obtain a load sufficient to generate tensile stresses of 1 MPa when the crosssectional area of the ridge is ~3 km 2. If the effective tensile failure limit at thetime of the tiger stripes’ formation is only 100 kPa, then the required ridge crosssection is just ~0.3 km 2. Estimating the height (≲150 m) and width (≲2 km perside) of the present-day ridges from Cassini imagery12,27 and accounting forthe deflection (~120 m) discussed above, the ridge cross-sections may only be
0.2–0.3 km 2 at present. This suggests that our proposed mechanism requires either
that the bending-induced failures occurred only after the effective tensile failure
limit was reduced to ~100 kPa (perhaps due to pervasive weakening following
the formation of the first fissure), or that the ridges were initially larger and haveexperienced some erosion or relaxation (perhaps due to localized heating) since
their formation.
Subsequent fractures. Once the bending stresses are great enough to cause
tensile failure at the surface, a new fracture is initiated (Fig. 1b). Whether or
not the new fracture can penetrate the ice shell depends on how the stressfield evolves during crack propagation. Although this is a dynamic process, we
can gain some insight by considering the hypothetical static scenario in which
the crack is arrested after propagating a distance y from the surface. At this point,
the load, which has not changed on the timescale of crack propagation, would
now have to be supported by the partially fractured plate. That is, focusing on thelocation of the fracture, where x = x m, the bending moment that must be supported
remains fixed at Mcrit, given by equation (19). The fibre stresses in the remaining
unbroken part of the lithosphere would have to support the same bending moment
but with a reduced effective elastic-layer thickness. The bending stress at thecrack tip is thus
6Mcrit(Te − y ) 2=

σcritTe2(Te − y ) 2(29)

where Te is the initial elastic-layer thickness such that the effective remaining elastic
thickness is always Te − y. The bending stresses at the crack tip are therefore always
tensile and increasing in magnitude as the crack tip propagates downward throughthe ice shell (dashed red lines in Supplementary Fig. 7).
NATURE ASTRONOMY | www.nature.com/natureastronomy

σcritTe2(Te − y ) 2− ρice gy(30)

The gradient as a function of depth is then given by
2σ T 2d
σtip (y ) = crit e 3 − ρice gdy(Te − y )

Teg <

(26)

The accumulated load gives rise to tensile stresses which are maximum at x = x m
and given by equation (18). Supplementary Fig. 6 shows this maximum tensilestress as a function of time and accumulation rate, with contours indicating two
examples of tensile-failure limits. Assuming a tensile-failure limit of 1 MPa andthe same elastic properties used in Supplementary Fig. 3, failure occurs whenV0 = ~ 3.1 × 108 N m−1. Assuming an accumulation rate of 10 kg s−1 per 100 km
length of ridge, the accumulated load would be sufficient to initiate tensile
failure after ~875 kyr. As a point of reference, the current total rate of erupted
solids has been estimated at [51 ± 18] kg s−1, of which 9% is estimated to escapefrom Enceladus while the remaining 91% falls ballistically back to the surface25.
The rate of material accumulation per ridge is of course somewhat smaller anddepends on how the erupted material is distributed among the tiger stripes,
whose rates of activity vary. The eruption rates also vary with time overvarious timescales43,44.
The accumulated load may be related to the cross-sectional area, A, of theridge byσtip,bending (y ) =σtip (y ) =(31)such that σtip (y ) increases monotonically with y as long as

where m is the mass of accumulated material and L is the total length of the ridges.
Given an accumulation rate of ṁ , the load as a function of time is given by
V0 (t ) =This effect is opposed by the background compressive stresses that exist due
to overburden pressure, which increases linearly with depth (dotted blue lines in
Supplementary Fig. 7), such that the net stress at the crack tip is

2σcritρice(32)

Provided that equation (32) is satisfied, the tensile stresses due to bending
always exceed the compressive stresses due to overburden pressure, such that the
net stresses at the crack tip are always tensile, and increasingly so as the crack
penetrates deeper through the lithosphere (Supplementary Fig. 7). Hence, once
initiated in this way, crack propagation cannot be arrested at any point within thelithosphere. Whereas the condition described by equation (32) is readily satisfied
for Enceladus due to its thin ice shell and small gravity, it may not be satisfied forlarger bodies or when a thicker elastic layer is present. For example, this condition
is not satisfied for Europa unless the elastic-layer thickness is at the low end of theestimated range29—less than about 1.5 km (Fig. 2). For Ganymede, Callisto andTitan, the gravity is so strong that the compressive stresses overwhelm the bending
stresses for plausible shell thicknesses29,45–48, precluding rapid crack propagation.
Hence, because this situation is permitted only for bodies with low surface gravity,or when the ice shell is extremely thin, Enceladus may be unique among ocean
worlds for its ability to develop through-going fractures due to bending stresses(positive or green region of Fig. 2).
What happens when the crack tip reaches the ductile part of the ice? The
lower part of the ice shell is ductile on long timescales, and therefore does notparticipate in the gradually accumulated bending stresses that support the edge
loading. If the new fracture is propagating sufficiently rapidly, however, the entire
ice shell behaves elastically (the Maxwell time is at least a few hours, even for the
low-viscosity ice near the base of the ice shell). That is, as the fracture propagates
downward, stresses may build within the otherwise ductile regions of the ice shell,
helping to support the load, at least transiently. Determining how this changes
the evolution of the stress profile during crack propagation will be the subject of
future work.

Data availability

All required data are available in the published literature as indicated.

Code availability

The computer code required to carry out the calculations discussed herein isavailable on request from the corresponding author.

Received: 8 June 2019; Accepted: 29 October 2019;
Published: xx xx xxxx

References

1. Porco, C. C. et al. Cassini observes the active south pole of Enceladus. Science
311, 1393–1401 (2006).
2. Porco, C., DiNino, D. & Nimmo, F. How the geysers, tidal stresses, andthermal emission across the south polar terrain of Enceladus are related.
Astron. J. 148, 45 (2014).
3. Ingersoll, A. P. & Nakajima, M. Controlled boiling on Enceladus. 2. Model ofthe liquid-filled cracks. Icarus 272, 319–326 (2016).
4. Kite, E. S. & Rubin, A. M. Sustained eruptions on Enceladus explained byturbulent dissipation in tiger stripes. Proc. Natl Acad. Sci. USA 113,
3972–3975 (2016).
5. Manga, M. & Wang, C. Y. Pressurized oceans and the eruption of liquid wateron Europa and Enceladus. Geophys. Res. Lett. 34, 1–5 (2007).
6. Rudolph, M. L. & Manga, M. Fracture penetration in planetary ice shells.
Icarus 199, 536–541 (2009).
7. Hemingway, D. J. & Mittal, T. Enceladus’s ice shell structure as a window oninternal heat production. Icarus 332, 111–131 (2019).
8. Postberg, F., Schmidt, J., Hillier, J., Kempf, S. & Srama, R. A salt-water
reservoir as the source of a compositionally stratified plume on Enceladus.
Nature 474, 620–622 (2011).
9. Iess, L. et al. The gravity field and interior structure of Enceladus. Science
344, 78–80 (2014).

LETTERS
10. Thomas, P. et al. Enceladus’s measured physical libration requires a global
subsurface ocean. Icarus 264, 37–47 (2016).
11. Spencer, R. et al. in Enceladus and the Icy Moons of Saturn (eds Schenk, P. M.,
Clark, R. N., Howett, C. J. A. & Verbiscer, A. J.) 163–174 (Univ. Arizona
Press, 2018).
12. Patthoff, D. A. & Kattenhorn, S. A. A fracture history on Enceladus provides
evidence for a global ocean. Geophys. Res. Lett. 38, 1–6 (2011).
13. Bland, M. T., McKinnon, W. B. & Schenk, P. M. Constraining the heat flux
between Enceladus’ tiger stripes: numerical modeling of funiscular plains
formation. Icarus 260, 232–245 (2015).
14. Yin, A. & Pappalardo, R. T. Gravitational spreading, bookshelf faulting, andtectonic evolution of the south polar terrain of Saturn’s moon Enceladus.
Icarus 260, 409–439 (2015).
15. Yin, A., Zuza, A. V. & Pappalardo, R. T. Mechanics of evenly spaced
strike-slip faults and its implications for the formation of tiger-stripe fractures
on Saturn’s moon Enceladus. Icarus 266, 204–216 (2016).
16. Nakajima, M. & Ingersoll, A. P. Controlled boiling on Enceladus. 1. Model ofthe vapor-driven jets. Icarus 272, 309–318 (2016).
17. Nimmo, F., Barr, A. C., Běhounková, M. & Mckinnon, W. B. in Enceladus andthe Icy Moons of Saturn (eds Schenk, P. M., Clark, R. N., Howett, C. J. A. &Verbiscer, A. J.) 79–94 (Univ. Arizona Press, 2018).
18. Kamata, S. & Nimmo, F. Interior thermal state of Enceladus inferred from theviscoelastic state of the ice shell. Icarus 284, 387–393 (2017).
19. Čadek, O., Běhounková, M., Tobie, G. & Choblet, G. Viscoelastic relaxationof Enceladus’s ice shell. Icarus 291, 31–35 (2017).
20. Choblet, G. et al. Powering prolonged hydrothermal activity inside Enceladus.
Nat. Astron. 1, 841–487 (2017).
21. Tajeddine, R. et al. True polar wander of Enceladus from topographic data.
Icarus 295, 46–60 (2017).
22. Nimmo, F. & Pappalardo, R. T. Diapir-induced reorientation of Saturn’s moon
Enceladus. Nature 441, 614–6 (2006).
23. Roberts, J. H. & Stickle, A. M. Break the world’s shell: an impact on
Enceladus: bringing the ocean to the surface. In Lunar and Planetary Science
XLVIII 1955 (LPI, 2017).
24. Nimmo, F., Spencer, J. R., Pappalardo, R. T. & Mullen, M. E. Shear heating as
the origin of the plumes and heat flux on Enceladus. Nature 447, 289–91 (2007).
25. Ingersoll, A. P. & Ewald, S. P. Total particulate mass in Enceladus plumes
and mass of Saturn’s E ring inferred from Cassini ISS images. Icarus 216,
492–506 (2011).
26. Degruyter, W. & Manga, M. Cryoclastic origin of particles on the surface ofEnceladus. Geophys. Res. Lett. 38, L16201 (2011).
27. Crow-Willard, E. N. & Pappalardo, R. T. Structural mapping of Enceladus
and implications for formation of tectonized regions. J. Geophys. Res. Planets
120, 928–950 (2015).
28. Turcotte, D. L. & Schubert, G. Geodynamics (Cambridge Univ. Press, 1982).
29. Billings, S. E. & Kattenhorn, S. A. The great thickness debate: ice shell
thickness models for Europa and comparisons with estimates based on
flexure at ridges. Icarus 177, 397–412 (2005).
30. tenBrink, U. Volcano spacing and plate rigidity. Geology 19, 397–400 (1991).
31. Hieronymus, C. F. & Bercovici, D. Discrete alternating hotspot islands formed
by interaction of magma transport and lithospheric flexure. Nature 397,
604–607 (1999).
32. Nimmo, F. Non-Newtonian topographic relaxation on Europa. Icarus 168,
205–208 (2004).
33. Nimmo, F., Pappalardo, R. T. & Giese, B. Effective elastic thickness and heat
flux estimates on Ganymede. Geophys. Res. Lett. 29, 1158 (2002).
34. Kalousová, K., Soucek, O. & Cadek, O. Deformation of an elastic shell withvariable thickness: a comparison of different methods. Geophys. J. Int. 190,
726–744 (2012).
35. Crouch, S. L. & Starfield, A. M. Boundary Element Methods in Solid
Mechanics (George Allen and Unwin, 1983).
36. Qin, R., Buck, W. R. & Germanovich, L. Comment on “Mechanics of tidally
driven fractures in Europa’s ice shell” by S. Lee, R. T. Pappalardo, and N. C.
Makris [2005. Icarus 177, 367–379]. Icarus 189, 595–597 (2007).
37. Dombard, A. J., Patterson, G. W., Lederer, A. P. & Prockter, L. M. Flanking
fractures and the formation of double ridges on Europa. Icarus 223,
74–81 (2013).

NATURE ASTRONOMY
38. Bland, M. T., Beyer, R. A. & Showman, A. P. Unstable extension of Enceladus’
lithosphere. Icarus 192, 92–105 (2007).
39. Giese, B. et al. Enceladus: an estimate of heat flux and lithosphericthickness from flexurally supported topography. Geophys. Res. Lett. 35,
1–5 (2008).
40. Hammond, N. P., Barr, A. C., Cooper, R. F., Caswell, T. E. & Hirth, G.
Experimental constraints on the fatigue of icy satellite lithospheres by tidal
forces. J. Geophys. Res. Planets 123, 1–15 (2018).
41. Petrenko, V. F. & Whitworth, R. W. Physics of Ice (Oxford Univ. Press, 1999).
42. Hurford, T. A., Helfenstein, P., Hoppa, G. V., Greenberg, R. & Bills, B. G.
Eruptions arising from tidally controlled periodic openings of rifts on
Enceladus. Nature 447, 292–294 (2007).
43. Nimmo, F., Porco, C. & Mitchell, C. Tidally modulated eruptions
on Enceladus: Cassini ISS observations and models. Astron. J. 148,
46 (2014).
44. Ingersoll, A. P. & Ewald, S. P. Decadal timescale variability of the Enceladus
plumes inferred from Cassini images. Icarus 282, 260–275 (2017).
45. Spohn, T. & Schubert, G. Oceans in the icy Galilean satellites of Jupiter?
Icarus 161, 456–467 (2003).
46. Béghin, C., Sotin, C. & Hamelin, M. Titan’s native ocean revealed beneathsome 45 km of ice by a Schumann-like resonance. C. R. Geosci. 342,
425–433 (2010).
47. Hemingway, D. J., Nimmo, F., Zebker, H. & Iess, L. A rigid and weathered ice
shell on Titan. Nature 500, 550–552 (2013).
48. Vance, S., Bouffard, M., Choukroun, M. & Sotin, C. Ganymede’s internal
structure including thermodynamics of magnesium sulfate oceans in contact
with ice. Planet. Space Sci. 96, 62–70 (2014).

Acknowledgements

This work was made possible by the NASA/ESA Cassini mission to Saturn and, inparticular, the work of the Imaging Science Subsystem team. D.J.H. was funded inpart by the Miller Institute for Basic Research in Science at the University of California
Berkeley and in part by the Carnegie Institution for Science in Washington DC. D.J.H.
and M.M. acknowledge support from the Center for Integrative Planetary Science (CIPS)
at the University of California Berkeley. M.M. was supported in part by NASA Solar
System Workings grant 80NSSC19K0557. M.L.R. was supported in part by NSF DMS1624776. We thank the CIDER working group, supported by NSF EAR-1135452, forearly discussions that contributed to parts of this work. We thank M. Bland, R. Citron,
J. Jordan, S. Kattenhorn, E. Kite and T. Mittal for helpful discussions.

Author contributions

M.M. developed the analytical equations for the secular cooling-induced tangential
stresses in the ice shell. M.L.R. developed the analytical and boundary element models
for crack penetration for the first fracture and carried out the related calculations.
M.L.R. computed the relationship between turbulent dissipation in the fissure and thecrack opening angle. D.J.H. proposed the mechanism of forming subsequent parallel
fractures due to bending stresses and carried out the related calculations. D.J.H. drafted
the manuscript with input from M.M. and M.L.R. All authors discussed the analysis andreviewed and commented on the manuscript.

Competing interests

The authors declare no competing interests.

Additional information

Supplementary information is available for this paper at https://doi.org/10.1038/
s41550-019-0958-x.
Correspondence and requests for materials should be addressed to D.J.H.
Peer review information Nature Astronomy thanks Gaël Choblet and An Yin for theircontribution to the peer review of this work.
Reprints and permissions information is available at www.nature.com/reprints.
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims inpublished maps and institutional affiliations.
© The Author(s), under exclusive licence to Springer Nature Limited 2019

NATURE ASTRONOMY | www.nature.com/natureastronomySupplementary Information for “Cascading parallel fractures on
Enceladus”
Douglas J. Hemingway*1,2 , Maxwell L. Rudolph3 , Michael Manga2
November 18, 20191

Department of Terrestrial Magnetism, Carnegie Institution for Science, Washington, DC, 20015, 2 Departmentof Earth & Planetary Science, University of California Berkeley, Berkeley, CA, 94720, 3 Department of Earth& Planetary Sciences, University of California Davis, Davis, CA, 95616. *Corresponding author: dhemingway@carnegiescience.edu.

Tangential stress and ocean overpressure (Pa)

10 8

10 7

10 6

10

5

10

4

10 3

10 2
10 0

10 1

10 2

10 3

10 4

Thickness of water that freezes (m)Supplementary Figure 1: Tensile stress accumulation due to ocean freezing — Accumulation of excess oceanpressure (dotted blue line) and the resulting tangential stresses in the ice shell (solid blue line) as a functionof the thickness of water that has frozen in the presence of an ice shell that is capable of supporting global
tensile stresses. Also shown is a nominal tensile failure limit for cold ice (dashed red line).

Elastic layer thickness (km)14
12Crack Reaches Ocean(1 MPa)10
8

Crack Reaches Ocean(3 MPa)6
4

Crack Arrested

2
0

2

4

6

8

10

12

14

16

18

20

Ice shell thickness (km)Supplementary Figure 2: Initial fracture penetration — Conditions determining whether or not a fracture
initiated at the surface can penetrate entirely through the ice shell, with the dark and light grey regions
corresponding to applied stresses of 1 MPa or 3 MPa, respectively, and delineated with the solid black lines
given by equation (8), assuming g = 0.113 m/s2 and ⇢i = 930 kg/m3 . The symbols represent the resultsof the numerical analysis, with blue symbols indicating that the crack reaches the ocean and grey symbols
indicating that the crack is arrested; circles and triangles correspond to applied stresses of 1 MPa or 3 MPa,
respectively.

-2

maximumtensile stress

xm

-1

1
2
3
4
5
6

edge loading

Plate Deflection (km)

0

7
0

50

100

150

200

250

Distance from fracture (km)Supplementary Figure 3: Elastic plate deflection and internal stresses — Assumed parameters for this example are E = 9 GPa, ⌫ = 0.25, ⇢w = 1020 kg/m3 , and g = 0.113 m/s2 , for which the elastic thickness required
to deliver xm = 35 km, given by equations (1) and (2), is Te ⇡ 5.22 km. Deflection is computed according
to equation (12) and is shown without any vertical exaggeration (though note that the illustration’s aspectratio is approximately 30:1). Fibre stresses are computed with equation (16) and are illustrated with shading, where warm colours correspond to tensile (positive) stresses and cool colours to compressive (negative)
stresses. Compressive stresses due to overburden pressure are not shown. In the example illustrated, thedeflection at the edge of the plate reaches a maximum value of ⇡ 120 m when the maximum fibre stressreaches 1 MPa.
35

9

60

45

30

60

55

40

50

Elastic Thickness (km)

8

50

25

7
6

35

40

40

20
25

5
4

45

30

35

30

30

15

20

25

20

3
2

15

10

20
15

10

1

10

10

5

5

1

2

Distance to Maximum Bending Stress, x m (km)

10

5

4

8

16

Young's Modulus for Ice, E (GPa)Supplementary Figure 4: Position of maximum tensile stress — Distance from the fracture to the positionof maximum tensile stress, computed via equation (15), as a function of the ice shell’s elastic properties:
namely the Young’s modulus (E, shown on a logarithmic scale) and the e↵ective elastic layer thickness (Te ).
A fixed Poisson’s ratio of ⌫ = 0.25 is assumed. As a point of reference, a nominal Young’s modulus of
E = 9 GPa is indicated with a vertical dashed white line, illustrating that the spacing of 35 km corresponds
to an elastic thickness of Te ⇡ 5.2 km, indicated by the horizontal dashed white line.

a

θ

b

Supplementary Figure 5: Melt-back wedge development — (a) Development of ice shell temperature structurein the vicinity of a liquid filled fissure. The red dots represent the evaluation points in the radial basis function
calculation. (b) Turbulent dissipation as a function of the fissure’s opening angle, ✓. The dashed black line
indicates the approximate dissipation required to keep the fissure from freezing shut. For opening angles
larger than ⇠ 10 3 degrees, the dissipation is not sufficient to keep the fissure open.
5

50
1

4

35

3.5

30

3
1

25

2.5

20

2

15

1.5

1

10

1

1

5

0.5

0.1

0.1

0
0

0.2

Maximum Fiber Stress (MPa)

40

0.1

Accumulation Rate (kg/s)

4.5

0.1

45

0.4

0.6

0.1

0.8

0

1

Time (Myrs)

Supplementary Figure 6: Stress accumulation due to ridge loading — Accumulated tensile stresses at xm =35 km as a function of time and the rate of material accumulation in the flanking ridge, whose length isassumed to be 100 km. The dashed white contours illustrate failure envelopes assuming tensile failure limitsof 100 kPa or 1 MPa.

ab

Enceladus

0

0.5

1

1

1.5

1.5

Crack Tip Depth (km)Crack Tip Depth (km)

0.5

2
2.5
3
3.5

2
2.5
3
3.5

4

4

4.5

4.5

5
-4

Ganymede

0

5
-3

-2

-1

0

1

Crack Tip Stress (MPa)

2

3

4

-4

-3

-2

-1

0

1

2

3

4

Crack Tip Stress (MPa)Supplementary Figure 7: Crack tip stress resulting from bending — Stresses at the crack tip as a functionof crack tip depth for the cases of (a) Enceladus and (b) Ganymede. Positive stresses are tensile, negative
stresses are compressive. The dotted blue line represents the compressive stresses due to lithostatic pressure.
The dashed red line represents the bending stresses given by equation (29). The solid black line represents
the net stress given by equation (30). Whereas for Enceladus, the crack tip stress is always increasing as thecrack propagates downward, the crack tip stress initially decreases for Ganymede, even when we assume an
equally thin ice shell. The surface gravity is g = 0.113 m/s2 for Enceladus and g = 1.428 m/s2 for Ganymede.
For both examples, we assume Te = 5.2 km, crit = 1 MPa, and ⇢ice = 930 kg/m3 .