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.

Open with DEXTER

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.
Open with DEXTER

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).

© ESO, 2012