VOLATILE TRANSPORT INSIDE SUPER-EARTHS BY ENTRAPMENT IN THE WATER-ICE MATRIX

The discoveries and characterization of planetary systems orbiting other stars has entered an exciting period where we are starting to gain access to observing the atmospheres of planets that are essentially solid in nature—high-density rocky or icy planets of 1–10 Earth masses. These planets, collectively called super-Earths, have been discovered in relatively large quantities, though only a handful have measured radii and masses so far (Carter et al. 2012 and references therein). The mean densities derived for these exoplanets reveal a range of possible bulk compositions ranging from rocky with high iron content (e.g., Kepler-10b; Batalha et al. 2011) to mini-Neptunes with high H/He fraction (e.g., Kepler-11d,e; Lissauer et al. 2011). One of these planets orbiting a nearby M-dwarf star, GJ1214b (Charbonneau et al. 2009), has been accessible to spectroscopic studies of its atmosphere with inferences to its composition (e.g., Bean et al. 2011; Berta et al. 2012 and references therein).

At intermediate densities, many of the super-Earths should represent a new type of planet, unknown in our solar system, composed of a rocky core and a water envelope that exceeds 10% by mass. We call them water planets; this paper studies the transport of volatiles inside them. Possible examples of water planets could be Kepler-11b, 18b, or 20b.

Massive water planets were introduced conceptually by Kuchner (2003) and Léger et al. (2004). Fu et al. (2010) developed models of the interior dynamics of the water layers of water planets and concluded that materials released from the silicate-iron core should reach the surface, despite the high pressures at the core–mantle boundary and phase transitions. Materials would traverse the water mantle, composed of high-pressure ice phases, and reach the surface on timescales of 0.1–100 Myr. These are convective overturn timescales: what are the actual materials (e.g., gases) that could be transported to the surface and affect atmospheric observables needs to be determined by a detailed study. Here, we attempt a first step in this direction by considering entrapment of volatiles under the extremely high pressures inside water super-Earths.

The entrapment in water ice of volatile gases, which are mostly hydrophobic species, occurs in structures known as clathrate hydrates: the guest molecules occupy cages of water molecules stabilized predominantly by repulsive interactions between them. Typical clathrate hydrates include gases such as methane, carbon dioxide, oxygen, and nitrogen. Volatile transport by clathrates has been studied in the context of Titan and the icy satellites in our solar system (Tobie et al. 2006; Halevy & Stewart 2008; Schubert et al. 2010; Lunine et al. 2010; Sohl et al. 2010). Because clathrates are of practical terrestrial, as well as of astrophysical importance, much effort has been devoted to measuring their properties in the laboratory (see, e.g., Fortes & Choukroun 2010 and references therein). Most experiments are carried out at pressures up to an order of 100 MPa, which is high enough to provide useful information for modeling Titan-sized bodies.

Experiments at higher pressure show that methane clathrates undergo a transition to a more compact form called filled ice at pressures of around 2 GPa (Loveday et al. 2001b) that can survive at even higher pressures above 86 GPa (Hirai et al. 2006). For comparison, in a typical water super-Earth Fu et al. (2010) estimate that the pressures at the boundary between the ice mantle and the silicate core will be of the order of 100 GPa while the temperature will be of the order of 10³ K.

The basic theory of clathrates that was developed by van der Waals & Platteeuw (1959) was later applied by Lunine & Stevenson (1985) to situations of astrophysical interest. Below we suggest how this theory may be extended to higher pressures and temperatures. We use this theory to estimate the stability regime of filled ice, and discuss the implications of this for volatile transport in super-Earths.

Clathrates are crystals whose lattice structure forms cells that act as cages for foreign molecules. The empty cage-like structure is usually thermodynamically unstable, but the captured foreign molecules (i.e., guest molecules) help stabilize the clathrate crystal (van der Waals & Platteeuw 1959). Clathrate hydrate is essentially water ice. In this case, a framework of groups of four coordinated water molecules creates a cage-like lattice. The hydrogen-bonded water molecules are slightly distorted, however, from the tetrahedral angle of ordinary ice. There are two basic low-pressure (<1 GPa) forms of clathrate hydrates, referred to as structure I (SI) and structure II (SII). Each structure is composed of two types of cages, one small and other large. There is an additional hexagonal structure (SH) whose thermodynamic stability regime is roughly in the 1–2 GPa range (Hirai et al. 2001).

Filled water-ice structures also have the capability of trapping volatiles within their hydrogen-bonded network, except that instead of the clathrate cages the volatiles are entrapped in channels within the water ice (Loveday et al. 2001b). All cage clathrates have cages whose diameters are much larger than the channels which connect them, while the channels occupied by the guest molecules in filled ice are of constant and smaller diameter along their length (see Figure 1). Their densities also differ significantly: for methane clathrate SI, the water–methane ratio is 5.75:1 while for methane filled ice it is 2:1. The increased guest–host and guest–guest intermolecular interaction may help explain the increased stability of the filled ice against high-pressure compression and decomposition (Machida et al. 2007). Diamond anvil experiments, using methane as the guest molecule, suggest that such structures may survive pressures beyond 86 GPa (Hirai et al. 2006).

Zoom In Zoom Out Reset image size

