Structure and Dynamics of Cold Water Super-Earths

Structure and Dynamics of Cold Water Super-Earths: The Case of Occluded CH and its Outgassing

A. Levi Dept. of Geophysics & Planetary Science, Tel Aviv University, Tel Aviv, Israel 69978 Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA D. Sasselov Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA M. Podolak Dept. of Geophysics & Planetary Science, Tel Aviv University, Tel Aviv, Israel 69978
planets and satellites: atmospheres, Planets and satellites: composition, Planets and satellites: interiors, planets and satellites: oceans, planets and satellites: surfaces, planets and satellites: tectonics


In this work we study the transport of methane in the external water envelopes surrounding water-rich super-Earths. We investigate the influence of methane on the thermodynamics and mechanics of the water mantle. We find that including methane in the water matrix introduces a new phase (filled ice) resulting in hotter planetary interiors. This effect renders the super-ionic and reticulating phases accessible to the lower ice mantle of relatively low mass planets ( M) lacking a H/He atmosphere. We model the thermal and structural profile of the planetary crust and discuss five possible crustal regimes, which depend on the surface temperature and heat flux. We demonstrate the planetary crust can be conductive throughout or partly confined to the dissociation curve of methane clathrate hydrate. The formation of methane clathrate in the subsurface is shown to inhibit the formation of a subterranean ocean. This effect results in increased stresses on the lithosphere making modes of ice plate tectonics possible. The dynamic character of the tectonic plates is analysed and the ability of this tectonic mode to cool the planet estimated. The icy tectonic plates are found to be faster than those on a silicate super-Earth. A mid-layer of low viscosity is found to exist between the lithosphere and lower mantle. Its existence results in a large difference between ice mantle overturn time scales and resurfacing time scales. Resurfacing time scales are found to be  Ma for fast plates and  Ma for sluggish plates, depending on the viscosity profile and ice mass fraction. Melting beneath spreading centres is required in order to account for the planetary radiogenic heating. The melt fraction is quantified for the various tectonic solutions explored, ranging from a few percent for the fast and thin plates to total melting of the upwelled material for the thick and sluggish plates. Ice mantle dynamics is found to be important for assessing the composition of the atmosphere. We propose a mechanism for methane release into the atmosphere, where freshly exposed reservoirs of methane clathrate hydrate at the ridge dissociate under surface conditions. We formulate the relation between the outgassing flux and the tectonic mode dynamical characteristics. We give numerical estimates for the global outgassing rate of methane into the atmosphere. We find, for example, that for a planet outgassing can release to  molec s of methane to the atmosphere. We suggest a qualitative explanation for how the same outgassing mechanism may result in either a stable or a runaway volatile release, depending on the specifics of a given planet. Finally, we integrate the global outgassing rate for a few cases and quantify how the surface atmospheric pressure of methane evolves over time. We find that methane is likely an important constituent of water planets’ atmospheres.

1 Introduction

Planets that are intermediate in size and mass between Earth and Neptune are very common (Mayor et al., 2011; Howard et al., 2012; Fressin et al., 2013). They show a wide range of densities, with some of them indicating the existence of large water-rich interiors (water mantles), as in the case of exoplanet Kepler-68b (Gilliland et al., 2013). Such a planet would have a mass of several Earths, with about half of that mass in water - a water-rich mantle surrounding a core of heavy elements (rocky-iron composition). Such planets are unknown in our Solar System and the extremely high pressures in their mantles pose challenges to theoretical models. Interior modelling has treated the water as pure (see Zeng & Sasselov, 2013, and references therein). while this may be sufficient to reproduce the planetary radius for a given planet mass, the neglect of treating volatiles emerging from the core and mixed in with the water is a serious impediment to understanding the atmospheres such planets might have.

Interpreting the atmospheres of super-Earths, their structure and composition as gleaned from spectroscopic observations, is now upon us with the high-quality data for the nearby exoplanet GJ1214b (see Kreidberg et al., 2013, and references therein). Modelling such atmospheres in the habitable zone, as is the case of exoplanet Kepler-78f (Kaltenegger et al., 2013), is crucial in determining the outer boundary of the zone and whether a particular exoplanet is potentially habitable. The success of such modelling depends heavily on the gas fluxes emerging from the planet’s interior (as well as returning, in most cases). The volatile transport inside massive planets that are rich in water is poorly studied and understood, and to that end we have initiated a study of the entrapment of gases in the water-ice matrix (Levi et al., 2013).

Methane and other abundant gases trapped inside water super-Earth planets will be occluded in the solid-state water matrix while traversing the deep water mantle between the silicate-rich planet core and the planet’s surface. Chemical occlusion of methane takes the form of methane hydrate, also known as methane clathrate, or filled ice (at very high pressure). Filled ice, a novel occlusion structure (Loveday et al., 2001b), is particularly relevant for water super-Earths because of the high pressures reached in their water mantles. We studied the relevance of filled ice to super-Earths in Levi et al. (2013); this paper is a natural continuation of that work into understanding the phase transitions and dynamics near the planet’s surface for a complete physical picture of the outgassing process.

The paper consists of sections and appendices. In section we discuss the planetary crust and its dynamics by following the common path of describing the thermodynamics (of clathrates), their rheology, and the convective stability; we end by describing five different structural regimes of the crust. In section we describe briefly the deeper planetary structure. In section we discuss the thermal profile of the (water) mantle, and in section we use all these results to discuss the tectonics of water planets. We summarize our findings by describing resulting methane outgassing in section . In section we discuss a number of details - possible caveats and improvements of the model developed here, and finally we summarize our findings in section .

2 The Planetary Crust

In this section we present a thermal model for a planetary crust which is composed of a structure I (SI) methane clathrate hydrate. A full solution of the near surface heat transfer problem would involve modeling convection cells, which, in turn, requires a detailed simultaneous treatment of the thermodynamic and hydrodynamic behavior of the system. To avoid this complexity, we rely on a parameterization of the problem using dimensionless parameters, such as the Nusselt and Rayleigh numbers (Holman, 1968), and scaling techniques. Such an approach is generally in good agreement with more complex numerical modeling. In addition, the physical insight gained in this way can be extremely valuable (Solomatov, 1995).

Perhaps the most important dimensionless parameter describing this region is the Rayleigh number, , a measure of the importance of buoyancy in driving convection, which is given by:


Here is the acceleration of gravity, is the volume thermal expansivity, is the convective length scale, is the temperature difference driving the convection, is the thermal diffusivity, and is the kinematic viscosity.

2.1 Clathrate Thermodynamics

We first assess the numerical values of the different thermodynamic parameters for a SI methane clathrate. Clathrates, though crystalline in structure, have thermal conductivities that are closer to those of amorphous solids (see, e.g. Tse et al., 1997; Krivchikov et al., 2005a, and the references therein). This is due to the ”rattling” effect of the guest molecule, whose rattling frequencies may overlap the host network acoustic frequencies, resulting in a large guest-host thermal interaction. The thermal conductivity of methane clathrate is experimentally found to obey the following rules (Krivchikov et al., 2006):


Where is in erg s cm K. In Fig. 1 we show the thermal conductivity of methane clathrate hydrate as compared with the thermal conductivity of water ice Ih, the latter given by Slack (1980).

Figure 1: Thermal conductivity of methane hydrate as deduced from experiment (Krivchikov et al., 2006). For comparison we also show the thermal conductivity of ice Ih (Slack, 1980).

We describe the density, , of the methane hydrate by a third order Birch-Murnaghan EOS. Davidson (1983) gives a value of  g cm for the bulk mass density of an SI methane hydrate under low pressure. We use this as the zero pressure density. We take a bulk modulus of  GPa and a bulk modulus pressure derivative of based on the experimental data of Shimizu et al. (2002).

From Handa (1986) the heat capacity for SI methane hydrate is represented by:


The thermal diffusivity is defined as (Hobbs, 2010):


Using the experimental data mentioned above we derive the thermal diffusivity for SI methane hydrate (see Fig. 2).

Figure 2: Thermal diffusivity of SI methane hydrate (left panel) and ratio of volume expansivity of SI hydrate to water ice (right panel). SI hydrate data from (Hester et al., 2007) and water ice Ih data from Feistel & Wagner (2006).

The linear thermal expansivity, for SI and SII (structure II) hydrates, was determined experimentally by Hester et al. (2007) for different guest species. These authors have shown that it is the guest molecule contribution to the anharmonic part of the crystal intermolecular potential that enhances the thermal expansivity of hydrates above that of hexagonal water ice. The linear thermal expansivity was further found to be weakly dependent on the guest species, unless an excessively large guest is considered, and a small dependency on the hydrate structure was suggested. The linear expansivity, , suggested for a SI hydrate crystal is (Hester et al., 2007):

