Issue 
A&A
Volume 542, June 2012



Article Number  A82  
Number of page(s)  24  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201118176  
Published online  13 June 2012 
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 onedimensional grid with N_{C} layers (see Fig. A.1). We define r_{i} as the position of the top of each layer from r_{1} = 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 P_{i} and temperature T_{i} are defined at the center of each cell i and are represented by the coefficient L_{i}. The coefficients of heat conduction and gas diffusion G_{i} are defined at the center and at the edge of each cell i. They are here represented by the coefficients L_{i} at the center and l_{i} at the edge of each cell i. The parameters r_{i} and Δr_{i} 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 T_{j} and pressure P_{j} of each cell j are calculated from the previous values T_{i} and P_{i} by interpolation. and are respectively the volume fractions of cell i and i + 1 that give the volume V_{j} 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 predictorcorrector 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 predictorcorrector 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 nonlinear 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 CranckNicolson scheme which is a semiimplicit 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^{4} (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 semiimplicit 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 tridiagonal 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.
Greek symbols.
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).
Latin symbols.
© ESO, 2012
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.