The ability of clathrate structures to incorporate relatively large amounts of volatile species under conditions where the pure volatiles are highly unstable in the condensed phase is a primary characteristic that makes them so important to solar system and exoplanet studies. Enclathrating volatiles in icy planetesimals results in their inclusion in planetary icy mantles. As the icy mantle grows its inner layers become pressurized, and at pressures higher than about 2 GPa even SH clathrate hydrates are destabilized. Experiments show that the rate of pressure increase controls the resulting water phase: if the pressure is raised in steps of 0.2–0.5 GPa every 1–2 hr, the methane clathrate will decompose to ice VII and pure methane ice. However, if the sample of methane clathrate hydrate is kept at more than 3 GPa for 12 hr the molecules have enough time to rearrange themselves into a filled ice structure that can sustain very high pressures without decomposing (Loveday et al. 2001a). Thus, it appears that the filled ice phase will be an important component of a planetary ice mantle, especially inside water super-Earths.

In this paper, we begin by revisiting the thermodynamic stability theory of clathrates and extend it to higher pressures and temperatures. We consider only methane as the guest molecule here. In Section 3, we consider the case of methane filled ice. We conclude with numerical results and brief application to super-Earths.

The basic theory of clathrate hydrates views their stability in terms of two phases: the ordinary water phase (e.g., ice Ih, liquid water, etc.—the α phase) and the empty clathrate-hydrate structure (the β phase) + a gas of guest molecules. The chemical potentials of the two phases are equal on the thermodynamic equilibrium boundary between them. For clathrates, this may be written as (van der Waals & Platteeuw 1959)

$$\begin{equation} \mu _Q^\alpha = \mu _Q^\beta +kT\sum _i\nu _i\ln \left(1-\sum _Ky_{Ki}\right) \end{equation} \tag{ 1 }$$

$$\begin{equation} y_{Ki}=\frac{f_KC_{Ki}}{1+\sum _Jf_JC_{Ji}}. \end{equation} \tag{ 2 }$$

Here, the subindex Q refers to the molecule that makes up the host molecular network, which in our case is H₂O. The second term on the right-hand side of Equation (1) represents the contribution of the guest molecules to the clathrate-hydrate chemical potential. T is the temperature, k is Boltzmann's constant, ν_i is the ratio between the number of i type cages to water molecules per cubic unit crystal, and finally y_Ki is the probability that a guest molecule of type K occupies a clathrate cage of type i. This last function is given in Equation (2) in terms of the volatile fugacity (f_K) and its Langmuir constant (C_Ki). We will give further information on the Langmuir constant below and refer the reader to van der Waals & Platteeuw (1959) for a more detailed derivation.

Since most high-pressure experimental information is for methane-filled water ices, we restrict ourselves to the case of a single species of guest molecule model and omit the summation over K. The equilibrium equations may therefore be written as

$$\begin{equation} \frac{\mu _{{\rm H_2O}}^\alpha }{kT} = \frac{\mu _{{\rm H_2O}}^\beta }{kT}+\sum _i\nu _i\ln\, (1-y_{{\rm CH_4},i}), \end{equation} \tag{ 3 }$$

$$\begin{equation} y_{{\rm CH_4},i}=\frac{f_{{\rm CH_4}}C_{{\rm CH_4},i}}{1+f_{{\rm CH_4}}C_{{\rm CH_4},i}}. \end{equation} \tag{ 4 }$$

Recalling that both the fugacity and the Langmuir constant are functions of temperature and pressure, we differentiate and combine the last set of equations to obtain

$$\begin{eqnarray} &&\frac{\partial }{\partial T}\left(\frac{\mu _{{\rm H_2O}}^\alpha }{kT}-\frac{\mu _{{\rm H_2O}}^\beta }{kT}\right)_PdT+\frac{1}{kT}\frac{\partial }{\partial P}\left(\mu _{{\rm H_2O}}^\alpha -\mu _{{\rm H_2O}}^\beta \right)_TdP\nonumber\\ &&\quad=-\sum _i\nu _i\frac{1}{1+f_{{\rm CH_4}}C_{{\rm CH_4},i}}\left[\frac{\partial }{\partial P}(f_{{\rm CH_4}}C_{{\rm CH_4},i})_TdP\right.\nonumber\\ &&\qquad\left.+\, \frac{\partial }{\partial T}(f_{{\rm CH_4}}C_{{\rm CH_4},i})_PdT\right]. \end{eqnarray} \tag{ 5 }$$

We introduce the thermodynamic relations:

$$\begin{eqnarray} &&\frac{\partial }{\partial T}\left(\frac{\mu _{{\rm H_2O}}}{kT}\right)_P\equiv -\frac{H_{{\rm H_2O}}}{kT^2}, \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} &&\frac{1}{kT}\frac{\partial }{\partial P}(\mu _{{\rm H_2O}})_T\equiv \frac{V_{{\rm H_2O}}}{kT}, \end{eqnarray} \tag{ 7 }$$

where $H_{{\rm H_2O}}$ and $V_{{\rm H_2O}}$ represent the enthalpy and volume per water molecule, respectively. Inserting these thermodynamic relations into Equation (5) and solving for dP/dT gives