where is in Kelvins. The volume thermal expansivity, , is three times the linear expansivity. For its dependency on pressure we adopt the model of Fei et al. (1993) to finally give:


where is the linear thermal expansivity, is the bulk modulus and is its pressure derivative. We assign a value of unity to (see Ashcroft & Mermin, 1976). In Fig. 2 we also plot the volume thermal expansivity ratio between that for a SI clathrate and that for water ice Ih. From the plot it is clear that the ratio is well approximated by , deviating from this value by no more than for temperatures above  K.

2.2 Rheology

In order to construct the near surface thermal profile for a water planet we need to evaluate the kinematic viscosity for a clathrate layer. Unfortunately, little work has been done on the viscous behavior of clathrates both experimentally and theoretically (the exceptions are Stern et al. (1996) and Durham et al. (2003)). Experimentally it is not known whether the rheology of the clathrate is guest-dependent or not, and whether there is a rheological difference between SI and SII clathrates. We argue that the clathrate viscosity should depend on the type of guest molecule and clathrate structure. Because convective instability is highly dependent on the value of viscosity, we develop the expected relationships in detail below.

The dynamic viscosity relates the deviatoric stress tensor to the strain rate tensor in the following way:


Experiment and theory indicate that both diffusion and dislocation creep yield the following relation between stress and strain rate (Schubert et al., 2001):


Where is the pre-exponential factor, is the shear modulus, is the grain size, b is the Burgers vector, is the second invariant of the deviatoric stress tensor, and are an activation energy and volume respectively, and are the pressure and temperature and is Boltzmann’s constant.

Quite generally a strain rate equation of this type is based on the idea that solid state creep is intimately connected with the thermal creation of crystal imperfections and their thermal ability to migrate under applied stress, hence the Boltzmann factor. The parameters in Eq.(7) may vary for different thermodynamic regimes and for different applied stresses, as diffusional processes of different physical nature become active and dominate the solid state creep. Commonly, a single diffusion mechanism will dominate under given pressure, temperature and stress conditions (Durham et al., 1997).

The kinematic viscosity is defined as:


Durham et al. (2003) found for methane clathrate hydrate SI, for a confining pressure of  MPa and stresses of order  MPa, the following values for the parameters in eq.(8):  J mol,  cm mol, , and  MPa s, where the subindex stands for methane clathrate hydrate SI. These parameters were derived from an experiment on a laboratory made clathrate sample. Deciding whether these are transferable to a natural setting requires some thought of the dependency of the rheology on grain size. Due to the short time scale of an experiment one expects a laboratory sample to be composed of grains smaller than those composing a naturally formed sample that has time to ripen. Durham et al. (2003) report that their sample was composed of methane clathrate grains in the size range of m, an order of magnitude smaller than the grain size composing naturally formed methane clathrate bulk (see discussion below on methane clathrate grain sizes). When quantifying convective instability we will estimate the viscosity at temperatures higher than of the melting temperature (i.e. dissociation temperature) of methane clathrate. This temperature criterion is also maintained in the experiment of Durham et al. (2003). This high temperature regime suggests the parameters above represent a viscosity whose rate-controlling step is dominated by dislocation climb which is fairly insensitive to grain size (Kohlstedt, 2007). Therefore the transferability of the experimental parameters to our larger grain size case may be considered permissible. Although the fact that and not as expected from a creep solely dominated by dislocation climb (Kohlstedt, 2007) hints that other possible creep mechanisms may have also been at work during the experiment of Durham et al. (2003). One such possibility is grain boundary sliding which is probably enhanced due to the small grain sizes of the laboratory sample. Such a creep mechanism indeed encourages (Kohlstedt, 2007). If that is the case then the creep measured in Durham et al. (2003) is grain size sensitive and not easily transferable to a sample composed of much larger grain sizes. One must consider though that the larger grains composing a naturally formed sample of methane clathrate will make grain boundary sliding less efficient. In this case if grain boundary sliding is indeed folded in the experimentally derived parameters of Durham et al. (1997) then the viscosity they represent is a lower bound on the naturally forming methane clathrate dislocation viscosity whose larger grains make it stronger. An example of such strengthening due to increased grain size was measured in clinopyroxene (see Kohlstedt, 2007, and references therein).

Another point that must be considered is that methane hydrate survives to pressures as high as  GPa (Sloan, 1998), where the viscosity parameters of Durham et al. (2003) may no longer be applicable. In order to make extrapolations of the viscosity to higher pressures, we utilize a scheme proposed by Weertman (1970).

Low pressure experiments have shown that one may write the following:


where is the melting temperature and is a dimensionless constant which depends on the crystal structure (see Weertman & Weertman, 1975, and references therein). Weertman then further proposed the following extension to higher pressures:


This last transformation requires the melting curve to contain the information of how the activation energy and volume change with pressure. This has some support as, for a given stress, contours of constant viscosity which are functions of pressure and temperature, do seem to correspond fairly well to the contour of the melting curve. This is seen to be true for water ice (Durham et al., 1997) and for other substances as well (Poirier, 1985). Indeed this method has been utilized by Borch & Green (1987) in extrapolating the viscosity data for olivine to conditions in the Earth’s upper mantle and by Spohn & Schubert (2003) for pure water ice crusts in the Galilean satellites.

Inserting Eq.(10) into Eq.(8) yields for the viscosity:


The melting (i.e. dissociation) curve for clathrates is both guest molecule and crystal structure dependent, and this dependence enters into the viscous behavior of a given clathrate. We need to solve for the thermodynamic stability regime for a SI methane hydrate in order to extrapolate its viscosity. A full discussion of how to derive a clathrate hydrate thermal stability field and how to extrapolate it to high pressure ( GPa) is beyond the scope of this paper. For an in-depth explanation of thermal stability calculations we refer the reader to the works of van der Waals & Platteeuw (1959) and Sloan (1998).

Solving for the case of a methane clathrate hydrate we find the pressure at the first quadruple point to be  bar. By first quadruple point we mean the point where a clathrate hydrate transforms from being in equilibrium with ice Ih to being in equilibrium with liquid water. Thus the four phases: clathrate, liquid water, water ice Ih and methane vapour coexist. We further find that the dissociation (i.e. melting) curve beyond this pressure is well represented by the following polynomial:


where . For clarity we have plotted this melting curve, together with the three phase hydrate-ice Ih-vapour curve, on top of a phase diagram for pure water (see Fig. 3).

Figure 3: Clathrate hydrate SI dissociation curve (thick solid curve, blue in the on-line version) and clathrate hydrate structure-H (SH) dissociation curve (blue dashed curve in the on-line version). Filled circles represent available experimental data for the dissociation curve. Dashed line (red in the on-line version) is the ice-Ih melting curve. Dashed-dotted lines (green, cyan and magenta in the on-line version) are the melting curves for ice III, V and VI respectively. Arrow is pointing to the position of the quadruple point, where ice-Ih, methane clathrate hydrate, liquid water and methane vapour coexist. The stability field for methane clathrate hydrate is to the left of its dissociation curve.

For pressures below the first quadruple point we set the melting temperature equal to that of pure water ice for purposes of viscosity estimation. With the aid of Eq. (12) and Eq. (10) we solve for , taking into account that the experiment of Durham et al. (2003) was conducted at a confining pressure of  MPa. We find . The non-Newtonian viscosity, as expressed in Eq.(11) with , is not an intrinsic quantity of a crystal but rather is dependent on external conditions. In other words it depends on the applied deviatoric stress. Estimating the importance of the non-Newtonian dislocation creep in the planetary crust, therefore, requires an estimation of the second invariant of the deviatoric stress, . Golitsyn (1979) has shown that in steady state convection, where the rate of the work done by buoyancy exactly equals the rate of energy dissipation via friction, one has the following relation:


where the integral on the LHS is the rate of kinetic energy dissipation in a convecting cell, is the heat flux entering the cell and is the cell surface through which the heat enters the cell. In deriving Eq. (13) it is assumed that the density is constant (i.e. a Boussinesq fluid). In the Boussinesq approximation the background adiabatic temperature is constant (Schubert et al., 2001) and the temperature increase, , is confined to the thermal boundary layer, , thus:


For small viscosity contrasts (SVC) this will be true for both upper and lower boundary layers. For the stagnant lid regime (SL) is the cold boundary layer (Solomatov, 1995).

