© ESO, 2012

1. Introduction

Comets are currently supposed to be the most primitive objects in the solar system. Their chemical composition suggest that cometary material is formed at low temperature in colder regions of the protoplanetary disk or in the interstellar medium (ISM), where most volatile molecules can condensate. The chemical composition and water ice structure included in cometary material are mainly affected by temperature and molecular composition of the surrounding environment both during their formation and their thermodynamical evolution in the protoplanetary disk and solar system after their incorporation in comets. Experimental and theoretical studies have shown that cometary material could be formed of pure crystalline or amorphous water ice, clathrate hydrates, or a mixture of these structures of water ice, depending of the location of formation of the comet in the solar system. The amorphous water ice structure could have been formed in the solar nebula that gave birth to the solar system, or in the ISM. This structure of ice could have been preserved in colder region of solar system (Heliocentric distance R_h ≥ 12 UA, Kouchi et al. 1994) before their incorporation in comets. In the frame of this hypothesis, the water ice structure of comets is currently assumed to be initially amorphous in numerical models (Espinasse et al. 1991; Mekler & Podolak 1994; Enzian et al. 1997; Orosei et al. 1999; Capria et al. 2003; Prialnik et al. 2004; Mousis et al. 2005; Huebner et al. 2006). The amorphous ice structure becomes crystalline in subsurface layers as the comet nucleus is progressively and repeatedly heated by turning around the Sun. Up to now only these two icy structures are considered by current models of cometary nuclei. However, models of ice formation in the protoplanetary disk (Lewis 1972; Gautier et al. 2001a,b; Iro et al. 2003; Hersant et al. 2004; Marboeuf et al. 2008; Mousis et al. 2008, 2009, 2010) show that comets formed beyond 15 UA could be fully made up of crystalline water ice and clathrate hydrates. Moreover, as shown by Marboeuf et al. (2010, 2011), this structure of water ice could form within all cometary nuclei, whatever the initial water ice structure considered in comets (amorphous or pure crystalline).

The clathrates are crystalline solids composed of water and gas. The water molecules structure is organized in the form of cages which are stabilized by the inclusion of gas molecules. Each cage contains a single gas molecule trapped thanks to van der Waals interactions. The formation and decomposition of this structure of ice inside cometary nuclei would add sinks and sources of volatile molecules and of energy. In addition, the thermodynamic properties of clathrate differs from those of the pure crystalline and amorphous water ice structures and can strongly change the thermal behavior of cometary nuclei. As a result, their formation and decomposition inside the nucleus can generate an outgassing profile of volatile molecules at the surface that differ markedly from that expected from models without formation of clathrate since this structure selectively retains and releases the gases trapped inside. The existence of clathrate inside comets has been discussed since a long time (Delsemme & Swings 1952; Klinger et al. 1986; Schmitt & Klinger 1987; Smoluchowski 1988; Mousis et al. 2000; Prialnik et al. 2004; Huebner 2008; Marboeuf et al. 2010, 2011). The authors proposed that the presence of clathrate in extraterrestrial ices might account for their anomalous volatile molecules retention and release (Blake et al. 1991; Iro et al. 2003). Thus, the existence of this structure of ice has often been invoked in order to account for the appearance of various molecules at particular heliocentric distances before or after perihelion of a comet (Smoluchowski 1988). The possible formation/dissociation of clathrate in cometary nuclei and their implication for the thermal evolution and the degassing of these objects have been poorly studied in numerical models. Only a couple of models of comets including such a structure of ice have already been developed (Houpis et al. 1985; Flammer et al. 1998), but unfortunately the hypotheses about the nucleus composition as well as the physics on which they are based are incorrect: they considered that the icy matrix of the cometary nucleus is initially entirely composed of clathrate and that their dissociation occurs only during the sublimation of H₂O ice at the nucleus surface.

The present work aims to present a model of cometary nucleus, which takes into account all structural water ices and phase changes during the thermal evolution of the comet. All initial water ice structures (amorphous, pure crystalline, clathrate hydrates or a mixture of these structures of water ice) can be taken into account in this model following assumptions considered on the formation location of the comet in the solar system. So, all phases changes of water ice (amorphous  →  pure crystalline, amorphous  →  clathrate hydrates and pure crystalline  ↔  clathrate hydrates) are taken into account following the thermodynamical evolution of the comet around the sun and initial physical assumptions considered in the model. In addition to physical processes that are also taken into account by most other models, we added the physics of formation/dissociation of clathrate within the porous network of the nucleus icy matrix. All thermodynamic changes and gas trapping/release induced by clathrate formation/dissociation are included in this cometary nucleus model. The goal of this model is to be able to interpret the outgassing observations of comets, and in particular the future detailed observations that will be made by the Rosetta mission in order to constrain the chemical composition and the water ice structure in the interior of the target comet 67P/Churyumov-Gerasimenko (hereafter 67P/C-G). This article is organized as follows: in a first step, we present the model and describe the physical processes taken into account. In a second step, we discuss the physical assumptions and the thermodynamic parameters adopted. In a third step, we present the outgassing profiles as a function of heliocentric distance of volatile molecules for a model with the same orbit parameters as 67P/C-G and for different initial water ice structures.

2. Description of the physical nucleus model

Fig. 1 Schematic view of the interior of a homogeneous comet nucleus composed of icy grains. The matrix of the nucleus can be formed of icy grains with amorphous, crystalline water ice structures, clathrate hydrates or a mixture of these icy grains and water ice structures.

The original comet nucleus model we use is the one developed by Marboeuf et al. (2008). This model considers a sphere composed of a porous predefined mixture of water ice, volatile molecules (in gas and solid states) and dust grains in specified proportions. It describes heat transmission, latent heat exchanges, amorphous-to-crystalline ice transition, sublimation/recondensation of volatiles in the nucleus, gas diffusion, as well as gas and dust release and mantle formation at the surface. The comet matter modeled here is composed of refractory grains with an icy mantle composed of water and some other volatile molecules (see Fig. 1). Water ice can be initially amorphous or pure crystalline, depending on the formation location of the body in the solar system. When the ice is amorphous, a fraction of the other volatile molecules can be trapped in the amorphous icy matrix, while the remaining is condensed in the pores as pure ices. When heated, the fraction of volatiles condensed in the pores sublimates first, and then the other fraction trapped within the ice matrix is released during the transition from amorphous to crystalline water ice. The gas released diffuses through the porous matrix and eventually escapes to space at the surface of the object.

In addition to these physical processes that are also taken into account by most other models, our model is able to take into account the clathrate structure, or a mixture of amorphous, pure crystalline water ice and clathrate hydrate, as the initial water ice structure in comets depending on assumptions on the formation of cometary material in the ISM or/and protoplanetary disk, and the formation location of the comet in the solar system. Moreover, we added the physics of clathrate formation and dissociation within the icy porous matrix of our comet nucleus model whatever the initial ice structure considered in the model¹. Hence, all phase changes of water ice structure are taken into account in this model. All these processes are presented below. The clathrates are ice-like solids formed from a nonstoichiometric mixture of water and low-molecular mass gases. The water molecules structure is organized in the form of cages which are stabilized by the inclusion of gas molecules. Each cage contains a single gas molecule trapped thanks to van der Waals interactions but it is possible to trap several different molecules in a structure composed of several cages. There are mainly two types of structure of clathrates whose formation depends on the size of the trapped molecules (Handa 1986a; Sloan & Fleyfel 1992; Sloan 1998; Rydzy et al. 2007). The formation/dissociation of clathrate hydrates occurs by the direct gas-ice interaction between gas phase molecules and crystalline water ice or clathrate hydrate, respectively. The formation of cages becomes possible, and the clathrate is stabilized, when the gas pressure is greater than the equilibrium pressure of the clathrate of this gas. At lower pressure, the cages dissociate. Thus, depending on the sign and amplitude of the difference between the gas pressure in the pores and the clathrate equilibrium pressure (temperature dependent), the volatile molecules can be trapped in, or released by, clathrate at rates determined by the kinetic laws of their formation or dissociation. The possible ranges of formation/dissociation rates are determined from laboratory works. We consider that gas and crystalline water ice or clathrate hydrate can always interact together (see footnote ¹). All thermal changes induced by clathrate formation/decomposition are also included in this model, such as the exchange of latent heat during their formation/dissociation and their change in thermal conductivity. Along the orbital evolution of comets, clathrate formation provides a third internal source of energy in addition to the ones of crystallization and gas recondensation and a second sink for volatile molecules (aside gas recondensation). On the other side clathrate dissociation supplies a third source of gas into the porous matrix (with pure volatile ices sublimation and gas released during crystallization) as well as a second sink of energy, in addition to ices sublimation. The transfers of mass and heat in the nucleus that take into account these physical effects are described by the heat diffusion (Sect. 2.1 and Eq. (1)) and gas diffusion (Sect. 2.2 and Eq. (22)) equations given below. In Appendix A, we describe the numerical model of the nucleus to solve the heat and mass equations of conservation. To help the reader to find his/her way through the equations and tables, a list of the principal symbols used in this paper is given in Appendix B. Moreover, in order to facilitate reading of the paper, we have chosen to insert a bullet • in front of paragraphs that define and explain the main physical parameters and assumptions related to the conservation Eqs. (1) and (22).

2.1. Equations of conservation of energy

The heat diffusion through the nucleus is described by the energy conservation equation: (1)with (2)

• T is the temperature (K), t the time (s), r the distance (m) from the center of the nucleus and ρ^l the density (kg m^-3) of the solid component l of the comet nucleus (l = d for dust, cl for clathrate or i for pure ices of elements x).

• c^l is the specific heat capacity (J kg^-1 K^-1) of the solid component l of the comet nucleus (l = d for dust, cl for clathrate or i for pure ices of elements x).

The specific heat capacity of clathrate is assumed to be the sum of the specific heat capacity of pure water ice and the one of the gas molecules trapped in the clathrate structure as described (Handa 1986a; Handa & Tse 1986) by: (3)where is the mass density (kg m^-3) of water molecules that form the clathrate hydrate structure and the specific heat capacity (J kg^-1 K^-1) of water ice. is the mass density (kg m^-3) of volatile molecules x that are trapped in the clathrate and their gas mass heat capacity (J kg^-1 K^-1) at constant volume. Equation (3) assumes that the heat capacity of the empty hydrate lattice is essentially equal to that of water crystalline ice (Handa & Tse 1986). The contribution of the enclathrated guest to the clathrate heat capacity is then given by the constant volume heat capacity of the guest molecule (see Sect. 3.3 for details).