$$\begin{equation} \frac{dP}{dT}=\frac{\frac{1}{kT^2}\left(H_{{\rm H_2O}}^\alpha -H_{{\rm H_2O}}^\beta \right)-\sum _i\nu _i\frac{1}{1+fC_i}\left(\frac{\partial f}{\partial T}C_i+f\frac{\partial C_i}{\partial T}\right)}{\frac{1}{kT}\left(V_{{\rm H_2O}}^\alpha -V_{{\rm H_2O}}^\beta \right)+\sum _i\nu _i\frac{1}{1+fC_i}\left(\frac{\partial f}{\partial P}C_i+f\frac{\partial C_i}{\partial P}\right)}. \end{equation} \tag{ 8 }$$

Here, we have omitted the index CH₄ for convenience. For a homogeneous substance (i.e., composed only of water molecules), the chemical potential equals the Gibbs free energy per particle, $\tilde{G}$, which, combined with the assumption of constant temperature, gives

$$\begin{equation} d\mu =d\tilde{G}=VdP=kTd\ln P. \end{equation} \tag{ 9 }$$

Here, all extensive parameters are per particle. For the final equality on the right-hand side, we have used the equation of state for an ideal gas so that the last relation between pressure and Gibbs free energy is only valid for this case. When the gas is not ideal, we can retain the functional form if we replace the pressure with the fugacity, f, which acts as an effective pressure function to correct for the effect of intermolecular interactions and which obeys the following relation (Smith & Van Ness 1975):

$$\begin{equation} d\mu =d\tilde{G}=VdP\equiv kTd\ln f. \end{equation} \tag{ 10 }$$

From the last relation, we have

$$\begin{equation} \left(\frac{\partial f}{\partial P}\right)_T=\frac{Vf}{kT}, \end{equation} \tag{ 11 }$$

where V is the volume per methane molecule. For high pressures, we can solve Equations (1) and (2) for the case of an SI methane clathrate hydrate numerically. We find that fC_i ∼ 10³ ≫ 1. A similar result was found by Lunine & Stevenson (1985). Inserting this numerical result together with Equation (11) into Equation (8) yields

$$\begin{eqnarray} &&\frac{dP}{dT}=\frac{\frac{1}{kT^2}\left(H_{{\rm H_2O}}^\alpha -H_{{\rm H_2O}}^\beta \right)-\sum _i\nu _i\left(\frac{\partial }{\partial T}\ln f+\frac{\partial }{\partial T}\ln C_i\right)}{\frac{1}{kT}\left(V_{{\rm H_2O}}^\alpha -V_{{\rm H_2O}}^\beta \right)+\sum _i\nu _i\left(\frac{V}{kT} +\frac{\partial }{\partial P}\ln C_i\right)}.\nonumber\\ \end{eqnarray} \tag{ 12 }$$

The definition of the Langmuir constant according to van der Waals & Platteeuw (1959) is

$$\begin{equation} C_{Ki}\equiv \frac{h_{Ki}}{kT\zeta _{K}(T)}, \end{equation} \tag{ 13 }$$

where h_Ki is the single cell canonical partition function for a K type guest molecule in an i type cage. ζ_K is the quantum number density function for a K type molecule in an ideal gas and is independent of pressure. The cell partition function depends on pressure both via the cell dimension and through the form of the guest–host potential, so

$$\begin{equation} \frac{\partial }{\partial P}\left(\ln C_i\right)_T=\frac{\partial }{\partial P}\left(\ln h_i\right)_T. \end{equation} \tag{ 14 }$$

Here again we have omitted the index K. We are left with estimating the derivative of the single cell partition function, for which we give the following quasi-classical form:

$$\begin{eqnarray} &&\frac{\partial }{\partial P}\left(\ln h_i\right)_T\nonumber\\ &&\quad=\frac{\partial }{\partial P}\ln \left\lbrace \frac{1}{\hbar ^3}\int d^3 {\boldsymbol{r}}\int d^3{\boldsymbol{p}}e^{-1/kT[\epsilon _{{\rm rot},i}+\epsilon _{{\rm vib},i}+\epsilon _{{\rm trans},i}+W_i]} \right\rbrace. \nonumber\\ \end{eqnarray} \tag{ 15 }$$

Here, we have divided the Hamiltonian of the guest molecule into its separate kinetic and potential contributions. As a first-order approximation, we assume that the kinetic degrees of freedom of the entrapped molecule are unaffected by its inclusion in the water network. This is a common approximation in clathrate modeling, and experiment shows that in solidified form methane molecules rotate freely as in an ideal gas (Hazen et al. 1980). In this approximation, the kinetic degrees of freedom will contribute a function that is only a function of temperature, and the logarithm will cancel upon differentiation with respect to pressure, thus yielding a simplified form:

$$\begin{equation} \frac{\partial }{\partial P}\left(\ln h_i\right)_T\approx \frac{\partial }{\partial P}\ln \int e^{-W_i/kT}d^3\textbf {{r}}. \end{equation} \tag{ 16 }$$

In clathrate hydrates, the cages entrapping the volatiles are assumed spherical (e.g., Sloan 1998; McKoy & Sinanoǧlu 1963) so that the integration is over a spherical cage. In filled ices, the cages are actually cylindrical channels within the water-ice structure (Loveday et al. 2001b) so that we may write:

