© ESO, 2010

1. Introduction

For several years, Jupiter-family comets (hereafter JFCs) have been the targets of spacecraft missions whose aims are to determine the comets’ composition, structure, and physical properties. The Stardust mission has collected dust originating from Comet 81P/Wild 2 in January 2004 (Tsou et al. 2004), and the Rosetta mission is currently flying towards Comet 67P/Churyumov Gerasimenko for a rendezvous in August 2014 (Schulz et al. 2004). Moreover, Comet 46P/Wirtanen, which was the initial target of the Rosetta mission (Schulz & Schwehm 1996), is again considered for a rendezvous in 2021 with the PriME (Primitive Material Explorer) mission, which is currently proposed to NASA as a Discovery mission (Cochran et al. 2010).

Marboeuf et al. (2010) recently proposed that CO and CO₂ clathrate hydrates (hereafter clathrates) could exist within short period comets¹ (hereafter SPCs). Clathrates are non-stoichiometric crystalline compounds, with hydrogen-bonded H₂O molecules forming cages in which gas molecules are trapped individually by van der Waals interactions. At a given temperature, the stability of these cages is ensured as long as the gas phase pressure equals or exceeds the dissociation pressure of the considered clathrate. Using a thermodynamic model of JFCs including a mixture of H₂O, CO, and CO₂ ices, Marboeuf et al. (2010) showed that the stability zone of clathrates in these comets is restricted to a small part of the nucleus. They also found that the stability zone of clathrates could almost disappear at perihelion in Halley-type comets.

Using a mixture of gases, Mousis & Schmitt (2008) showed that the resulting clathrate is also a mixture of volatiles whose equilibrium pressure is somewhere between that of the individual gases. In the present work, we consider a larger and more representative set of volatiles (CO, CO₂, CH₄, and H₂S) that might exist in the nuclei compared to the one (CO and CO₂) initially adopted by Marboeuf et al. (2010). We then investigate the formation conditions of clathrates within the two JFCs 67P/Churyumov-Gerasimenko (hereafter 67PCG) and 46P/Wirtanen (hereafter 46PW) and show that irrespective of the latitude of the considered area, a permanent layer of multiple guest clathrates is likely to form in the subsurfaces of the two comets, and in a broader manner, in the subsurfaces of all SPCs. This suggests new directions of investigation for comet nucleus modelling which will be addressed in a forthcoming paper.

2. Nucleus model

The nucleus model employed in this work is based on the one-dimensional model developed by Marboeuf (2008) and Marboeuf et al. (2009). The model considers an initially homogeneous sphere composed of a predefined porous mixture of ices and dust in specified proportions. It describes heat transmission, gas diffusion, sublimation/recondensation of volatiles within the nucleus, water ice phase transition, dust release, and mantle formation. Water ice can be initially amorphous or crystalline, depending on the formation location of the body in the primitive nebula. When the ice is initially amorphous, a fraction of the volatiles is trapped in the water ice matrix and the remaining fraction forms pure condensates in the pores. When heated, the fraction of volatiles condensed in the pores sublimates first, and then the other fraction trapped within the matrix is released during the transition from amorphous to crystalline water ice. The released gas diffuses throughout the porous matrix by tortuous capillary tubes and eventually escapes from the object. The crystallization occurs on a temperature-dependant timescale $\mathit{\tau }\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{=}\mathrm{9.54}\mathrm{\times }{\mathrm{10}}^{-14}{\mathrm{e}}^{\frac{\mathrm{5370}}{\mathit{T}}}$, where T is the temperature (K) (Schmitt et al. 1989). Hence crystallization becomes efficient and effectively occurs at temperatures ~110 K.

These physical processes are described by the equation of diffusion of heat and mass. The heat diffusion throughout the nucleus is described by the energy conservation equation:

$\mathit{\rho c}\frac{\mathit{\partial T}}{\mathit{\partial t}}\mathrm{=}\mathrm{\nabla }\mathit{.}\left(\mathit{K}\frac{\mathit{\partial T}}{\mathit{\partial r}}\right)\mathrm{-}\sum _{\mathit{x}}{\mathit{H}}_{\mathit{x}}\mathrm{(}{\mathit{Q}}_{\mathit{x}}\mathrm{+}{\mathit{Q}}_{\mathit{x}}^{\mathrm{cr}}\mathrm{)}\mathrm{+}{\mathit{Y}}_{\mathrm{cr}} \mathrm{(}{\mathrm{J}\hspace{0.17em}\mathrm{m}}^{-3}\hspace{0.17em}{\mathrm{s}}^{-1}\mathrm{)}$(1)

