The persistence of oceans on Earth-like planets: insights from the deep-water cycle
aa r X i v : . [ a s t r o - ph . E P ] J a n Draft version January 6, 2015
Preprint typeset using L A TEX style emulateapj v. 5/2/11
THE PERSISTENCE OF OCEANS ON EARTH-LIKE PLANETS: INSIGHTS FROM THE DEEP-WATERCYCLE
Laura Schaefer and Dimitar Sasselov
Harvard-Smithsonian Center for Astrophysics, 60 Garden St., Cambridge, MA 02138
Draft version January 6, 2015
ABSTRACTIn this paper we present a series of models for the deep water cycle on super-Earths experiencingplate tectonics. The deep water cycle can be modeled through parameterized convection modelscoupled with a volatile recycling model. The convection of the silicate mantle is linked to the volatilecycle through the water-dependent viscosity. Important differences in surface water content are foundfor different parameterizations of convection. Surface oceans are smaller and more persistent for singlelayer convection, rather than convection by boundary layer instability. Smaller planets have initiallylarger oceans but also return that water to the mantle more rapidly than larger planets. Super-Earths may therefore be less habitable in their early years than smaller planets, but their habitability(assuming stable surface conditions), will persist much longer. INTRODUCTION
In the post-Kepler era, the number of known planetsof Earth and super-Earth size continues to grow. TheTESS mission, which is a planned survey of the bright-est and closest stars, is expected to detect thousands ofsmall planets (Ricker et al. 2015). The TESS planets,by virtue of their closeness and the brightness of theirstars, will be much easier to follow up with ground-basedphotometric and spectroscopic techniques than the plan-ets of similar size detected by Kepler. Currently, themini-Neptune GJ 1214 is the only planet in the super-Earth mass range (1 − M ⊕ ) for which atmosphericmeasurements have been made (e.g., Wilson et al. 2014;Fraine et al. 2013; Berta et al. 2012). With future tele-scope resources such as GMT and JWST, we will be ableto characterize even more super-Earths, some of whichmay be in the habitable zones (HZ) of their stars. Theclassical habitable zone is defined as the orbital region inwhich an Earth-like planet can maintain liquid water onits surface given a variety of atmospheric compositions.However, the Earth’s atmosphere, although it permitssurface water to remain liquid, does not control the sup-ply of water on the surface. For that, we must look to themantle. In extending this to super-Earths, it is presentlyunclear what the effect of mass may be in controlling thesupply of water at the surface.Water on the Earth is found in two primary reser-voirs: (1) the surface oceans and (2) the silicate man-tle. Water, either in molecular form or as an H + orOH − group, can dissolve in silicate minerals. This istrue even of the so-called nominally anhydrous minerals(NAMs), which do not have hydrogen in their chemicalformula. Furthermore, the solubility of water in silicateminerals increases with pressure (Kohlstedt et al. 1996;Hirschmann et al. 2005).Geochemical studies of the materials that make upEarth’s mantle have long suggested that it may hold any-where from 1 to 10 times as much water as is found in thesurface oceans (e.g., Table 1 of Bounama et al. 2001),although more recent work seems to discount larger val- [email protected] ues in favor of values ranging from 0.5 to 2.5 oceanmasses (OMs) of water (Hirschmann & Kohlstedt 2012).The presence of water in these minerals strongly influ-ences their material properties, such as melting temper-atures, rheology, phase changes, electrical conductivity,etc. (see e.g., Hirth & Kohlstedt 1996; Karato 1990;Litasov & Ohtani 2007). Studies of the effect of wateron silicate rheology show that the presence of even minoramounts of water can reduce the viscosity by several or-ders of magnitude (Karato & Wu 1993). This can havea significant effect on the material flow of the mantle.In fact, it has been shown that the abundance of wa-ter on the Earth’s surface is controlled by the deep wa-ter/silicate cycle, which is tied to plate tectonics. Theabundance of water at the surface is a balance betweenthe rate of return via volcanic outgassing from the mantleat mid-ocean ridges (MORs) and the rate of loss to themantle through so-called ingassing, which is the return ofwater into the deep mantle by the sinking, or subducting,of water-rich oceanic seafloor. Much of this water is re-leased immediately back to the surface through shallow,water-induced volcanism. However, a small but signif-icant fraction of the water contained in the subductingoceanic slabs can be transported to deeper levels of themantle. For a detailed review of the exchange mecha-nisms between the surface and mantle, see Hirschmann(2006).The deep water cycle on Earth has been stud-ied through the use of parameterized convection mod-els incorporating a water-dependent viscosity (e.g.McGovern & Schubert 1989; Bounama et al. 2001;Sandu et al. 2011; Crowley et al. 2011). Parameter-ized convection models are simplified 1D models usinga parameterization derived from more complicated andexpensive 2D and 3D models of convection for differ-ent systems (e.g. spherical vs. plane-parallel, heatedfrom within vs. below, etc.). The parameterized mod-els rely primarily on two dimensionless parameters: theRayleigh number, which describes whether a system isunstable to convection, and the Nusselt number, whichcompares the convective and conductive heat flows. TheRayleigh number depends on the viscosity of the system, Schaefer & Sasselovwhich itself is dependent on temperature, pressure andwater fugacity. The abundance of water in the mantle,which helps determine the viscosity, evolves along withthe mantle temperature.Here we will use a parameterized convection model tostudy the deep water cycles of super-Earths. In thecurrent era of exoplanet studies, we are still search-ing for Earth-like planets in Earth-like orbits aroundSun-like stars, but what we have found and what wecan soon characterize, are super-Earths ( ∼ − R ⊕ , ∼ − M ⊕ ). Although super-Earths may have cool enoughsurfaces to retain liquid water, tectonics and mantleconvection plays an important role in controlling itsabundance at the surface. The question of whetherthese 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;Noack & Breuer 2013; Stamenkovi´c & Breuer 2014),although the issue has not been settled. Here, we as-sume plate tectonics as a starting point in order to applya similar water cycle model. The question this work thenaddresses is whether the different pressure regimes insidesuper-Earths affect the deep water cycle, and whethersurface oceans are tectonically sustainable on these plan-ets. This has important implications for the potentialhabitability of super-Earths.In this paper, we address these questions using a pa-rameterized thermal evolution model coupled with a wa-ter cycle model. We give a full description of the modelin Section 2. In Section 3, we describe results for mod-els which have either single layer convection or bound-ary layer convection for a number of super-Earth mod-els. We also describe the response of the models to rea-sonable variations in the parameter values. Section 4discusses implications for the persistence of oceans onsuper-Earths. METHODS
Mantle convection in the terrestrial planets can becontrolled either by the stability of the whole man-tle layer (single layer convection, e.g. Schubert et al.(2001)) or by the stability of two boundary layers: acold boundary layer at the surface and a hot bound-ary layer at the interface of the silicate mantle with themetallic core (see e.g. Turcotte & Schubert 2002). Heatfrom either secular cooling or decay of radioactive ma-terials (or both) is transferred by conduction throughthe boundary layers. The interior region is convectiveand thus approximately adiabatic. See Figure 1 for aschematic representation of the thermal profile. Singlelayer convection models typically neglect the heat fluxfrom the core, which is small, for simplicity and areentirely heated from within by radioactive decay. Nu-merical simulations are used to determine scaling lawsrelating the heat flux and mean temperature to the de-gree of convection. These simulations have been donefor a variety of different boundary conditions (free slip,no slip, etc.), and for different material properties (iso-viscous, highly temperature-dependent viscosity, etc.)(e.g., Honda & Iwase 1996; Deschamps & Sotin 2000;Dumoulin et al. 2005). Choosing the proper scaling re-lations is therefore important for a successful model.The parameter which is of the most impor-tance in determining the convective behavior is typi- R ad i u s R c T s R p T l T c δ c δ u solidus TemperatureT u T p Figure 1.
Schematic of temperature-depth profile. See text fordetails. cally the characteristic viscosity. Discussion in the lit-erature is often contradictory on where to define thecharacteristic viscosity: as an average value (see e.g.McGovern & Schubert 1989; Tajika & Matsui 1992;Sandu et al. 2011), or at the base of a boundary layer(see e.g. Deschamps & Sotin 2000; Dumoulin et al.2005; Stamenkovi´c et al. 2012). Here we will try bothtypes of models, starting with the single-layer modelstypically used with volatile evolution models, which useaverage viscosities. We will then look at a boundary layermodel. We explore both models here to see the effect onthe behavior of water, and note that we are not attempt-ing to reproduce the Earth, but explore two valid, butdifferent, models.The water cycle model is coupled to the thermal modelthrough the viscosity, which is dependent on the waterabundance in the mantle (see e.g. McGovern & Schubert1989; Sandu et al. 2011). Models of the Earth haveshown that the viscosity and mantle temperature createa feedback loop (Schubert et al. 2001; Crowley et al.2011). When the mantle is warm, convection is vigor-ous and the mantle cools quickly. As the temperaturedrops, the viscosity increases, causing the convection tobecome sluggish. Sluggish convection means that lessheat is removed from the mantle, causing it to heat up.The viscosity increases, and the cycle repeats. This cyclehas been shown to be enhanced by the presence of wa-ter (McGovern & Schubert 1989; Crowley et al. 2011).Crowley et al. (2011) describe how the water and tem-perature feedbacks interact for temperature and water-dependent viscosities. Cowan & Abbot (2014) studiedthe effect of sea floor pressure on a steady state modelof the deep water cycle, without considering the effect ofplanetary thermal evolution.In the following sections, we first describe selection ofthe super-Earth model parameters. We then describethe parameterized thermal model, followed by the watercycle model. Section 3 will discuss the results for thesemodels.
Super-Earth Models ceans on Super-Earths 3
Table 1
Super-Earth Model Parameters.Mass( M ⊕ ) R p ( R ⊕ ) R core ( R ⊕ ) h ρ m i (kg m − ) h g i (m s − )1 1 0.547 4480 9.82 1.21 0.649 4960 13.43 1.34 0.717 5470 16.45 1.54 0.814 5970 20.7 Note . — Core mass fraction is held fixed at 0.3259, andscaling relations from Valencia et al. (2006) are used to deter-mine R p and R core . See text for details. We use the scaling relations of Valencia et al. (2006)to calculate the planetary parameters for planets of 1, 2,3, and 5 Earth masses (1 M ⊕ = 5 . × kg). Largerplanets are not considered here because Noack & Breuer(2013) find that the peak likelihood for plate tectonicsoccurs for planets between 1 − M ⊕ . The scaling lawsof Valencia et al. (2006) assume a constant core massfraction of 0.3259. The planetary and core radii thenscale by R i ∼ R i, ⊕ ( M p /M ⊕ ) a i , where i = c (core) or p (planet), a c = 0.247, a p = 0.27, and values for theEarth are indicated by ⊕ . The average mantle density h ρ m i is calculated from the mantle mass and volume.The average gravitational acceleration h g i is found from GM p /R p . Values for these parameters are given in Table1. We take a constant water mantle mass fraction of1 . × − for the nominal models. This is equivalentto 4 ocean-masses (OM) of water for the Earth, where 1OM is equal to 1 . × kg H O. We will later explorethe effect of variable water on the results.
Thermal Evolution Model
Following models of Earth’s deep water cycle(e.g. McGovern & Schubert 1989; Schubert et al. 2001;Sandu et al. 2011), we will first consider models heatedonly from within (i.e., heat flux q c = 0). These modelsare considered single-layer convection because the wholemantle convects. There is a conductive upper thermalboundary layer that governs heat loss from the surface.For most of the lifetime of the Earth, such models havebeen shown to give good fits to the observed mantle vis-cosity and heat flux (Schubert et al. 2001). The ther-mal evolution model requires solution of the mantle heattransfer equation: ρ m C p V m d h T m i dt = − A s q s + A c q c + V m Q ( t ) (1)where ρ m is the mantle density, C p is the mantle heatcapacity, V m is the mantle volume, A s , A c are the sur-face area’s of the planet and core, q s and q c are the con-ducted heat fluxes through the surface and core-mantleboundary (CMB), respectively, and Q ( t ) is the heat pro-duced by radioactive decay within the mantle. In thismodel, there is no heat flux from the core, so the sec-ond term on the right side vanishes. The temperaturemodeled in any parameterized convection model is some-what ambiguous. Here, we will follow the convention ofMcGovern & Schubert (1989) and take the temperatureto be the spherically-averaged mantle temperature. Wecan relate this averaged temperature to the temperatureof the mantle adiabat extrapolated to the surface (i.e. the mantle potential temperature T p ) through the equa-tion: h T m i = ǫ m T p = 3 R p − R c Z R p R c r T ( r ) dr (2)See Tajika & Matsui (1992) for a derivation of this equa-tion. We approximate the adiabat as: T ( r ) = T p + T p αgc p ∆ r (3)Values of ǫ m for the different planet masses are given inTable 5.In the mantle, heat is generated by decay of radioactiveelements. The heat flux from radioactive decay is domi-nated by U, U, Th, and K. The heat producedis calculated from the equation: Q ( t ) = ρ m X C i H i exp [ λ i (4 . × − t )] (4)where C i is the mantle concentration of the element bymass, H i is the heat production per unit mass, and λ i isthe decay constant. The decay constants, heat produc-tion rates, and abundances relative to total uranium aretaken from Turcotte & Schubert (2002). In this paper,we assume that all super-Earths have the same ratios ofradioactive elements as the present day Earth. The nom-inal bulk silicate Earth (i.e., primitive mantle) contains ∼
21 ppb U (McDonough & Sun 1995).Using boundary layer theory, the heat flux out of themantle is given by: q m = k ( T u − T s ) δ u (5)where k is the mantle conductivity, δ u is the boundarylayer thickness, and the T u is the temperature at the baseof the boundary layer (see Figure 1). T u is calculatedfrom h T m i using the adiabatic temperature profile of themantle.The thickness of the upper boundary layer is given bythe global Rayleigh number: δ u = Z (cid:18) Ra cr Ra (cid:19) β = (cid:18) κη ( T, P ) Ra crit gαρ m ∆ T (cid:19) β (6)where Z is the thickness of the mantle, Ra is theRayleigh number of the whole mantle, Ra crit is the crit-ical Rayleigh number for convective instability, κ is themantle thermal diffusivity, η ( T, P ) is the characteristicmantle viscosity, and α is the thermal expansivity. ∆ T is the temperature drop across the mantle minus the adi-abatic temperature change. We use the mantle viscositycalculated with h T m i and h P i , which is characteristic ofthe whole mantle layer. This choice dictates the behav-ior of our model, as will be discussed later. The viscosityparameterization is described in Section 2.4.Two other parameters are derived from the thermalmodel. The areal spreading rate is the rate at whichnew ocean crust is being created. McGovern & Schubert(1989) parameterized the areal spreading rate using thecurrent volume of the oceans and the present day heatflux, which are unconstrained for exoplanets. Instead, wefollow Sandu et al. (2011), who relate the areal spread-ing rate to the convective velocity u c ) and the length of Schaefer & Sasselov Table 2
Nominal Thermal Model Parameters.param. units Upper Mantle Lower Mantle Core Ra crit . Ra . — α × − K − k W m − K − κ m s − − − — C p J kg − K − ρ kg m − h ρ m i a a See Table 1
Table 3
Water Cycle Parameters.param. name value χ r regassing efficiency 0.03 χ d degassing efficiency 0.02 f bas hydrated basalt fraction 0.03 ρ bas basalt density (kg m − ) 3000 L ridge Mid-ocean ridge length 1.5 × (2 πR p ) K Solidus depression constant (K wt% − γ ) 43 γ Solidus depression coeffiient 0.75 θ Melt fraction exponent 1.5 D H O Silicate/melt partition coefficient 0.01 the spreading centers ( L ridge ) where ocean crust is cre-ated: S = 2 L ridge u c (7)The convective velocity is determined by the convec-tive layer overturn time from boundary layer theory(Schubert et al. 2001, ch. 8) using the equation: u c = 5 . κ ( R p − R c ) δ u (8)We parameterize the length of the mid-ocean ridges as1.5 times the planetary circumference. This parameteri-zation is chosen to give the present day mid-ocean ridge(MOR) length on the Earth of ∼ L ridge values in Section 4.Values for the parameters used in the thermal evolu-tion model are given in Table 2. The value of Ra crit forthe mantle is taken from Schubert et al. (2001). We useconstant values for the heat capacity, thermal conductiv-ity, thermal expansivity and thermal diffusivity, althoughthese parameters are all known to be pressure-dependent. Volatile Evolution Model
Volatile evolution models harken back toMcGovern & Schubert (1989), repeated with vari-ations by many others. The volatile evolution modelinvolves calculation of outgassing and ingassing rates forwater based on mantle melting and surface hydration.The rate of change of the water abundance in the mantleis given by combining the ingassing and outgassingrates: dM H O,m dt = r ingas − r outgas (9)This equation is solved simultaneously with the heattransfer equation at each time. In the simplest form ofMcGovern & Schubert (1989), the outgassing rate is pa- rameterized as: r outgas = ρ m,v d m S (10)where ρ m,v is the density of volatiles in the mantle,, d m isthe depth of melting, and S is the areal spreading rate ofthe mid-ocean ridges. In McGovern & Schubert (1989), d m is kept at a constant value of 100 km. The ingassingrate is parameterized as: r ingas = f bas ρ bas d bas Sχ r (11)where f bas is the mass fraction of volatiles in the hy-drated basalt layer, ρ bas is the density of basalt, d bas isthe average thickness of the basalt (held constant at 5km) and χ r is an efficiency factor reflecting the incom-plete transport of the water in the hydrated basalt layerinto the deep mantle. With all parameters except S and ρ m,v held constant in equations (10) and (11), we foundthat all planets necessarily reached a steady state (i.e., dM H O,m dt = 0), where the mantle water abundance isgiven by setting r ingas equal to r outgas , and solving: ρ m,v = M H O,m V m = f bas ρ bas d bas Sχ r d m S (12) M H O,m = f bas ρ bas d bas χ r V m d m (13)While a case may be made that the Earth is in a steadystate, there is little reason to suppose that this is anecessary condition for all exoplanets experiencing platetectonics. We therefore follow here the volatile evolu-tion model of Sandu et al. (2011), described briefly be-low. In this model, the depths of melting and hydrationare not held constant, but vary based on local tempera-tures. We found that in these models, steady-state wasrarely achieved. Parameters used in the volatile evolu-tion model are given in Table 3. The outgassing rate isgiven by: r outgas = ρ m h F melt ih X melt i D melt Sχ d (14)where h F melt i and h X melt i are the average fraction ofmelting and average abundance of water in the melt overthe melt layer thickness, D melt , S is the areal spreadingrate and χ d is the degassing efficiency, which accounts forincomplete transport of water to the surface. WhereasMcGovern & Schubert (1989) used a constant value forthe melt layer thickness, Sandu et al. (2011) used themantle thermal profile and the peridotite solidus to deter-mine the melt layer thickness. The thermal profile usedis composed of the conductive upper thermal boundarylayer and the upper mantle adiabat, and the intersectionof this profile with the hydrated solidus curve for peri-dotite determines where melt forms (see Fig. 1). Waterdissolved in silicates lowers their solidus (the temperatureat which partial melting begins), and the water partitionspreferentially into the melt. Using the parameterizationof Katz et al. (2003) for wet melting, the solidus depres-sion is given by: T sol,wet = T sol,dry − ∆ T H O = T sol,dry − KX γmelt (15)where K and γ are empirically determined constants forperidotite (see Table 3). The melt fraction and waterabundance in the melt are determined where the mantleceans on Super-Earths 5thermal profile is above the wet solidus temperature andare given by: F melt = (cid:20) T − T sol,wet T liq,dry − T sol,dry (cid:21) θ (16) X melt = X H O,m D H O + F melt (1 − D H O ) (17)where X H O,m is the water mass fraction in the mantle, D H O is the silicate/melt partition coefficient, and θ isan empirically determined exponent. The dry liquidusand solidus equations are taken from Zhang & Herzberg(1994) and Hirschmann et al. (2009), respectively. Themelt fraction and water fraction are averaged over themelt zone thickness at each timestep for use in equation(9).The return of water from the surface back into themantle through ocean plate subduction is described by: r ingas = f bas ρ bas D hydr Sχ r (18)where f bas is the mass fraction of water in a hydratedbasalt layer of thickness D hydr , ρ bas is the density ofbasalt, and χ r is the regassing efficiency, which accountsfor imperfect return of water to the mantle. The hy-drated layer thickness is measured from the surface, downto the depth at which the temperature reaches the stabil-ity boundary of serpentinite (i.e., serpentine is not stableat lower depths). The upper thermal boundary layer hasa linear conductive temperature profile (see eq. 5), wherethe temperature change with depth is ( T u − T s ) /δ u . Thedepth to the serpentine stability temperature is thereforegiven by: D hydr = δ u T serp − T s T u − T s = k T serp − T s q m (19)where T serp is the highest temperature of serpentine sta-bility at pressures below ∼ ∼
973 K(Ulmer & Trommsdorff 1995). The hydrated layer isnecessarily smaller than the thermal boundary layer, andis further restricted to hold no more water than is presentat the surface in a given instant in order to maintain wa-ter mass-balance. We use the values for χ r and χ d de-termined by Sandu et al. (2011).A final word about parameterized volatile evolutionmodels: As noted by a reviewer, the equations describedabove assume that water transported into the mantle isinstantaneously transported throughout the mantle andavailable for outgassing. In the real world, of course, themantle is not homogeneous and the spreading zone meltcenters will likely become dehydrated before they can bereplenished by advection from subducted water. There-fore outgassing rates calculated here are upper limits onthe true values. Viscosity
Water dissolved in silicate minerals such as olivine re-duces the mineral’s strength (e.g., Chopra & Paterson1981), and therefore its viscosity. In experimental liter-ature, the water dependence of the viscosity is parame-terized via the water fugacity f H O , which depends ontemperature and pressure. The rheological or constitu-tive law relating stress τ and strain rate ˙ ǫ for olivine is Table 4
Viscosity Parameters.parameter Olivine (wet) Perovskite(dry) η (Pa s) 1 . × × r E a (kJ mole − ) 335 300 V a (cm mole − ) 4.0 2.5 R gas (J mole − K − ) 8.314 then given by:˙ ǫ = A crp τ n d − p lnf rH O exp (cid:18) − E a + P V a R gas T (cid:19) (20)where n = 1 for diffusion creep, A crp is a material param-eter, d is the grain size in microns, E a is the activationenergy, V a is the activation volume, and R gas is the idealgas constant. This equation shows that the response tostress depends on the pressure and temperature, as wellas the water fugacity ( f H O ) in a non-linear way. Theviscosity is then derived from the constitutive law: η eff = τ ǫ (21)We combine the constant parameters A crp and the grainsize into a normalization factor eta to arrive at the ef-fective viscosity. The form of the effective viscosity isthen: η eff = η f − rH O exp (cid:20) E a R gas ( 1 T − T ref )+1 R gas ( P V a T − P ref V a T ref ) (cid:21) (22)We normalize so that η ( T ref = 1600, X H O = 500 ppm, P ref = 0)= 10 Pa s, which gives reasonable viscositiesfor the Earth. Parameters are given in Table 4, and aretaken primarily from Hirth & Kohlstedt (2003) for wetdiffusion in olivine. These authors measured values of r of ∼ . − . r = 0 − .
33) (Fei et al.2013; Girard et al. 2013). However, more work needsto be done to reconcile these experiments with the poly-crystalline experiments. In this paper we use results fromthe earlier studies for the nominal models, but we alsoexplore the effect of different r values on our findings.To convert from mass fraction water in the mantle tofugacity f H O , we use the formulation of Li et al. (2008),which is an empirical relationship between the water fu-gacity and the concentration of water in olivine ( C OH ,atomic H/10 Si): lnf H O = c + c lnC OH + c ln C OH + c ln C OH (23)where c = − . , c = 4 . , c = − . c = 0 . ∼ − C OH and f H O , so we apply it tothe whole temperature range considered here for lackof other available data. It is straightforward to convert Schaefer & Sasselov Table 5
Initial Temperature Parameters.Mass( M ⊕ ) ǫ m h T m,i i (K) T c,i (K)1 1.19 3000 38102 1.32 3330 46303 1.44 3630 53505 1.64 4130 6640 from the mantle water mass fraction to the concentrationin C OH using the molecular weights of water and olivine. Initial parameters
The initial temperature parameters for each of thesuper-Earth models are given in Table 5. Previ-ous parameterized convection models have shown thatinitial temperatures do not significantly affect thethermal evolution beyond a few hundred Myr (seee.g. Schubert et al. 1980; McGovern & Schubert 1989).Therefore, the present models are all started with thesame initial mantle potential temperature T p , which isthe temperature of the mantle adiabat extrapolated tothe surface. The initial mantle potential temperature isset to 2520 K for all models, which is equivalent to anaverage mantle temperature of ∼ RESULTS
For our nominal model, we focus on the evolution oftemperature and water abundances in both the mantleand on the surface. In the following sub-section, we in-troduce a second model which includes core-cooling anduses a different characteristic mantle viscosity for the up-per mantle that produces significantly different results.
Nominal Model - Single Layer Convection
Figures 2 and 3 show results for the nominal modelusing parameters from Table 3 and a mantle waterabundance of ∼ ⊕ ). Figure 2shows the evolution of the mantle potential temperature,viscosity and water abundance for the nominal modelusing two different viscosity parameterizations. For thecurves shown in red, the activation volume is set to 0, sothe viscosity is pressure-independent. The blue curvesinclude a non-zero activation volume (see Table 4), sothe viscosity depends on pressure as well as temperatureand water fugacity. The viscosity is calculated with theaverage mantle temperature h T m i and the mid-mantlepressure.For the pressure-independent viscosity, Fig. 2a seemsto indicate that the large planets cool off more rapidlythan the smaller planets. The 5 M ⊕ planet has a finalpotential temperature ∼
240 K lower than that of the1 M ⊕ planet. This is counter to standard intuition:large planets should cool off more slowly than small T p ( K ) v i sc o s i t y ( P a s ) −3 time (Myrs) X w a t e r ⊕ ⊕ ⊕ ⊕ Figure 2.
Nominal model with single layer convection. (a) man-tle potential temperature, (b) mantle viscosity, and (c) mantle wa-ter mass fraction. Red lines represent calculations done with apressure-independent viscosity, blue lines represent the use of apressure-dependent viscosity. Line styles indicate planet mass ac-cording to the legend. ceans on Super-Earths 7 w a t e r ( O M ) ⊕ ⊕ ⊕ ⊕ Figure 3.
Surface water abundance for the nominal models. Thewater abundance is equivalent to 4 ocean masses of water in theEarth-sized planet. planets due to their smaller surface-area/volume ratios.However, note that these are potential temperatures(i.e, the mantle temperature extrapolated to the surfacealong an adiabat). The average mantle temperature ofthe 5 M ⊕ planet is in fact hotter than the 1 M ⊕ , but byonly ∼
80 K. The adiabatic gradient is much steeper forthe large planets due to dependence on gravity, so thenear-surface temperatures are therefore much lower. Asseen in Fig. 2b, the hotter average mantle temperaturesof the large planets results in lower mantle viscositiesand therefore more rapid cooling.The lower near-surface temperatures of larger planetsreduces the degree of melting, which strongly affectstheir water cycles, shown in Fig. 2c. There is a sharpdecrease in mantle water abundance in the first 2 Gyrdue to rapid early outgassing for all planets. As mantletemperature drops, near-surface temperatures dropbelow the solidus temperature, which halts melting andoutgassing. This is followed by a rapid ingassing period.Ingassing is limited by the mass of the water at thesurface and the thickness of the hydrated surface layer,which is regulated by the surface heat flux. All but thesmallest planet have ingassed nearly all of their waterby ∼ M ⊕ planet is ∼ M ⊕ planet. The slow cooling of the planets can beattributed to the substantially larger viscosities in thismodel, due primarily to the pressure dependence of theviscosity (see Fig. 2b). The changes in water abundanceshown in Fig. 2c are much more gradual than withpressure-independent viscosities. There is a rapid earlyoutgassing phase only for the smallest planet, which has the lowest mantle viscosity, but this outgassing ismuch less complete than for the pressure-independentcase. The larger planets show steady but very gradualoutgassing, as their temperatures increase and theirviscosities drop. In fact, the 5 M ⊕ planet has delayedonset of outgassing by ∼ mol − , which is derived fromexperimental data. However, the activation volumefor other silicates has been shown to decrease withpressure (Stamenkovi´c et al. 2011). In their thermalmodels, Stamenkovi´c et al. (2012) use activation vol-umes for pressures at the planet’s CMB (2.5 cm − mol − for Earth). Therefore, a slightly lower activationvolume may be more appropriate, particularly forthe larger planets. Intermediate values of V a givemantle temperatures intermediate to those shown inFig. 2a for the pressure-dependent viscosities. Themantle water abundance drops less precipitously inearly times, and ingassing is much slower than for thepressure-independent case, but more complete than forthe pressure-dependent model shown in Fig. 2b. Thepressure-dependence of the viscosity is therefore animportant parameter to include in the models. Thevalue of activation volume will affect the final volumesof the surface oceans and the mantle temperature.Figure 3 shows the surface abundance of water for thenominal model for each of the super-Earths. This figureis the surface corollary to Fig. 2c. We show only thepressure-dependent models for clarity. Note that whilethe relative abundance of water is the same for eachplanet, the total planetary water mass is 4 (8, 12, 20)ocean masses for the 1 (2, 3, 5 M ⊕ ) planet. The time-dependent behaviors of the planets do not scale simplywith planet mass. The 1 M ⊕ planet has a much moresignificant outgassing phase at early times, followed byingassing. The larger planets show much more gradualoutgassing, with delayed onset of outgassing for the 5 M ⊕ planet. Due to the delayed outgassing, the 5 M ⊕ planet has less surface water than the 3 M ⊕ planetfor most of their lifetimes. However, the 3 and 5 M ⊕ planets have roughly equal surface water abundances at10 Gyr. For smaller values of the activation volume, thesurface oceans will be substantially larger. Boundary Layer Convection
Volatile evolution models have typically assumedsingle-layer convection, which has been shown to workwell for the Earth. However, thermal evolution modelsthat neglect volatile evolution often assume that mantleconvection is controlled by boundary layer instability, inwhich convective instability is determined by the localconditions at the boundary layers rather than the man-tle as a whole. We will now describe how these mod-els differ from the nominal model described above. Forthe boundary layer models, we assume mixed heating,with a non-zero heat flux from the core (taken here to beboth solid and isothermal) and a conductive lower bound-ary layer. Here we are less interested in reproducing the Schaefer & SasselovEarth than in exploring the behavior of the models fordifferent planet masses. When using the same physically-plausible parameters here, we show that the two types ofmodels give fundamentally different results, due to thechoice of characteristic mantle viscosity.For the boundary layer models, a second heat transferequation is needed for the core: ρ c C p,c V c dT c dt = − A c q c (24)where variables are as in eq. (1), but defined for the corerather than the mantle. The core is isothermal, and sois characterized by a single temperature T c . We neglectradioactive decay in the core, as well as latent heat dueto possible core freeze-out. The heat flux out of the coreis given by: q c = k ( T c − T l ) δ c (25)where δ c is the lower thermal boundary layer, T c is thecore temperature, and T l is the mantle temperature atthe top of the lower thermal boundary layer (see Fig. 1).We use a local Rayleigh number to define the thickness ofthe lower boundary layer. The boundary layer thicknessis determined by setting it equal to its critical thickness,the point at which it becomes unstable to convection. δ c = (cid:18) κη ( T c , P cmb ) Ra crit,l gαρ m ( T c − T l ) (cid:19) (26)where η ( T c , P cmb ) is the viscosity of the lower bound-ary layer. For the lower boundary layer, we use a per-ovskite rheology. The choice of rheologies will be dis-cussed more below. The critical Rayleigh number of thelower boundary layer is given by Ra crit,l = 0 . Ra . (Deschamps & Sotin 2000), where Ra is the globalRayleigh number used in the definition of δ u . The globalRayleigh number is defined here as in eq. (6), except forthe characteristic viscosity. For the characteristic vis-cosity, we use the value defined by the temperature andpressure at the base of the upper thermal boundary layer.This has a signifcant effect on the outcome of the modelsas shown below.The viscosity of olivine is used for the upper thermalboundary layer, but olivine is not a stable phase at thepressures and temperatures of the lower thermal bound-aries. We therefore use the viscosity law for perovskitederived by Stamenkovi´c et al. (2011) for the lower ther-mal boundary. Parameter values are given in Table 4.The lower mantle of super-Earths likely consists of per-ovskite transitioning to post-perovskite for larger planets(see e.g. Valencia et al. 2006; Tackley et al. 2013). Noexperimental data exists on the water-dependence of theviscosity of perovskite, so the value of r is set to zero.Stamenkovi´c et al. (2011) also give the dependence ofthe activation volume as a function of pressure. In thenominal models, we use a value of 2.5 cm mol − , whichis the value at the Earth’s CMB.We use a slightly smaller value of α for the lowerthermal boundary layer (see Table 2), which is takenfrom Tackley et al. (2013) for the Earth’s CMB.Stamenkovi´c & Breuer (2014) found that scaling α withplanetary mass was as important as scaling ρ m to thecalculation of planetary temperatures. However, we find little variability in the results for surface oceans whenholding α constant.Results for the boundary layer convection model arediscussed in the following subsection. Afterwards, wediscuss how parameters affect the two different models. Results - Boundary Layer Convection Model
Results are shown in Figure 4. Core temperatures arenormalized to their initial values, because they vary by ∼ ⊕ planet’s core cools by ∼ M ⊕ planet’s core cools only ∼ M ⊕ planet has a peak surface waterabundance of ∼ ∼ ∼ M ⊕ planet has the largest remaining surface waterabundance, with ∼ ∼ Dependence on parameters
Many of the parameters used in the thermal andvolatile evolution models are poorly constrained for theEarth, much less super-Earths. In the following section,we exam the effect of varying several of these parame-ters on the results of both the nominal and the bound-ary layer convection models. We refer the reader toStamenkovi´c & Breuer (2014) for an analysis of the ef-fect of β , α m , κ m , and ρ m on thermal evolution models.ceans on Super-Earths 9 T p ( K ) ⊕ ⊕ ⊕ ⊕ T c / T c , i n i t ( K ) v i sc o s i t y ( P a s ) w a t e r ( O M ) ⊕ ⊕ ⊕ ⊕ Figure 4.
Results for nominal models with boundary layer convection heating. (a) mantle potential temperature, (b) core temperaturenormalized to initial core temperature, (c) thermal boundary layer viscosities (upper = blue, lower = red), (d) surface water abundance inocean masses. Line styles indicate planet mass according to the legend.
Fugacity coefficient
Figure 5 explores the effect of the fugacity coefficient r on the results of the previous models for planets of 1 and5 M ⊕ . Other planet masses are not shown for clarity.This figure shows r values of 0.7, 1.0 (nominal), and1.2. Hirth & Kohlstedt (2003) give values of 0 . − . r = 1 . r values. Therefore thevariances reflect the absolute changes in viscosity.For the single layer convection model (Fig. 5a), theinitial outgassing phase is significantly stronger for r =1.2 for the 1 M ⊕ planet, but the final abundanceis nearly the same as the nominal value. In contrast,outgassing is slightly delayed for r = 0.7 but the finalwater abundance is ∼ ⊕ planet, the larger r value increases the amount of wateroutgassed, whereas the lower r value both furtherdelays outgassing and limits the total amount of water outgassed significantly.For the boundary layer convection model (Fig. 5b),the higher r value increases the amount of outgassingfor both the 1 and 5 M ⊕ planets and shifts the peak ofoutgassing to slightly earlier times. However, ingassingalso occurs more rapidly, so the oceans persist onlyuntil ∼ r value reduces the amountof outgassing for both planets, and causes the surfacewater to persist for longer. For the 1 M ⊕ planet and r = 0.7, the water abundance is relatively constant overthe planet’s lifetime. There is about 0.4 OM of waterremaining on the surface after 10 Gyr. Total water abundance
Figure 6 shows the surface water abundances usingdifferent initial mantle water abundances. The surfacewater abundances appear to be a fairly straightforwardfunction of water abundance. Models with largerwater abundance have larger surface inventories. Theoutgassing occurs earlier for the single layer convection0 Schaefer & Sasselov w a t e r ( O M ) r = 1.0r = 1.2r = 0.7 w a t e r ( O M ) Figure 5.
Surface water abundances for (a) the nominal single layer convection model and (b) boundary layer convection models usingdifferent values of the fugacity exponent r . Blue lines are for 1 M ⊕ , red lines for 5 M ⊕ . w a t e r ( O M ) w a t e r ( O M ) Figure 6.
Surface water abundances for (a) the nominal single layer convection model and (b) boundary layer convection models fordifferent initial water abundances. Total water abundances are equivalent to 2, 4 and 6 ocean masses of water for the Earth-mass planet.Abundances for the larger planets are the same in terms of mantle mass fraction of water. Colors indicate planet mass (blue 1 M ⊕ , red5 M ⊕ ), and line styles indicate water abundance. models (Fig. 6a), and oceans persist later for theboundary layer convection models (Fig. 6b). Theeffect of varying initial water abundances producesresults similar to changes in the fugacity coefficientshown in Fig. 5. However, one major difference tonote is that although both larger water abundancesand larger fugacity coefficients increase outgassing,the persistence of the oceans differs. Surface oceanspersist longer (til ∼ ⊕ planet) for en-hanced water, whereas increasing the fugacity coefficientshortened the ocean lifetime. Another thing to noteis that the surface water abundances for the boundarylayer model with 2 (10) OM of water are nearly identical. Initial mantle temperature
Although models for the Earth show limited sensitiv-ity to initial temperature (e.g. McGovern & Schubert(1989), Tajika & Matsui (1992)), this appears notto be the case for the super-Earth models. Figure7 compares results for the single layer model for thenominal starting mantle temperature of 2520 K, 2000K, and 3000 K for the 1 (blue) and 5 (red) M ⊕ planets.Mantle temperatures for the smaller planets convergewithin ∼ M ⊕ planetremains persistently hotter throughout its lifetime. Thehotter initial temperatures effect the water cycle for allof the planets. Although the 1 M ⊕ planet converges tonearly the same temperatures, the initially hotter planetoutgasses 3 × more water within the first 500 Myr. Thewater is gradually ingassed over the planet’s lifetime,but the final abundance of water remains slightly largerceans on Super-Earths 11 w a t e r ( O M ) Figure 7.
Surface water for the single layer convection model.Solid lines show the nominal model (also shown in Fig. 2 and3), compared with models with either a higher (3000 K) or lower(2000 K) starting mantle potential temperature. Line styles indi-cate temperature according to the legend. Colors indicate planetmass (blue 1 M ⊕ , red 5 M ⊕ ). than for the colder starting planet. The 5 M ⊕ planetbegins outgassing much more rapidly, and continuessteadily outgassing for its lifetime. It ends with ∼ M ⊕ planet does not begin degassing until ∼ Mid-ocean ridge length
While we use a parameterization that producesthe Earth’s present day mid-ocean ridge length, wenote that the length of the ridges on Earth may havechanged throughout time. However, results for bothmodels change only slighlty with different values for L ridge . For the single layer model, the temperatureswith variable L ridge do not change, and surface waterabundances only slightly decrease. The effect on thewater abundance for the boundary layer models ismoderate, as shown in Figure 8. We show resultsfor models using lower values of L ridge , which wouldcorrespond with slower plate growth. The effect on w a t e r ( O M )
150 %100 %50 %
Figure 8.
Surface water abundances for the boundary layer con-vection model. Solid lines show the nominal model (also shownin Fig. 4), compared with models using either L ridge = 100%(dash) or 50% (dash-dot), respectively of the planetary circum-ference. Colors indicate planet mass (blue 1 M ⊕ , red 5 M ⊕ ). Thenominal model uses L ridge = 150% of the planetary circumference. temperature is minor for both the 1 and 5 M ⊕ planet,so we do not show it here. For the surface waterabundances, the lower L ridge values reduce the peaksurface water abundance and significantly prolong theingassing of the surface water after the initial outgassingphase. The 5 M ⊕ planet has ∼ M ⊕ planet has ∼ L ridge = 50%. Other parameters
Another parameter that can significantly affect themodel results is the viscosity parameterization. Wehave chosen here to normalize the viscosity to 10 Pas for a reference state of 1600 K, 0 Pa, and 500 ppmwater. This gives a reasonable viscosity for the Earth’smantle. However, there are other choices that couldbe made. Sandu et al. (2011) tuned their model tomatch the present day Earth’s viscosity and heat flux,by normalizing their viscosity to 2 . × Pa s at 2300K and 500 ppm of water. Using this reference valuefor the mid-mantle pressure, we get η = 2 . × Pa s, roughly an order of magnitude lower than thevalue used here. For the single layer model, using thisnormalization factor results in mantle temperaturescooler by ∼ M ⊕ planet heatup within the first 2 Gyr, rather than cooling. The 5 M ⊕ planet reaches a peak temperature of ∼ ∼ ∼
100 K. The peak surface waterabundances do not change from the nominal results,but the ingassing of the water back to the mantle takesa longer amount of time. The 1 M ⊕ has substantialsurface water until ∼ M ⊕ planet’ssurface water persists to ∼ ∼ DISCUSSION OF RESULTS
Stagnant lid regime
Many terrestrial planets, such as Mars and Venus, donot experience plate tectonics, but are in a stagnant-lid regime. For these planets, a thick lithosphere insu-lates the convecting portion of the mantle from the sur-face. Stagnant lids develop when planets are too cool tomaintain convection across the entire silicate layer. Itis likely that most planets can transition between platetectonics and the stagnant or sluggish lid regimes (Sleep2000). However, the mechanisms by which this transi-tion occurs are poorly understood. We can speculate,however, on how the transition would effect the modelsdescribed. The formation of a stagnant lid insulates themantle and allows the temperature to increase, whichwould lower the viscosity and increase convective vigor.Regassing would halt in this regime, as there is no trans-port mechanism for water to return to the mantle (seee.g. Morschhauser et al. 2011). For all of the planets inthe boundary layer model, the formation of a stagnantlid would likely extend the lifetime of the surface oceans,by halting regassing. The planets could also heat upenough to re-initiate mantle melting and degassing. Fu-ture work will look at the effect of such a transition onthe persistence of surface oceans.
Steady state versus thermal evolution
As discussed in section 2.3, using the volatile evolu-tion parameterization of McGovern & Schubert (1989),planets were found to achieve a volatile steady state.For a planet with a large global water abundance, mostof the water will be found at the surface of the planetfor the planet’s lifetime. For planets with global waterabundances smaller than the steady state value, allwater will be trapped within the mantle. Note thatCowan & Abbot (2014) considered a steady-state solu-tion to determine the water mass fraction necessary fora planet to be considered a water planet ( >
90% surfaceocean coverage). However, in the models consideredhere, steady state was never reached. In fact, the resultsdiscussed above show widely divergent outcomes for thesame planet based upon the characteristic upper mantleviscosity chosen. Once the upper mantle temperaturedrops below the peridotite solidus, the melt layer disap-pears and the volatile cycle effectively ceases, althoughslow ingassing may continue until the surface becomesdepleted of all water.
Ocean depths on Super-Earths
We have focused in the discussion of results onthe surface water abundances of the different modelsbecause this is a potentially observable parameter, andone that has implications for the possible habitabilityof super-Earths. Concerns have been raised about thehabitability of, or more importantly the ability of lifeto begin on, planets with global ocean coverage. Somefraction of continental surface, which provides a shallowwater environment, may be necessary for life to beginand for the evolution of complex life. Cowan & Abbot(2014) derive a maximum ocean depth that separatesplanets with continents from totally water-coveredplanets using a crustal buoyancy model. They find thatmaximum ocean depth scales with surface gravity as d max ∼ . g/g ⊕ ) − km. Applying this to the resultsfor the models above gives us the minimum surfacearea coverage on each planet. We show these resultsfor both the single layer and boundary layer modelsin Figure 9. The single layer convection planets haveslightly less than 50% areal coverage by oceans. The 5 M ⊕ planet has minimal ocean areal fraction until ∼ ⊕ planets have arealcoverages greater than 1, which indicates that they willhave global oceans. These planets are not likely to haveexposed continents. However, the 1 and 2 M ⊕ planetshave significant continental area for both convectionmodes, which indicates that these planets may be bestplaces to search for life. Additional planetary processes
The models presented here are not comprehensive pa-rameterizations of processes that can affect the waterbudget of a planet’s surface. One particularly impor-tant process is the loss of water from a planet’s at-mosphere (Wordsworth & Pierrehumbert 2013). Ourmodel is agnostic to the form that water takes at thesurface (i.e., water, ice, steam, etc.), other than byconstraining the surface temperature. However, thesurface temperature will vary with stellar type, or-bital distance, atmospheric composition, and age. AsWordsworth & Pierrehumbert (2013) show, loss of wa-ter vapor depends not only on the surface tempera-ture, but also on the CO abundance of the atmo-sphere. Significant loss of water vapor on hot CO -richplanets would reduce the amount of regassing, so themantle would become more dehydrated over time. Wealso neglect continental crust formation and weathering(Honing et al. 2014). Continental crust effectively actsas an insulating barrier, but also serves as a sink forradioactive elements which are extracted from the man-tle. Weathering of continents and sedimentation rates,which are enhanced by living organisms, can affect theregassing rate into the oceans. Honing et al. (2014) sug-gest that planets without life will evolve to have moresurface water (therefore lower continental coverage), anddryer mantles than planets with life. Weathering alsoacts as climate control by stabilizing atmospheric CO ,which affects the retention of water in the atmosphereas described above (Abbot et al. 2012). Much furtherwork needs to be done to truly understand the feedbacksceans on Super-Earths 13 s u r f a c e a r ea f r a c t i on s u r f a c e a r ea f r a c t i on ⊕ ⊕ ⊕ ⊕ Figure 9.
Minimum surface area coverage for the nominal models with a) single layer convection and b) boundary layer convection. Theminimum surface area coverage is determined assuming the maximum ocean depth based on the scaling of Cowan & Abbot (2014). Planetswith greater than 90% ocean coverage are considered water planets. that will contribute to the persistence of oceans on super-Earths. SUMMARY
We explored two different scaling parameterizationsfor plate tectonics planets using either single layerconvection or boundary layer convection. Mantle tem-peratures are significantly hotter for the first model, andthe surface water abundance lower, but more persistent.For many different parameters, smaller planets willhave initially larger surface oceans, but lose them morerapidly than larger planets. Larger planets show delayedoutgassing, which may compromise them as locations onwhich life can originate. The boundary layer convectionmodel cools very rapidly due to vigorous convectiondriven by both an upper and a lower thermal boundarylayer. The surface water abundances on the massiveplanets in the boundary layer model are extremelyhigh and suggest that these planets will have limited,if any, continental coverage. However, these oceanspersist for less than half of the planet’s lifetime. Uponcomplete ingassing of the oceans, the massive planetsare effectively tectonically dead, and therefore unlikelyto be habitable.Observations of rocky exoplanets in the 1 - 5 M ⊕ range are already producing very accurate planetradii and mass determinations (Dressing et al. 2014).Many of these exoplanets have ages, determined fromasteroseismic ages of their host stars. In the era ofJWST and large ground-based telescopes, some of thenearest exoplanets will have atmospheric spectroscopycapable of distinguishing different geochemical regimes,like some of the extremes described here. Understandingthe general features of the deep water cycle acrossrocky planets of different mass and structure might giveus unique windows into their interior through its sin-gular effect on atmospheric and surface water abundance.We thank Li Zeng for helpful discussions and an anony-mous referee for a detailed review that greatly enhanced the quality of this paper. REFERENCESAbbot, D.S., Cowan, N.B., & Ciesla, F.J. 2012, ApJ, 756, 178.Berta, Z. K., and 9 coauthors. 2012. ApJ, 747, 35.Bounama, C., Franck, S., von Bloh, W. 2001, HESS, 5, 569.Breuer, D., Labrosse, S., Spohn, T. 2010, SSRv 152, 449.Chopra, P.N., Paterson, M.S. 1981, Tectp. 78, 453.Cowan, N., Abbot, D. 2014, ApJ, 781, 27.Crowley, J. W., Gerault, M., O’Connell, R. J. 2011, E&PSL, 310,380.Deming, D., and 11 coauthors. 2009, PASP 121, 952.Deschamps, F., Sotin, C. 2000, Geophys. J. Int. 143, 204.Dressing, C., et al. 2014, in preparation.Dumoulin, C., Doin, M.P., Arcay, D., Fleitout, L. 2005, Geophys.J. Int. 160, 344.Fei, H., Wiedenbeck, M., Yamazaki, D., Katsura, T. 2013, Nature498, 213.Fraine, J.D. and 9 coauthors. 2013, ApJ, 765, 127.Frank, E.A., Meyer, B.S., Mojzsis, S.J. 2014. Icarus, 243, 274.Girard, J., Chen, J., Raterron, P., Holyoke, C. W. III. 2013,PEPI, 216, 12.Hirschmann, M.M. 2006, AREPS, 34, 629.Hirschmann, M.M., Kohlstedt, D. 2012, PhT, 65, 40.Hirschmann, M.M., Aubaud, C., Withers, A.C. 2005, E&PSL236, 167.Hirschmann, M. M., Tenner, T., Aubaud, C., Withers, A. C.2009, PEPI, 176, 54.Hirth, G., Kohlstedt, D.L. 1996, E&PSL 144, 93.Hirth, G., Kohlstedt, D.L. 2003,Honda, S., Iwase, Y. 1996, E&PSL, 139, 133.Honing, D., Hansen-Goos, H., Airo, A., Spohn, T. 2014, P&SS,98, 5.Karato, S. 1990, Nature. 347, 272.Karato, S. 2012, Icarus, 212, 14.Karato, S., Wu, P. 1993, Sci, 260, 771.Katz, R.F., Spiegelman, M., Langmuir, C.H. 2003, GGG, 4, 1073Kohlstedt, D. L., Keppler, H., Rubie, D.C. 1996, CoMP, 123, 345.Korenaga, J. 2010a, JGR, 115, B11405.Korenaga, J. 2010b, ApJ, 725, L43.Li, Z. A., Lee, C.A., Peslier, A.H., Lenardic, A., Mackwell, S.J.2008, JGR, 113, B09210.Litasov, K.D., Ohtani, E. 2007, GSA Special paper. 421, 115.McDonough, W.F., Sun, S.S. 1995, ChGeo, 120, 223.McGovern, P. J., Schubert, G. 1989, E&PSL, 96, 27.McNamara, A.K., van Keken, P. E. 2000, GGG, 1, 1027.Morschhauser, A., Grott, M., Breuer, D. 2011. Icarus, 212, 541.
Noack, L., Breuer, D. 2013, P&SS, 98, 41.O’Neill, C., Lenardic, A. 2007. JGR, 34, L19204.Ricker, G.R., Winn, J.N., Vanderspek, R., et al. 2015, J. Astron.Telesc. Instrum. Syst., 1, 014003.Sandu, C., Lenardic, A., McGovern, P. 2011, JGR, 116, B12404.Schubert, G., Turcotte, D.L., Olson, P. 2001, Mantle Convectionin the Earth and Planets. (New York, NY: Cambridge Univ.Press.)Schubert, G., Stevenson, D., Cassen, P. 1980. JGR, 85, 2531.Schubert, G., Turcotte, D.L., Olson, P. 2001. in MantleConvection in the Earth and Planets. (New York, NY: Camb.Univ. Press) Ch. 13.Sleep, N.H. 2000. JGR, 105, 17563.Stamenkovi´c, V., Breuer, D. 2014. Icarus, 234, 174. Stamenkovi´c, V., Breuer, D., Spohn, T. 2011, Icarus 216, 572.Stamenkovi´c, V., Noack, L., Breuer, D., Spohn, T. 2012, ApJ,748, 41.Tackley, P.J., Ammann, M., Brodholt, J.P., Dobson, D.P.,Valencia, D. 2013, Icarus, 225, 50.Tajika, E., Matsui, T. 1992, E&PSL, 113, 251.Turcotte, D.L., Schubert, G. 2002, Geodynamics. (New York, NY:Cambridge Univ. Press).Ulmer, P., Trommsdorff, V. 1995, 268, 858.Valencia, D., O’Connell, R.J., Sasselov, D. 2006, Icarus 181, 545.Valencia, D., O’Connell, R.J., Sasselov, D. 2007, ApJL, 670, L45.Wilson, P. A. and 13 coauthors. 2014, MNRAS, 438, 2395.Wordsworth, R.D., Pierrehumbert, R.T. 2013. ApJ, 778, 154.Zhang, J., Herzberg, C. 1994, JGRB, 99. 17729. −3 time (Myrs) X w a t e r v i sc o s i t y ( P a s ) v i sc o s i t y ( P a s ))