$$\begin{eqnarray} &&\frac{\partial }{\partial P}\ln \int e^{-W_i/kT}d^3\boldsymbol {r}\nonumber\\ &&\quad=\frac{\partial }{\partial P}\ln \int _{z1}^{z2}dz\int _0^{2\pi }d\phi \int _0^a e^{-W_i(r,\phi,z)/kT}rdr, \nonumber\\ \end{eqnarray} \tag{ 17 }$$

where we assume the cylindrical water-made channel has a radius a and since intermolecular potentials fall rapidly with increasing distance we limit the integration along the z-coordinate to the finite values z₁ and z₂.

Both the limits of integration and the potential interaction of the methane molecule with its surroundings, W_i, may depend on the pressure. Raman spectra of methane-filled water ice shows increases in the attractive potential between methane and its water network host with pressure increases (Machida et al. 2007). This was also suggested by Hirai et al. (2006). The rearrangement of the water network, from cage clathrate to the filled ice structure, reduces molecular distances by ∼0.5 × 10⁻⁸ cm. It is also suggested, from intermolecular distance considerations, that the guest–host Lennard-Jones potential interaction estimated for clathrates is enhanced by weak hydrogen bonds between the guest and the host in the filled ice structure (Loveday et al. 2001b). Since the guest–host potential energy changes considerably with pressure, a good first approximation will be to consider only the change of W_i with pressure, therefore allowing us to insert the derivative with respect to pressure into the integrand, giving

$$\begin{equation} \frac{\partial }{\partial P}\left(\ln C_i\right)_T=-\frac{1}{kT}\frac{\int _{z1}^{z2}\int _0^{2\pi }\int _0^are^{-W_i/kT}\left(\frac{\partial W_i}{\partial P}\right)_Td\phi dzdr}{\int _{z1}^{z2}\int _0^{2\pi }\int _0^are^{-W_i/kT}d\phi dzdr}. \end{equation} \tag{ 18 }$$

Finally, taking a spatial average for the partial derivative appearing in the numerator, we have

$$\begin{equation} \frac{\partial }{\partial P}\left(\ln C_i\right)_T=-\frac{1}{kT}\left\langle \left(\frac{\partial W_i}{\partial P}\right)_T\right\rangle. \end{equation} \tag{ 19 }$$

Inserting this last relation into Equation (12) yields

$$\begin{eqnarray} &&\frac{dP}{dT}=\frac{\frac{1}{T}\left(H_{{\rm H_2O}}^\beta -H_{{\rm H_2O}}^\alpha \right)+kT\sum _i\nu _i\left(\frac{\partial }{\partial T}\ln f+\frac{\partial }{\partial T}\ln C_i\right)}{V_{{\rm H_2O}}^\beta -V_{{\rm H_2O}}^\alpha +\sum _i\nu _i\left\langle \left(\frac{\partial W_i}{\partial P}\right)_T\right\rangle -\sum _i\nu _iV}.\nonumber\\ \end{eqnarray} \tag{ 20 }$$

Since the pressure varies by many orders of magnitude and the temperature does not, and since pressure has a dominant effect in determining the correction for the system non-ideality, we assume that

$$\begin{eqnarray*} &&\frac{\partial C_i}{\partial T}\gg \frac{\partial f}{\partial T}. \end{eqnarray*}$$

We may further approximate the following:

$$\begin{eqnarray} &&kT\sum _i\nu _i\frac{\partial }{\partial T}\ln C_i=kT\sum _i\nu _i\frac{\partial }{\partial T}\ln \left(\frac{1}{kT}\int d^3\boldsymbol {r}e^{-W_i/kT}\right)\nonumber\\ &&=kT\!\sum _i\nu _i\!\left[\!\frac{\int d^3\boldsymbol {r}\frac{W_i}{kT^2}e^{-W_i/kT}}{\int d^3\boldsymbol {r}e^{-W_i/kT}}-\frac{1}{T}\!\right]\!=\!\sum _i\nu _i\left[\!\frac{1}{T}\left\langle W_i\right\rangle -k\!\right], \nonumber\\ \end{eqnarray} \tag{ 21 }$$

where in the last step we have again averaged the intermolecular interaction of the methane molecule spatially over its surroundings. Inserting the last formula into the relation for dP/dT gives a Clausius–Clapeyron type equation of the form:

$$\begin{equation} \frac{dP}{dT}=\frac{\frac{1}{T}\left(H_{{\rm H_2O}}^\beta -H_{{\rm H_2O}}^\alpha \right)+\sum _i\nu _i\left[\frac{1}{T}\left\langle W_i\right\rangle -k\right]}{V_{{\rm H_2O}}^\beta -V_{{\rm H_2O}}^\alpha +\sum _i\nu _i\left\langle \left(\frac{\partial W_i}{\partial P}\right)_T\right\rangle -\sum _i\nu _iV}. \end{equation} \tag{ 22 }$$

This modified Clausius–Clapeyron type equation is a generalization of the equation given in Lunine & Stevenson (1985) in order to explain the behavior of clathrates. We will use this formalism to obtain the thermodynamic field of stability for the filled water ice as well. For this purpose, we need to evaluate the different terms appearing in this equation, which we do in the next section.