with $\mathit{\rho c}\mathrm{=}{\sum }_{\mathit{l}}{\mathit{\rho }}^{\mathit{l}}{\mathit{c}}^{\mathit{l}}\mathrm{=}{\sum }_{\mathit{x}}{\mathit{\rho }}_{\mathit{x}}^{\mathit{i}}{\mathit{c}}_{\mathit{x}}^{\mathit{i}}\mathrm{+}{\mathit{\rho }}^{\mathit{d}}{\mathit{c}}^{\mathit{d}} \mathrm{(}{\mathrm{J}\hspace{0.17em}\mathrm{m}}^{-3}\hspace{0.17em}{\mathrm{K}}^{-1}\mathrm{)}$. Here t is the time (s), K the heat conduction coefficient (J s^-1 m^-1 K^-1) of the porous matrix, whose formula is based on geometric assumptions (Espinasse et al. 1993; Orosei et al. 1999), r the distance (m) from the centre of the nucleus, ρ^l the density (kg m^-3) and c^l the specific heat capacity (J kg^-1 K^-1) of the solid component l in the comet nucleus (l =  d for dust or i for pure ices of elements x). H_x is the molar latent heat of sublimation of ice x (J mol^-1) and ${\mathit{Q}}_{\mathit{x}}^{\mathrm{cr}}$ (mol m^-3 s^-1) and Y_cr are the rate of moles of gas x and the power per unit volume released during the crystallization process of amorphous water ice respectively (Espinasse et al. 1991; Orosei et al. 1999). Q_x is the rate of volatile molecule x (mol m^-3 s^-1) that sublimates/condenses in the pores of the matrix, and is given by the gas diffusion equation that describes the diffusion of gas through the porous matrix for each molecule x:

$\frac{\mathit{\partial }{\mathit{\rho }}_{\mathit{x}}^{\mathrm{g}}}{\mathit{\partial t}}\mathrm{=}\mathrm{\nabla }\mathit{.}\left({\mathit{G}}_{\mathit{x}}\frac{\mathit{\partial }{\mathit{P}}_{\mathit{x}}}{\mathit{\partial r}}\right)\mathrm{+}{\mathit{Q}}_{\mathit{x}}^{\mathrm{cr}}\mathrm{+}{\mathit{Q}}_{\mathit{x}} \mathrm{(}{\mathrm{mol}\hspace{0.17em}\mathrm{m}}^{-3}\hspace{0.17em}{\mathrm{s}}^{-1}\mathrm{)}\mathit{,}$(2)

where ${\mathit{\rho }}_{\mathit{x}}^{\mathrm{g}}$ is the molar density of gas x (mol m^-3), P_x its partial pressure (Pa) and G_x its diffusion coefficient (mol m^-1 s^-1 Pa^-1). The flux of gas diffusing through the porous matrix can be either a free molecular flow or a viscous flow, because the mechanism governing the diffusion of volatiles is selected according to the mean free path of molecules (Prialnik et al. 2004).

At the surface, the local temperature is given by a thermal balance between the adsorbed solar energy, the thermal emission, the heat diffusion towards the interior and the energy of sublimation of the existing ices (Orosei et al. 1999). Note that at each time step the porosity and the radius of pores are recomputed for each layer, taking into account the density variations of the solid phase.

3. Thermodynamic parameters and initial composition

At the beginning of the computation, the initial water ice of our nucleus is assumed to be amorphous (Taylor 1992; Kouchi et al. 1994). The JFCs come probably from the Kuiper Belt (Lowry et al. 2008), which is populated by objects likely made from amorphous ice. Kouchi et al. (1994) showed that temperatures lower than ~110 K, corresponding to heliocentric distances greater than ~12 AU in the solar nebula would preserve the amorphous ice initially formed in the ISM. Below ~12 AU, the disk temperature exceeded ~110 K, resulting in the crystallization of the amorphous grains embedded in that part of the nebula. In addition to water ice, the model takes into account the solid phases of CO, CO₂, CH₄, and H₂S. These molecules are among the most abundant volatile species (production rates relative to water greater than 1%) observed in cometary nuclei (Bockelée-Morvan et al. 2004) and are susceptible to form clathrates at low-pressure conditions (Lunine & Stevenson 1985). Other molecules such as H₂CO, CH₃OH, and NH₃, which are also abundant in comets, have not been considered in this study because there are no experimental data concerning the equilibrium of H₂CO and CH₃OH clathrates (Fray & Schmitt 2009; Marboeuf et al. 2008; Mousis et al. 2009), while NH₃ does not form clathrates (Lunine & Stevenson 1985). Table 1 gives the X/H₂O (J_X) mole fractions relative to water (with X = CO, CO₂, CH₄ or H₂S) of the volatiles initially present in the nucleus, either condensed in the porous network (first number), or trapped in the amorphous matrix (second number, in parenthesis). The latter depends on the environment temperature, the molecule (equilibrium pressure, size, polarizability) and its initial abundance in the molecular cloud (Kouchi et al. 1994; Bar-Nun et al. 2007), which leads to very large differences among the trapping efficiencies of various gases (Bar-Nun et al. 2007). Schmitt et al. (1989) showed that amorphous ice can trap other volatiles only up to 8% in mole of water. Unfortunately, there are currently no experimental data that give the relative proportions of volatiles trapped in amorphous ice and in porous network. We therefore choose plausible arbitrary initial concentrations, considering the initial abundances and equilibrium pressure of volatiles. Changing the concentrations by a factor of a few only speeds up or delays the progression of the crystalization and sublimation fronts, but not the actual temperature of the sublimation interface and the vapor pressure in the pores. In other words, this does not change the overall effects described here. The values of J_X are consistent with the observations in cometary comae of molecular species that are directly released from the nucleus (Bockelée-Morvan et al. 2004).

