The persistence of oceans on Earth-like planets: insights from the deep-water cycle

0.33 MB·English

by  Laura Schaefer

Draftversion January6,2015
PreprinttypesetusingLATEXstyleemulateapjv.5/2/11
THE PERSISTENCE OF OCEANS ON EARTH-LIKE PLANETS: INSIGHTS FROM THE DEEP-WATER
CYCLE
Laura Schaefer and Dimitar Sasselov
Harvard-SmithsonianCenterforAstrophysics,60GardenSt.,Cambridge,MA02138
Draft version January 6, 2015
ABSTRACT
In this paper we present a series of models for the deep water cycle on super-Earths experiencing
5
plate tectonics. The deep water cycle can be modeled through parameterized convection models
1
0 coupled with a volatile recycling model. The convectionof the silicate mantle is linked to the volatile
2 cyclethroughthewater-dependentviscosity. Importantdifferencesinsurfacewatercontentarefound
fordifferentparameterizationsofconvection. Surfaceoceansaresmallerandmorepersistentforsingle
n
layer convection, rather than convection by boundary layer instability. Smaller planets have initially
a
larger oceans but also return that water to the mantle more rapidly than larger planets. Super-
J
Earthsmaythereforebe less habitableintheir earlyyearsthansmallerplanets,buttheir habitability
5
(assuming stable surface conditions), will persist much longer.

]
P
1. INTRODUCTION ues in favor of values ranging from 0.5 to 2.5 ocean
E
masses(OMs)ofwater(Hirschmann & Kohlstedt 2012).
. In the post-Kepler era, the number of known planets
h The presence of water in these minerals strongly influ-
of Earth and super-Earth size continues to grow. The
p ences their material properties, such as melting temper-
TESS mission, which is a planned survey of the bright-
- atures, rheology, phase changes, electrical conductivity,
o est and closest stars, is expected to detect thousands of
etc. (see e.g., Hirth & Kohlstedt 1996; Karato 1990;
r small planets (Ricker et al. 2015). The TESS planets,
t by virtue of their closeness and the brightness of their Litasov & Ohtani 2007). Studies of the effect of water
s onsilicaterheologyshowthatthepresenceofevenminor
a stars,willbemucheasiertofollowupwithground-based
amounts of water can reduce the viscosity by severalor-
[ photometricandspectroscopictechniquesthantheplan-
ders of magnitude (Karato & Wu 1993). This can have
ets of similar size detected by Kepler. Currently, the
1 a significant effect on the material flow of the mantle.
mini-Neptune GJ 1214 is the only planet in the super-
v In fact, it has been shown that the abundance of wa-
5 Earth mass range (1 − 10M⊕) for which atmospheric ter on the Earth’s surface is controlled by the deep wa-
measurementshavebeenmade(e.g.,Wilson et al. 2014;
3 ter/silicate cycle, which is tied to plate tectonics. The
Fraine et al. 2013;Berta et al. 2012). Withfuturetele-
7
abundance of water at the surface is a balance between
0 scoperesourcessuchasGMTandJWST,wewillbeable
therateofreturnviavolcanicoutgassingfromthemantle
0 to characterize even more super-Earths, some of which
at mid-ocean ridges (MORs) and the rate of loss to the
. may be in the habitable zones (HZ) of their stars. The
1 mantlethroughso-calledingassing,whichisthereturnof
classicalhabitablezoneisdefinedastheorbitalregionin
0 waterintothedeepmantlebythesinking,orsubducting,
which an Earth-likeplanet can maintain liquid water on
5 of water-rich oceanic seafloor. Much of this water is re-
its surface given a variety of atmospheric compositions.
1 leased immediately back to the surface through shallow,
However, the Earth’s atmosphere, although it permits
: water-induced volcanism. However, a small but signif-
v surface waterto remainliquid, does not controlthe sup-
icant fraction of the water contained in the subducting
i plyofwateronthesurface. Forthat,wemustlooktothe
X oceanic slabs can be transported to deeper levels of the
mantle. Inextendingthistosuper-Earths,itispresently
r unclearwhattheeffectofmassmaybeincontrollingthe mantle. For a detailed review of the exchange mecha-
a nisms between the surface and mantle, see Hirschmann
supply of water at the surface.
(2006).
Water on the Earth is found in two primary reser-
The deep water cycle on Earth has been stud-
voirs: (1) the surface oceans and (2) the silicate man-
tle. Water, either in molecular form or as an H+ or ied through the use of parameterized convection mod-
OH− group, can dissolve in silicate minerals. This is els incorporating a water-dependent viscosity (e.g.
McGovern & Schubert 1989; Bounama et al. 2001;
true even of the so-called nominally anhydrous minerals
Sandu et al. 2011; Crowley et al. 2011). Parameter-
(NAMs), which do not have hydrogen in their chemical
ized convection models are simplified 1D models using
formula. Furthermore, the solubility of water in silicate
a parameterization derived from more complicated and
mineralsincreaseswith pressure(Kohlstedt et al. 1996;
expensive 2D and 3D models of convection for differ-
Hirschmann et al. 2005).
ent systems (e.g. spherical vs. plane-parallel, heated
Geochemical studies of the materials that make up
from within vs. below, etc.). The parameterized mod-
Earth’smantlehavelongsuggestedthatitmayholdany-
els rely primarily on two dimensionless parameters: the
wherefrom1to10timesasmuchwaterasisfoundinthe
Rayleigh number, which describes whether a system is
surface oceans (e.g., Table 1 of Bounama et al. 2001),
unstable to convection, and the Nusselt number, which
although more recent work seems to discount larger val-
comparesthe convectiveandconductive heat flows. The
[email protected] Rayleighnumberdependsontheviscosityofthesystem,
2 Schaefer & Sasselov
which itself is dependent on temperature, pressure and
water fugacity. The abundance of water in the mantle, R
p δ
which helps determine the viscosity, evolves along with u
the mantle temperature.
Here we will use a parameterized convection model to solidus
study the deep water cycles of super-Earths. In the
current era of exoplanet studies, we are still search-
ing for Earth-like planets in Earth-like orbits around
Sun-like stars, but what we have found and what we s
u
can soon characterize, are super-Earths (∼ 1−2R⊕,∼ adi
R
1−10M⊕). Althoughsuper-Earthsmayhavecoolenough
surfaces to retain liquid water, tectonics and mantle
convection plays an important role in controlling its
abundance at the surface. The question of whether
these planets will have plate tectonics has been dis-
cussed previously in the literature (e.g. Valencia et al.
2007; O’Neill & Lenardic 2007; Korenaga 2010b; δ
c
Noack & Breuer 2013; Stamenkovi´c & Breuer 2014), R T TT T T
although the issue has not been settled. Here, we as- c s p u l c
Temperature
sumeplatetectonicsasastartingpointinordertoapply
Figure 1. Schematic of temperature-depth profile. See text for
asimilarwatercyclemodel. Thequestionthisworkthen
details.
addressesiswhetherthedifferentpressureregimesinside
super-Earths affect the deep water cycle, and whether
cally the characteristic viscosity. Discussion in the lit-
surfaceoceansaretectonicallysustainableontheseplan-
erature is often contradictory on where to define the
ets. This has important implications for the potential
characteristic viscosity: as an average value (see e.g.
habitability of super-Earths.
McGovern & Schubert 1989; Tajika & Matsui 1992;
In this paper, we address these questions using a pa-
Sandu et al. 2011), or at the base of a boundary layer
rameterizedthermalevolutionmodelcoupledwithawa-
(see e.g. Deschamps & Sotin 2000; Dumoulin et al.
ter cycle model. We give a full description of the model
2005; Stamenkovi´c et al. 2012). Here we will try both
in Section 2. In Section 3, we describe results for mod-
types of models, starting with the single-layer models
els which have either single layer convection or bound-
typically used with volatile evolution models, which use
ary layer convection for a number of super-Earth mod-
averageviscosities. Wewillthenlookataboundarylayer
els. We also describe the response of the models to rea-
model. We explore both models here to see the effect on
sonable variations in the parameter values. Section 4
thebehaviorofwater,andnotethatwearenotattempt-
discusses implications for the persistence of oceans on
ing to reproduce the Earth, but explore two valid, but
super-Earths.
different, models.
2. METHODS Thewatercyclemodeliscoupledtothethermalmodel
through the viscosity, which is dependent on the water
Mantle convection in the terrestrial planets can be
abundanceinthemantle(seee.g.McGovern & Schubert
controlled either by the stability of the whole man-
1989; Sandu et al. 2011). Models of the Earth have
tle layer (single layer convection, e.g. Schubert et al.
shown that the viscosity and mantle temperature create
(2001)) or by the stability of two boundary layers: a
a feedback loop (Schubert et al. 2001; Crowley et al.
cold boundary layer at the surface and a hot bound-
2011). When the mantle is warm, convection is vigor-
ary layer at the interface of the silicate mantle with the
ous and the mantle cools quickly. As the temperature
metallic core(see e.g.Turcotte & Schubert 2002). Heat
drops, the viscosity increases, causing the convection to
from either secular cooling or decay of radioactive ma-
become sluggish. Sluggish convection means that less
terials (or both) is transferred by conduction through
heat is removed from the mantle, causing it to heat up.
the boundary layers. The interior region is convective
Theviscosityincreases,andthecyclerepeats. Thiscycle
and thus approximately adiabatic. See Figure 1 for a
has been shown to be enhanced by the presence of wa-
schematic representation of the thermal profile. Single
ter (McGovern & Schubert 1989; Crowley et al. 2011).
layer convection models typically neglect the heat flux
Crowley et al. (2011) describe how the water and tem-
from the core, which is small, for simplicity and are
perature feedbacks interact for temperature and water-
entirely heated from within by radioactive decay. Nu-
dependent viscosities. Cowan & Abbot (2014) studied
merical simulations are used to determine scaling laws
the effect of sea floor pressure on a steady state model
relating the heat flux and mean temperature to the de-
ofthe deepwatercycle, withoutconsideringthe effectof
gree of convection. These simulations have been done
planetary thermal evolution.
for a variety of different boundary conditions (free slip,
In the following sections, we first describe selection of
no slip, etc.), and for different material properties (iso-
the super-Earth model parameters. We then describe
viscous, highly temperature-dependent viscosity, etc.)
the parameterizedthermal model, followedby the water
(e.g., Honda & Iwase 1996; Deschamps & Sotin 2000;
cycle model. Section 3 will discuss the results for these
Dumoulin et al. 2005). Choosing the proper scaling re-
models.
lations is therefore important for a successful model.
The parameter which is of the most impor-
2.1. Super-Earth Models
tance in determining the convective behavior is typi-
Oceans on Super-Earths 3
the mantle potential temperature T ) throughthe equa-
p
Table 1 tion:
Super-EarthModelParameters.
3 Rp
hT i=ǫ T = r2T(r)dr (2)
Mass Rp Rcore hρmi hgi m m p R3−R3 Z
(M⊕) (R⊕) (R⊕) (kgm−3) (ms−2) p c Rc
1 1 0.547 4480 9.8 SeeTajika & Matsui (1992)foraderivationofthisequa-
2 1.21 0.649 4960 13.4
tion. We approximate the adiabat as:
3 1.34 0.717 5470 16.4
5 1.54 0.814 5970 20.7 αg
T(r)=T +T ∆r (3)
Note. — Core mass fraction is held fixed at 0.3259, and p pc
p
scalingrelationsfromValenciaetal. (2006)areusedtodeter-
mineRp andRcore. Seetextfordetails. Values of ǫm for the different planet masses are given in
Table 5.
We use the scaling relations of Valencia et al. (2006) Inthemantle,heatisgeneratedbydecayofradioactive
to calculate the planetaryparametersfor planets of1, 2, elements. The heat flux from radioactive decay is domi-
3, and 5 Earth masses (1M⊕ = 5.97×1024 kg). Larger natedby238U,235U,232Th,and40K.Theheatproduced
planetsarenotconsideredherebecauseNoack & Breuer is calculated from the equation:
(2013) find that the peak likelihood for plate tectonics
occurs for planets between 1−5M⊕. The scaling laws Q(t)=ρm CiHiexp[λi(4.6×109−t)] (4)
of Valencia et al. (2006) assume a constant core mass X
fraction of 0.3259. The planetary and core radii then where Ci is the mantle concentration of the element by
scale by Ri ∼ Ri,⊕(Mp/M⊕)ai, where i = c (core) or mass, Hi is the heat production per unit mass, and λi is
p (planet), a = 0.247, a = 0.27, and values for the the decay constant. The decay constants, heat produc-
c p
Earth are indicated by ⊕. The average mantle density tion rates, and abundances relative to total uranium are
hρ i is calculated from the mantle mass and volume. taken from Turcotte & Schubert (2002). In this paper,
m
The average gravitationalacceleration hgi is found from we assume that all super-Earths have the same ratios of
GM /R2. ValuesfortheseparametersaregiveninTable radioactiveelementsasthepresentdayEarth. Thenom-
p p inal bulk silicate Earth (i.e., primitive mantle) contains
1. We take a constant water mantle mass fraction of
1.4×10−3 for the nominal models. This is equivalent ∼ 21 ppb U (McDonough & Sun 1995).
Using boundary layer theory, the heat flux out of the
to 4 ocean-masses(OM) of water for the Earth, where 1
mantle is given by:
OMisequalto1.39×1021kgH O.Wewilllaterexplore
2
the effect of variable water on the results. (T −T )
u s
q =k (5)
m
δ
2.2. Thermal Evolution Model u
where k is the mantle conductivity, δ is the boundary
Following models of Earth’s deep water cycle u
layerthickness,andtheT isthetemperatureatthebase
(e.g.McGovern & Schubert 1989;Schubert et al. 2001; u
of the boundary layer (see Figure 1). T is calculated
Sandu et al. 2011), we will first consider models heated u
fromhT iusing the adiabatic temperatureprofile ofthe
only from within (i.e., heat flux q = 0). These models m
c mantle.
are considered single-layerconvectionbecause the whole
The thickness of the upper boundary layer is given by
mantle convects. There is a conductive upper thermal
the global Rayleigh number:
boundary layer that governs heat loss from the surface.
For most of the lifetime of the Earth, such models have β β
Ra κη(T,P)Ra
cr crit
been shown to give good fits to the observedmantle vis- δ =Z = (6)
u
(cid:18) Ra (cid:19) (cid:18) gαρ ∆T (cid:19)
cosityandheatflux (Schubert et al. 2001). The ther- m
malevolutionmodelrequiressolutionofthemantleheat where Z is the thickness of the mantle, Ra is the
transfer equation:
Rayleighnumber of the whole mantle, Ra is the crit-
crit
ical Rayleigh number for convective instability, κ is the
dhT i
m
ρ C V =−A q +A q +V Q(t) (1) mantle thermal diffusivity, η(T,P) is the characteristic
m p m s s c c m
dt
mantle viscosity, and α is the thermal expansivity. ∆T
where ρ is the mantle density, C is the mantle heat isthetemperaturedropacrossthemantleminustheadi-
m p
capacity, V is the mantle volume, A , A are the sur- abatic temperature change. We use the mantle viscosity
m s c
face area’s of the planet and core, q and q are the con- calculated with hT i and hPi, which is characteristic of
s c m
ducted heat fluxes through the surface and core-mantle the whole mantle layer. This choice dictates the behav-
boundary(CMB),respectively,andQ(t)isthe heatpro- iorofourmodel,aswillbediscussedlater. Theviscosity
duced by radioactive decay within the mantle. In this parameterizationis described in Section 2.4.
model, there is no heat flux from the core, so the sec- Two other parameters are derived from the thermal
ond term on the right side vanishes. The temperature model. The areal spreading rate is the rate at which
modeledinanyparameterizedconvectionmodelissome- newoceancrustisbeingcreated. McGovern & Schubert
what ambiguous. Here, we will follow the convention of (1989) parameterized the areal spreading rate using the
McGovern & Schubert (1989)andtakethetemperature current volume of the oceans and the present day heat
to be the spherically-averaged mantle temperature. We flux,whichareunconstrainedforexoplanets. Instead,we
canrelatethis averagedtemperature to the temperature follow Sandu et al. (2011), who relate the arealspread-
of the mantle adiabat extrapolated to the surface (i.e. ing rate to the convective velocity u ) and the length of
c
4 Schaefer & Sasselov
rameterized as:
Table 2
NominalThermalModelParameters. r =ρ d S (10)
outgas m,v m
param. units UpperMantle LowerMantle Core whereρm,v isthedensityofvolatilesinthemantle,,dm is
thedepthofmelting,andS isthearealspreadingrateof
αRacrit ×10−5 K−1 11200 0.28R1a0.21 —— the mid-ocean ridges. In McGovern & Schubert (1989),
k Wm−1 K−1 4.2 4.2 — dm is kept at a constant value of 100 km. The ingassing
κ m2 s−1 10−6 10−6 — rate is parameterized as:
Cp Jkg−1 K−1 1200 1200 840
ρ kgm−3 3300 hρmia 8400 ringas =fbasρbasdbasSχr (11)
a SeeTable1 where f is the mass fraction of volatiles in the hy-
bas
drated basalt layer, ρ is the density of basalt, d is
bas bas
the average thickness of the basalt (held constant at 5
Table 3 km) and χ is an efficiency factor reflecting the incom-
r
Water CycleParameters.
plete transportof the water in the hydrated basalt layer
into the deep mantle. With all parametersexcept S and
param. name value
ρ held constant in equations (10) and (11), we found
m,v
χr regassingefficiency 0.03 that all planets necessarily reached a steady state (i.e.,
fχbdas hdyegdarsastiendgbeaffisaclitenfcrayction 00..0032 dMHd2tO,m = 0), where the mantle water abundance is
ρbas basaltdensity(kgm−3) 3000 given by setting ringas equal to routgas, and solving:
Lridge Mid-oceanridgelength 1.5×(2πRp)
Kγ SSoolliidduussddeepprreessssiioonnccooneffistiaenntt(Kwt%−γ) 04.735 ρm,v = MH2O,m = fbasρbasdbasSχr (12)
V d S
θ Meltfractionexponent 1.5 m m
DH2O Silicate/meltpartitioncoefficient 0.01 f ρ d χ V
bas bas bas r m
M = (13)
H2O,m d
m
the spreading centers (L ) where ocean crust is cre-
ridge
While a case may be made that the Earth is in a steady
ated:
state, there is little reason to suppose that this is a
S =2L u (7)
ridge c necessarycondition for all exoplanets experiencing plate
The convective velocity is determined by the convec- tectonics. We therefore follow here the volatile evolu-
tive layer overturn time from boundary layer theory tion model of Sandu et al. (2011), described briefly be-
(Schubert et al. 2001, ch. 8) using the equation: low. In this model, the depths of melting and hydration
are not held constant, but vary based on local tempera-
5.38κ(Rp−Rc) tures. We found that in these models, steady-state was
u = (8)
c δ2 rarely achieved. Parameters used in the volatile evolu-
u
tion model are given in Table 3. The outgassing rate is
We parameterize the length of the mid-ocean ridges as
given by:
1.5 times the planetary circumference. This parameteri-
zation is chosen to give the present day mid-ocean ridge r =ρ hF ihX iD Sχ (14)
outgas m melt melt melt d
(MOR)lengthontheEarthof∼60,000km. Wedescribe
where hF i and hX i are the average fraction of
results using smaller L values in Section 4. melt melt
ridge
meltingandaverageabundanceofwaterinthemeltover
Values for the parameters used in the thermal evolu-
the melt layer thickness, D , S is the areal spreading
tion model are given in Table 2. The value of Ra for melt
crit
rateandχ isthedegassingefficiency,whichaccountsfor
themantleistakenfromSchubert et al. (2001). Weuse d
incomplete transport of water to the surface. Whereas
constantvaluesfortheheatcapacity,thermalconductiv-
McGovern & Schubert (1989) used a constant value for
ity,thermalexpansivityandthermaldiffusivity,although
the melt layer thickness, Sandu et al. (2011) used the
theseparametersareallknowntobepressure-dependent.
mantlethermalprofileandtheperidotitesolidustodeter-
2.3. Volatile Evolution Model mine the melt layer thickness. The thermal profile used
is composed of the conductive upper thermal boundary
Volatile evolution models harken back to layerandthe uppermantleadiabat,andthe intersection
McGovern & Schubert (1989), repeated with vari- of this profile with the hydrated solidus curve for peri-
ations by many others. The volatile evolution model dotite determines where melt forms (see Fig. 1). Water
involvescalculationof outgassingandingassingratesfor dissolvedinsilicateslowerstheirsolidus(thetemperature
water based on mantle melting and surface hydration. atwhichpartialmeltingbegins),andthewaterpartitions
Therateofchangeofthewaterabundanceinthemantle preferentially into the melt. Using the parameterization
is given by combining the ingassing and outgassing ofKatz et al. (2003)forwetmelting,thesolidusdepres-
rates: sion is given by:
dM
Hd2tO,m =ringas−routgas (9) Tsol,wet =Tsol,dry−∆TH2O =Tsol,dry−KXmγelt (15)
This equation is solved simultaneously with the heat where K and γ are empirically determined constants for
transfer equation at each time. In the simplest form of peridotite (see Table 3). The melt fraction and water
McGovern & Schubert (1989),theoutgassingrateispa- abundance in the melt are determined where the mantle
Oceans on Super-Earths 5
thermalprofileis abovethe wetsolidustemperatureand
are given by: Table 4
ViscosityParameters.
θ
T −T
sol,wet
Fmelt = (16) parameter Olivine(wet) Perovskite(dry)
(cid:20)T −T (cid:21)
liq,dry sol,dry η0 (Pas) 1.08×1025 1×1021
X r 1.0 ...
Xmelt = H2O,m (17) Ea (kJmole−1) 335 300
DH2O+Fmelt(1−DH2O) Va (cm3 mole−1) 4.0 2.5
where X is the water mass fraction in the mantle, Rgas (Jmole−1 K−1) 8.314
H2O,m
DH2O is the silicate/melt partition coefficient, and θ is then given by:
an empirically determined exponent. The dry liquidus
andsolidusequationsaretakenfromZhang & Herzberg ǫ˙=A τnd−plnfr exp −Ea+PVa (20)
(1994) and Hirschmann et al. (2009), respectively. The crp H2O (cid:18) R T (cid:19)
gas
melt fraction and water fraction are averaged over the
wheren=1fordiffusioncreep,A isamaterialparam-
melt zone thickness at each timestep for use in equation crp
eter, d is the grain size in microns, E is the activation
(9). a
energy,V istheactivationvolume,andR istheideal
The return of water from the surface back into the a gas
gas constant. This equation shows that the response to
mantle through ocean plate subduction is described by:
stress depends on the pressure and temperature, as well
r =f ρ D Sχ (18) as the water fugacity (f ) in a non-linear way. The
ingas bas bas hydr r H2O
viscosity is then derived from the constitutive law:
where f is the mass fraction of water in a hydrated
bas
τ
basalt layer of thickness Dhydr, ρbas is the density of ηeff = (21)
basalt,and χ is the regassingefficiency, whichaccounts 2ǫ˙
r
for imperfect return of water to the mantle. The hy- We combine the constantparametersA andthe grain
crp
dratedlayerthicknessismeasuredfromthesurface,down size into a normalization factor eta to arrive at the ef-
0
tothedepthatwhichthetemperaturereachesthestabil- fective viscosity. The form of the effective viscosity is
ityboundaryofserpentinite(i.e.,serpentineisnotstable then:
atlowerdepths). Theupperthermalboundarylayerhas
E 1 1
alinearconductivetemperatureprofile(seeeq. 5),where η =η f−r exp a ( − )+
the temperaturechangewithdepth is(Tu−Ts)/δu. The eff 0 H2O (cid:20)Rgas T Tref
(22)
depthtotheserpentinestabilitytemperatureistherefore
1 PV P V
a ref a
given by: ( − )
R T T (cid:21)
gas ref
T −T T −T
serp s serp s
Dhydr =δu T −T =k q (19) We normalize so that η(Tref = 1600, XH2O= 500 ppm,
u s m P = 0)= 1021 Pa s, which gives reasonable viscosities
ref
whereT is thehighesttemperatureofserpentinesta- for the Earth. Parameters are given in Table 4, and are
serp
bility at pressures below ∼ 3 GPa and is ∼ 973 K taken primarily from Hirth & Kohlstedt (2003) for wet
(Ulmer & Trommsdorff 1995). The hydrated layer is diffusion in olivine. These authors measured values of r
necessarilysmallerthanthethermalboundarylayer,and of ∼0.7−1.2 on polycrystalline olivine samples. Recent
isfurtherrestrictedtoholdnomorewaterthanispresent workonSidiffusioninolivineandsingle-crystaldeforma-
atthesurfaceinagiveninstantinordertomaintainwa- tion experiments suggests that the dependence on water
ter mass-balance. We use the values for χ and χ de- abundancemaybe muchlower(r =0−0.33)(Fei et al.
r d
termined by Sandu et al. (2011). 2013; Girard et al. 2013). However, more work needs
A final word about parameterized volatile evolution to be done to reconcile these experiments with the poly-
models: As noted bya reviewer,the equationsdescribed crystallineexperiments. Inthispaperweuseresultsfrom
above assume that water transported into the mantle is the earlier studies for the nominal models, but we also
instantaneously transported throughout the mantle and explore the effect of different r values on our findings.
availablefor outgassing. Inthe realworld,ofcourse,the To convert from mass fraction water in the mantle to
mantle is not homogeneous and the spreading zone melt fugacityf ,weusetheformulationofLi et al. (2008),
H2O
centerswilllikelybecome dehydratedbeforethey canbe which is an empirical relationship between the water fu-
replenished by advection from subducted water. There- gacity and the concentration of water in olivine (C ,
OH
fore outgassing rates calculated here are upper limits on atomic H/106 Si):
the true values.
lnf =c +c lnC +c ln2C +c ln3C (23)
H2O 0 1 OH 2 OH 3 OH
2.4. Viscosity
where c = −7.9859,c = 4.3559,c = −0.5742, and
0 1 2
Water dissolved in silicate minerals such as olivine re- c = 0.0227. This relationship is derived from olivine
3
duces the mineral’s strength (e.g., Chopra & Paterson solubility data between ∼ 1373 − 1600 K and is only
1981), and therefore its viscosity. In experimental liter- strictly valid within that temperature range. However,
ature, the water dependence of the viscosity is parame- therelationshipvariesonlymarginallywithtemperature
terized via the water fugacity f , which depends on for a wide range of C and f , so we apply it to
H2O OH H2O
temperature and pressure. The rheological or constitu- the whole temperature range considered here for lack
tive law relating stress τ and strain rate ǫ˙ for olivine is of other available data. It is straightforward to convert
6 Schaefer & Sasselov
3000
Table 5
InitialTemperatureParameters.
Mass ǫm hTm,ii Tc,i 2500
(M⊕) (K) (K)
1 1.19 3000 3810
2 1.32 3330 4630
3 1.44 3630 5350 2000
5 1.64 4130 6640 K)
fromthemantlewatermassfractiontotheconcentration T (p
inC usingthemolecularweightsofwaterandolivine. 1500
OH
2.5. Initial parameters
1000
The initial temperature parameters for each of the
super-Earth models are given in Table 5. Previ-
ous parameterized convection models have shown that
500
initial temperatures do not significantly affect the 0 2000 4000 6000 8000 10000
thermal evolution beyond a few hundred Myr (see time (Myrs)
e.g.Schubert et al. 1980;McGovern & Schubert 1989).
24
Therefore, the present models are all started with the
same initial mantle potential temperature T , which is
p
the temperature of the mantle adiabat extrapolated to
the surface. The initial mantle potential temperature is
set to 2520 K for all models, which is equivalent to an
average mantle temperature of ∼ 3000 K for the Earth. 22
The averagemantle temperatures are calculated accord- s)
a
ingly for all super-Earthmodels. In models described in P
a later section, which include core evolution, we assume y (
an initial temperature contrast across the CMB of 100 sit
o
K. Therefore, the initial core temperature is set equal sc
to the temperature extrapolated along the mantle adia- vi20
bat to the CMB plus 100 K. We describe the effect of
different initial temperatures in Section 4.
3. RESULTS
For our nominal model, we focus on the evolution of
temperature and water abundances in both the mantle 18
0 2000 4000 6000 8000 10000
and on the surface. In the following sub-section, we in- time (Myrs)
troduce a second model which includes core-cooling and
usesadifferentcharacteristicmantleviscosityfortheup- x 10−3
1.4
per mantle that produces significantly different results.
3.1. Nominal Model - Single Layer Convection 1.2
Figures 2 and 3 show results for the nominal model
using parameters from Table 3 and a mantle water 1
abundance of ∼1400 ppm (4 OM for 1 M⊕). Figure 2
showstheevolutionofthemantlepotentialtemperature,
0.8
viscosity and water abundance for the nominal model er
at
using two different viscosity parameterizations. For the w
X
curvesshowninred,the activationvolumeis setto 0,so 0.6
the viscosity is pressure-independent. The blue curves
include a non-zero activation volume (see Table 4), so 0.4
the viscositydepends onpressureaswellas temperature 1 M⊕
and water fugacity. The viscosity is calculated with the 2 M⊕
average mantle temperature hTmi and the mid-mantle 0.2 3 M⊕
pressure. 5 M
⊕
For the pressure-independent viscosity, Fig. 2a seems 0
0 2000 4000 6000 8000 10000
to indicate that the large planets cool off more rapidly
time (Myrs)
than the smaller planets. The 5 M⊕ planet has a final
potential temperature ∼240 K lower than that of the Figure 2. Nominalmodelwithsinglelayerconvection. (a)man-
1 M⊕ planet. This is counter to standard intuition: ttleerpmotaesnstifarlactteimonp.eraRteudre,li(nbe)smreapnrtelseevnitsccoaslictuy,laatniodn(sc)domnaentwleitwhaa-
large planets should cool off more slowly than small
pressure-independent viscosity, blue lines represent the use of a
pressure-dependent viscosity. Line styles indicate planet mass ac-
cordingtothelegend.
Oceans on Super-Earths 7
1.8 has the lowest mantle viscosity, but this outgassing is
1 M much less complete than for the pressure-independent
⊕
1.6 case. The larger planets show steady but very gradual
2 M
⊕
outgassing, as their temperatures increase and their
1.4 3 M⊕ viscosities drop. In fact, the 5 M⊕ planet has delayed
1.2 5 M⊕ onset of outgassing by ∼1 Gyr, due to a very large
thermal boundary and low surface heat flux, both of
OM) 1 which can be attributed to the large initial mantle
ater (0.8 visTchoseitya.ctivation volume that we use here for the
w olivine viscosity is 4 cm3 mol−1, which is derived from
0.6 experimental data. However, the activation volume
for other silicates has been shown to decrease with
0.4
pressure (Stamenkovi´c et al. 2011). In their thermal
models, Stamenkovi´c et al. (2012) use activation vol-
0.2
umes for pressures at the planet’s CMB (2.5 cm−3
0 mol−1 for Earth). Therefore, a slightly lower activation
0 2000 4000 6000 8000 10000
volume may be more appropriate, particularly for
time (Myrs)
the larger planets. Intermediate values of V give
a
Figure 3. Surfacewaterabundanceforthenominalmodels. The mantle temperatures intermediate to those shown in
water abundance is equivalent to 4 ocean masses of water in the
Fig. 2a for the pressure-dependent viscosities. The
Earth-sizedplanet.
mantle water abundance drops less precipitously in
planets due to their smaller surface-area/volume ratios. early times, and ingassing is much slower than for the
However, note that these are potential temperatures pressure-independent case, but more complete than for
(i.e, the mantle temperature extrapolated to the surface the pressure-dependent model shown in Fig. 2b. The
along an adiabat). The average mantle temperature of pressure-dependence of the viscosity is therefore an
the 5 M⊕ planet is in fact hotter than the 1 M⊕, but by important parameter to include in the models. The
only ∼80 K. The adiabatic gradient is much steeper for value of activation volume will affect the final volumes
the large planets due to dependence on gravity, so the of the surface oceans and the mantle temperature.
near-surface temperatures are therefore much lower. As Figure 3 shows the surface abundance of water for the
seen in Fig. 2b, the hotter averagemantle temperatures nominal model for each of the super-Earths. This figure
of the large planets results in lower mantle viscosities is the surface corollary to Fig. 2c. We show only the
and therefore more rapid cooling. pressure-dependent models for clarity. Note that while
The lower near-surface temperatures of larger planets the relative abundance of water is the same for each
reduces the degree of melting, which strongly affects planet, the total planetary water mass is 4 (8, 12, 20)
their water cycles, shown in Fig. 2c. There is a sharp ocean masses for the 1 (2, 3, 5 M⊕) planet. The time-
decrease in mantle water abundance in the first 2 Gyr dependent behaviors of the planets do not scale simply
due to rapid early outgassing for all planets. As mantle with planet mass. The 1 M⊕ planet has a much more
temperature drops, near-surface temperatures drop significant outgassing phase at early times, followed by
below the solidus temperature, which halts melting and ingassing. The larger planets show much more gradual
outgassing. This is followed by a rapid ingassing period. outgassing, with delayed onset of outgassing for the 5
Ingassing is limited by the mass of the water at the M⊕ planet. Due to the delayed outgassing, the 5 M⊕
surface and the thickness of the hydrated surface layer, planet has less surface water than the 3 M⊕ planet
which is regulated by the surface heat flux. All but the for most of their lifetimes. However, the 3 and 5 M⊕
smallest planet have ingassed nearly all of their water planets have roughly equal surface water abundances at
by ∼2 Gyr. The deep water cycles of these planets have 10 Gyr. For smaller values of the activation volume, the
effectively ceased. However, it should be noted that for surface oceans will be substantially larger.
all planets, this leaves a small residual surface reservoir
of water that effectively cannot be lost to the mantle.
Forthepressure-dependentviscositycalculations(blue
3.2. Boundary Layer Convection
curves) both the average and potential temperatures of
the smaller planets cool more quickly than the larger Volatile evolution models have typically assumed
planets (see Fig. 2a). In fact, the largest planets single-layer convection, which has been shown to work
initially heat up substantially before beginning to cool, well for the Earth. However, thermal evolution models
so near-surface melting persists for much longer than in that neglect volatile evolution often assume that mantle
the pressure-independent case. The final temperature convectionis controlledby boundary layer instability, in
of the 5 M⊕ planet is ∼1000 K higher than that of the which convective instability is determined by the local
1 M⊕ planet. The slow cooling of the planets can be conditions at the boundary layers rather than the man-
attributed to the substantially larger viscosities in this tle as a whole. We will now describe how these mod-
model, due primarily to the pressure dependence of the els differ from the nominal model described above. For
viscosity(see Fig. 2b). The changesinwater abundance the boundary layer models, we assume mixed heating,
shown in Fig. 2c are much more gradual than with withanon-zeroheatfluxfromthecore(takenheretobe
pressure-independent viscosities. There is a rapid early bothsolidandisothermal)andaconductivelowerbound-
outgassing phase only for the smallest planet, which ary layer. Here we are less interested in reproducing the
8 Schaefer & Sasselov
Earth than in exploring the behavior of the models for little variability in the results for surface oceans when
differentplanetmasses. Whenusingthesamephysically- holding α constant.
plausibleparametershere,weshowthatthetwotypesof Results for the boundary layer convection model are
models give fundamentally different results, due to the discussed in the following subsection. Afterwards, we
choice of characteristic mantle viscosity. discuss how parameters affect the two different models.
For the boundary layer models, a second heat transfer
equation is needed for the core: 3.2.1. Results - Boundary Layer Convection Model
dT Results are shown in Figure 4. Core temperatures are
c
ρ C V =−A q (24)
c p,c c c c normalized to their initial values, because they vary by
dt
∼3000 K. The upper boundary layer is calculated using
wherevariablesareasineq. (1),butdefinedforthecore the pressure-dependent olivine viscosity using temper-
rather than the mantle. The core is isothermal, and so ature and pressure at the base of the boundary layer.
is characterized by a single temperature Tc. We neglect The lower boundary is calculated using the pressure-
radioactive decay in the core, as well as latent heat due dependent perovskite viscosity from Stamenkovi´c et al.
to possible corefreeze-out. The heat flux out ofthe core (2011) using temperatures and pressures at the CMB.
is given by: This viscosity does not depend on the mantle water
(Tc−Tl) abundance.
q =k (25)
c δ Mantle potential temperatures decrease more rapidly
c
and significantly here than for the nominal model. All
where δ is the lower thermal boundary layer, T is the
c c planets have nearly identical potential temperature
core temperature, and T is the mantle temperature at
l evolution, in contrast to the pressure-dependent results
thetopofthelowerthermalboundarylayer(seeFig. 1).
shown in Fig. 2a. The upper thermal boundary layer
WeusealocalRayleighnumbertodefinethethicknessof
thickness is comparable for all planets, so the upper
the lower boundary layer. The boundary layer thickness
mantle viscosity used here is only weakly dependent on
is determined by setting it equalto its criticalthickness,
pressure(bluecurves,Fig. 4c). Incontrast,theviscosity
the point at which it becomes unstable to convection.
of the lower thermal boundary layer is strongly depen-
1 dent on pressure, and therefore the core temperature
κη(T ,P )Ra 3
δ = c cmb crit,l (26) (Fig. 4b) is highly dependent on planet mass. The 1
c
(cid:18) gαρm(Tc−Tl) (cid:19) M⊕ planet’s core cools by ∼35%, whereas the 5 M⊕
planet’s core cools only ∼5%.
where η(T ,P ) is the viscosity of the lower bound-
c cmb Figure 4d shows the surface water inventories for
ary layer. For the lower boundary layer, we use a per-
the boundary layer model, with 4 (8, 12, 20) OM of
ovskite rheology. The choice of rheologies will be dis-
initial water. In comparison to Figure 3, it is obvious
cussed more below. The critical Rayleigh number of the
that significantly more water is outgassed from the
lower boundary layer is given by Ra = 0.28Ra0.21
(Deschamps & Sotin 2000), where cRriat,lis the global mantle here. The 5 M⊕ planet has a peak surface water
abundance of ∼8 OM, in comparison to ∼1.6 in Fig. 3.
Rayleighnumberusedinthe definitionofδ . Theglobal
u However, the residence time is much shorter. All of the
Rayleighnumber is defined here as in eq. (6), except for
planet’s lose most of their surface water inventory back
the characteristic viscosity. For the characteristic vis-
into the mantle by ∼4.5 Gyr. After complete ingassing,
cosity, we use the value defined by the temperature and
pressureatthebaseoftheupperthermalboundarylayer. the 1 M⊕ planethas the largestremaining surfacewater
abundance, with ∼0.2 OM at 10 Gyr. All planets have
Thishasasignifcanteffectontheoutcomeofthemodels
significantly lower surface water abundances at 10 Gyr
as shown below.
for the boundary layer models than for the nominal
The viscosity of olivine is used for the upper thermal
models shown in Figure 3. For the boundary layer
boundary layer, but olivine is not a stable phase at the
models, we find little difference between the pressure-
pressures and temperatures of the lower thermal bound-
dependent and pressure-independent viscosities, except
aries. We therefore use the viscosity law for perovskite
for the cooling of the core. The peak surface water
derivedbyStamenkovi´c et al. (2011)forthe lowerther-
abundances remain essentially unchanged, but the pres-
mal boundary. Parameter values are given in Table 4.
ence of surface water does persist for ∼0.5 Gyr longer
The lower mantle of super-Earths likely consists of per-
in the pressure-independent case. This is not surprising,
ovskitetransitioningtopost-perovskiteforlargerplanets
since the pressures of the upper thermal boundary layer
(see e.g. Valencia et al. 2006; Tackley et al. 2013). No
areverylowanddonotsignficantlyaffecttheviscosities.
experimentaldataexists onthe water-dependenceofthe
viscosity of perovskite, so the value of r is set to zero.
Stamenkovi´c et al. (2011) also give the dependence of
3.3. Dependence on parameters
the activation volume as a function of pressure. In the
nominalmodels, we use a value of 2.5 cm3 mol−1, which Many of the parameters used in the thermal and
is the value at the Earth’s CMB. volatile evolution models are poorly constrained for the
We use a slightly smaller value of α for the lower Earth, much less super-Earths. In the following section,
thermal boundary layer (see Table 2), which is taken we exam the effect of varying several of these parame-
from Tackley et al. (2013) for the Earth’s CMB. ters on the results of both the nominal and the bound-
Stamenkovi´c & Breuer (2014)foundthatscalingαwith ary layer convection models. We refer the reader to
planetary mass was as important as scaling ρ to the Stamenkovi´c & Breuer (2014) for an analysis of the ef-
m
calculationof planetary temperatures. However,we find fect of β, α , κ , and ρ on thermal evolution models.
m m m
Oceans on Super-Earths 9
2600 1
1 M
⊕
0.95
2 M
2400 ⊕
3 M
⊕ 0.9
5 M
⊕
2200
0.85
K)
T (K)p2000 T (c,init 0.8
/c
T 0.75
1800
0.7
1600
0.65
1400
0 2000 4000 6000 8000 10000 0 2000 4000 6000 8000 10000
time (Myrs) time (Myrs)
32 8
1 M
⊕
30 2 M
⊕
3 M
28 6 ⊕
5 M
⊕
s)26
a M)
P
sity (24 er (O4
o at
visc22 w
20 2
18
16 0
0 2000 4000 6000 8000 10000 0 2000 4000 6000 8000 10000
time (Myrs) time (Myrs)
Figure 4. Results for nominal models with boundary layer convection heating. (a) mantle potential temperature, (b) core temperature
normalizedtoinitialcoretemperature,(c)thermalboundarylayerviscosities(upper=blue,lower=red),(d)surfacewaterabundancein
oceanmasses. Linestylesindicateplanetmassaccordingtothelegend.
3.3.1. Fugacity coefficient outgassed significantly.
For the boundary layer convection model (Fig. 5b),
Figure5explorestheeffectofthefugacitycoefficientr the higher r value increases the amount of outgassing
ontheresultsofthepreviousmodelsforplanetsof1and for both the 1 and 5 M⊕ planets and shifts the peak of
5 M⊕. Other planet masses are not shown for clarity. outgassing to slightly earlier times. However, ingassing
This figure shows r values of 0.7, 1.0 (nominal), and also occurs more rapidly, so the oceans persist only
1.2. Hirth & Kohlstedt (2003) give values of 0.7−1.0 until ∼4 Gyr. The lower r value reduces the amount
for wet diffusion creep of olivine and r = 1.2 for wet of outgassing for both planets, and causes the surface
dislocation creep of olivine, which is why we chose these water to persist for longer. For the 1 M⊕ planet and r
values. It should be noted that the viscosities were not = 0.7, the water abundance is relatively constant over
renormalized with the change of r values. Therefore the the planet’s lifetime. There is about 0.4 OM of water
variances reflect the absolute changes in viscosity. remaining on the surface after 10 Gyr.
For the single layer convection model (Fig. 5a), the
initial outgassing phase is significantly stronger for
3.3.2. Total water abundance
r=1.2 for the 1 M⊕ planet, but the final abundance
is nearly the same as the nominal value. In contrast, Figure 6 shows the surface water abundances using
outgassing is slightly delayed for r = 0.7 but the final different initial mantle water abundances. The surface
water abundance is ∼0.2 OM larger. For the 5 M⊕ water abundances appear to be a fairly straightforward
planet, the larger r value increases the amount of water function of water abundance. Models with larger
outgassed, whereas the lower r value both further water abundance have larger surface inventories. The
delays outgassing and limits the total amount of water outgassing occurs earlier for the single layer convection
10 Schaefer & Sasselov
2.0
r = 1.0 12
r = 1.2
1.6 r = 0.7
10
8
M)1.2 M)
O O
er ( er ( 6
at at
w0.8 w
4
0.4
2
0 0
0 2000 4000 6000 8000 10000 0 2000 4000 6000 8000 10000
time (Myrs) time (Myrs)
Figure 5. Surface water abundances for (a) the nominal singlelayer convection model and (b) boundary layer convection models using
differentvaluesofthefugacityexponent r. Bluelinesarefor1M⊕,redlinesfor5M⊕.
3.5 18
2 OM
16 4 OM
3 6 OM
14
2.5
12
OM) 2 OM)10
ater (1.5 ater ( 8
w w
6
1
4
0.5
2
0 0
0 2000 4000 6000 8000 10000 0 2000 4000 6000 8000 10000
time (Myrs) time (Myrs)
Figure 6. Surface water abundances for (a) the nominal single layer convection model and (b) boundary layer convection models for
differentinitialwater abundances. Total water abundances areequivalent to 2,4and6ocean massesofwater forthe Earth-massplanet.
Abundances for the larger planets are the same in terms of mantle mass fraction of water. Colors indicate planet mass (blue 1 M⊕, red
5M⊕),andlinestylesindicatewaterabundance.
models (Fig. 6a), and oceans persist later for the Although models for the Earth show limited sensitiv-
boundary layer convection models (Fig. 6b). The ity to initial temperature (e.g. McGovern & Schubert
effect of varying initial water abundances produces (1989), Tajika & Matsui (1992)), this appears not
results similar to changes in the fugacity coefficient to be the case for the super-Earth models. Figure
shown in Fig. 5. However, one major difference to 7 compares results for the single layer model for the
note is that although both larger water abundances nominal starting mantle temperature of 2520 K, 2000
and larger fugacity coefficients increase outgassing, K, and 3000 K for the 1 (blue) and 5 (red) M⊕ planets.
the persistence of the oceans differs. Surface oceans Mantle temperatures for the smaller planets converge
persist longer (til ∼8 Gyr for the 5 M⊕ planet) for en- within ∼ 4 Gyr on the nominal results (Fig 2a) but
hancedwater,whereasincreasingthefugacitycoefficient for the hotter starting temperature the 5 M⊕ planet
shortened the ocean lifetime. Another thing to note remains persistently hotter throughout its lifetime. The
is that the surface water abundances for the boundary hotter initial temperatures effect the water cycle for all
layermodelwith2(10)OMofwaterarenearlyidentical. of the planets. Although the 1 M⊕ planet converges to
nearlythe sametemperatures,the initially hotterplanet
outgasses 3× more water within the first 500 Myr. The
3.3.3. Initial mantle temperature water is gradually ingassed over the planet’s lifetime,
but the final abundance of water remains slightly larger