3.1. Clathrate Hydrate

Starting with clathrate hydrates, let us examine the numerator of Equation (22). At low pressures and temperatures, the α phase is water-ice Ih whose enthalpy is taken to be equal to the enthalpy of the empty hydrate (β) phase (Lunine & Stevenson 1985). The potential of interaction W_i is attractive (negative) and of a Lennard-Jones type, so generally the numerator is negative. In the denominator, at low pressures, the gaseous volatile volume (V) is the dominant factor appearing with a minus sign. Therefore, the derivative dP/dT is positive. The potential of interaction hardly changes with increasing pressure in this low-pressure regime.

As we increase the temperature, the pressure increases until we reach the melting point for water ice Ih and the α phase now represents liquid water. Due to the enthalpy of fusion, the enthalpy difference appearing in the numerator is no longer negligible and causes a sharp increase in the absolute value of the numerator. This is manifested as a sharp increase in the derivative dP/dT. Every increase in temperature is now accompanied by a steeper increase in pressure. At room temperature at 5 MPa, methane gas already deviates from ideality enough so that we need to consider a second virial correction; at 10 MPa a third virial correction is required, and so on (Hirschfelder et al. 1966). This means that the volume a methane molecule occupies in the gas decreases with pressure. At high enough pressure, the volatile gas contracts enough so that the empty clathrate volume equals the liquid water volume plus the compressed methane volume and the derivative dP/dT diverges.

Any further increase in pressure brings about a situation where the volume occupied by a clathrate water molecule is larger than the sum of the volumes occupied by a water molecule in the liquid phase and a methane molecule in the gas phase weighted by the hydration number (ν). The result is that the derivative dP/dT becomes negative and the high-pressure stability limit of the clathrate hydrate is attained and the clathrate dissociates. This general type of clathrate behavior is shown in Figure 2, which we have derived by numerically solving the set of Equations (1) and (2).

Zoom In Zoom Out Reset image size

Now let us consider what happens for the case of filled water ice. This ice is formed in the laboratory at room temperature and at a pressure of ∼2 Gpa. The structure is an ice Ih water network which is distorted in such a way that the interconnecting channels widen to accommodate the methane molecules within the hydrogen-bonded network. This distortion from the usual structure of water ice Ih is manifested by a change in the tetrahedral angles (Loveday et al. 2001b). The β phase will now describe the filled ice structure while the α phase represents either water ice VII or fluid water, as appropriate.

To obtain the thermodynamic stability field for filled water ice, we shall require the temperature and pressure dependencies for both the different constituents' volumes and for the attractive potential between a methane molecule and its surroundings. In addition, the enthalpy difference between filled ice and fluid water ought be estimated for the case of stability with respect to fluid water.

3.2. High-pressure Equations of State

We start with the general relation:

$$\begin{equation} \frac{dV}{V}=\chi dT-\kappa dP, \end{equation} \tag{ 23 }$$

where χ and κ are the volumetric thermal expansivity and isothermal compressibility, respectively. We assume the bulk modulus, B, is linearly dependent on pressure and may be written as

$$\begin{equation} B\equiv \frac{1}{\kappa }=B_0+\tilde{B_0}P. \end{equation} \tag{ 24 }$$

Combining Equations (23) and (24) yields after integration:

$$\begin{eqnarray} V(T,P)&=&V(T_0,P_0)\left(\frac{B_0+\tilde{B_0}P}{B_0+\tilde{B_0}P_0}\right)^{-1/\tilde{B_0}}\nonumber\\ &&\times\exp \left(\int _{T_0}^T\chi (T,P) dT\right). \end{eqnarray} \tag{ 25 }$$

If we keep the temperature constant at T₀, we end up with the Birch–Murnaghan equation of state. By setting the reference temperature (T₀) to room temperature, we can assign to each solid constituent a proper value for its bulk modulus (B₀) and its pressure gradient ($\tilde{B_0}$) by using published room temperature hydrostatic compression experiments. We use the data published in Hemley et al. (1987) for water ice up to 128 GPa for assigning a bulk modulus for ice VII. For assigning a bulk modulus for solid methane, we use the data published in Hazen et al. (1980) for phase I and in Umemoto et al. (2002) for phases A and B up to 37 GPa. Solid methane transforms at high pressures to a hexagonal close packed structure (Bini & Pratesi 1997). This phase was difficult to account for as we found no volumetric data for it, rather we extrapolated from the phase B data in Umemoto et al. (2002) to higher pressures. For the filled ice structure, we used room temperature, experimentally deduced, unit cell volumes up to 42 GPa by Hirai et al. (2003).

For the thermal expansivity, we have adopted the approach of Fei et al. (1993), who determined the thermal expansivity of water ice VII by fitting their volumetric experimental data to an equation of state of the form given above (Equation (25)). These authors assumed the following dependence for χ on the pressure:

$$\begin{equation} \chi (T,P)=\chi (T,P_0)\left(\frac{B_0+\tilde{B_0}P}{B_0+\tilde{B_0}P_0}\right)^{-\eta }, \end{equation} \tag{ 26 }$$