The vapour pressures of the pure condensates are obtained from Fray & Schmitt (2009) and the amorphous to crystalline ice phase transition is postulated exothermic² and irreversible. Thermodynamic quantities characterizing other materials derive from Marboeuf et al. (2010). Table 1 summarizes the orbital parameters adopted in our study for comets 67PCG and 46PW.

The presence of several volatile compounds in the gas phase of the porous network can generate the formation of a multiple guest clathrate (hereafter MG clathrate) whose equilibrium pressure varies as a function of the gas phase composition and the temperature. The equilibrium pressure of the MG clathrate P_c is given by (Lipenkov and Istomin 2001; Hand et al. 2006)

${\mathit{P}}_{\mathrm{c}}\mathrm{=}{\left(\begin{array}{}\sum \end{array}\frac{{\mathit{y}}_{\mathit{i}}}{{\mathit{P}}_{\mathit{i}}^{\mathrm{c}}}\right)}^{-1}\mathit{,}$(3)

where y_i is the mole fraction of the volatile i in the gas phase and ${\mathit{P}}_{\mathit{i}}^{\mathrm{c}}$ the equilibrium pressure of the corresponding clathrate. Here we assume that the clathrate hydrate behaves as an ideally diluted solution and that the ratio of occupancies for cages is constant and the same for all guest molecules. The equilibrium pressures of the different single guest clathrates are modelled by $\log {\mathit{P}}_{\mathit{i}}^{\mathrm{c}}\mathrm{=}\mathit{A}\mathit{/}\mathit{T}\mathrm{+}\mathit{B}$, where A and B are constants (Miller 1961) whose values are given in Table 2 (with ${\mathit{P}}_{\mathit{i}}^{\mathrm{c}}$ in bars and T in K). When the gas pressure is higher than the equilibrium pressure of the MG clathrate, the ice and gas phases can combine to form the clathrate cages in the pores. Below the dissociation pressure, the cages become unstable, which leads to their dissociation, and only the ice and gas phase remain.

4. Results

We computed the thermodynamic evolution of comets 67PCG and 46PW as a function of time and for two angles of insolation. We find that the stability conditions of MG clathrates are permanently met in the subsurface, irrespective of the adopted angles of insolation and dust conductivities. Only the thickness of the clathrate layer is affected by these latter quantities. Figures 1 and 2 represent the evolution of the stratigraphies of comets 67PCG and 46PW over 70 yr and 80 yr timespans, respectively, for a latitude θ of 0° (i.e., the maximum insolation angle) and a dust conductivity of 4 W m^-1 K^-1. The two figures show that the stability zone of the MG clathrate extends from the base of the amorphous-to-crystalline water ice phase transition layer up to the surface of the nuclei. This region grows slowly with time and also always contains the sublimation interfaces of H₂S and CO₂ pure ices. In both comets, the gas phase, from which the MG clathrate can be formed, is dominated by CO (~95%) and CH₄ (~5%) and contains small fractions of H₂S (~0.1%) and CO₂ (~0.05%). When the comet approaches perihelion for the first time, the amorphous layers start to crystallize and CO, CO₂, CH₄, and H₂S released in the pores can be enclathrated by the crystalline water ice available on their surface because the gas pressure (a few hundreds of Pa) is greater than the clathrate equilibrium pressure (order of 1 Pa).

Fig. 1 Stratigraphy of the nucleus 67PCG as a function of time and for θ = 0°. The lines represent the surface and the minimum depths at which solid CO2 and CH4 (dashed lines), H2S and CO (bold solid lines) exist. The dotted area corresponds to the zone where the amorphous-to-crystalline water ice phase transition occurs and the dashed area to the zone of the MG clathrate stability.

Fig. 2 Same as in Fig. 1, but for 46PW.