Given that the creep velocity under the cold thermal boundary layer is of order , the strain rate is of order , and in terms of scales Eq.(13) may be written as:


From boundary layer theory (see Solomatov, 1995, and references therein):


This results in the following estimate for the second invariant of the deviatoric stress tensor:


Assuming for the thermal boundary layer (which is the planetary crust) a length scale of  km, a temperature difference of  K, a surface gravity of  cm s and the thermal properties mentioned above, we estimate to be of order  Pa.

In order to properly describe the rheological behaviour of the crust considering dislocation creep alone is not sufficient. In particular, in case the crust is acted upon by low deviatoric stresses diffusional creep may best estimate the crustal rheology. Properly estimating the crustal rheology requires a viscosity map spanning crustal conditions, stating which creep mechanism minimizes the viscosity for varying crustal stress, pressure and temperature conditions. To that aim we shall also estimate the Newtonian (diffusion) creep for methane clathrate hydrate SI.

The microscopic manifestation of a Newtonian creep mechanism is the diffusion of lattice vacancies and interstitial molecules. The solid state viscosity related to this physical mechanism was formulated and analysed by Herring (1950). A slightly different formulation, given by Weertman & Weertman (1975), for the same mechanism is:


Where is a dimensionless constant for which we adopt the numerical value of (Frost & Ashby, 1982), is the average crystal grain diameter, is the atomic volume, and is the creep diffusion coefficient. For the atomic volume we use the value of  cm molec for a water molecule, derived from hydrogen bond length.

A question now arises about the diffusion of lattice vacancies and interstitial molecules in methane hydrate: Peters et al. (2008) studied the diffusivity of methane in a SI hydrate. The diffusion is considered to be due to thermal jumping of a methane molecule from a cage it occupies to a neighbouring vacant cage. Three jumping paths are considered, one is from a small cage to a large cage via a five membered water ring (pentagon face of a cage), and two different paths from a large cage to a neighbouring vacant large cage once via a five membered water ring and once through a six membered water ring (hexagonal face). They find the methane molecule to be too large to jump thermally through water rings without causing massive distortion to the water lattice, therefore, the authors invoke a water vacancy (defect) between the occupied and vacant cages, so as to lower the thermal barrier to jumping. The ability of the methane molecule to diffuse between cages becomes dependent not only on the degree of cage occupancy but also on the probability of water vacancy formation. From their results we may derive the following form for the diffusion coefficient of methane in a SI hydrate:


Here is the fraction of unoccupied water cages. By solving for the thermodynamic stability regime, for a SI methane hydrate, we obtain the variation of with temperature (see Fig. 4).

Figure 4: The fraction of unoccupied cages in a SI methane hydrate and its dependence on temperature. The blue (green) curve represents the degree of unoccupancy in the small (large) cage of a SI hydrate. The red curve is an averaged curve which takes into consideration that for every cubic unit cell there are eight cages, of which two are small and six are large.

Liang & Kusalik (2011) proposed a mechanism for creating and migrating HO defects in the hydrate water lattice, defects that may help reduce the thermal barrier to guest molecule jumping between cages, as required by the model of Peters et al. (2008). Liang & Kusalik (2011) proposed, based on MD simulations, that some small fraction of the hydrate cages may actually become occupied by a water molecule, leaving a defect in the lattice. This defect may help a guest molecule to diffuse in the lattice. The authors also found that these interstitial water molecules represent the most mobile defect in the hydrate lattice. When a water molecule, that occupies a cage, gets too near to one of the cage boundaries, whether it is a five (pentagonal face) or six (hexagonal face) membered water ring, it creates hydrogen bonds with molecules forming the water ring, which results in a metastable structure. This metastable structure then collapses by emitting a water molecule back to the center of the cage it came from, or to a first or second neighboring empty cage. In this way they derived a coefficient for self diffusion of water molecules within the hydrate lattice. Using their results we calculate the diffusion coefficient to be:


The diffusion coefficient, , in Eq. 18 is actually a weighted average of the contributions of water and methane to the diffusion. This weighted diffusion for a SI hydrate may be written as:


In the last equation we adopt the weighing procedure for lattice diffusion in multicomponent solids (see Weertman & Weertman, 1975), for the case of clathrate hydrate solid solutions. The weights in the denominator consider a cubic crystal unit cell, composed of water molecules and cages that with high probability are fully occupied, with a single methane molecule per cage.

Estimating the Newtonian viscosity also requires a value for the average diameter of the clathrate crystal grains, . Laboratory experiments show that single crystal grains of synthetically produced hydrates have diameters in order of several tens of micrometers. Smaller grains have higher growth rates. The smaller grains could be a consequence of the short time scales to which the laboratory experiment is confined. An examination carried out by Klapp et al. (2007) of actual geological samples of hydrates retrieved from the Gulf of Mexico and from Hydrate Ridge, revealed that in natural samples the crystal grain size diameter was in the range of m. The fact that naturally occurring hydrates have grains an order of magnitude larger than their synthetic counterparts was explained by these authors to be due to an Ostwald ripening process. In the latter process minimization of the free energy causes big grains to grow on the expense of smaller ones over a geological time scale. As a final note, the proper viscosity, whether it be Newtonian or non-Newtonian, is chosen to be such that the viscosity is a minimum for the given conditions.

2.3 Convective Stability Analysis

We are now ready to evaluate the thermal profile, in the near surface layer of a water planet, using the values derived above. As in Fu et al. (2010) we define the crust of a planet as the domain where conduction is the dominant mechanism for the heat transport. If the crust is of a radial dimension , we may write for it:


where and are the temperatures at the crustal base and planetary surface respectively and is the surface heat flux. As in Fu et al. (2010) we shall leave and as independent variables. Assuming the flux is only due to radioactive decay in the planetary metallic and silicate interior, and that the power released per gram of silicates and metals equals that for Earth, one may write:


where  W m (Turcotte & Schubert, 2002), and are the mass of Earth and the studied planet respectively, is the mass fraction of silicates and metals in the studied planet, and are the planetary radii respectively and and are the appropriate surface accelerations of gravity.

Both scaling analysis and assuming the thermal boundary layer is on the verge of convective instability are independent and equivalent techniques (Solomatov, 1995). The termination of the planetary crust occurs at some deep sublayer whose Rayleigh number is maximal and equals a critical value (). By deriving the width of the sublayer, that maximizes its Rayleigh number, the transition between the small viscosity regime and the stagnant lid regime was found by Solomatov (1995) to occur when:


where the LHS is the ratio of viscosities across the cold boundary layer and is the deviatoric stress power [see Eq. (8)].

First we assume a Newtonian viscosity. The condition that the crustal layer be on the verge of convective instability may be written, with the aid of Eq. (1), as:


where the subindex means each parameter is assigned the appropriate value for methane clathrate hydrate SI. Further, assuming the crust is in hydrostatic equilibrium, one may write:


where and are the crust bottom pressure and planetary surface pressure respectively. For the small viscosity contrast (SVC) we estimate the viscosity at the mid-layer temperature (Solomatov, 1995), using Eqs. (26) and (22).


The condition for convective instability (eq.25) can then be written as:


This last equation is a univariant equation for the pressure at the crustal base. We further assume for the case of the small viscosity contrast. This value is appropriate for (SVC regime for Newtonian fluids) as was shown for various wave numbers (see Schubert et al., 2001, and references therein). It is not likely that will be so low as to turn this value excessively large because the thermal conductivity of clathrates is low.

For the stagnant lid regime ( for Newtonian fluids) the viscosity is estimated at the temperature of the bottom of the crust:


The critical Rayleigh number for the stagnant lid regime is (Schubert et al., 2001):


The relation between the viscosity contrast, , and the creep mechanism enthalpy of activation, , was shown to be (McKinnon, 1999):


Where is a characteristic adiabatic temperature in the convecting sub-layer. It is important to note here that in the stagnant lid regime the actual temperature difference across the stagnant lid () does not constitute the temperature difference which drives the convection. Rather, between the convection cell and the stagnant lid there exists a boundary layer in which the viscosity increases exponentially with decreasing depth while at the same time the strain rate decreases exponentially from its value at the convective cell to its negligible value at the stagnant lid. The temperature difference across this boundary layer (referred to as the rheological temperature difference, ) is what drives convection in the stagnant lid regime (Solomatov, 1995). It was shown by Solomatov (1995) that the rheological temperature difference obeys:


It is also roughly given by:


From Eqs. (31-33) one may obtain a second order polynomial for , whose roots are:


To choose the physical root we note that . Therefore for the scenario where we expect [see Eq. (24)], and this is satisfied by the root with the plus sign in front of the square root. The other root will go to unity. The last equation is used in conjunction with Eq. (30) for the cases where .

2.4 Results and Discussion on the Planetary Crust

For the Newtonian creep mechanism, discussed above, we find there are five possible structures for the near surface thermodynamic behavior of water planets.

Regime I - In this case the surface temperature is low and so is the surface heat flux (low metallic and silicate content for a given surface gravity). For a low enough surface temperature even a very small surface atmospheric pressure will suffice to stabilize methane clathrates on the planetary surface (the dissociation pressure of methane clathrate at  K is  kPa). In this regime the low surface heat flux will result in relatively small increases in temperature with depth even though the thermal conductivity of clathrates is very small. Therefore, in this regime we expect a conductive planetary crust composed of methane clathrate beginning from the surface and ending in the depth where convective instability is reached.

Regime II - In this regime the surface temperature is still low enough so that the corresponding clathrate dissociation pressure is so low that probable surface atmospheric pressures ( bar) will suffice to stabilize clathrates on the planetary surface. Contrary to regime I the surface heat flux is now high (a large mass fraction of metals and silicates for a given surface gravity). In this case the very low thermal conductivity of clathrate hydrates will result in a steep temperature increase with depth. The resulting conductive thermal profile, in this case, will reach the hydrate dissociation (i.e. melt) curve before becoming convectively unstable, and will try to intrude into the liquid water phase region. Liquid water, having a very low viscosity, will introduce convection resulting in an adiabatic profile whose gradient is much steeper than that of the hydrate dissociation curve. This will drive the system immediately back to the hydrate stability regime which will again try to penetrate the liquid water regime. The result is a planetary crust whose upper part has a conductive profile and its lower part is restricted to the hydrate dissociation curve until convective instability is reached. The part of the crust with the on-melt behavior is expected to have a small viscosity contrast due to the fact that the viscous topology follows the melt curve, as explained above. We shall refer to this on-melt layer as the dissociation boundary layer (DBL).

Regime III - This regime is the counterpart of regime I, except that the surface temperature is now high enough so that the appropriate hydrate dissociation pressure could be higher than the atmospheric surface pressure (in case  bar, the quadruple point pressure). In this case ice Ih rather then methane hydrate will be stable on the planetary surface. Since the first quadruple point for methane hydrate is  bar, then for a surface gravitational acceleration of  m s the depth at which clathrate hydrates will become stable is, at most, of the order of  m. Since in this regime (as in regime I) the surface heat flux is low, and since the thermal conductivity of ice Ih is much higher than that of hydrates, the temperature will increase slowly with depth reaching the stability field for methane hydrates at a depth of, at most,  m. As hydrates become stable their lower thermal conductivity will lead to a faster temperature increase with depth and the system will converge to the situation described in regime I.

Regime IV - In this regime, as in regime III, the upper part of the crust is a thin sheet of water ice Ih and only at a depth of, at most, an order of  m, is methane hydrate stabilized. In this regime the high surface heat flux means the near surface behavior will converge to that depicted in regime II.

Regime V - In this regime the crust is made of hexagonal ice, and the surface temperature and/or the surface heat flux are high enough so that even in a relatively thin sheet of ice Ih the temperature may rise fast enough with depth to reach melting before entering the stability field for methane clathrate hydrate, where the solution would have been stuck in a regime III or IV type behaviors. In this regime we expect to find a subterranean ocean. Since a liquid ocean has a higher mass density than methane clathrate hydrate, any sublayer of clathrate hydrates trying to form will travel upward due to buoyancy consequently experiencing depressurization and decomposition.

In Fig. (5) we plot the domains of the five crustal structures as a function of the planetary surface temperature and surface heat flux, expressed in silicate and metal mass fraction [see Eq. 23], for a surface gravitational acceleration of  m s. The solid lines represent the boundaries between the crustal regimes for a  bar surface atmospheric pressure while the dashed lines and the markers are for the  bar surface atmospheric pressure scenario.

Figure 5: The parameter space for each of the five near surface thermodynamic regimes described in subsection . Regimes I and III will have a conductive upper thermal boundary layer ending with initiation of convection. In regimes II and IV a layer confined to the SI hydrate melt curve will separate the conductive boundary from the convection cell. In regime V an ocean will be formed, a configuration where the less dense pure methane hydrate will be unstable. The vertical solid line (red in the on-line version) marks the stability for methane hydrate on the planetary surface assuming a  bar surface atmospheric pressure. Markers and dashed curves represent the change in the parameter space of each regime when assuming a  bar surface atmospheric pressure.

The transition from regime I to regime II and the transition from regime III to IV are hardly affected by the increase in the assumed surface atmospheric pressure, from to  bar. This is because these regime transitions occur deeper in the crust where the effects of surface pressure are negligible. The most pronounced change, due to the increase in surface atmospheric pressure, is the shift to higher temperature (from to  K) in the boundary between regime I and III and from II to IV. As we have already mentioned, in order to be in regime V, the thermal profile in the ice Ih crust (whose depth is at most  m and whose thermal conductivity is high, relative to that of hydrates) must reach the melting curve for liquid water before stabilizing methane hydrates. This is hard to accomplish in such a narrow ice Ih layer and requires high surface heat fluxes or very high initial surface temperatures. From Fig. (5) we see that for the  bar atmosphere such an effect becomes possible for a minimum surface temperature of  K. For the  bar atmosphere the minimum stands at  K, as the higher surface atmospheric pressure results in an even thinner ice Ih layer which must reach the liquid water melting curve before stabilizing hydrates. We point out that for the ice Ih viscosity we adopt the formalism and parameters given in Spohn & Schubert (2003). We now wish to venture deeper into the planet, beginning with its internal structure.

3 The Planetary Structure

We consider in our model a differentiated body composed of four distinct regions: An Fe core, a surrounding silicate mantle, and an outer water mantle which is itself divided into two regions. As we are interested in the details of the transport and cycling of CH in water planets we mainly focus on the fine structure of the water mantle where such transport exists. We divide the water mantle into a high pressure region composed of filled ice, and a lower pressure region where methane clathrate hydrate is stable. The filled ice is a water ice polymorph created under high pressure in the presence of methane (see Levi et al. (2013) - hereafter called paper I). With regards to the iron core and silicate mantle we will restrict ourselves to making simpler though adequate assumptions.

For the iron core we use the Vinet EOS with the parameters given in Seager et al. (2007), which, according to these authors, describes the phase of Fe up to a pressure of  GPa. This pressure is never exceeded in any of the planets we consider. The external water layer we assume ensures that the contribution from low pressure silicate phases is negligible. We therefore model the silicate mantle to be solely composed of the perovskite phase of MgSiO, where for the EOS we use the fit suggested in Seager et al. (2007), which smoothly connects a fourth-order Birch-Murnaghan EOS with the Thomas-Fermi-Dirac EOS.

For methane filled ice Ih we adopt a third order Birch-Murnaghan EOS with bulk modulus,  GPa, derived from Hirai et al. (2003). We also use this EOS for the methane clathrate hydrate, with  GPa and determined experimentally by Shimizu et al. (2002). In tables 1-3 we present various internal structure results for our , and planets. For each mass we assume both a and a water mass fraction. For the case of the planet we examined a wider range of ice mass fractions ranging from to as low as in order to better understand the dependence of mantle convection on this parameter.

Ice Fraction g P D P R P R
% (m s) (GPa) (km) (GPa) (km) (GPa) (km)
3 13.4 7 377 822 3899 300 7363
5 12.9 11 576 814 3874 300 7299
10 12.1 22 990 794 3808 301 7143
15 11.5 32 1370 771 3741 301 6988
20 10.9 41 1732 747 3672 300 6833
25 10.4 50 2082 722 3600 298 6675
30 10.1 59 2382 696 3525 295 6512
35 9.7 67 2753 669 3446 291 6346
40 9.4 75 3050 641 3363 287 6172
45 9.2 83 3351 612 3275 282 5991
50 8.9 90 3680 582 3182 276 5800
55 8.7 97 4011 550 3082 270 5598
60 8.4 104 4363 517 2973 262 5383

Note. – g-surface gravity, P-silicate and water mantle boundary pressure, D-water mantle depth, P-pressure at the center, R-iron core radius, P-iron core and silicate mantle boundary pressure, R-distance from center to water mantle.