When water ice is amorphous, a fraction of other volatile molecules is trapped inside. The specific heat capacity of this agglomerate becomes then: (4)where is the mass density (kg m^-3) of amorphous water ice, its specific heat capacity (J kg^-1 K^-1) and the mass density (kg m^-3) of the volatile molecules trapped in amorphous ice.

• is the molar latent heat of sublimation of ice x (J mol^-1) and (mol m^-3 s^-1) the rate of moles of gas x per unit volume released (≥0) by amorphous ice during the process of crystallization. Note that the energy needed to release the volatile molecules initially trapped in amorphous ice remains unknown and we assume in this work that it is the same as the latent heat of sublimation (see Eq. (1)). represents the rate of volatile molecule x (mol m^-3 s^-1) that sublimates/condenses (≥0/ ≤ 0) in the pores of the matrix. Its expression is given by the inversion of the gas diffusion Eq. (22) given in Sect. 2.2. This term is present only if the volatile x is condensed, or condenses, in the layer as pure ice otherwise it takes the value 0.

• Y^cr is the power per unit volume released during the crystallization process of amorphous water ice, and is described by (Espinasse et al. 1991): (5)where H^cr is the latent heat of crystallization (J mol^-1), M_H2O the molar mass of water (kg mol^-1), and τ^cr the time (s) of crystallization of the amorphous water ice, as given by Schmitt et al. (1989): (6)

• Y^chs is the power per unit volume exchanged between the gas phase of molecules, which diffuse in the porous network, and the solid matrix described by: (7) is the change of mass density of the gas x per unit time (kg m^-3 s^-1) and ΔT is the difference of temperature of the gas that diffuse through the matrix from the previous deeper layer at the temperature T_i + 1 (i is the index of the layer, see Appendix A) to current layers at the temperature T_i.

• Y^cl is the power per unit volume released/taken during the formation/dissociation of cages (H₂O molecules structure in addition of trapped volatile molecules) of clathrates described by: (8) is the number of moles of gas x trapped/released (≤0/ ≥0) per units volume and time (mol m^-3 s^-1) described in Sect. 2.3 and (J mol^-1) is the enthalpy of formation/dissociation of clathrate per mole of gas x trapped/released. Although is only function of the volatile molecule x trapped/released in/from the clathrate structure, this term takes also into account the energy to form/dissociate the cages of H₂O molecules which trap/release the volatile molecule x. Here, is assumed to be constant with temperature (see Sect. 3.3 for details).

Note that, in the process of clathrate formation by crystallization of an amorphous ice mixture, the energy taken or released remains unknown. We assume, from Eq. (1), that the energy taken or released during this process is equal to the latent heat of formation of clathrate of gas x minus the molar latent heat of sublimation of ice since .

• K^m is the heat conduction coefficient (J s^-1 m^-1 K^-1) of the porous matrix. This parameter is linked to the porosity and to the contact of icy grains with other. It influences directly, with the porosity, the thermal inertia of the nucleus. Groussin et al. (2007) estimate a low thermal inertia (≤50 W K^-1 m^-2  ) for the comet 9P/Tempel 1, implying a low heat conductivity (≤3 × 10^-3 W K^-1 m^-1). Davidsson et al. (2009) re-evalued IR spectra of comet 9P/Tempel 1 obtained by the Deep Impact spacecraft and find a thermal inertia generally high (1000−3000 W K^-1 m^-2 ), although it may be substantially lower (40−380 W K^-1 m^-2 ) in specific areas. Using the assumed thermal inertia () given by Davidsson et al. (2009), the high thermal inertia induces a thermal conductivity between 1 and 10 W K^-1 m^-1 while the low thermal inertia induces values between 10^-2 and 1. Due to the uncertainty of the thermal conductivity of the porous icy matrix, the model is able to take into account either the Russel’s formula (1935) expressed below (see Espinasse et al. 1993), and which gives heat conductivity of about 1 W K^-1 m^-1, or the Hertz factor h expressed below (see Kossacki et al. 1999; Davidsson & Skorov 2002; Prialnik et al. 2004; Huebner et al. 2006), used to correct the effective area of the matrix material through which heat flows, and which can impose heat conductivity between 10^-3 and 1 W K^-1 m^-1.

Russel’s formula (1935) (see Espinasse et al. 1993):

(9)where Ψ is the porosity of the matrix and K^p the radiative conductivity across the pores described by Squyres et al. (1985): (10)where ϵ is the infrared surface emissivity of the nucleus, σ the Stefan-Boltzmann constant and r_p the radius of the pores.

K^s is the conductivity of the solid phase of the components (dust and ices) described by: (11)where f^l is the volume fraction of the component l (with l = i for the solid phase of the icy components x, d for the dust grains and cl for the clathrate) in the solid matrix and k^l its conductivity. Note that the term f^clk^cl of the conductivity of clathrate includes the water structure and the gases trapped inside. In this study, the thermal conductivity of the clathrate structure is set independent of the quantity and composition of gases trapped inside because we are missing experimental data to constrain these dependencies (See Sect. 3.3). f^l is given by: (12)where ρ^l is the mass density of the element l by nucleus unit volume and ρ^b,   l the bulk density of the non-porous solid phase of the component l. For clathrates, f^cl gives: (13)with . and are the masses of H₂O and gas x molecules that form the clathrate per unit volume of nucleus. ρ^b,   cl is defined (Tonnet 2007; Thiam 2008) by: (14) is the bulk mass density of water molecules in clathrate (kg m^-3) and the mean molar mass of gases trapped in the clathrate (kg mol^-1) defined as:

(15) is the molar mass of the gas x and the molar fraction of gas x trapped in the clathrate. n_hyd is the hydrate number defined as the molecular ratio in the clathrate:

(16)where and are respectively the number and the mass density (kg m^-3) of gas molecules that are trapped in the clathrate structure, and the number of water molecules forming the bulk clathrate structure. This term gives the number of water molecules for each gas molecule trapped (see Sect. 3.3 for its value).

Hertz factor formula:

If the area of contact between adjoining grains is small, the resultant conductivity of the medium will be reduced even further by a Hertz factor h (Davidsson & Skorov 2002; Prialnik et al. 2004): (17)K^s is the conductivity of the solid phase of the components (dust and ices) described by Eq. (11). h is expressed by considering two spheres of radius R that are pressed together and have a contact area of radius r_c, so (Kossacki et al. 1999). Its value can vary between 10^-3 and 1 (see Davidsson & Skorov 2002; Huebner et al. 2006; Volkov & Lukyanov 2008).

2.1.1. Conservation of energy at the surface of the nucleus

At the surface, the local temperature is given by a thermal balance between the solar energy absorbed by the surface on one side and on the other side, from left to right, its thermal emission, the heat diffusion towards the interior and the energy of sublimation of the ice species existing at the surface: (18)where C_s is the solar constant (W m^-2), A_b the bolometric Bond albedo, R_h the heliocentric distance (UA), and ξ the solar zenith distance calculated as (see Sekanina 1979; Fanale & Salvail 1984; Prialnik 2004; Gortsas et al. 2011): (19)θ is the latitude on the comet, t the time since the beginning of the computation (s), t₀ the initial time of computation (s), and where P_r is the nucleus rotation period of the comet (s). θ_s is the cometocentric latitude of the sub-solar point given by: (20)δ is the obliquity, φ the argument of the sub-solar meridian at perihelion. The true anomaly γ of the comet is calculated by using Kepler’s equations.

Note that using a 1D model, the insertion of the term cos   θ   cos(ω   (t − t₀)) in the left hand side of Eq. (18) assumes that the lateral transfers of energy are negligible in the comet. Also, using this formula, we study the thermodynamic and outgassing behavior of the comet only on one point on it surface.

ϵ is the infrared surface emissivity, σ the Stefan-Boltzmann constant (W m^-2 K^-4), the surface fraction covered by the pure ice species x (including H₂O) and the one covered by water molecules that form the clathrate structure. is the molar fraction of the volatile x (relative to water) initially trapped in the clathrate structure. ϕ_x is the free sublimation rate of the species x (mol m^-2 s^-1) given by the expression of Delsemme & Miller (1971): (21)where M_x is the molar mass of the corresponding gas specie, its vapor sublimation pressure (Pa), and R the perfect gas constant (J mol^-1 K^-1). Here, the last term of Eq. (18) assumes that the water molecules still present in the clathrate structure at the surface of the nucleus (if not dissociated before) are subjected to the same rate of sublimation as for pure crystalline ice.

2.2. Equation of conservation of mass

For each molecule x the diffusion of gas through the matrix pores is described by the mass conservation equation:

(22)

• is the the net sources of gas x released in the pores during water ice crystallization: (23)where is the initial x/H₂O mole fraction of the gas x trapped in amorphous ice.

• and are respectively the pure ice sublimation/condensation (see part 2.1) and clathrate dissociation/formation (see part 2.3) rates. is the mass density of gas x.