As shown in Figs. 1 and 2, the thickness of the stability zone of the MG clathrate slightly oscillates with time as a function of the heliocentric distance, but never vanishes. In particular, when approaching the Sun, about the top 10 m of the clathrate layer destabilize and, inversely, become stable again as the nucleus cools. The destabilization of these top layers is caused by an increase of the local temperature at perihelion, which makes the MG clathrate equilibrium pressure rises above the gas pressure in the pores. When the comet moves away from perihelion, the layer cooling decreases the equilibrium pressure of the MG clathrate below the gas pressure in the pores, favouring again the possible formation of MG clathrate layers close to the surface, provided that there is enough available crystalline water ice. Note that at each perihelion passage, at least 1 to 2 m of new clathrate layer would be destabilized for the first time in subsurface, i.e. the gas pressure becomes lower than the clathrate equilibrium pressure. As an example, a region located in the middle of this oscillation zone has a gas pressure above that of clathrate equilibrium during more than three years for each orbit of Comet 67P/Churyumov-Gerasimenko.

It is interesting to note that the top clathrate layers, which are closer to the nucleus surface, are heated and destabilized when the comets approaches perihelion before the amorphous ice will be transformed into crystalline ice. We also performed tests with lower abundances of CH₄ (4 times lower) and H₂S (20 times lower) to the lowest values in Bockelée-Morvan et al. (2004), and our results show that the stability area of clathrates varies only by a few percent.

5. Discussion

Our calculations show that assuming a plausible set of volatile abundances and irrespective of the adopted thermodynamic parameters of the nuclei, JFCs possess a region where clathrates would be stable, extending from the amorphous-to-crystalline water ice phase transition interface up to the surface. Additional calculations lead to the same results for Halley-type comets and then for all SPCs. The main improvement over Marboeuf et al. (2010) is the use of two additional volatiles (CH₄ and H₂S), which increases the stability of MG clathrates. This allows the trapping of CO and CO₂ at higher temperatures than if they were alone. Thus, clathrates containing CO and CO₂ are permanently stable in much larger regions than shown in Marboeuf et al. (2010).

The permanent stability of clathrates in the subsurface of JFCs and SPCs raises two important questions. The first concerns the kinetics of clathrate formation. Do the various gases spend enough time inside the pores before they are outgassed at the surface to actually be trapped in clathrate cages? A rough estimate of this possibility for SPCs can be given with the following argument. The rate of formation/dissociation of clathrates in their stable/instable regions is mainly the function of the kinetics constant λ^cl, the ice surface to volume ratio A_s, and the difference between the gas pressure and the clathrate equilibrium pressure ΔP, and can be writen as λ^clA_sΔP (Kim et al. 1987; Englezos et al. 1987; Schmitt 1986; Sun & Mohanty 2006). On average, one has λ^cl ~ 5 × 10^-12 mol m^-2 Pa^-1 s^-1 (Englezos et al. 1987), A_s ~ 2 × 10⁴ in our model and ΔP ≃ 10−100 Pa in most of the region of stability of clathrates. Integrating the formation rate over the thickness of the clathrate stability layer (~40 m as seen in Figs. 1 and 2) we get a formation rate per unit surface of 4−40 × 10^-5 mol m^-2 s^-1. This is to be compared to the typical gas production rate at the surface of the nucleus of 0.8−4 × 10^-4 mol m^-2 s^-1. Hence the clathrate formation rate is similar to the classical degasing. At the same time, the amount of volatiles present in gas form per unit area of nucleus surface is ~4 mol m^-2. Hence the gaseous reservoir is only capable of supplying material for clathration for a few hours before there is a strong competition between clathration and degasing. Unfortunately, a precise evaluation of this effect can only be done with a fully self-consistent model.

This leads to the second question concerning the thermodynamics of the reactions. The present work is based on a comparison between the equilibrium pressure of clathrates and that of the gaseous phase in the porous network. In our model of a comet nucleus no clathrate forms (and thus decomposes) within the porous network. The influence of clathrate formation/dissociation on the internal temperature of the nucleus, pressure of the gaseous phase and gas production of SPCs is therefore not accounted for. We expect the temperature and gas pressure to be different in a self-consistent model, but we reckon this should not modify our main conclusions on the stability regions because the gas pressure is much higher than the clathrate equilibrium pressure. As for gas production, again only a fully self-consistent model can answer this question. Developping such a model accounting for the kinetics and thermodynamics of the formation/dissociation of clathrates is the subject of a forthcoming paper.

References

All Tables

All Figures

	Fig. 1 Stratigraphy of the nucleus 67PCG as a function of time and for θ = 0°. The lines represent the surface and the minimum depths at which solid CO2 and CH4 (dashed lines), H2S and CO (bold solid lines) exist. The dotted area corresponds to the zone where the amorphous-to-crystalline water ice phase transition occurs and the dashed area to the zone of the MG clathrate stability.
In the text

	Fig. 2 Same as in Fig. 1, but for 46PW.
In the text