Table 1: Internal Structure for 2M Planet
Ice Fraction g P D P R P R
% (m s) (GPa) (km) (GPa) (km) (GPa) (km)
25 16.3 119 2520 1786 4537 716 8566
50 14.2 228 4472 1442 4012 679 7402

Note. – g-surface gravity, P-silicate and water mantle boundary pressure, D-water mantle depth, P-pressure at the center, R-iron core radius, P-iron core and silicate mantle boundary pressure, R-distance from center to water mantle.

Table 2: Internal Structure for 5M Planet
Ice Fraction g P D P R P R
% (m s) (GPa) (km) (GPa) (km) (GPa) (km)
25 23.4 246 2900 3763 5330 1483 10180
50 20.9 493 5088 3055 4713 1436 8739

Note. – g-surface gravity, P-silicate and water mantle boundary pressure, D-water mantle depth, P-pressure at the center, R-iron core radius, P-iron core and silicate mantle boundary pressure, R-distance from center to water mantle.

Table 3: Internal Structure for 10M Planet

4 The Mantle Thermal Profile

For the purposes of this section we define the planetary water mantle as the water ice layer bounded from above by the crust and from below by the silicate mantle. Starting from the planetary crust and making our way deeper into the planet we need to make a distinction between crustal regimes I/III and II/IV. In regimes I and III the surface heat flux is low enough so that the crust becomes unstable with respect to convection within the methane hydrate stability field. In these scenarios a convective cell, still in the methane hydrate layer, lies immediately underneath the conductive crust. As explained above, in regimes II and IV, underneath the conductive crust lies a layer confined to the hydrate dissociation (i.e melt) curve, which separates the crust from the underlying convective cell. We shall refer to this layer as the dissociation boundary layer (DBL) whose radial dimension, , we constrain by assuming it is on the verge of convective instability (see Fu et al., 2010, and references therein). Since the DBL follows the dissociation curve we do not expect a large viscosity contrast across its length.

Assuming the DBL is on the verge of convective instability brings up the question of the proper formalism for the convective threshold calculation. In Fig. 6 we present a map of the Newtonian and non-Newtonian viscosities of methane hydrate as a function of temperature for a confining pressure of  bar and the parameters given in section . The lower viscosity of the two general mechanisms (i.e diffusion versus dislocation creep) is the proper viscosity (Schubert et al., 2001). It is clear from Fig. 6 that for the viscosity parameters adopted and for the temperatures expected at the DBL, dislocation creep will be slightly more efficient than diffusion creep, resulting in a non-linear stability problem. Tough, as discussed above, the parameters adopted for describing dislocation creep in methane hydrate could actually represent a lower bound on the dislocation viscosity due to the grain size difference between experimental and naturally forming clathrate samples. Also, for the viscosity parameters adopted and at the probable temperatures of the DBL both diffusion creep and dislocation creep do not show many orders of magnitude difference. In addition, the non-linear dislocation creep requires finite disturbances whose existence is not certain. For these reasons we shall treat the DBL stability problem with the formalism of the linear diffusion creep viscosity.

Assuming the DBL is on the verge of convective instability yields the following relation [see Eq. 25]:


where we take advantage of the fact that the thermal profile in the DBL follows the hydrate dissociation curve. Therefore, the temperature at its bottom is the hydrate dissociation temperature at the pressure prevailing at the base of the DBL layer, . The temperature at the base of the crust is . If the dissociation boundary layer is, relatively, narrow (which will prove to be the case) then Eq. 26 may be used to yield a univariant equation for the DBL base pressure:


where is the pressure at the base of the crust. Due to the mild viscosity contrast that is expected across the DBL (because it follows the methane hydrate melt curve) we assume a value of for its critical Rayleigh number (see Turcotte & Schubert, 2002). The viscosity is approximated by its value at the average temperature of the DBL layer.

Figure 6: Diffusion creep (solid blue curve) and dislocation creep (dash and dash-dot curves) dependence on temperature for a solid of CH clathrate hydrate and a reference pressure of  bar. Up to a temperature of K diffusion yields the lower viscosity while at higher temperatures dislocation creep yields the lower viscosity. Non-Newtonian viscosity is stress dependent and the viscosity curves shown are confined between stresses of - Pa.

For our choices for the ice mass fractions we find that all the planets we investigated have a dissociation boundary layer. The radial dimension of the DBL is listed, for six of our studied planets, in table 4. is approximately  km for all six bodies. The radial dimension of the crust varies from a few hundred meters up to a kilometer for the six bodies listed. Changing the composition from to water, for a given body, widens the crust by a factor of about two. The lower mass planet has a slightly thicker crust than the more massive planet.

Parameters (km) (km)
M=2M  25 HO 0.61 1.36
M=2M  50 HO 1.13 1.52
M=5M  25 HO 0.39 1.11
M=5M  50 HO 0.71 1.22
M=10M  25 HO 0.27 0.94
M=10M  50 HO 0.48 1.02

Note. – For six characteristic planets, out of our studied water planets, we give the radial dimension of the planetary crust () and the radial dimension of the dissociation boundary layer ().

Table 4: Crust and DBL data

Regardless of whether the convecting cell underlies the conductive crust or the DBL it will follow a thermodynamic adiabatic profile into the abyss. The adiabatic temperature gradient is (Schubert et al., 2001):


The heat capacity measurements of Handa (1986) [see Eq. 3] are indeed appropriate for the planetary crust as these measurements were conducted at a low confining pressure (approximately  MPa). The low pressure increases the probability that the sample studied contained pores filled with methane gas. The higher confining pressure (about  MPa) in the experiments of Waite et al. (2007) indicate a higher compaction and a lower probability for the existence of gas filled pores, which is more appropriate for the higher pressures in the convecting cell. We therefore assume for the clathrate hydrate section of the convecting cell the following heat capacity:


From the base of the DBL we follow the adiabat in the convecting clathrate hydrate cell to find that in all the six planets we studied a point of intersection is reached where the adiabat tries to cross the SH methane hydrate dissociation curve into the liquid water regime. For all our planets this happens at a pressure of between  GPa. To understand the thermal behavior beyond this point one has to calculate the adibatic profile in liquid water. The temperature dependence of the volume thermal expansivity can be calculated from the tables of Kell (1975). We calculate its pressure dependence using a scheme similar to what we used for the clathrate expansivity:


where , , and . For the bulk modulus of liquid water and its derivative we assume numerical values of  GPa and respectively (Manghnani et al., 1999). For we again suggest a value of unity as for hydrates. The reference pressure here is  bar.

For the heat capacity of liquid water we adopt without change the formulation in Waite et al. (2007). To the liquid water bulk mass density dependence on temperature (see Waite et al., 2007, and references therein) we add a pressure dependency, yielding the following form:


Calculating the adiabat in the liquid water, beyond the point of intersection just mentioned, we find its gradient () to be even larger than the gradient of the SH methane hydrate dissociation curve. In that case, as in the DBL, the penetration into the liquid water regime drives the thermal profile back to the SH methane hydrate stability field which in turn will try to re-penetrate the liquid regime. We are once more in a situation where the thermal profile is confined to the methane hydrate dissociation curve, only now it is the dissociation curve for the SH methane hydrate. This implies that, under our assumptions, our studied planets do not have a liquid subterranean ocean.

The SH methane hydrate dissociation curve is based on the experimental data points of Dyadin et al. (1997) which are given up to a pressure of  GPa. Extrapolating these experimental data points we find the SH hydrate dissociation curve crosses the pure water ice VI melting curve at  GPa and  K. Going even deeper into the planet the SH methane hydrate dissociation curve is now set by the chemical potential equality of the SH methane hydrate and water ice VI. Unfortunately, no experimental data points are available for this regime of the methane hydrate dissociation curve. Beyond a pressure of  GPa the SH methane hydrate will transform to a filled ice-Ih structure (Loveday et al., 2001a) on which we elaborate below. Therefore, the range of uncertainty in the location of the SH hydrate dissociation curve spans a pressure difference of only about  GPa. Since the pressure range from  GPa to the introduction of the filled ice-Ih structure is relatively narrow, it will not make a major impact on our results if we assume for it an adiabatic profile of hydrates or of water ice VI.

We expect that an adiabatic profile with SH methane hydrate characteristics prevails in the pressure range from  GPa to  GPa because the filled ice-Ih structure is able to maintain within it more methane per water molecules than the SH methane hydrate. It is unlikely that along the adiabat somewhere between and  GPa the SH hydrate would dissociate to pure water ice VI and solid methane, only to incorporate methane with even a greater efficiency within the water structure due to a small increase in pressure of the order of  GPa. On that ground we tentatively assume that a direct transition from clathrate hydrate to filled-ice occurs not only at room temperature, where it is seen experimentally, but also up to  K (the adiabatic temperature at  GPa for our planets).