where χ(T, P₀) is taken to be a linear function of the temperature and η is well fitted with a numerical value of 0.9. For the thermal expansivity of solid methane, we adopt the formalism of Equation (26) for its dependency on pressure and set the expansivity value at P₀ to be 10⁻³ K⁻¹, according to experimental data given by Heberlein & Adams (1970). The thermal expansivity of filled ice is not known. As filled ice is a combination of a hydrogen-bonded network and methane molecules between which there is van der Waals type attraction, we suggest its thermal expansivity to be intermediate between that for ice VII and for solid methane.

For the equation of state of fluid water at high pressure (up to 100 GPa), we adopt the formalism derived using molecular dynamic simulations by Belonoshko & Saxena (1991). Although this formalism is inherently dependent on the type of model used for the water molecules' intermolecular potential it does show good agreement with recent experimental data for high-pressure water fluid density (Goncharov et al. 2009).

We shall now turn to evaluate the energy of interaction of a methane molecule with its surroundings in the water filled ice structure.

3.3. The Interaction Potential W_i

In order to build the thermodynamic stability field for methane-filled water ice, we need to approximate how the potential of interaction, of the methane molecule with its surroundings, depends on the temperature and the pressure. As was shown by Raghavendra & Arunan (2008), a hydrogen bond may form between the water electron poor hydrogen and the center of the methane tetrahedral face which is electron-rich. It was further shown by these authors that the bond energy (E_bond) for such an interaction is −6.3 kJ mol⁻¹. For comparison, a methane–methane van der Waals potential well is about −1.2 kJ mol⁻¹ (Hirschfelder et al. 1966). It was also suggested by Loveday et al. (2001b) that weak hydrogen bonds are formed in the filled water-ice structure between the water network and the dissolved methane molecules.

Let us consider a simple approximate model where, with increasing pressure, more of the tetrahedral faces per methane molecule create such hydrogen bonds with the water network. An increase in temperature will have the opposite effect. That the number of hydrogen bonds per molecule increases with pressure and decreases with temperature is known for water structures (e.g., Kalinichev & Bass 1994; Pattanayak et al. 2011). We may therefore write for the spatially averaged potential of the interaction of a methane molecule with its surroundings:

$$\begin{equation} \left\langle W_i\right\rangle = n(T,P)E_{{\rm bond}}, \end{equation} \tag{ 27 }$$

where n(T, P) is the number of hydrogen bonds between a methane molecule and the water network at a given pressure and temperature, and is bounded between 0 and 4, where in the upper limit all of the methane four tetrahedral faces are hydrogen bonded to the water network.

We normalize n to give it a probability interpretation (i.e., what is the probability a bond will form at a given pressure–temperature point). We further assume a division into a temperature-dependent and pressure-dependent probabilities, thus

$$\begin{equation} n(T,P)\equiv n_1(T)n_2(P). \end{equation} \tag{ 28 }$$

For the temperature-dependent probability, we assume a Boltzmann type form of

$$\begin{equation} n_1(T)=\frac{1}{1+\exp (-4|E_{{\rm bond}}|/kT)}. \end{equation} \tag{ 29 }$$

To obtain the form for the pressure-dependent probability (n₂), we use the fact that a formation of a hydrogen bond is accompanied by a substantial penetration of the hydrogen bonding molecules into each other. For the case of the hydrogen bonding between water and methane, the combined van der Waals radius of 2.90 × 10⁻⁸ cm, before the bonding, reduces to 2.47 × 10⁻⁸ cm, after bonding (Raghavendra & Arunan 2008). This interpenetration may therefore account for a 40% reduction in the crystal volume. This is considerable enough so that we may relate this interpenetration (to first approximation) to the solid compressibility, κ. At low pressures, only a few hydrogen bonds are formed and the solid is easily compressed, as many bonds are still ready to be formed. As the pressure increases more bonds form per molecule, and it becomes more difficult to further compress the solid. Hypothetically, at a high enough pressure, all bonds per molecule are already formed and it is no longer possible to further compress the solid via the route of hydrogen bonding molecular interpenetration. Adopting this model to first approximation, we can write for the molecular volume, V, the following:

$$\begin{equation} V=n_{{\rm HB}}V_{{\rm HB}}+n_{{\rm nHB}}V_{{\rm nHB}}, \end{equation} \tag{ 30 }$$

where n_HB and n_nHB are the probabilities whether a methane molecule is fully hydrogen bonded to the water network or not hydrogen bonded at all, respectively, and V_HB and V_nHB are the molecular volumes associated with these two molecular situations, respectively. If the number of hydrogen bonds indeed determines, to first approximation, the crystal volume, we may say

$$\begin{eqnarray} \kappa &\equiv& {-}\frac{1}{V}\left(\frac{\partial V}{\partial P}\right)_T=\frac{-1}{n_{{\rm HB}}V_{{\rm HB}}+n_{{\rm nHB}}V_{{\rm nHB}}}\nonumber\\ &&\times\left(V_{{\rm HB}}\frac{dn_{{\rm HB}}}{dP}+V_{{\rm nHB}}\frac{dn_{{\rm nHB}}}{dP}\right). \end{eqnarray} \tag{ 31 }$$

Since n_HB + n_nHB = 1, an integration of Equation (31) yields