• Φ_x its molar flow given by: (24)where P_x is the partial pressure of gas x (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 a free molecular (Knudsen) flow (see 1. below), a visquous flow (see 2. below), or a mixture of them (see 3. below). The mechanism governing the diffusion of volatile molecules is determined by the Knudsen number Kⁿ of the gas mixture: (25)r_p is the average radius of the pores (m) and λ the mean free path of molecules from kinetic theory: (26)k_B is the Boltzmann constant (J K^-1), P_t the total pressure of the gas (Pa) and d the mean diameter of gas molecules defined as: (27)where d_x is the diameter of molecule x (m) and P_x its partial pressure (Pa).

• 1.

  When Kⁿ is greater than 1, molecule-wall collisions predominate over molecule-molecule collisions (Knudsen 1909; Kast & Hohenthanner 2000). The gradient of the partial pressure is the driving force and each specie can diffuse independently of others. The gas diffusion through the porous matrix is then described by Knudsen (free-molecule) flow and the coefficient of diffusion G_x of the molecule x in a single cylindrical pore is (Mekler et al. 1990): (28)with τ the tortuosity that is defined as the ratio of the length of a pore to the distance between its ends. Further on, we adopt the appropriate value for the tortuosity of an unconsolidated medium (Carman 1956; Mekler et al. 1990; Kossacki & Szutowicz 2006).

• 2.

  When Kⁿ is less than 10^-2, the mean free path of molecules is much lower than the diameter of pores: molecule-molecule collisions predominate over molecule-wall collisions. The gas diffusion is then viscous and the coefficient of diffusion G_x for a single cylindrical pore becomes (see Espinasse et al. 1991): (29)with η the dynamic viscosity of gas derived from kinetic theory of gases as: (30)ρ_t is the total gas density (kg m^-3) and v^t the thermal velocity of molecules: (31)where M_gas is the molar mass of the gas (kg mol^-1).

• 3.

  For values of Kⁿ between 10^-2 and 1 (transition region between the free-molecule flow and the viscous flow), the mean free path of molecules is of the size of the pore diameter. In this region, there are molecule-molecule collisions as well as molecule-wall collisions: mass transfer occurs due to Knudsen flow and viscous flow. The total flux can then be described as the sum of the fluxes due to viscous flow and due to Knudsen flow (Knudsen 1909; Kast & Hohenthanner 2000). Fanale & Salvail (1987) used, for the transition regime of CO₂ in porous material, the following equation: (32)Coefficients 0.9 and 0.5 are issue from an adaptation of an equation presented by Scheidegger (1974). These coefficients are valid only for CO₂ since they may be different for other chemical species (Bouziani & Fanale 1998). However, they have been used by Espinasse et al. (1991) for the transition regimes of CO, CO₂ and H₂O. In this work, we will use it regardless of the chemical specie that diffuse in the porous matrix in this transition region since no experimental data concerning the flow of other species has been reported in the literature.

In the case of a parallel network arrangement of cylindrical pores², the total molar flux of gas x perpendicular to the surface can be written as: (33)where N_p is the number of cylindrical pores per unit of surface (m^-2) defined as: (34)We finally obtain: (35)Note that the gas is always supposed to be perfect so we can link the pressure P to the gas density following the ideal gas equation of state .

All ices in a cell of the matrix are at the same temperature. Similarly, all gases are in thermal equilibrium at the temperature of the corresponding cell. When the volatile species x is condensed in the pores, we impose equilibrium between the gas and solid phases. Hence the partial gas pressure equals the vapor saturation pressure: . Otherwise the value of the pressure is obtained by solving the gas diffusion equation (Eq. (22)).

2.3. Term of formation/dissociation of clathrates Q

The presence of several volatile molecules in the gas phase of the porous network can generate the formation of multiple guest (hereafter MG) clathrates whose equilibrium pressure varies as a fonction of the gas phase composition and the temperature. The equilibrium pressure of the MG clathrate is given by (Lipenkov & Istomin 2001; Hand et al. 2006): (36)where is the mole fraction of the volatile x in the gas phase and the equilibrium pressure of the corresponding single guest clathrate. This equation is only valid when the fractions of volatile molecules in the gas phase and trapped in clathrate are identical, but we will use it regardless, in the absence of a current better and well established model. Indeed some statistical thermodynamic models aims at determining the clathrate equilibrium pressure and the composition of the gas trapped inside MG clathrates in equilibrium with the gas phase (van der Waals & Platteeuw 1959; Parrish & Prausnitz 1972; Lunine & Stevenson 1985; Thomas et al. 2009) but they depend on several parameters that are poorly constrained or unknown for some molecules. So we decided here, as a first step, to use the equilibrium pressure of clathrates given by Eq. (36) whatever the difference in composition between the gas phase and the gas in MG clathrates.

When the total gas pressure P_t in the pores is higher than the equilibrium pressure of the MG clathrate , the water ice and the gas phase can combine to form cages of clathrate at the surface of the pores. Below the equilibrium pressure, the cages become unstable leading to their dissociation into crystalline water ice and to the release of the trapped gases.

The molar amount of volatile molecules trapped or released by unit of nucleus volume and per second during the formation/dissociation (≤0/ ≥0) of clathrate is given by the term : (37)f^l is the volume ratio (see Eq. (12)) of water ice (l = i when formation of clathrate occurs) or of clathrate in the solid matrix (l = cl when dissociation of clathrate occurs), is the number of moles of gas trapped (k = f for formation of clathrate) or released (k = d for dissociation of clathrate) by the clathrate per unit volume of solid water ice (formation) or of clathrate (dissociation) and per second described by (Sun & Mohanty 2006; Clarke & Bisnoi 2004; Englezos et al. 1987; Kim et al. 1987; Schmitt 1986): (38) (mol m^-2 Pa^-1 s^-1) is the kinetic parameter of formation (k = f) or dissociation (k = d) of clathrate, P_t the total gas pressure, the equilibrium pressure of MG clathrate given by Eq. (36) and A_s (m^-1) the reaction surface area per unit of volume of solid reactant defined as the total surface area of the pores S_p divided by the volume of solid: (39)Note that, when water ice undergoes crystallization, we assume that the rate of clathrate formation is the one of crystallization. This assumption is justified as the mobility of the water molecules, the limiting factor of clathrate formation, is high during crystallization. It is only taken into account if the equilibrium pressure of MG clathrate is lower than the total gas pressure P_t and if the volatile molecules needed to form cages are in sufficient amount. The number of moles of gas that form clathrate per unit of volume and per second given previously by Eq. (37) becomes then: (40)The number of moles of molecule x trapped/released (≤0/ ≥0) in/by clathrate per unit of nucleus volume and per second is defined as: (41)where is the molar fraction of the molecule x trapped in the case of formation of MG clathrate (l = f) or released in the case of dissociation of MG clathrate (l = d). The calculation of this parameter being somewhat complex and uncertain at low temperature and low pressure, we decided to use here a simple law that gives with a good approximation (see Sect. 3.3 for details) the mole fraction of gas trapped during the formation of cages: (42)with . We assume in this model that this equation remains valid whatever the total gas pressure P_t of the gas phase in the porous network provided it is greater to the MG clathrate equilibrium pressure . So, the mole fraction of gas trapped doesn’t directly depend of the total gas pressure P_t but on the gas composition and on the equilibrium pressure curves of single guest clathrates of molecules x present in the gas phase.

When the clathrate equilibrium pressure is greater than the total gas pressure, the gas trapped in the cages can be released and its composition is then given by: (43)with . is the molar fraction, averaged over the layer, of the volatile x trapped in the clathrate structure (see Sect. 3.3 for details).

In order to help the reader to follow kinetic law adopted and power exchanged during water ice transitions, a resume is given in Fig. 2.

Fig. 2 Schematic view resuming the kinetic law adopted and the power released or taken by water ice structure during its transition.

2.4. Porosity and pore radius changes

At the end of each time step Δt, we calculate the variation of mass density of species x in each layer using (l = i for solid phase and cl for clathrates of elements x): (44)For a more realistic physical representation of the nucleus, the porosity and the radius of the pores are also recomputed for each layer: where ρ^l and ρ^b,l are the mass density of the solid phase and bulk density of the component l, and Ψⁱ and the initial porosity and pores radius.

2.5. Dust ejection and mantle formation at the surface of the nucleus

The radius a of dust grains in the cometary nucleus is given by a power law size distribution (Rickman et al. 1990): (47)β is the power law of the size distribution and N₀ a normalization factor. The grains are initially encased in H₂O ice and can be freed from it by its sublimation at the surface. At this time, the grains can either be ejected from the nucleus or accumulate at its surface thus forming a dust mantle covering the icy layers. A critical radius gives the largest dust grain a^∗ that can be ejected from the comet. It is computed by comparing the sum of gas drag and centrifugal force with the gravitational attraction of the nucleus (see Orosei et al. 1999): (48)where Φ_x is the flux of the gas x at the surface of the comet (kg m^-2 s^-1) and V_x its velocity (m s^-1). ρ^b,d is the bulk mass density of the non porous dust grain (kg m^-3), G_c is the gravitational constant (m³ kg^-1 s^-2), M_n is the comet nucleus mass (kg) and R_n its radius (m). Dust grains whose radius is smaller than r^∗ are immediately lost to space, while the larger grains stay at the surface and contribute to the formation of a dust mantle. In this model, no cohesive forces between grains (see description in Huebner et al. 2006) are taken into account. The exact value of a cohesive force would remain very uncertain although it is possible that a cohesive energy between particles dominates compared to gravitational energy of the nucleus (see Huebner et al. 2006). This choice is at least consistent with the assumption that the dust flux is proportional to the gas flux in models and observations (Jewitt et al. 1999; Prialnik et al. 2004; Rosenberg & Prialnik 2009).

3. Discussion about the physical assumptions and thermodynamic parameters adopted

3.1. Volatile molecules considered in the model

At the beginning of the computation the cometary nucleus has a homogeneous composition made of ices and dust. The initial ice phase of our model can be any type of water ice (amorphous, crystalline, clathrate or a mixture of these structures) mixed with CO, CO₂, CH₄ and H₂S, either as pure ice phases, or trapped in amorphous ice or in clathrate. These four molecules are the most abundant volatile species (production rates relative to water greater than 1%) observed in cometary comas (Bockelée-Morvan et al. 2004) whose data on the equilibrium pressure of single guest clathrates exist. The other equally abundant molecules H₂CO, CH₃OH and NH₃ are not considered in this model for the following reasons. H₂CO is not mostly produced in the nucleus of comets but is rather the result of a distributed source in the coma, which could be provided by the photo and thermal degradation of dust (Fray et al. 2004; Fray et al. 2006; Cottin & Fray 2008). Moreover, no data is available on its stability curves either in the form of clathrate or as a pure condensate (Fray & Schmitt 2009) and thus this molecule cannot be considered in the mixture of ices of our comet nucleus model. In addition, to our best knowledge, no experimental data concerning the stability curve of the CH₃OH clathrate has been reported in the literature (Marboeuf et al. 2008) and the conditions under which it forms stoichiometric hydrates or clathrate (Blake et al. 1991; Notesco & Bar-Nun 2000) is still unclear. Finally, NH₃ doesn’t form clathrate, but rather stoichiometric hemihydrates (2NH₃ − H₂O) and/or monohydrates (NH₃ − H₂O) under some conditions (Lewis 1972; Bertie & Shehata 1984; Lunine & Stevenson 1987; Kargel 1998; Moore et al. 2007). The volatile molecules considered in the model can be condensed as pure ices in the porous nucleus network and/or trapped in the amorphous water ice or in the clathrate hydrate structure. Their initial distribution between these states strongly depends on the environment temperature, the molecule (equilibrium pressure, size, polarizability) and their initial abundance in the molecular cloud (Kouchi et al. 1994; Bar-Nun et al. 2007) or in the protoplanetary disk, leading to a very large diversity among trapping efficiencies of various gases in amorphous water ice (Bar-Nun et al. 2007) and clathrate hydrates. However there are some constraints on these trapping efficiencies. Schmitt et al. (1989) showed that amorphous ice can firmly trap other volatile molecules inside its structure only up to a total of 8% (in mole, relative to water) until crystallization occurs (Schmitt et al. 1992). The clathrate hydrate structures can trap until n_hyd volatile molecules (hereafter n_hyd = 6, see Sect. 3.3). However, up to now no general agreement exists on the respective amounts of trapped gas (in amorphous water ice or clathrate hydrate) and gas condensed as pure ices in the initial comet nucleus material. A sensitivity study of the model on different initial scenarios of repartition between species trapped and condensed should thus be conducted (Paper II).

3.2. Thermodynamic parameters

The thermodynamic parameters of the materials are given in Table 1. The modeling of the thermodynamical processes depends on several hypotheses that we will explain here in more details. We also provide the assumptions we made on the values of their parameters.

Assumptions on the phase transition from pure amorphous to crystalline water ice:

the phase transition from pure amorphous to crystalline water ice is exothermic and irreversible but, as it has been shown experimentally by Kouchi & Sirono (2001), the crystallization of amorphous mixtures made of water and some other volatile molecules can become endothermic. This is fully included in our model as the overall energy liberated, or absorbed, during the crystallization process is assumed to be the crystallization energy of the water ice minus the sublimation energy of each of the trapped species: (49)This assumes that (1) the crystallization energy of the amorphous water structure is the same for mixed H₂O-dominated ices as for pure water ice and that (2) the energy required to expel the guest molecules from the amorphous water lattice is the same as the sublimation energy from its pure solid. Both these energies are possibly larger for mixed ices, and should partly compensate in the sum, but without experimental knowledge of their values and dependences they cannot be taken into account in our model.

Densities, saturation pressure and enthalpy of sublimation of volatile ices:

the density and saturation pressure of all the pure volatile ices taken into account in the model are given in Table 1. Expressions of saturation pressure are taken from the critical compilation by Fray & Schmitt (2009). Enthalpy of sublimation of volatile ices are determined using the Clausius-Clapeyron equation given by: (50)This equation supposes that the gas is considered as perfect and that the molar volume of the corresponding solid ice phase can be neglected against that of the gas phase.

Heat capacities of volatile ices:

the heat capacities of CO, CO₂ and H₂S ices and their variation with temperature are obtained from fits of the data from Clayton & Giauque (1932), Giauque & Egan (1937) and Giauque & Blue (1936) respectively. CO has two crystalline phases (Clayton & Giauque 1932) with a discontinuity of its heat capacity at 61.55 K (temperature of the phase transition). The heat capacity of CH₄ ice is assumed to be the same as the one of CO because of the lack of data. This discontinuity for solid phases of CO and CH₄ does not generate numerical instability in the model, since these volatile ices remain minor elements in relation to water ice and dust in comets. The heat capacity of crystalline water ice comes from Giauque & Stout (1936). The one of amorphous water ice is assumed to be about 111 J kg^-1 K^-1 higher than the values for crystalline water ice below the glass transition temperature, this increase being probably due in part to the slightly weaker hydrogen bonding in low-density amorphous ice (Handa & Klug 1988).

Heat conductivities of volatile ices:

the heat conductivity of CO ice is obtained by fiting data from Stachowiak et al. (1998) between 20 and 37 K. We assume here that the extrapolation of this law is valid between 37 and 60 K. Thoses of CO₂ and CH₄ ices are obtained by fitting the data from Koloskova et al. (1974) and Manzhelii & krupskii (1968). The heat conductivity of H₂S is assumed to be the same as CO because of the lack of data. The heat conductivity of amorphous ice containing trapped molecules is assumed to be the same as that of pure amorphous ice. The heat conductivity of dust grains is from Ellsworth & schubert (1983). For a dust mantle at the surface of the nucleus, laboratory experiments (see Grn et al. 1993; and Huebner et al. 2006 for a review), showed that the thermal conductivity of dust layers composed by silicate particles is very low (about 10^-2 to 10^-3 W m^-1 K^-1). Following these results, we assigned a value of 0.01 W m^-1 K^-1 to the thermal conductivity of the dust mantle.

3.3. Modeling parameters of clathrates

The physical parameters for modeling the formation of clathrates are given in Tables 2 and 3. There are mainly two types of structure of clathrate whose cages number and size differ and depend on the size of the trapped molecules. The volatile molecules we consider (CO, CH₄, CO₂, H₂S) can be trapped in clathrates, and all create single guest clathrates of structure I individually (Davidson et al. 1987, Dartois 2011; Handa 1986b; Anderson 2003; Miller 1961). So, we consider in this model that only the structure I is formed inside the comet nucleus during the formation of MG clathrates. The unit cell of structure I is cubic, of volume V^cl equal to 1.741 × 10^-27 m³ and contains 46 water molecules (). This structure can trap up to 8 gas molecules respectively in 6 large and 2 small cavities. The ideal hydrate number n_hyd (see Eq. (16)) in the case of a maximum occupancy of the cages is then equal to 46/8 = 5.75. In reality, the occupancy rate of cages depends on the thermodynamic conditions (temperature, gas pressure), and on the size, shape and abundance of the gas molecules trapped. If only the large cages are filled, the hydration number is 46/6 = 7.67. However, the occupation of the cages by molecules depends on their size but also on their relative abundance for mixed clathrates. It is currently impossible to determine the occupancy of the cages in the thermodynamic conditions existing within cometary nuclei and we therefore fix the hydrate number to 6, closer to real values (Handa 1986a,b; Sloan 1998; Sun & Mohanty 2006), as the nominal value of the number of hydration. This assumption is in good agreement with Avlonitis et al. (2005) who showed that the hydration number decreases slowly with temperature if both types of cavities are occupied by gas.

3.3.1. Heat capacity of clathrates

The heat capacity of the clathrate hydrate structure is assumed, in this model (see Eq. (3)), equal to the sum of the heat capacity of pure crystalline water ice and of the constant volume heat capacity of the guest molecule (the contribution of the enclathrated guest to the crystalline water ice structure). This model predicted successfully the heat capacity of clathrates of small molecules such as Ar, Kr, N₂, O₂, CO and CH₄ (Parsonage & Staveley 1984; Handa 1986c; Handa & Tse 1986; Avlonitis 1994). However, for larger molecules and multicomponent systems trapped in clathrate structure, Eq. (3) is no more accurate since the heat capacity of the empty hydrate lattice differs from that of crystalline water ice (Avlonitis 1994). Calculations show that the heat capacity of the hydrate lattice alone can be either higher or lower than the heat capacity of crystalline water ice, depending on the guest molecule. Similarly, the partial molar heat capacity of the enclathrated gases can be higher or lower than the corresponding experimental ideal gas heat capacity. These differences depend on the size of the guest relative to the two clathrate cavities, the hydrate number (see Eq. (16)), and the temperature (Avlonitis 1994). Since there is not yet enough experimental data for a better estimation of the heat capacity of MG clathrates, we choose to keep the law given by Eq. (3) for MG clathrates, irrespective of the nature of the guest molecules trapped in the clathrate structure.

3.3.2. Equilibrium pressure of single and multiple guest clathrates

The equilibrium pressure curves of single guest clathrates of molecules x (with x =  CO, CO₂, CH₄ or H₂S) required to calculate the equilibrium pressure of the MG clathrate given by Eq. (36), are fits of experimental data by a semi-empirical law: (51)where A and B are constants (Miller 1961), is in Pa and T in K. Their values are given in Table 2.

3.3.3. Enthalpy of formation and dissociation of clathrates

The enthalpy of the clathrate formation reaction H per mole of encaged gas x, is the heat released when one mole of gas and n_hyd moles of crystalline water ice are converted into clathrate at the prevailing conditions of thermodynamic equilibrium. It ’s effective value can strongly influence the thermal behavior and the physical evolution of the nucleus. Experimental measurements of the enthalpy of formation/dissociation are not easily obtained and such data are very scarce (Avlonitis 2005). Only the enthalpy of formation of clathrate hydrates of structure I as Xe, CH₄, C₂H₆ (Handa 1986a,b), and clathrate hydrates of structure II as tetrahydrofuran, ethylene oxide (Leaist et al. 1982), Kr, C₃H₈ and iC₄H₁₀ (Handa 1986a,b) have been measured, but in a limited temperature range for each molecule, all around or above 200 K. Avlonitis (2005) developed and implemented a method for predicting enthalpies of clathrate formation and demonstrated that a major energetic component in hydrate formation is the heat of enclathration, defined as the residual enthalpy of the enclathrated gas, and that the enthalpy of transformation of crystalline ice into the empty clathrate structure is comparatively small. So, the enthalpy of formation/dissociation of single guest clathrates below the ice point is mainly a property of the gas itself: it depends solely on the nature of the enclathrated molecule. The enthalpy of formation of single guest clathrates is then independent of the clathrate structure, the hydration number n_hyd, the gas distribution between cavities, the temperature and the pressure (Avlonitis 2005). His model corroborates in every case the experimental data of Handa (1986a,b) and Anderson (2003) as well as the results of Yoon et al. (2003), which were derived by application of the Clausius-Clapeyron equation (See Eq. (50)). This equation is a good approximation to calculate the formation/dissociation heat of single guest clathrates and show fair agreements with experimental data (Handa 1986a,b). Roberts et al. (1940), Barrer & Edge (1967), and Skovborg & Rasmussen (1994) derived analytically this Clausius-Clapeyron equation for simple gas hydrates. Furthermore, the values of enthalpy of clathrate dissociation obtained from the equilibrium pressure data of CH₄ and Xe clathrates measured by Fray et al. (2010), using the Clausius-Clapeyron equation, fit very well the data from Handa (1986a,b), with a relative difference of only 0.2% and 0.9% respectively. Given, the simplicity of the Clausius-Clapeyron equation and the good correlation with experimental values, we can use this equation to obtain approximate values of the enthalpy of dissociation of single guest clathrate hydrates. Moreover, the model of Avlonitis (2005) showed that the values of enthalpy of all of the seven simple gas hydrates (N₂, CH₄, CO₂, Xe, C₂H₆, C₃H₈ and iC₄H₁₀) he studied in the range 180−270 K are so slightly decreasing with temperature that they may be considered approximately constant. However, Anderson (2004), by using the Clausius-Clapeyron equation, has shown a decrease of about 10% of the enthalpy of dissociation of the CH₄ clathrate at 150 K compared to its value at 200−220 K. It is also known that vary with T as follows (Makogon & Sloan 1994): (52)with the heat capacity difference between the ice and clathrate structures at temperature T: (53)where and are the heat capacities of the clathrate of gas x and of pure water ice, respectively, and the heat capacity of the gas x. However, vary with temperature and the guest molecules in the clathrate structure. These variations depend on the size of the guest relative to the cavity, the hydrate number n_hyd (see Eq. (14)), and the temperature (Avlonitis 1994). Unfortunately, there is no experimental data allowing to an estimate and hence in the 100−200 K range temperature for the guest molecules CO, CO₂ and H₂S used in this paper. However, the heat capacities of clathrate hydrates have been measured over a wide temperature range, moslty from 85 K to 270 K, for clathrate hydrates of structure I as Xe, CH₄, C₂H₆ (Handa 1986a,b), and for structure II as tetrahydrofuran (Leaist et al. 1982; Handa et al. 1984; Yamamuro et al. 1988a), ethylene oxide (Leaist et al. 1982; Yamauro et al. 1990), propylene oxide, 1.3-dioxolane, 2.5-dihydrofuran, 1.3-dioxane (Handa 1985), acetone (Kuratomi et al. 1991), Ar (Yamamuro et al. 1988b), Kr, C₃H₈ and iC₄H₁₀ (Handa 1986a,b). Using the gas heat capacities of CH₄ (McDowell & Kruse 1963), C₃H₈ and C₂H₆ (Chao et al. 1973), tetrahydrofuran (Chao et al. 1986) and assuming the values of R for the noble gases Xe and Kr, we obtain a ratio that doesn’t exceed 5% between 85 and 270 K for all these single guest clathrates. Given the small variation of H with temperature, we assume hereafter that it remains constant in the model (see Table 3). Moreover, this assumption is justified since the simple law adopted for the equilibrium pressure of single guest clathrates (Eq. (51)) assumes implicitly that H provided by using the Clausius-Clapeyron equation (Eq. (50)) remains constant over its valid temperature range.

An additional difficulty remains with the calculation of the enthalpy of MG clathrates since it depends on the type and abundance of the different molecules trapped inside. Recently, Rydzy et al. (2007) studied the influence of incorporation of guest molecules such as C₂H₆, C₃H₈ and CO₂ in methane clathrate on its thermal behavior. They concluded that the enthalpy of formation of MG clathrates increases with increasing size of the guest molecules or more precisely with increasing guest-to-cavity size ratios for the guests in the large cages (and hence the cavity occupation of molecules). This is contradictory to the hypothesis proposed by Sloan & Fleyfel (1992) that the heat of dissociation should be, within size constraints, independent of the type and concentration of the mixed guest gas molecules. In addition, Hachikubo et al. (2008) showed that the dissociation enthalpy of MG clathrates composed of CH₄ and C₂H₆ increases with ethane concentration. To conclude, the value of the enthalpy of formation/dissociation of MG clathrates remains difficult to predict as it depends on temperature, molecular abundances and guest to cavity size ratios. We thus use in our model the simple relation given by Eq. (8) that provides the energy released/taken per unit volume of clathrate formed/dissociated as a linear combination of the single guest clathrate values, following the abundance and type of gas molecule trapped/released in/from the clathrate structure. The values of the enthalpy of formation of the single guest clathrates listed in Table 3 are mostly provided by the Clausius-Clapeyron equation.

3.3.4. Thermal conductivity of clathrates

The thermal conductivity of clathrates differs drastically from that of crystalline water ice both in magnitude and temperature dependence (Krivchikov et al. 2007) because their cages are very poor thermal conductors. In contrast to the crystalline structure of H₂O ice, in which the thermal conductivity falls with increasing temperature following a T^-1 dependence (for T > 100 K), the thermal conductivity of clathrate hydrates mostly increases slightly with increasing temperature (Tse & White 1988). This is a behavior very similar to the ones of the thermal conductivity of amorphous solids (Tse & White 1988; Krivchikov et al. 2005a,b) but with smaller values than for amorphous H₂O ice. There are limited data on the conductivity of clathrates at low temperatures. Their thermal conductivity has been measured on gas hydrates of structure I as CH₄ (Krivchikov et al. 2005a), Xe (Handa & Cook 1987; Krivchikov et al. 2006) and ethylene oxide (Cook & Laubitz 1983), and gas hydrates of structure II as tetrahydrofuran (Ross & Andersson 1982; Tse & White 1988; Andersson & Suga 1996; Krivchikov et al. 2005b), dioxolane (Andersson & Ross 1983; Ahmad & Phillips 1987) and cyclobutanone (Andersson & Ross 1983) in a wide range of temperatures. It is very challenging to measure the bulk thermal conductivity of a homogeneous continuous solid (Krivchikov et al. 2007). This bulk value is thus deduced from the measured effective thermal conductivity of powders using simple models of heat transfer (parallel layers model, Jagjiwanram model, Adler model, Rayleigh model, geometric model or Maxwell model; see Krivchikov et al. 2007, for a brief description) in inhomogeneous media disregarding the thermal resistance of the boundary between the two media. This results in a dispersion of the thermal conductivities values obtained for each guest molecule, but by less than one order of magnitude. Moreover, the thermal conductivity of clathrates depends on the guest molecules trapped inside and decreases when its size increases (Andersson & Ross 1983). However, at high temperature (170 K), the thermal conductivities of clathrates of CH₄, Xe and tetrahydrofuran differ by only 20% (Krivchikov et al. 2007). The temperature dependences of the thermal conductivity of structure I of Xe and CH₄ are similar (Krivchikov et al. 2006) but are different from the ones of structure II of tetrahydrofuran, dioxolane and cyclobutanone. The temperature dependence of the thermal conductivity of Xe and CH₄ is divided into four (I−IV) distinct temperature regimes. The temperature dependences in the intervals I−II (below 54 K) and IV (greater than 95 K) are similar to what is observed in tetrahydrofuran and dioxolane clathrate hydrates of structure II (Krivchikov et al. 2006). The thermal conductivity grows with temperature and the curve has a shape typical for amorphous substances. In the intervals I−II (below 54 K), the temperature dependence of clathrate hydrate is only weakly dependent on the type of clathrate structure (Krivchikov et al. 2006). In these intervals the thermal conductivity is practically independent on the nature of the guest molecules. In the interval IV (above 95 K), the thermal conductivity is approximately proportional to T. In the interval III (54−95 K) the thermal conductivity of the clathrates of Xe and CH₄ exhibits an anomalous behavior: the thermal conductivity decreases by more than 50% as the temperature increases for Xe but is less pronounced for CH₄. For temperatures greater than 80 K (mainly in the interval IV), Krivchikov et al. (2005a,b) fited the thermal conductivity of the single guest clathrate of CH₄ with a good accuracy by a linear dependence with T (law given in Table 3). Moreover, the thermal conductivity of the CH₄ hydrate coincides within the experimental error with the thermal conductivity of the dioxolane clathrate (Krivchikov et al. 2005a,b). At present there is no theoretical model permitting a quantitative description and prediction of the thermal conductivity of clathrate hydrates in a wide range of temperatures (Krivchikov et al. 2006). Thus, in the absence of more data on the thermal conductivity of single and MG clathrates, we adopt the temperature-dependent thermal conductivity of the methane clathrate given in Table 3 and use it independently of the clathrate composition. The thermal conductivity of methane clathrate between 10 and 80 K is obtained by fitting the data from Krivchikov et al. (2005a,b, 2007). Above 80 K, it is approximated by a linear temperature dependence (Krivchikov et al. 2005a,b).

3.3.5. Kinetics of formation and dissociation of clathrates

The kinetics parameters, , , of formation/dissociation of clathrates from/to crystalline water ice and gas are poorly constrained at low temperature. Conventionally, two main stages of the clathrate formation process are distinguished (Kuhs et al. 2006). The first stage, relatively rapid, corresponds to the formation of a clathrate film over the crystalline ice surfaces (stage I), and the second, wich dominates when the whole ice grain is covered by a clathrate layer, is the growth of the shell of clathrate phase around ice grains (stage II). This last stage includes two steps: (1) the clathration reaction itself at the inner (ice-clathrate) as well as external (clathrate-gas) interfaces and (2) gas and water mass transports (diffusion) through the clathrate layer. In this paper, we assume that the first stage mainly dominates the kinetics of formation of the cages at the surface of the pores because the ice matrix is very porous and the grains diameter very small. It means that we consider that the maximum thickness of the clathrate shell that can be formed during this first stage is always larger than the radius of our micron sized grains. The second, much slower, stage of grain growth by diffusion of water molecules and gas through the clathrate shell is then neglected. For the first stage, Schmitt (1986) determined values between 3 × 10^-13 and 10^-12 mol m^-2 Pa^-1 s^-1 for the formation rate of CO₂ clathrate from CO₂ gas and crystalline water ice at temperatures around 200 K. Staykova (2004) and Kuhs et al. (2006) determined values of the kinetics parameter of methane clathrate formation between 245 and 272 K starting from hydrogenated and deuterated ices. They showed a temperature dependence with values ranging between about 9 × 10^-14 mol m^-2 Pa^-1 s^-1 at 245 K and 5 × 10^-13 mol m^-2 Pa^-1 s^-1 at 263 K. These authors and others (Wang et al. 2002; Staykova et al. 2003; Genov & Kuhs 2003) assumed that the value of the kinetics parameter of formation/dissociation of clathrates follows an Arrhenius-type function of temperature: (54)where is a constant kinetics parameter and E_a the activation energy of the process. Kuhs et al. (2006) calculated E_a and find a value (92.8 kJ mol^-1) two times greater than the ones of Staykova et al. (2003) and Genov et al. (2004) for the formation of CH₄ and CO₂ clathrates respectively. In the case of CO₂-clathrate, Genov et al. (2004) estimated the activation energies to be 5.5 kJ mol^-1 at low temperatures and 31.5 kJ mol^-1 above 220 K, indicating that water molecule mobility plays a considerable role in the clathration reaction. Furthermore, Falenty et al. (2011) find a value of 72.5 kJ mol^-1 for CO₂ and temperatures between 185 and 195 K. Thus, the physical parameters that may affect the value of are multiple: temperature, activity of the ice surface (mobility of water molecules), type of gas molecule, thickness of clathrate formed, thermal history, ... (Schmitt 1986). In the absence of more detailed knowledge of the kinetics laws at low temperature, we assume in this work that the kinetics parameter of clathrate formation is constant (no temperature dependence). We also adopt the minimum value measured by Schmitt (1986) and called hereafter the “nominal value”, (3 × 10^-13 mol m^-2 Pa^-1 s^-1). Concerning their decomposition, Schmitt (1986) has experimentally shown that the dissociation of clathrates can be much faster than their formation (up to 30 times faster). However, under some conditions a few monolayers of pure crystalline water ice can form above the clathrate structure, leading to its isolation from the gas phase in the porous network and preventing its dissociation for same time even in the presence of a gas pressure well below the equilibrium pressure (metastability). We thus assume in this model, when the gas pressure is lower than the one of MG clathrate, a nominal kinetics parameter for dissociation equal to 30 times the one adopted for the formation of clathrate. However, the influences of the values of and (and of a possible “decomposition latency time”) on the formation/dissociation of clathrates will be studied in next paper. Note that the values of the kinetics parameters used to determine the number of moles of gas trapped or released by the clathrate (Eq. (38)) are valid only in the nucleus layers where water ice is crystallized and not completely converted into clathrate, whatever their temperature. When clathrate formation occurs in layers where amorphous ice is undergoing crystallization, then the crystallisation rate of ice is used (see Eq. (40)).

3.3.6. Composition of volatile molecules both trapped and released in/by clathrate

The main difficulty with the modeling of clathrates is to accurately determine the abundance of molecules trapped inside the clathrate structure () in equilibrium with their abundance in the gas phase. For this, it is necessary to know the fractionation factors between the gas and solid phases. Many statistical thermodynamic models (Van der Waals & Platteeuw 1959; Parrish & Prausnitz 1972; Lunine & Stevenson 1985; Sun & Duan 2005; Anderson 2007; Thomas et al. 2009) predict the molecular abundances trapped in the cages of clathrates formed from a gas phase of known composition. They are mostly based on the Van der waals & Platteeuw (1959) statistical model, with some variations. However, these models depend on several parameters that are not well known yet for all molecules. Moreover, these models are fairly well constrained near the triple point of ice but were never validated in the 100−200 K temperature range because there is no experimental data to our knowledge which gives both the composition of guest molecules and that of the gas phase in this range of temperature. So we decided, as a first step, to calculate the molecular abundances trapped during the formation of clathrates by a simple law given by Eq. (42). This equation is determined from the analysis of the data from Kang et al. (2001) and Rydzy et al. (2007). It shows fair agreements (within few %) at the triple point of ice (Kang et al. 2001) and below (Rydzy et al. 2007) with both the molecular abundances in the gas phase and trapped in the clathrate. Secondly, we compared the molecular abundances obtained from this equation with those from the model of Thomas et al. (2009), which is derived from a statistical model of Van der Waals & Platteeuw (1959). The results show that the relative difference between the abundances of molecules trapped in cages calculated by Eq. (42) and the ones determined by the model of Thomas et al. (2009) is about or less than 1% for the most abundant species trapped, which is sufficient for the purposes of our calculations. However, this relative difference increases strongly when the abundance of the molecules trapped becomes low. But in such cases, larger uncertainty also remains about the results obtained by this type of model. Since the calculation of this fractionation parameter is somewhat complex and uncertain at low temperature and low pressure with such statistical models, we decided to use in this paper the law we determined, given by Eq. (42), that provides a good approximation of the main molecular abundances trapped in cages and also strongly simplifies the calculation.

An other difficulty remains with the molar fraction of volatile molecules x released (, see Eq. (43)) in the gas phase of the porous network when the clathrate structure dissociates (for total pressure P_t of the gas phase lower than the MG clathrate equilibrium pressure ). Unfortunately, there is no experimental or theoretical data that allows us to predict the behavior of MG clathrate during its dissociation in the presence of a gas phase out of equilibrium with clathrate, and hence the molar fraction of volatile molecules released. We thus assumed in this model that the molar fraction of the volatile x released by the clathrate structure is the average over the layer of the volatile x trapped in the clathrate structure (Eq. (43)). This assumption can be justified as follows: the clathrate structure is stable as long as the volatile molecules trapped remain inside. If some of these molecules are released, they will destabilize the clathrate structure and the resulting dissociation of the cages will expel the remainder of the volatile molecules, whatever the total gas pressure P_t () and the molar fraction of the volatile molecules in the gas phase of the porous network. However the assumption that the mole fraction released is simply the average of all the volatile molecules trapped in clathrate in a given nucleus layer is not exact because MG clathrate with variable composition may form by layers depending on the gas composition (if ). The resultant onion-like clathrate grains will not equilibrate between layers, or extremely slowly, due to the extremely slow diffusivity of molecules inside the clathrate structure. Ideally it would thus be necessary to manage the formation/decomposition of each clathrate layer at the grain scale inside each nucleus layer but in a first step the averaging hypothesis allows us to strongly simplify the calculation in the model.

4. Outgassing behavior of four different models of a comet

In this section, we study the outgassing profile of volatile molecules of the comet 67P/C-G, for four models made respectively fully of crystalline ice, amorphous ice, clathrate or a mixture of the three. The goal is to show that the outgassing behavior of the comet strongly depend on the initial water ice structure and that the model is able to constrain the chemical composition and the water ice structure from the analysis of the temporal profiles of gas production of 67P/C-G that will be observed by the Rosetta mission.

4.1. orbital and physical parameters adopted

The orbital and physical parameters adopted in this paper for the nucleus are the ones of the comet 67P/C-G given in Table 4. At the beginning of the computation, the cometary nucleus is assumed to have a homogeneous composition made of ices and dust. In this section, we consider four models whose initial ice phase is assumed to be composed respectively of crystalline water ice, amorphous water ice, clathrate hydrates or a mixture of the three (called hereafter respectively crystalline, amorphous, clathrate and mixed models), all mixed with pure solid CO, CO₂, CH₄ and H₂S. In this study, it is assumed that no clathrate formation/dissociation occurs in the amorphous and crystalline models. Table 4 gives the initial x/H₂O (J_x) mole fractions of the species x (x =  CO, CO₂, CH₄ or H₂S) condensed either as pure ices in the porous network or trapped in the water ice following its structure: amorphous or clathrate³. The values of J_X (sum of the three states) of all models are consistent with the observations in cometary comae of molecular species that are released directly from the nucleus (Bockelée-Morvan et al. 2004).

Whatever the model, we adopt only one nominal global initial composition of the nucleus with about 14% of CO, 5% of CO₂, 2% of CH₄ and H₂S relative to water molecules. As discussed in the Sect. 3.1, the abundance of trapped species in amorphous water ice depends on many parameters. As no reliable experimental data exists on the composition of mixed gas trapped in amorphous ice, we therefore choose, for the amorphous model, one nominal arbitrary set of plausible distribution of the volatile molecules between the states “trapped inside amorphous water ice” (up to a total of 8% in mole, Schmitt et al. 1992) and “condensed as pure ice” (thus segregated from water ice) in the porous network, considering the equilibrium pressure of species at very low temperature (hereafter 30 K) and our current knowledge on trapping processes. For the crystalline model, the volatile molecules are only in the state “condensed as pure ice”. For the clathrate model, the distribution between the states “trapped inside clathrate structure” (up to about 17% in mole) and “condensed as pure ice” come from Mousis et al. (2010). CO, CH₄ and H₂S represent respectively 79%, 10% and 11% of the molecules trapped in the clathrate structure. The “mixed” model is arbitrarily composed of 33% of each of the three water ice structures presented before.

The value of the dust/ice ratio J_dust is assumed to be equal to 1 in all models, about the value indicated for 1P/Halley by the Giotto-DIDSY measurements (McDonnell et al. 1987) and prescribed by Greenberg’s (1982) interstellar dust model (Tancredi et al. 1994). In this work, we assume a power of order −3.5 for the size distribution of dust grains (McDonnell et al. 1986; Huebner et al. 2006) with a cut-off at a radius of 1 cm (Prialnik 1997). The initial temperature is assumed to be equal to 30 K. Davidsson & Gutierrez (2005) give a density of 100−600 kg m^-3 for the comet 67P/Churyumov-Gerasimenko, leading to porosities greater than 60%. Lamy et al. (2007) estimated the mean nucleus density at  ≈ 370 kg m^-3. Note that these values should be taken with caution, since most of these estimates are based on many assumptions and are generally regarded as very uncertain (Richardson et al. 2007; Weissman et al. 2004). Here, we have chosen to test a porosity of 70% leading to a density of about 430 kg m^-3.

An important parameter that remains unknown is the thermal inertia I, i.e. the thermal conductivity K^m by using the relation (), of the porous icy matrix in comets (see Sect. 2.1). A low thermal inertia could limit the depletion of volatile molecules and the differentiation of the nucleus by thermal diffusion only to close subsurface layers. The comet would remain almost intact until the sublimation of water ice and other solid volatile molecules occurs near the surface during its perihelion passage. Conversely, a large thermal inertia could induce a great and deep differentiation of the nucleus. The temporal outgassing profile of the volatile molecules of the comet would be thereby quite different between these two extremes thermal inertia scenarios since the interfaces of sublimation of ices and the surface positions would be different. Due to the uncertainty in the thermal conductivity of the porous icy matrix, we test heat conductivities of about 1 and 0.01 W K^-1 m^-1 using Russel formula, i.e. high thermal inertia, and Hertz factor (hereafter value of 10^-2), i.e. low thermal inertia, given in Sect. 2.1 of the paper. Using the Russel formula, this leads to a thermal inertia of the nucleus of about 350 W K^-1 m^-2  for clathrate and mixed models, and about 500 W K^-1 m^-2  for crystalline and amorphous models. Using the Hertz factor, this leads to thermal inertia of 40 W K^-1 m^-2  for clathrate and mixed models, and 60 W K^-1 m^-2  for crystalline and amorphous models.

4.2. results

Figure 3 represents the thermodynamic evolution of the nucleus of 67P/C-G as a function of distance to the sun, during one revolution and after about 50 years of revolution around the sun, for the crystalline model and for the two types of thermal conductivity. The lines represent the nucleus surface and the minimum depths at which solid CO₂, H₂S, CH₄ and CO exist (hereafter called interfaces of sublimation of solid ices). Below, the gas phase of a given molecule is in equilibrium with its pure condensed phase. Above, only the gas phase exists in the porous network.

Fig. 3 Physical differentiation of 67P/C-G as a function of distance to the sun, during one revolution and after about 50 years of revolution around the sun, for crystalline model and for the two types of thermal inertia. The vertical dashed line represents the perihelion passage. Scale of the two figures is different.

The lower thermal inertia (Hertz factor used) induces a limited physical differentiation of the nucleus within only about 3 m of depth. At each perihelion passage, the ablation of the surface (about 6 m) reaches the interfaces of the volatile ices phases, leading to a strong sublimation of volatile molecules. After each passage, the comet is like new. For greater thermal inertia (Russel formula used), the depth of the physical differentiation is about 80 m. The ablation of the surface, about 2 m per revolution, never reaches the interfaces of the volatile ices in the interior of the comet. So, with such a high thermal inertia, the interfaces of sublimation of the volatile ices are only slightly affected by the ablation of the surface of the nucleus. The nucleus remains then differentiated in subsurface layers.

Figure 4 represents the outgassing profiles of CO, CO₂, CH₄ and H₂S, in molecules per second and per unit of surface, of the nucleus as a function of distance to the sun, during one revolution after about 50 years of revolution around the sun, for the four models and for the two types of thermal conductivity. As the water ice production doesn’t change as a function of the initial state of the water ice, we plot its curve only for the crystalline model with high thermal inertia.

Let’s start the description by the models with a high thermal conductivity given by Russell’s formula. For the crystalline model, the maximum productions in the outgassing profiles of CO₂ and H₂S are shifted compared to perihelion (at about 1.3 AU). The outgassing profiles of CO and CH₄ remain approximately constant whatever the position of the comet around the sun. For all the molecules, their outgassing comes from the sublimation of their pure ice phase. The deeper the interface of sublimation of a volatile molecule inside the nucleus is, more shifted relative to perihelion its outgassing peak is.

For the amorphous model, the crystallization of amorphous H₂O ice only slightly changes the production of CO and CH₄ during the perihelion passage. The main difference with the crystalline model, i.e. the global decrease in gas production, especially of CO, comes from the distribution of the volatile molecules between the “trapped” and “condensed” states in the icy matrix.

For the clathrate model, the outgassing profiles of all molecules (except CO₂ since this molecule is not trapped in the clathrate structure at low temparature) are changed compared to the others models. The maximum of outgassing of CO, CH₄ and H₂S occurs at perihelion passage, with some strong fluctuations of the outgassing (hereafter spikes). The spikes represent the release, close to the surface, of these molecules initially trapped in the clathrate structure. They are related to the day/night changes of insolation, and hence to sublimation of H₂O, at the surface of the nucleus.

For the mixed model, the outgassing profiles of molecules are as the ones of the crystalline and clathrate models together, with a thermal inertia close to this last one. As in the clathrate model, the maximum outgassings of CO, CH₄ and H₂S molecules occur at perihelion passage.

For models using a low thermal inertia determined by the Hertz factor, the outgassing profiles of the volatile molecules change compared to the above models with high thermal inertia. In most cases, the outgassing of molecules increases strongly just before and during perihelion passage since the ablation of the surface reaches the interfaces of sublimation of the volatile ices. The global production of all molecules is increased compared to models with greater thermal conductivity. For the crystalline model, the maximum of outgassing of CO₂ and H₂S occurs at perihelion passage. For these two molecules, which have their interface close to the surface, their outgassing fluctuates in phase with H₂O due to the day/night insolation variation. The maximum of the outgassing profiles of CO and CH₄ are shifted compared to the ones of CO₂ and H₂S since the depth of their interface of sublimation are deeper in the nucleus (see Fig. 3).

Fig. 4 Outgassing profiles, in molecule per second and per unit of surface, of 67P/C-G as a function of the distance to the sun, during one revolution and after about 50 years of revolution around the sun, for the four models and for the two types of thermal inertia (high on the left and low on the right). The vertical dashed line represents the perihelion passage.

For the amorphous model, the crystallization of amorphous H₂O ice fully changes the production of all the molecules during the perihelion passage. The outgassing of molecules increases significantly during perihelion passage since the ablation of the surface reaches both the interfaces of crystallization of amorphous water ice and of sublimation of volatile ices. For all the molecules, the outgassing fluctuates with the day/night insolation variation, but with a variable amplitude reflecting their abundance close to the surface.

For the mixed and clathrate models, the outgassing of all molecules (including CO₂) is concentrated during the perihelion passage with fluctuations due to the day/night insolation variation. The differences between the outgassing profiles of these two models come from the different distributions of molecules between the “trapped” and “condensed” states in the nucleus.

To summarize, the outgassing profiles of the volatile molecules is function of the distribution of these molecules between the “trapped” and “condensed” states, the nature of the trapping (amorphous ice or clathrate structure), the thermal conductivity of the icy matrix and hence the depth of the interfaces of phase change of the volatile ices. Knowing the outgassing profiles of several molecules from a comet nucleus, we expect that this model, by inversion, should be able to determine its water ice structure and the abundances of volatile molecules between the “trapped” and “condensed” states. In the case of the 67P/C-G comet, the target of the Rosetta mission, which will be monitored all along its orbit such an inversion of the actual and initial states and compositions of the nucleus is very promising.

5. Conclusion

The chemical composition and water ice structure of the cometary nuclei remain currently unknown. For a better insight into these witnesses of the origins of comets, we have developed a model that takes into account all water ice structures (amorphous, crystalline, clathrate hydrate or any mixture of these structures of ice) and all their phase changes (amorphous water ice to pure crystalline water ice, amorphous water ice to clathrate hydrate, pure crystalline water ice to clathrate hydrate and vice-versa), as well as sublimation/condensation of volatile molecules in the porous network. This model describes the heat transfer, latent heat exchanges, gas diffusion, and the gas releases and trapping by crystallization and clathrate formation/dissociation processes. Taking into account all these physico-chemical processes, this model has been applied to the comet 67P/Churyumov-Gerasimenko with four initial different ices structures: crystalline ice, amorphous ice, clathrate and a mixture of the three. Results showed that the outgassing profile of volatile molecules from the nucleus is mainly function of the structure of water ice, the distribution of the volatile molecules between the “trapped” and “condensed” states in the icy matrix, and the thermal inertia of the icy matrix. We expect that this model will be able to constrain the chemical composition and the water ice structure in cometary nuclei from the outgassing profiles of volatile molecules observed during comet activity, and especially those of the 67P/Churyumov-Gerasimenko, the target comet of the Rosetta mission.

In the process of building this model, we did a complete review of all physical processes and thermodynamic properties of ices below the freezing point of water. We have identified several gaps to be filled by experimentals studies, essentially linked with the lacks of data in the thermodynamic properties of ices and clathrates, the pressure equilibrium conditions of clathrates, and the trapping and release processes of volatile molecules in the water ice structures (amorphous water ice and clathrate hydrate). these limitations are one of the reasons why we have considered only four volatile molecules (CO, CO₂, CH₄ and H₂S) in addition to H₂O although more than twenty volatile molecules have been identified in comets (Bockelée-Morvan et al. 2004). Since all these four molecules create single guest clathrates of structure I, only this structure of clathrate hydrate has been taken into account in the model. The presence of molecules in comets such as C₂H₂, C₂H₆ (observed in comets) or as Ar, Kr, N₂,... (not yet observed) could lead to the formation of structure II or a mixture of structure I and II as shown in Rydzy et al. (2007). In these conditions, the formation and decomposition of clathrate hydrate during the physico-chemical evolution of comets could then be more complex in terms of energy and trapping/release of volatile molecules than the simpler case presented in this model. Moreover, the main difficulty, taking into account all phases changes in the model, remains the correct estimate of the abundances of volatile molecules trapped and released by amorphous water ice and clathrate hydrate, at and far from equilibrium pressure conditions, as well as the fractionation between the gas phase and the clathrate compositions. The trapping of volatile molecules during amorphous water ice condensation and clathrate hydrate structure formation, and their release during the amorphous-to-crystalline water ice transition and the dissociation of the clathrate structure, out of equilibrium pressure conditions, should be more deeply studied experimentally. This will allow to better constrain the models of cometary nuclei which currently adopt arbitrary initial abundances of volatile molecules and poorly constrained processes for their release in the porous network. Currently, the remaining uncertainties about the behavior of water ice and the trapping/release processes of molecules inside the nucleus are problematic, without experimental and theoretical data, for the interpretation of the outgassing profiles of the volatile molecules that will be observed by the Rosetta spacecraft. In order to assess the uncertainties about the behavior of the different water ice structures, especially during the various phase transitions, this model will be submitted to a sensitivity study on the influences of the physical parameters that are poorly constrained. We will devote a particular attention to issues of the kinetic rates of formation/dissociation of clathrates and of the degree of release, or not (metastability), of the volatile molecules trapped in water ice (amorphous or clathrate hydrate). This study will be presented separately in a next paper whose the objective is to determine whether the formation of clathrate is physically viable, thermodynamically and kinetically in the 67P/Churyumov-Gerasimenko comet without clathrate structure at the origin.

Acknowledgments

This work has been supported by the French Centre National d’Études Spatiales (CNES) and by the University Joseph Fourier – Grenoble 1 through post-doctoral fellowships (U.M.). All the computations presented in this paper were performed at the Service Commun de Calcul Intensif de l’Observatoire de Grenoble (SCCI) and using the UTINAM Cluster thanks to the grant of Franche-Comté region.

References

Online material

Appendix A: The numerical scheme

In this section, we describe the numerical comet nucleus model used to solve the heat and mass diffusion equations and developed originally by Marboeuf et al. (2008). In order to ensure perfect conservation of mass and energy in the model, we use finite volume method (Patankar 1980) for the discretisation of Eqs. (1) and (22). For illustrative purpose, we detail below the space and time integration of Eq. (1) which is the most complex equation.

Our numerical model of comet nucleus is a one-dimensional grid with N_C layers (see Fig. A.1). We define r_i as the position of the top of each layer from r₁ = R_nucleus at the surface to r_Nc + 1 at the center of the nucleus.

The layers are spherical shells and the temperature and other state variables are assumed constant inside each layer, a common approximation in discretization schemes.

We note Δr_i the thickness of each layer i: (A.1)The initial thickness of the layers increases with depth following a geometric progression Δr_i = aΔr_i − 1, where a is a constant factor slightly larger than 1.

The variable L in Fig. A.1 represents either the temperature T, the pressure P, the coefficient of diffusion of mass G or the coefficient of diffusion of heat K^m that we calculate at the center of each layer: this corresponds to the physical conditions in the center of each cell.

Fig. A.1 Schematic view of the numerical nucleus. Pressure Pi and temperature Ti are defined at the center of each cell i and are represented by the coefficient Li. The coefficients of heat conduction and gas diffusion Gi are defined at the center and at the edge of each cell i. They are here represented by the coefficients Li at the center and li at the edge of each cell i. The parameters ri and Δri are the distance from the center of the nucleus and the thickness of each cell i.

When a chemical species begins to sublimate, i.e. the temperature is such that the vapor saturation pressure is equal to the partial gas pressure of the species x in the pores, the sublimation front of the species divides the nucleus into two zones. In the first one, from the surface down to the interface, species x is only present in the gas form and simply obeys Eq. (22) without the source/well terms. In the second zone, from the interface down to the center of the nucleus, the gas x is in equilibrium with its ice, its pressure is the vapor saturation pressure, a function of temperature, and all the other variables are adjusted by the model to account for gas flow due to gradients.

The two Eqs. (1) and (22) are solved in this structure following the method of Orosei et al. (1999). In order to increase the numerical stability, both the fluxes of heat and mass are solved together in the resolution of Eq. (1) for the layers where chemical species are condensed. For a more realist treatment of the gas diffusion through the porous matrix between the interface and the surface, each chemical species can diffuse following its own time step of diffusion.

When the total time of gas diffusion t_tot, between the interface of species x and the surface, is smaller than f_tg Δt (f_tg is an accuracy parameter that we impose and Δt is the time step for integration of Eq. (1)), the gas diffusion of species x is assumed to be in a steady state and the right term of the mass conservation Eq. (22) is taken equal to zero. In the following simulations, we use f_tg = 0.1.

A.1. Spatial integration

The finite volume method, described first by Patankar (1980), consists in integrating the equations of conservation (1) and (22) on each volume of control. The main advantage of this method is that the divergence term present in both equations of conservation is integrated analytically, leading to a simple algebraic difference equation of terms at the surface of the control volume, which can be easily solved numerically.

Hereafter, we develop successively the different terms of the spatial integration of the conservation energy Eq. (1), which is the most interesting equation:

A.1.1. Spatial integration of the first terms of Eq. (1)

Term of the left hand side of Eq. (1). Performing the spatial integration of the left hand side term of Eq. (1) over the control volume yields: (A.2)with and i the index of the layer. Here we make use of the underlying assumption in discretization schemes that the temperature in constant inside a cell. One can also consider the temperature used in the model as an average of the real temperature over the cell.

First term of the right hand side of Eq. (1). In the following, we will always use spherical coordinates, assuming spherical symmetry, so that only the radial component of the equations remains.

In order that two cells see the same flux on both sides of the interface of the volumes, we choose to express the coefficients of diffusion G_i and conduction at the surface of each volume of control, i.e. each radial cell, and denote g_i and their value at the limit between cells i − 1 and i. Figure A.1 shows the variable L_i that represents the coefficients of diffusion G_i and conduction at the center of the cell i and its corresponding variable l_i on each interface for the coefficients of diffusion g_i and conduction . In this way, each layer possesses a set of coefficients of diffusion (G,K^m) defined in the center of the cell and an other set of coefficients (g,k^m) defined at its surface. The coefficients l_i can be written as functions of the coefficients L_i − 1 and L_i located above and below the interface of the layers. We consider that theses coefficients form a parallel network defined by equation: (A.3)The result of replacing the coefficients L_i by l_i in the divergence terms of Eqs. (1) and (22) is the exact conservation of the flux of mass and energy through the interface. Note that the interpolation in divergence terms of the coefficients at the boundary of the cells must be applied to all, but only, the coefficients of the gradient terms (see below Eqs. (A.9) and (A.10)).

The spatial integration of the divergence becomes then: (A.4)where the gradient of temperature is written as: (A.5)

Second term of the right hand side of Eq. (1). The term Q_x only exists in Eqs. (1) and (22) if the chemical species x is condensed in the pores. It is derived from Eq. (22) and the integration of this term becomes: (A.6)with (mol m^-3 s^-1). Here we have used in place of P_x because the species x exists in condensed form and hence the gas is in thermodynamic equilibrium at pressure .

A.1.2. Spatial integration of the Eq. (A.6) (source/well term )

We expand below the different terms of the spatial integration of Eq. (A.6):

First term of the right hand side of Eq. (A.6): (A.7)where T is assumed to be constant during the gas pressure variation. Finally, because is only a function of T, we can write: (A.8)Second term of the right hand side of Eq. (A.6): since we know analytically, we express this gradient in terms of the temperature gradient . Similarly to what we did for the heat and gas diffusion coefficients, using Eq. (A.3), we interpolate the quantity , that appears in the divergence term, and which is known only at the centers of the control volumes, to its value θ_i,x at the edge of the cells, in order to ensure proper conservation.

The radial component of the integration gives then for each volatile molecule:(A.9)\arraycolsep1.75ptwhere is the gradient of vapor saturation pressure of species x.

The spatial integration over the control volume of the remaining terms Y^cr, Y^chs , and Y^cl is straightforward. from Eqs. (1) and (A.6) vanishes.

A.1.3. Final spatial integration of the energy Eq. (1)

Finally, as a result of the spatial integration of Eq. (1), we obtain the conservation of the energy at a given time step: (A.10)with (A.11)and (A.12)where θ_i,x and k_i are defined at the interface of the layers i and i − 1, as explained above.

A.2. Temporal integration

Fig. A.2 Reconstruction of a cell after removal of the interface. The new temperature Tj and pressure Pj of each cell j are calculated from the previous values Ti and Pi by interpolation. and are respectively the volume fractions of cell i and i + 1 that give the volume Vj of the new cell j.

Equations (1) and (22) are nonlinear and the best method to solve them would be to use an iterative fully implicit method for each spatial step in order to find the exact value of T and P for each time step. However, such a method is very time consuming. We choose then to follow Espinasse et al. (1991) in using the predictor-corrector method and imposing a slow evolution with time. The different choices of integration methods for both the space and time derivative equations of conservation, in the finite difference framework, have been already discussed in Huebner et al. (2006) and Prialnik et al. (2004). Here we recall only the principles of the integration scheme and give more details. The idea of the predictor-corrector scheme is to first estimate the values of T and P at the middle of the time step Δt for each cell and time step (predictor). The knowledge of P^Δt/2 and T^Δt/2 permits us to calculate the values of the non-linear coefficients (G,g), (K,k) and () at the middle of the current time step. These last values are then used to solve Eqs. (1) and (22) with the current time step Δt at the corrector. Now, we give, as an example, the time integration of Eq. (A.10) on a time step nΔt, with for the predictor and 1 for the corrector: (A.13)To solve this equation, we need to know the temporal evolution of the temperature gradient on the time step nΔt. To this end, we replace it by a linear combination of the known and unknown temperature gradients balanced by a parameter f that takes values between 0 and 1 according to the type of resolution scheme we choose: (A.14)As a result of the time integration, we obtain the final equation: (A.15)with m = 0 for the predictor and for the corrector.

For the predictor we choose to take a fully implicit scheme (f = 1), while for the corrector, we use a Cranck-Nicolson scheme which is a semi-implicit method (f = 1/2). Both numerical schemes are totally mathematically stable. The predictor scheme is accurate at the first level order in time step Δt and predicts a good value of the variables as long as the time step is short⁴ (Huebner et al. 2006; Prialnik et al. 2004). The corrector is accurate to second order in Δt. This last scheme is restricted by a boundary condition on the time step if anyone wants fair physical solutions (Patankar 1980): . Hence, mixing the fully implicit and the semi-implicit schemes improve the stability and accuracy of the overall scheme.

After expanding each term in Eq. (A.15), we obtain a tridiagonal matrix which takes the form: (A.16)For each layer i we can solve Eq. (A.16) with a tri-diagonal matrix algorithm (TDMA, Patankar 1980). Following the same scheme, Eq. (22) for species x is transformed into a tridiagonal system, with pressure P_x in place of T.

Remember that the sublimation interface of species x separates the nucleus into two different regions. Below the interface, the full Eq. (1) is solved with all source terms, while the pressure in Eq. (22) is the vapor saturation pressure , a function of T. The source term Q_x in Eq. (22) is then explicitly computed from T and and its gradient. Above the interface, the source term Q_x in Eq. (22) vanishes and hence the full Eq. (A.6) disappears from Eq. (A.15). Then Eq. (22) is solved between the surface and the sublimation interface of species x, while Eq. (1) is solved from the surface to the center of the nucleus, with the change of terms mentioned above at the sublimation interface.

To apply the TDMA, we need to impose boundary conditions at each end of the integration space, both for Eq. (1) and Eq. (22). At the surface, the temperature is given by Eq. (18) (Dirichlet boundary condition). At the center of the nucleus, we impose the flux of heat equal to zero (Neumann boundary condition): (A.17)If species x is present in condensed form everywhere in the nucleus, then we never solve Eq. (22), and even at the surface. If, on the contrary, there exists a sublimation interface for species x, then Eq. (22) is solved between the surface of the nucleus and the sublimation interface. Here, we impose a Dirichlet boundary condition at both ends. At the surface, we impose

and at the sublimation interface, we have

At the end of the corrector, when all the variables are recalculated, the interface front of sublimation of each specie is moved towards the center of the nucleus. In order to avoid instabilities during the computation, each interface is moved continuously towards the center following the prescription of Orosei et al. (1999). The size and the number of cells are then recalculated at each time step (see Fig. A.2). We do this so that the size of the cells do not change too much with time at a given location, therefore avoiding abrupt changes in the gradients of P_x and T and ensuring stability of the computation. In order to ensure the exact conservation of energy and mass in the nucleus, the new temperature T_j and pressure P_j of each new cell j are calculated from the previous values by interpolation.

For illustration, let us consider a new cell j that overlaps with previous cells i and i + 1 (see Fig. A.2). The new temperature T_j is computed so that the internal energy of the cell does not change, i.e. the energy change of the fraction of mass from cell i changing from T_i to T_j is exactly the opposite of that due to the fraction of mass from cell i + 1 changing from T_i + 1 to T_j, throught their respective heat capacities following the equation: (A.18)where and are respectively the volume fractions of cell i and i + 1 that give the volume V_j of the new cell j.

Similarly the total mass of each constituent (dust, water, molecules in gas and condensed forms) is strictly conserved by recomputing the condensed mass density and keeping either the same amount of gas when there is no condensed molecule or forcing the pressure to the vapor saturation pressure at the temperature of the cell.

Appendix B: Symbols used in the paper

In this section, we present tables with principal symbols used in the model (Tables B.1 and B.2).

All Tables

All Figures

	Fig. 1 Schematic view of the interior of a homogeneous comet nucleus composed of icy grains. The matrix of the nucleus can be formed of icy grains with amorphous, crystalline water ice structures, clathrate hydrates or a mixture of these icy grains and water ice structures.
In the text

	Fig. 2 Schematic view resuming the kinetic law adopted and the power released or taken by water ice structure during its transition.
In the text

	Fig. 3 Physical differentiation of 67P/C-G as a function of distance to the sun, during one revolution and after about 50 years of revolution around the sun, for crystalline model and for the two types of thermal inertia. The vertical dashed line represents the perihelion passage. Scale of the two figures is different.
In the text

	Fig. 4 Outgassing profiles, in molecule per second and per unit of surface, of 67P/C-G as a function of the distance to the sun, during one revolution and after about 50 years of revolution around the sun, for the four models and for the two types of thermal inertia (high on the left and low on the right). The vertical dashed line represents the perihelion passage.
In the text

Fig. A.1 Schematic view of the numerical nucleus. Pressure Pi and temperature Ti are defined at the center of each cell i and are represented by the coefficient Li. The coefficients of heat conduction and gas diffusion Gi are defined at the center and at the edge of each cell i. They are here represented by the coefficients Li at the center and li at the edge of each cell i. The parameters ri and Δri are the distance from the center of the nucleus and the thickness of each cell i.

In the text

	Fig. A.2 Reconstruction of a cell after removal of the interface. The new temperature Tj and pressure Pj of each cell j are calculated from the previous values Ti and Pi by interpolation. and are respectively the volume fractions of cell i and i + 1 that give the volume Vj of the new cell j.
In the text