At a pressure of about  GPa the methane clathrate hydrate will transform into a filled ice-Ih structure. An informative depiction of the filled ice-Ih crystal structure may be found in Loveday et al. (2001b). In the filled ice-Ih structure the methane molecules occupy the widened channels of the filled ice lattice instead of the quasi-spherical cages they occupy in classic clathrate hydrate crystals. For more information on filled ice-Ih we refer the reader to paper I where we have estimated the filled ice-Ih thermodynamic stability field, a point to which we shall return after obtaining the water mantle thermal profile.

The introduction of the classical clathrate hydrate to filled ice phase transition raises the important issue of phase change induced partitioning of the convective cell. Such a partitioning may result in higher temperatures inside the planet, as the partitioning introduces conductive boundary layers between the various convective cells. We have already mentioned above the pronounced low thermal conductivity of methane hydrate. In appendix A we derive the clathrate hydrate to filled ice-Ih phase transition curve and show that the mantle convective cell is not likely to partition due to this phase transformation.

We therefore follow the adiabat in the SH methane hydrate layer until the transition to the filled ice phase, where we continue along an adiabat for the latter phase. Estimating the adiabat in the filled ice layer requires knowledge of its equation of state, for which we adopt the formalism derived in paper I. We also require the volume thermal expansivity for filled ice, which is experimentally unknown. As we explain in paper I, it is expected to be intermediate to the values for water ice VII and pure solid methane. By analogy with clathrate hydrates (see Fig. 2), which also represent a methane-water solid solution, we assign filled ice a volume thermal expansivity twice the value determined for water ice VII (Fei et al., 1993).

The heat capacity of the filled-ice mantle is taken to be a linear combination of the values for water ice VII and pure solid methane, weighted according to their abundances in the crystal, and respectively. The heat capacity for water ice VII is taken from Fei et al. (1993) and the heat capacity for a homogeneous system of methane is taken from Chase (1998). We point out that the data used here for the heat capacity of methane is from experiments on the gaseous phase. In other words we assume the entrapment of methane in the water ice lattice does not restrain the degrees of freedom of the methane molecules. Though this assumption is probably correct for low pressure (see Sloan, 1998, and references therein) it will gradually lose validity as the pressure increases, when venturing deeper into the planet. At high pressure the methane molecule may partly lose its ability to rotate freely and its heat capacity will decrease. Therefore, the adiabatic temperature across the ice mantle will rise [see eq.37]. However, we do not expect the resulting uncertainty in the heat capacity to have a large influence on our results. Rather we find the uncertainty associated with the probable range of values for the thermal expansivity of filled ice to dominate the overall uncertainty in the filled ice layer adiabat.

Following the adiabatic thermal profile, in the methane filled-water ice Ih layer deeper into the planet, we assume it terminates at a boundary layer which connects the water mantle with the silicate-metal interior. This boundary layer at the bottom of the filled ice mantle is henceforth referred to as the BBL. As for the case of the DBL, we also assume the BBL is on the verge of convective instability. Before formulating the BBL’s appropriate scaling Rayleigh number we first wish to estimate its kinematic viscosity ().

As we have already discussed above, a methane molecule will find it hard to diffuse through a clathrate structure water ring without causing local deformation of the surrounding water lattice. The filled-ice water lattice, which is far more compressed than the water lattice of cage clathrates, will impose even greater limitations on the ability of methane to diffuse. We therefore assume that diffusion creep is subdued in the methane filled water ice mantle. In conjunction with the high stress acting in the deep mantle, the viscosity associated with non-Newtonian mechanisms should be the dominant creep mechanism.

The non-Newtonian viscosities of high pressure water ice poly-morphs, such as ice VII and X, are not known experimentally. This is also the case for highly pressurized solid solutions such as filled ice. We can remedy this lack of data, to some modest extent, by making use of the algorithm given in subsection , where the melting curve is assumed to correctly describe the dependency of the viscosity activation enthalpy on pressure and temperature. The caveat here is that the factor [see eq.9] is assumed to be a constant which depends on the crystal structure, whereas in reality it is also a function of pressure. Therefore, while over limited pressure ranges may be assumed constant, over large pressure ranges extending over the entire ice mantle, should be allowed to vary with pressure. Since the physical basis for is not yet well formulated, the method of the homologous temperature is somewhat lacking in its ability to predict viscosity for cases where no experimental data exists. This said, we shall adopt the non-Newtonian viscosity for the highest pressure water ice poly-morph whose viscosity was determined, i.e. water ice VI, and try to adjust its activation energy and volume to better suit the characteristics of methane filled ice.

We assume that the viscosity of filled ice, which comprises the BBL, has characteristics analogous to those of methane clathrate hydrate. This assumption stems from the general point of view that both clathrates and filled ice are solutions to the basic problem of methane-water solid solubility and therefore probably share similar characteristics. In addition, the inclusion of methane in both crystal structures introduces voids and represents a similar impurity inserted into the water lattices, which tends to increase the viscosity (Durham et al., 1992).

It is interesting to note that while the activation energy for cage clathrates ( kJ mol) is lower than the activation energy of water ice VI ( kJ mol), the activation volume for cage clathrates ( cm mol) is higher than the activation volume of water ice VI ( cm mol) (Durham et al., 1997, 2003). A tentative explanation for these values would be that the introduction of methane into the water lattice introduces weaker methane-water bonds and some distortion of the water lattice, resulting in a reduced activation energy. In addition, the gliding of crystal planes one along the other requires breaking followed by re-connection of molecular bonds in the new location. As molecular bonds break the local molecules tend to expand resulting in an activation volume. Raghavendra & Arunan (2008) have shown that the non-bonded radius of the oxygen atom in pure water is  cm, while the bonded radius is  cm. This represents a volume change of  cm mol, remarkably close to the experimentally determined activation volume for pure water ice VI.

Raghavendra & Arunan (2008) further give the penetration distance of methane into water due to the formation of a weak hydrogen bond, and that can be compared with methane’s van der Waals radius (see paper I) to yield a volume difference of  cm mol. For the case of filled ice we adopt an activation volume which is weighted according to the abundance of the two different constituents:


Following O’Connell (1977) the activation volume is assumed to decrease with pressure as a lattice vacancy, with an effective bulk modulus () given by:


Where is the Poisson ratio estimated to be (Fu et al., 2010). and are methane filled ice bulk modulus and its pressure derivative, which experimentally are found to be  GPa and (see paper I), respectively.

For the activation energy of filled ice, also expected to be lower than that for water ice VI, we simply adopt the value from cage clathrates of  kJ mol. All the other parameters are adopted from water ice VI (Durham et al., 1997), giving for the filled ice kinematic viscosity the following form:


where the index refers the parameter to the filled ice structure. For we adopt the value of and for we assume  MPa s, both from ice VI viscosity measurements (Durham et al., 1997).

Now that we have an estimate for the viscosity at the BBL, we can analyze its stability with respect to convection. Here we need consider the complication of the viscosity being non-Newtonian. This effect was analyzed by Solomatov (1995), whose scaling for the non-Newtonian Rayleigh number, in combination with the viscosity formalism as expressed in Eq. (43), gives the following Rayleigh instability criterion for the BBL:


where the subscript refers the parameter to its value at the bottom boundary layer.

The temperature difference across the BBL is here represented by the temperature at the water/silicate boundary () and the temperature along the filled ice adiabat at the pressure prevailing in the outer boundary of the BBL, . The length scale of the BBL, , and the heat flux at the BBL, , are related through the relation:


Using the last relation to eliminate the temperature difference across the BBL, yields:


Further, assuming the BBL is narrow enough so that it may be represented using a constant density and acceleration of gravity [see Eq. 26], the last equation transforms to:


where is the pressure at the water/silicate boundary.

Assuming the water mantle has neither heat sources nor sinks we may scale the heat flux at the surface to the water/silicate boundary:


The acceleration of gravity at the BBL obeys:


Which, in combination with Eq. (47), yields a univariant equation for the pressure at the outer boundary of the BBL:


The thermal conductivity at the BBL, , is estimated in appendix C. When solving Eq. (50) we set the temperature to the average temperature at the BBL:

Figure 7: Bottom boundary layer, BBL, parameters for the planet for a varying ice mass fraction. The solid (blue) curve is the length scale of the BBL, in km. The dashed (red) curve is the temperature difference across the BBL in kelvins (normalized by ) and the dashed-dotted curve (green) is the estimated temperature on the transition to the silicate mantle in kelvins (normalized by ).

In Fig. 7 we show the variation of as a function of the ice mass fraction for the planet. It is interesting that the bottom boundary length scale has a maximum for an ice mass fraction between to . The situation is more complicated for the more massive planets. Quite generally, we find that the higher thermal expansivity of the filled ice mantle compared to a pure water ice mantle results in hotter planetary interiors, i.e. less steep adiabatic profiles.

Figure 8: Thick (red) curves represent the thermal profile in the clathrate hydrate and filled-ice mantle of the planet with a ice mass fraction (left panel) and a planet with and ice mass fractions (right panel). Thin (green) curves represent the melting curve for pure molecular water ice (see Lin et al., 2004; Goncharov et al., 2009). Thin (blue) curve represents the dissociation boundary for SI methane clathrate hydrate , SH methane clathrate hydrate and methane filled ice Ih, in order of increasing pressure, respectively. The dashed (red) curve is the adiabat for the case where filled ice is replaced with water ice VII.

In Fig. 8 we show the thermal profile (thick red curve) in the ice mantle for the planet, with ice mass fraction, against our estimated water-methane phase diagram (blue curve) and the melting curve of pure molecular water (green curve). The dashed (red) curve is the adiabat for the case where filled-ice is replaced with water ice VII. For the water ice VII mantle the interior temperatures are shown to be lower. In addition, as shown in the figure, the thermal profile in the mantle of the planet crosses the estimated stability field for methane filled ice (filled ice is not stable to the right of the blue curve). The extension of the mantle adiabat to the right of the filled ice stability boundary is therefore not appropriate. A proper extension of the adiabat beyond this point of intersection would require knowledge of the thermodynamic properties of the mixtures that exist beyond the stability of the methane filled ice. For the case of the planet the general behavior of the icy mantle thermal profile is similar to that of the planet. We find a point of intersection with the methane filled ice stability boundary to occur for both the and planets, regardless of whether the ice mass fraction is or .

For these more massive water planets, a question now arises, of what lies beyond the filled-ice stability regime. According to Benedetti et al. (1999), at high temperatures the C-H bond may break, resulting in the dissociation of the methane molecules. At high pressure condensation of the freed carbon atoms may ensue. It is of particular interest to compare our model, which we confine to the molecular solid regime, with the phase diagram for synthetic Uranus, derived experimentally by Chau et al. (2011). The introduction of carbon atoms to a water surrounding, at pressure above  GPa and temperature beyond  K, introduces a super-ionic phase whose extent of stability is narrower than the corresponding phase for a pure water system. This is because the introduction of carbon atoms increases diffusivity among the oxygen atoms, resulting in destruction of the super-ionic phase. At even higher temperatures (- K), depending on pressure, a reticulating phase is introduced, where, methane dissociates, releasing excess hydrogen. The relatively long lifetime of the C-C bond results in the formation of dense carbon clusters that should tend to segregate and sink (Chau et al., 2011).

It is possible that these high temperature and high pressure phases are present in the lower part of the icy mantle of our and planets, underlying an upper mantle composed of methane filled-ice acting as a thermal insulator to keep the interior warm. If that is indeed the case, then the creation of carbon clusters and their segregation will limit the ability of carbon expelled from the silicate interior from reaching the filled ice upper mantle, thus hindering its convection upward to the surface and the atmosphere.

It is tempting to generalize these arguments and say that for the more massive planets we describe, where the interior temperatures and pressures expected are higher, the lower part of the ice mantle may indeed be in the reticulating phase. For a less massive body the decrease in the expected interior temperature and pressure may result in a lower ice mantle in the super-ionic regime, leading to different consequences for carbon transport. For low mass planets (as we will show for our planets) the much lower interior temperatures and pressures could lead to an icy mantle which is entirely in the molecular crystal regime of filled-ice. For these low mass planets, the transport of methane expelled from the interior may be entirely due to its incorporation in the filled ice phase.

Due to the difficulties just mentioned from this point onward we will continue with emphasis on the planet alone. We have already derived the thickness of the BBL for the planet above. The temperature difference across the BBL, which is a conductive layer, obeys:


In Fig. 7 we give the temperature difference across the BBL and the expected temperature at the ice/silicate boundary for the planet for various ice mass fractions. As expected, as increases so does the temperature difference across the BBL. The temperature at the transition to the silicate mantle increases monotonically with increasing ice mass fraction, even though the thickness of the BBL decreases. This is because the total scale of the ice mantle increases with the ice mass fraction.

In the right hand panel of Fig. 8 we give the thermal profiles in the icy mantles of the and water mass fraction, planet. The water mass fraction scenario is to the left of the filled ice dissociation curve (blue), and therefore its entire ice mantle is probably composed of filled ice molecular solid. In the case of the and water mass fraction planet, a point of intersection with the filled-ice stability curve exists, as for the more massive planets. Though, contrary to the case of the more massive planets, the uncertainty in determining the exact location of the filled-ice dissociation curve is great enough so that we cannot rule out the possibility that its entire mantle is composed of filled-ice as well. After deriving probable thermal profiles in the icy envelopes of our water planets we wish to estimate the surface outgassing flux of methane into the atmosphere. Determining the outgassing mechanism is intimately linked to the geophysical behavior of the lithosphere, and therefore depends on the active tectonic mode. This issue is addressed in the following section.

5 Tectonics in Water Planets

It has been suggested that due to their larger masses super-Earths are even more likely than Earth to establish plate tectonics (Valencia et al., 2007). This stems from the hypotheses that the lithospheres of more massive super-Earths will be thinner and experience higher applied stresses, therefore having a greater ability to deform. This is a vital condition for plate tectonics. On the other hand, O’Neill & Lenardic (2007) argue that the increased fault strength due to scaling up of the planetary mass will make stagnant lid more probable. The reason for this apparently contradictory behavior was recently shown to stem from the fact that tectonic modes may have multiple solutions for the same parameter space. This was shown both analytically (Crowley & O’Connell, 2012) and numerically (Lenardic & Crowley, 2012).

The multiple solution nature of tectonics tells us that listing a planet’s parameters (e.g. mass, composition, viscosity, etc.) does not guarantee a unique tectonic mode. Rather, the geologic and climatic history must also be taken into account (Lenardic & Crowley, 2012). Certainly such a detailed history is not known for any water planet. Therefore, in this section we try to map the characteristics of the different multiple tectonic mode solutions, bearing in mind that all the derived modes are possible since we presently lack the knowledge required to rule out particular modes. If the different modes result in sufficiently different atmospheric regimes, it may be possible to distinguish among the different possibilities observationally. A first step towards this goal will be addressed in the next section.

Although the theory for multiple tectonic modes was originally developed for rocky planets, we assume that similar forces are responsible for maintaining plate motion and subduction in planets with ice layers. We take into consideration the fact that our icy mantle will have a rheological profile with depth which is very different from that assumed for the Earth’s mantle.

Above, we mentioned the existence of two layers whose thermal profiles are confined to the melting curves of SI methane clathrate and SH methane clathrate. Being confined to the melting curve implies these layers have relatively low viscosities. Indeed, thermal profiles depicted in Fig. 8 reveal a mid-layer whose thermal profile is fairly close to the local melting (dissociation) curve, suggesting that a low viscosity layer exists between the planetary lithosphere and lower mantle. This corresponds to the asthenosphere in the Earth. Considering the effect of clathrates on the planetary thermal profile, such a low viscosity mid-layer may be the rule rather then the exception in these icy worlds.

The ratio of the asthenospheric to lower mantle dynamic viscosities () in our case is somewhat difficult to constrain, since, as discussed above, the viscosities are probably non-Newtonian. Therefore, the viscosity ratio depends on the applied stress profile with depth, which, in turn, depends on the vertical and horizontal velocity profiles along with the temperature and pressure profile. Using boundary layer theory [see Eq .17], we find that the shear stress second invariant,  MPa for the planet. In conjunction with our estimated non-Newtonian viscosity model for filled ice [see Eq. 43] this gives an average lower mantle viscosity ranging from  Pa s to  Pa s for the planet, when varying the ice mass fraction between to respectively.