$$\begin{equation} n_2(P)=n_{{\rm HB}}(P)=\frac{1-\exp \left(-\int _0^P\kappa dP\right)}{1-\frac{V_{{\rm HB}}}{V_{{\rm nHB}}}}. \end{equation} \tag{ 32 }$$

Inserting Equations (29) and (32) into Equation (27) gives, for the spatially averaged potential energy of a methane molecule with its surroundings, the following form:

$$\begin{equation} \left\langle W_i\right\rangle \approx \frac{4}{1+e^{-4|E_{{\rm bond}}|/kT}}(1-e^{-\kappa _0 P})E_{{\rm bond}}. \end{equation} \tag{ 33 }$$

Its derivative with respect to pressure will therefore be

$$\begin{equation} \left\langle \left(\frac{\partial W_i}{\partial P}\right)_T\right\rangle \approx \frac{4\kappa _0 E_{{\rm bond}}}{1+e^{-4|E_{{\rm bond}}|/kT}}e^{-\kappa _0 P}. \end{equation} \tag{ 34 }$$

Now that we have the approximate temperature and pressure dependencies for the terms appearing in Equation (22), we may integrate it numerically to obtain the stability regime for the filled ice structure.

3.4. Numerical Results and Discussion

The numerical integration of Equation (22) is shown in Figure 3. As we have mentioned above, we consider it reasonable to assume the thermal expansivity of the filled ice to be intermediate, between that for solid methane (a van der Waals solid) and water ice VII. From the laboratory data (see Section 3.2), we know that the thermal expansivity for solid methane is an order of magnitude larger than that for ice VII. We integrate Equation (22) assuming thermal expansivity for the filled ice two, three, and five times larger than the experimental value for ice VII. As a limiting case, we also solve assuming filled ice has a thermal expansivity equal to that of water ice VII. The four curves are given in the figure. It is appropriate to note here that Sloan (1998) gives for an SI clathrate hydrate (also a combination of methane and water) a thermal expansivity some five times larger than the thermal expansivity for ice Ih. More recent experiments confine the latter ratio to be between two and four (Hester et al. 2007).

Zoom In Zoom Out Reset image size

In addition, the parameter ν, the hydration number, is important to Equation (22), as defined earlier. Generally speaking, this parameter represents the ability to include volatiles in the water network. As the ratio between water and methane in the filled ice is 2:1 (Loveday et al. 2001b), then filled ice modeling requires a value of 1/2 for ν. Given that, we find the packing efficiency for filled ice to be greater than for the case of separation to pure solid constituents (i.e., water ice VII and pure solid methane). The contribution of the potential of interaction (in comparison to packing efficiency considerations) is found to be small and restricted to the lower pressure regime of the stability field.

Integrating Equation (22) and following the dissociation curve for filled ice with respect to ice VII, a point of intersection with the melting curve for water ice VII is reached. Such an intersection is commonly referred to as a quadruple point. Up to the quadruple point, the enthalpy difference between the α and β phases is relatively small, as both phases are solid. Continuing the integration beyond the quadruple point, the α phase will now represent fluid water. The enthalpy difference between the α and β phases therefore increases and becomes dominated by the enthalpy of fusion of filled ice. Unfortunately, the enthalpy of fusion of methane filled ice Ih is experimentally undetermined and we estimate its value using the enthalpy of fusion for pure solid water at high pressure.

Goncharov et al. (2009) experimentally deduced the melting curves for water ice VII and for super-ionic water. By further experimentally determining the volume difference between solid and melting along the melting curve, they were able to derive the enthalpy of fusion for pure water. Their reported error is approximately 50%. As explained in Goncharov et al. (2009), the enthalpy of fusion increases with pressure since the increase in pressure means the melting transition is between a molecular solid and a fluid whose molecules become ever more dissociated and ionized. At a pressure of about 47 GPa, a branching occurs in the melting curve for pure water due to the introduction of the super-ionic phase of water (Goncharov et al. 2005). The result is a reduction in the melting curve gradient and a sharp decrease of the enthalpy change upon melting (Goncharov et al. 2009). This behavior is clearly seen in Figure 4, where the dashed curve (red) reproduces the results from Goncharov et al. (2009) and the dash-dotted curve (blue) is a hypothetical enthalpy of fusion, for which we have arbitrarily reduced the point of branching in the melting curve to 27 GPa. The reason for the introduction of this hypothetical behavior is that the enthalpy of fusion we adopt is from experiments on a homogeneous water system and therefore it is only an approximation for our methane filled ice system. Molecular dynamic simulations by Iitaka & Ebisuzaki (2003) demonstrate that filled ice exhibits behaviors similar to those appearing in homogenous water systems, but at somewhat lower pressures. A branching of the melting curve is therefore probable in our system of methane mixed with water as well, though the branching may happen at a lower pressure than the 47 GPa of the pure water system. To test to what extent our results depend on the point of branching, we solve for the hypothetical enthalpy of fusion as well.

Zoom In Zoom Out Reset image size

In Figure 5, we solve Equation (22) where we assume for filled ice a thermal expansivity twice that for water ice VII and solve for each of the enthalpy of fusion scenarios depicted in Figure 4. In the same figure, we also show the solution of Equation (22) with the experimental enthalpy of fusion from Goncharov et al. (2009) but vary it globally by 50%, which is its experimental error. From the figure, we see that most of the variations in the filled ice Ih dissociation curve due to the changes examined in the enthalpy of fusion occur below 10 GPa.

Zoom In Zoom Out Reset image size

Combined with thermodynamic profiles for the interior of water exo-Earths, the thermodynamic stability field can help us decide whether filled ice structures can form at the core–mantle boundary and help convect methane and other hydrocarbons toward the surface.

For the interior structure of water planets with masses up to 10 M_⊕, Fu et al. (2010) consider a silicate-metal core, surrounded by a water-ice mantle. They find solid-state convection to prevail between the core–ice boundary and the conductive crust. Examining different surface temperatures and heat fluxes, the authors find it possible for a liquid water ocean to exist beneath an ice–Ih solid surface layer. Different compounds expelled from the silicate-metal core may be incorporated in the ice matrix and convect outward. For the 2 M_⊕ planet, Fu et al. (2010) find the silicate core–ice boundary pressure to range from approximately 50 to 90 GPa for water mass fractions between 25% and 50%, respectively. The expected temperatures are between 700 and 1000 K.

Our extension of the theory of van der Waals & Platteeuw (1959) suggests that the filled ice structure may be stable at such pressures and temperatures. In such a case, any CH₄ released from the core could be trapped in the lower part of the ice mantle. Experimentally, Loveday et al. (2001b) inferred a 2:1 water–methane ratio for the filled ice phase. As the mantle convects, the CH₄ would be carried to lower pressures where the filled ice will undergo a transition to a clathrate hydrate. Using the statistical model given by van der Waals & Platteeuw (1959), we can compute the probability, y_Ki that an SI clathrate-hydrate cage will be occupied. This is shown in Figure 6. It is clear that the occupancy is very close to full occupancy (probability of unity) for which there are 5.75 water molecules per methane molecule. For a structure H methane hydrate the water–methane ratio ranges from 4.25:1 to 3.40:1, depending on the number of methane molecules occupying its large cage (Koh 2002). Thus as the pressure decreases in an upwelling, excess CH₄ will be forced out of the water network. This may lead to the formation of a local methane reservoir. Such reservoirs may naturally occur on the transition from filled ice to structure H clathrate hydrate at about 2 GPa, and at the transition from structure H clathrate hydrate to SI clathrate hydrate at about 1 GPa.

Zoom In Zoom Out Reset image size

In A. Levi et al. (2013, in preparation), we show that the introduction of filled ice as a major constituent in an icy mantle has a large effect on the mantle thermodynamic profile. The probable higher thermal expansivity of filled ice compared with that for water ice VII results in a more moderate adiabatic thermal profile. The higher temperatures in the icy mantle, compared with those for a pure water mantle, create a physical route through which super-Earths (objects less massive than Uranus and Neptune) may develop lower mantles in the super-ionic and reticulating phases. Phases so far related with the interior of bodies whose mass is equivalent to that of the icy giants of our solar system.

The unique characteristics of methane clathrate hydrate, namely its low thermal conductivity and the topology of its melting curve, yield water planets with thin crusts (<1 km) and a tendency to form a layer closely confined to the local melting condition beneath it. In that respect, the geology and surface–atmosphere coupling in water super-Earths are quite different than simply assuming a water planet is a scaled-up version of an icy satellite, such as Titan. Icy satellites tend to form thick crusts, on the order of 100 km, and stabilize in a stagnant lid regime. This also means that chemical and geocycles in super-Earths are not simply resolved using analogies on icy moons.

The formation of a layer, in which conditions are close to melting, underneath the crust acts as a low viscosity layer that may decouple the convection cell and the crust, and represent a channel for relatively fast horizontal flow from upwelling to downwelling (Crowley & O'Connell 2012). Such a fast horizontal flow may preferentially drag methane reservoirs to areas of downwelling where it may again become incorporated into the water matrix, preventing it from reaching the surface and the atmosphere. The timescale for transport of methane to the atmosphere is therefore broken down to several stages: the filled ice mantle overturn, the rate of transport in the low viscosity layer underlying the crust, and the rate of collapse of the thin crust followed by exposure of fresh material.

Other substances trapped in filled ice structures may be dragged preferentially along with the convecting ice. If they reach the top of the convection cell and meet with the bottom of a subsurface ocean, some compounds may favor clathrate-hydrate formation. This will deplete those substances relative to others that are more easily transported in the liquid water and released to the atmosphere via cryovolcanism and surface ice tectonics. The volatile composition and abundance may also affect the ocean freezing temperature. This will change the thickness and temperature of the solid crust. These issues are developed in more detail in A. Levi et al. (2013, in preparation).

In this work, we have estimated the thermodynamic stability field for methane filled ice. From this field of stability, it is clear that the blanketing effect of a thick H/He atmosphere, resulting in high internal temperatures, will most likely destabilize any filled ice (see Figure 3 for the P–T conditions at the base of the H/He atmosphere surrounding Uranus and Neptune). Water super-Earths lacking a substantial H/He atmosphere are more favorable for the formation of filled ice, which in turn will have a large influence on them. Additional experimental and theoretical work is needed to further constrain the properties of methane filled ice in order to better understand this influence.

We wish to thank the anonymous referee for constructive comments.