In Fig. 9 we show the actual viscosity profile with depth in the filled ice mantle for the planet assuming a ice mass fraction. Fig. 9 also serves as a test for the filled ice viscosity model by comparing it to viscosities of high pressure water ice polymorphs and silicates under the same thermal conditions. We expect the viscosity of silicates to be much higher than filled ice, and that of filled ice to be close to that of pure water ice. In Fig. 9 the viscosity profiles are obtained by keeping the mantle thermal profile the same while varying the viscosities between that for olivine, both wet and dry (Karato & Wu, 1993), our filled ice viscosity model and the viscosity for water ice VI (Durham et al., 1997). Indeed our model for the viscosity of filled ice yields a viscosity intermediate between that for ice VI and that for olivine, though it is much closer to ice VI than to olivine.

Figure 9: We compare the kinematic viscosities along the adiabat, for the planet, assuming ice mass fraction, for four different material viscosities, using the same planetary conditions. Thick solid (blue) curve is assuming a filled ice mantle. Dashed (red) curves confine the kinematic viscosity assuming dry olivine (Karato & Wu, 1993). Dashed-dotted (green) curves confine the kinematic viscosity assuming wet olivine (Karato & Wu, 1993). Thin (red) curve is the kinematic viscosity assuming water ice VI (Durham et al., 1997). The x-axis is normalized pressure, therefore spans the entire mantle, where is the SH clathrate hydrate to filled-ice transition and is the transition from ice to the silicate-metal core.

For the asthenospheric viscosity we use the form suggested for clathrate hydrates [see Eq. 11]. Although the viscosity of SH clathrate hydrate is unknown, one can estimate it by replacing the melting curve for clathrate SI with that for clathrate SH in Eq. 11, while keeping all the other parameters unchanged. Clearly these approximations for the viscosities of the lower mantle and asthenosphere are somewhat crude and more experimental data is required. Due to the possible large errors in the viscosities assigned to the different layers an exact solution is currently beyond reach. However, the viscosity formulations we use indicate, at least qualitatively, that an asthenosphere exists and that its viscosity may easily be at least three orders of magnitude smaller than the viscosity of the lower mantle. The ratio of asthenospheric to lower mantle viscosity in our water planets may therefore be significantly lower than what is assumed for Earth. As a result this asthenosphere will be of great significance for the ability to develop plate tectonics.

A plate tectonic theory that has the ability to account for the effect of an asthenosphere was recently developed by Crowley & O’Connell (2012). It has the ability to quantify both an active and a sluggish lid and to qualitatively describe the transition to a stagnant lid scenario. Adopting their model we solve for the dynamics of the lithospheric plate, except that we use the composition for icy worlds rather then rocky planets. Due to the large uncertainties in the viscosities of our icy crystal structures we solve for ratios of asthenospheric to lower mantle viscosity ranging to . This is done by keeping the lower mantle viscosity constant at the values mentioned above while varying the asthenospheric viscosity accordingly. The purpose of this exercise is to obtain a general qualitative understanding of the way plate tectonics may behave in water worlds for different viscosity ratios. Below we briefly summarize the theory. Further details can be found in Crowley & O’Connell (2012).

The theory is based on scaling considerations with the exception that it involves the derivation of the horizontal velocity profile in the convection cell. The latter may be used to derive the maximal vertical flow velocity in the cell, . Two mechanical energy balance equations are formulated, one for the lithosphere alone and the other for the entire convection cell. The energy balance equation for the lithosphere is:


Here and are the lithospheric plate thermal expansivity and isobaric heat capacity respectively. The gravitational acceleration is . The thermal thickness of the lithosphere is , estimated by Crowley & O’Connell (2012) to obey the half space cooling model. This gives , where is the lithospheric thermal diffusivity and is a characteristic time scale. If is the length of the plate and is its horizontal speed, then . is the advective rate of heat transfer through a horizontal cross section of the lithosphere averaged along the lithospheric depth. is the flow horizontal pressure gradient estimated at the base of the plate. is the shear stress operating on the base of the plate due to coupling with the underlying flow (upper part of the asthenosphere). The net resistive stress is , which obeys:


where is an effective bending stress representing a weighting factor for the action of plate bending at subduction and its ability to dissipate plate kinetic energy. is the fault stress from friction with the overlying plate during subduction. It too is a weighting factor for the dissipative efficiency of this mechanism. The final weighting factor is , the normal stress associated with slab pull, this term weighs the ability of the pulling slab to generate kinetic energy in the plate.

In Eq. 53 the different terms represent the work of buoyancy () which is a source of kinetic energy for the lithosphere, the work due to the horizontal pressure gradient () which is basically the flow pressure difference between the ridge and subduction zone, the work of traction from the underlying asthenosphere (), and the work of the net resistive forces (). Finally the energy balance equation for the entire convection cell is:


Here the subindex replaces the subindex meaning that the parameter is now a representative average for the entire cell rather then for the lithosphere alone. Furthermore, the lithospheric scale () is replaced with the depth scale of the entire cell, . is a term representing the dissipation in the lower mantle from both the horizontal flow at mid-cell and vertical flow at the convection cell corners. Crowley & O’Connell (2012) solve for both Eqs. (53) and (55) simultaneously to obtain and .

We also solve for and for the case of the super-Earth for ice mass fractions of and . After producing several rheological profiles for different asthenospheric stresses we find it to be a good approximation to partition the lower mantle and the asthenosphere at the clathrate hydrate to filled ice phase transition at  GPa. This results in lower mantle depths of  km and  km, for the and ice mass fraction respectively. The sinking slab, which is still attached to the lithosphere, will apply a stress due to its negative buoyancy. We estimate the normal stress due to this slab pull to be:


where is the bulk mass density of clathrate hydrates, is assumed to be  K and the slab length, , is assumed to be  km. This gives  MPa and  MPa for the and ice mass fraction respectively.

Estimating the effective bending stress () is somewhat more complicated. According to Crowley & O’Connell (2012) and Conrad & Hager (1999) treating the lithosphere as a beam that experiences bending at subduction, and thus dissipation, results in the following parametrization for its effective bending stress:


Here is the dynamic viscosity appropriate for the lithosphere and is the radius of curvature of the bent lithosphere. To obtain the term on the far right hand side of the last equation one simply has to replace the lithospheric length scale, with its half space cooling model estimate. For a given planetary mass and ice mass fraction the terms in the expression on the far right hand side may be considered constant. This results in effective bending stresses of  Pa and  Pa for the and ice mass fractions respectively. These seemingly low stress values should not be surprising. By themselves they are not physically significant, rather the physical significance is in the energy dissipation term due to plate bending at subduction whose scale is . The low values for simply mean that relatively little energy is dissipated due to plate bending at subduction. This is mainly due to the dynamic viscosity difference between silicates and water ice. For Earth  Pa s while for a water world  Pa s. This means that if in a water world the bending stress is about  Pa then for a rocky world it is  MPa. Comparing these results with shows that while in a rocky planet the dissipation term due to plate bending at subduction can almost counteract the effect of slab pull, which is a kinetic energy source for the plate, it can hardly do so in a frozen water world. If this were the whole story then plate tectonics could be said to be more likely in water planets in comparison to rocky planets.

One may question the proper choice for the radius of curvature, , for which we assign a value an order of magnitude less than the depth scale of the whole icy mantle, . This is a reasonable estimate for Earth (Crowley & O’Connell, 2012), and there is no reason to assume is much smaller in super-Earths. Actually numerical models suggest that the flow system will not allow to decrease too much as the system favors minimizing the dissipation due to plate bending (Capitanio et al., 2009).

Another complication that may increase is the fact that the above theory for assumes the lithosphere is a flat sheet, whereas in reality it is a 2D surface on a 3D sphere. Much like a flat slice of pizza is more easily bent at the tip than a folded slice, so too is the subducting lithosphere that has to fold into itself during down-welling (Mahadevan et al., 2010). This effect may be important in increasing the dissipation at subduction but due to its purely geometrical nature it will have the same effect on either a water or a rocky planet. Therefore, if this folding were to increase the dissipation due to plate bending substantially it would more readily stop plate tectonics on Earth, for which the flat sheet assumption gives . Thus this 2D on 3D folding effect is probably not large enough to change the fact that in water planets .

The dissipation in the lithosphere due to friction with the overriding plate is proportional to the fault zone stress, . The numerical value of the latter is not known. Due to the exploratory nature of this section we vary its value between zero and  MPa. These end values characterize two scenarios, one with a positive and the other with a negative , representing an accumulative tendency of the plate bending the fault zone friction and the slab pull to either restrain or encourage plate motion, respectively.

sign (cm yr) (MPa) (km) (cm yr) (erg cm s)