%%==================================================================%% %% A S T R O N O M Y AND A S T R O P H Y S I C S %% %% %% %% %% Plain TeX Support Version 3.0 %% %% %% %% This is file P-AA.TEX the startup file of the macro package. %% %% Use PP-AA.FMT or CP-AA.FMT for formatting. %% %% %% %%==================================================================%% % % \refereelayout % a referee layout % \MAINTITLE={ Thermochemistry of cometary nuclei } % must be given % % % \FOOTNOTE{ ????? } \FOOTNOTE is valid only in % \MAINTITLE or in \SUBTITLE % \SUBTITLE={I. The Jupiter family case} % is optional % \AUTHOR={ G. Tancredi @1 @2, H. Rickman @1, J. M. Greenberg @3} % Firstname Name<@number><, Nextauthor@number> % items in "< >" are optional % \PRESADD{Present address} is optional % \INSTITUTE={@1 Astronomical Observatory, Box 515, {\sl S}-75120, Uppsala, Sweden @2 Dpto. Astronom\'{\i}a, Fac. Ciencias, Trist\'an Narvaja 1674, Montevideo, Uruguay @3 Laboratory Astrophysics Leiden, Box 9513, 2300 RA Leiden, The Netherlands} % Take care of the affiliations used with \AUTHOR % @1 Institute No.1<@2 next institute> % \ABSTRACT={ New experimental results related to the physical characteristics of the material components of the cometary nuclei as well as new ideas about several aspects of the modelling of the thermochemical process in the interior of this objects lead us to make a new attempt to analyse the physical evolution of Jupier family comets over time scales comparable to their lifetime. A new model is described in this paper where we present results concerning the evolution of Jupiter family comets and make comparisons with previous models. Our model of the cometary material includes a porous solid matrix and vapour filling the pores. As basic constituents of the solid matrix we consider three omnipresent species: dust, water ice (in two phases: amorphous and crystalline), and H$_2$O vapour. In addition to the above, we include one substance more volatile than H$_2$O, CO, initially trapped in the amorphous matrix. We improved on the earlier models by accounting for the state of near saturation attained by the vapour inside the nucleus, by including a separate treatment of an unsaturated surface layer and by explicitly including the erosional velocity of the surface. As far as physical parameters are concerned, our basic improvements on earlier models were: 1) the representation of this matrix as an aggregate of micron-sized core-mantle grains; 2) the adoption of a very low thermal conductivity of the amorphous ice mantles; and 3) a correct account of the energetics of gas release and the allowance for condensation of CO ice. We defined a ``standard'' nuclear model with the best guesses of the many unknown or poorly known parameters, and we ran it for 500 years in a typical Jupiter family orbit ($q=1.5$ AU, $Q=6$ AU), a time comparable to $\sim$ 10 \% of the their lifetime. The CO bursts (associated with crystallization spurts) are notorious in the first revolutions, but then, they gradually evolve from the sharp spikes to a much more subdued appearance. A set of variant models were run to explore the consequences of some of our assumptions. Variations of the following parameters were considered: dust to ice ratio, porosity and amorphous ice conductivity. We note the broad similarity between our standard and variant models. The ``standard'' model was also run in a capture scenario, where the comets first stay n a high-$q$ orbits and then into a low-$q$ one. For high-$q$ orbits , the rate of CO outgassing exceeds the perihelion H$_2$O outgassing rate by several orders of magnitude. Upon capture, the comet basically behaves in accordance with the burial depth of the crystallization zone independent of which previous orbital evolution has led to this state. Concluding on the behaviour of Jupiter family comets, we find that the complete crystallization of a sizeable nucleus with an initial radius of several km should take $\sim 10^4$ years. This means that Jupiter family comets with our assumed properties should still retain their CO, although in most cases buried deep below the nuclear surface. } % must be given % \KEYWORDS={ comets -- nuclei -- physical processes} % must be given % \THESAURUS={07(07.03.1; 02.04.2)} % must be given % \OFFPRINTS={ ????? } % is optional % \DATE={ ????? } % must be given at the last moment % % Please don't overrule this command, it will format your titlepage. \maketitle % % Now start your text using the Springer macros. See the % instructions in file P-AA.DOC or the examples in P-AA.DEM % \titlea{Introduction} \def\lta{~\raise.4ex\hbox{$<$}\llap{\lower.6ex\hbox{$\sim$}}~} The concept of cometary nuclei as porous objects initially containing an important fraction of gas-laden amorphous ice has recently gained wide acceptance (see Rickman 1991). The latent heat release of crystallization and the associated release of trapped gases have long been recognized as important processes influencing the evolution of cometary nuclei, secularly as well as on short time scales (Patashnick {\it et al.} 1974; Smoluchowski 1981; Bar-Nun {\it et al.} 1985; Prialnik and Bar-Nun 1987, 1988, 1990; Espinasse {\it et al.} 1989, 1991). However, numerical modelling of these processes over long periods of time and exploration of the unconstrained ranges of inaccurately known yet influential parameters still has not advanced very far. An obvious reason is that the physical scenario in question is a complicated one, expressed by a set of coupled, nonlinear partial differential equations where strongly variable time scales and nasty feedbacks are involved. The first long-term results for a comet in the orbit of P/Halley were published by Prialnik and Bar-Nun (1987, 1988). This model only included crystallization involving release of latent heat, supposed to occur instantaneously at a critical temperature $T_{cr}=137$ K. Indications were that the advance of the crystallization front toward the interior of the nucleus might occur in discontinuous spurts covering a typical depth range of 10--20 m, due to the self-sustaining nature of the exothermic process. The bursts would be halted as the front reaches the very cold amorphous ice assumed to characterize the whole nucleus initially -- the thermal imprint of the initial process and of the billions of years spent in the Oort cloud. After a long time of gradual heating from the surface, and the gradual approach of the surface to the phase transition level by the erosion caused by the sublimation of ice, the conditions for a new crystallization spurt would be met, and the cycle would repeat. It was also found, however, that the inclusion of a dust component of the nuclear material material, possibly forming an inert mantle on the surface, or the assumption of a large porosity would change the crystallization behaviour into a more continuous process, whereby the front would remain at a depth of $\sim 15$ m below the surface. Experimental results from the KOSI project (Kochan {\it et al.} 1989) indicated that at high enough temperatures ($T \ga 210$ K) the contribution of vapour flow with sublimation and recondensation in the interior of the nucleus becomes the dominant mode of energy transfer, as suggested by Smoluchowski (1982). Modelling of the crystallization behaviour of a cometary nucleus containing trapped gases, including the diffusion of the released vapours, was thus started by Espinasse and co-workers. Such was the situation reviewed by Rickman (1991) concerning the thermochemical evolution of comets. The models by Espinasse (1989) and Espinasse {\it et al.} (1989, 1991) dealt with pure ice nuclei moving in the orbits of P/Halley or P/Churyumov-Gerasimenko, containing 5 or 10~\% trapped gas (either CO or CO$_2$), and having a porosity in the range from 0.3 to 0.8. Crystallization was now modelled as a continuous process following an activation law as a function of temperature. These models were run for a few revolutions only, but for most combinations of parameters the crystallization front in the P/Churyumov-Gerasimenko case advanced so rapidly that a quick runaway was expected, whereby the whole nucleus would crystallize in a period of several centuries only. This result had the interesting implication that one would tentatively expect Jupiter family comets in general to lose their volatiles (CO or CO$_2$) shortly after capture, and thus it seemed likely that only those that have recently been captured for the first time from the trans-jovian reservoir would exhibit outgassing of such volatiles (cf. Rickman and Tancredi 1993). Another interesting consequence of crystallization with CO release turned out to be that bursts of CO outgassing might occur at practically any orbital position in the P/Churyumov-Gerasimenko case and far away from perihelion (after perihelion passage) in the P/Halley case. This was tentatively offered (Schmitt {\it et al.} 1991) as an explanation for the post-perihelion outburst of this comet at 14 AU heliocentric distance (West {\it et al.} 1991). Similar indications also came from the updated model of Prialnik and Bar-Nun (1990) including gas release and diffusion as well as continuous crystallization, i.e., that major bursts of outgassing would occur due to crystallization spurts at heliocentric distances well beyond the activity limit expected on the basis of H$_2$O sublimation. The question remained, however, what long-term evolution would ensue and whether the outburst behaviour would persist even after the $\ga 150$ revolutions spent by P/Halley in an orbit similar to the present one (Carusi {\it et al.} 1987). For Jupiter family comets, moreover, one should take into account the instability of the orbits with frequent jumps between different ranges of perihelion distance, since such sudden changes of insolation might influence the secular evolution of crystallization and CO outgassing in important ways. A new inducement to revisit the thermochemical modelling came with the experimental finding of extremely low thermal conductivity of amorphous ice by Kouchi {\it et al.} (1992). A decrease by several orders of magnitude with respect to the standard value used in all earlier models led to the expectation that quite different crystallization behaviour might ensue from using the lower conductivity. For instance, the long-term radiogenic heating of cometary nuclei had earlier been found to lead to quite modest temperatures (e.g., Yabushita and Wada 1988), but Haruyama {\it et al.} (1993) found with the conductivity in accordance with Kouchi {\it et al.}'s results that such internal heating could even lead to essentially complete crystallization. More work is needed in order to explore such scenarios, but the presence of CO outgassing in P/Halley and other comets indicates that at least some comets have escaped this fate. Meanwhile it is of great interest to explore the long-term crystallization behaviour of Jupiter family comets with special regard to the effect of lingering uncertainties over a number of physical parameters. For this purpose we need a simple, adaptable theory and an efficient numerical treatment, allowing the running of many models over extended periods of time. In the present paper we describe such a theory and the principles of our numerical procedure, and after checking our results by comparisons with earlier work, we perform a first set of simulations for a ``standard'' model and a number of variant models. We highlight the ways in which the behaviour varies with the choice of parameters and draw conclusions regarding features that appear robust. The model is described in Sect. 2, the results in Sect. 3, and the discussion and conclusions in Sect. 4. \titlea{The Model} \titleb{Components of the material} Our model of the cometary material includes a porous solid matrix and vapour filling the pores. The solid phase is not mobile, but the vapour flows down the pressure gradient. As basic constituents of the solid matrix we consider a ``dust'' component, made up of non-volatile substances, and a water ice component. The latter exhibits two phases: amorphous and crystalline. Thus our model does not include any restructuring of amorphous ice (i.e., annealing), nor does it take the cubic-hexagonal transition of crystalline ice into account. The amorphous and crystalline phases are considered to coexist, the finite transition rate being a continuous function of temperature. The dust and ice components are omnipresent in our model, i.e., we do not consider any ``dust mantle'' in the sense of a surface layer of finite thickness devoid of ice. A third omnipresent species is in principle H$_2$O vapour, although its saturated pressure at very low temperatures, say below 100 K, can be considered negligible. In addition to the above, we include one substance more volatile than H$_2$O, which may occur in three different forms. First, the amorphous ice holds a certain relative amount of the volatile as ``trapped gas''; secondly, upon crystallization, this is released into the pores as vapour; and thirdly, at low enough temperature, the vapour condenses into an additional ice component of the solid matrix. Our modelling procedure allows a free choice of this volatile substance, but in this paper we only consider CO. The basic physical parameters describing the material are the densities, i.e., the mass of each particular component contained in a unit volume of the cometary nucleus. Thus the dust density is $\rho_d$, the water ice density is $\rho_w^i$ and the CO ice density is $\rho_{CO}^i$. Since these are the three components of the solid matrix, the total density of this matrix is: $\rho_s = \rho_d + \rho_w^i + \rho_{CO}^i$. The porosity $p$ (fraction of the unit volume not occupied by solid material) is found, relating the above densities to the corresponding densities of compact materials ($\rho_{d,c}$, $\rho_{w,c}^i$, $\rho_{CO,c}^i$), by $$ p = 1 - { \rho_d \over \rho_{d,c}} - { \rho_w^i \over \rho_{w,c}^i} - { \rho_{CO}^i \over \rho_{CO,c}^i} \eqno(1)$$ Splitting the water ice density into different components, we have $$ \rho_w^i = \rho_{H_2O}^c + \rho_{H_2O}^a + \rho_{CO}^t \eqno(2)$$ where $\rho_{H_2O}^c$ and $\rho_{H_2O}^a$ are the densities of crystalline and amorphous H$_2$O ice, respectively, and $\rho_{CO}^t$ is the density of trapped CO forming part of the amorphous component. For simplicity we assume that condensation of water vapour in the interior of the nucleus leads to crystalline ice at all temperatures -- this should not imply any large errors, since the sublimation rate is extremely low at all temperatures where amorphous ice is stable (see below). We also assume the amount of trapped CO to be unaffected by the possible flow of CO vapour through the amorphous ice. This means that $\rho_{CO}^t$ is strictly related to $\rho_{H_2O}^a$ via $$ \rho_{CO}^t = f_t \ \rho_{H_2O}^a \eqno(3) $$ where $f_t$ is a constant factor to be chosen as a model parameter. In order to keep track of both water sublimation and crystallization we introduce the solid H$_2$O density $$ \rho_{H_2O}^i = \rho_{H_2O}^c + \rho_{H_2O}^a \eqno(4)$$ The dust/ice ratio of the solid matrix $$R_{di} = {\rho_d \over \rho_w^i + \rho_{CO}^i} \eqno(5)$$ is initially given a constant value throughout the nucleus, but its local values change with time due to sublimation, condensation and crystallization. Most of the above parameters vary with depth in the nucleus at any time. In particular, a critical issue is whether any CO ice is present or not. Therefore we introduce a step function $\delta_{CO}$ with the property $$ \delta_{CO} = \cases{1 &if $\rho_{CO}^i > 0$\cr 0 &if $\rho_{CO}^i = 0$\cr} \eqno(6)$$ The vapour densities, $\rho_w^v$ and $\rho_{CO}^v$, can be considered as either saturated or unsaturated. In the former case their values, $\rho_w^s(T)$ and $\rho_{CO}^s(T)$, are uniquely related to the temperature $T$. In the latter case, if ice is present, sublimation takes place tending to increase the vapour density. \titleb{Sublimation and vapour flow} According to kinetic theory, the sublimation rate is related to the rate of collisions between gas molecules and the surface of the icy matrix. Considering a plane surface where the adjacent gas has number density $n$ and is thermalized with an isotropic velocity distribution and mean speed $v_{th} = \sqrt{8 R_g T / \pi \mu}$ ($R_g = $ gas constant; $\mu = $ molar mass), the number of molecules impinging on unit surface area per unit time is $n v_{th} / 4$. We thus consider a ``mobility'' parameter $$ Q_x = \rho_x^v \ \sqrt{{R_g T \over 2 \pi \mu_x}} \eqno(7)$$ where the subscript $x$ denotes the species in question ($w$ or $CO$), such that the mass rate of collisions with the matrix is $S Q_x$, if we let $S$ denote the specific surface (surface area per unit volume) of the matrix. For a saturated vapour, by analogy with the above, we have $Q_x = Q_x^s(T)$, replacing $\rho_x^v$ in (7) by $\rho_x^s(T)$. Assuming, as usual, a sticking efficiency $s=1$, i.e., that every molecule colliding with the matrix sticks to it, $S Q_x$ is the mass condensation rate per unit volume and $S Q_x^s(T)$ is the mass evaporation rate. Thus the net mass sublimation rate, i.e., vapour mass production rate per unit volume, is $$ q_x = S \ \sqrt{{R_g T \over 2 \pi \mu_x}} \ \bigl[ \rho_x^s(T) -\rho_x^v \bigr] \eqno(8)$$ The specific surface can be computed as a function of porosity using a simple model for the structure of the pores. Following Mekler {\it et al.} (1990) and Prialnik (1992), we consider two models: one with interstices between randomly packed spheres of diameter $d_0$, and one with cylindrical capillaries whose circular cross-sections have diameter $d_0$. In these two cases we have $$~\hskip 1 true cm S = \cases{{6 (1 - p) \over d_0} & \hskip 47 true mm (9{\it a}) \cr {4 p \over d_0} & \hskip 47 true mm (9{\it b})\cr }$$ respectively. Note that $S$ decreases with $p$ in the packed spheres case and increases with $p$ in the capillaries case; as $p$ goes from 0 to 1, $S$ goes from $6/d_0$ to 0 in the packed spheres case and from 0 to $4/d_0$ in the capillaries case. According to (8), the sublimation time scale can be estimated as: $\tau_{sub} = \sqrt{2 \pi \mu / R_g T} / S$, and inserting (9a) or (9b) with $p \sim 1/2$, we get: $\tau_{sub} \sim d_0 \sqrt{2 \pi \mu / R_g T}$. The order of magnitude for $d_0 \sim 1 \mu$m and $T \sim 100$ K is: $\tau_{sub} \sim 10^{-8}$ s. This time scale is so much shorter than those of heat and gas diffusion that the vapour must be very close to saturation if ice is present, except in a surface boundary layer to be treated in Sect. 2.4. Consequently we assume complete saturation, i.e.: $\rho_w^v \equiv \rho_w^s(T)$ inside the boundary layer and $\rho_{CO}^v \equiv \rho_{CO}^s(T)$ wherever $\delta_{CO} = 1$. A gradient in the vapour mobility causes a flow of vapour through the pores in the direction of $-\nabla Q_x$, and we can write $$ {\bf J}_x = - L \ \nabla Q_x \eqno(10)$$ where ${\bf J}_x$ is the vapour flux and $L$ is the mean free path of the molecules. We consider a low-density regime with Knudsen flow, for which we have $$~\hskip 1 true cm L = \cases{{16 \over 3}{p^{3/2} d_0 \over (1 - p)^{1/3}}& \hskip 39 true mm (11{\it a})\cr {4 \over 3} {p d_0 \over \xi^2}& \hskip 39 true mm (11{\it b}) \cr}$$ corresponding to the two above-mentioned pore models (Mekler {\it et al.} 1990, Prialnik 1992). In (11$b$), representing the cylindrical capillary model, $\xi$ is the tortuosity -- a dimensionless parameter of order unity. Note that $L$ increases with $p$ in both pore models; as $p$ goes from 0 to 1, $L$ goes from 0 to $\infty$ in the packed spheres case and from 0 to $2d_0/3$ in the capillaries case for a typical value of $\xi=\sqrt{2}$. In the range of porosities to be considered here ($0.5 \le p \le 0.8$), $L$ goes from $2.4 d_0$ to $6.5 d_0$ and from $0.33 d_0$ to $0.53 d_0$ in the two cases, respectively. Since our assumption of the Knudsen regime, to be discussed in Sect. 4, implies that a molecule typically collides with the pore walls long before it encounters another molecule, we are justified in assuming that the different vapours flow independently of each other. Introducing for convenience $$ F_x = L \ \sqrt{{R \over 2 \pi \mu_x}} \eqno(12)$$ we get $$ {\bf J}_x = - F_x \ \nabla \left( \rho_x^v \sqrt{T} \right) \eqno(13)$$ \titleb{Evolutionary equations} Crystallization of amorphous ice is considered to proceed at a rate that varies with temperature according to an exponential activation law as found by Schmitt {\it et al.} (1989): $$ {\partial \rho_{H_2O}^a \over \partial t } = -\lambda(T) \ \rho_{H_2O}^a \eqno(14)$$ with $$\lambda(T) = 1.05 \times 10^{13} \,\exp (-5370/T) \ \ {\rm s}^{-1}\eqno(15)$$ As a consequence of crystallization, a release of latent heat occurs, to which we shall return when discussing the heat balance. A second consequence is the release of trapped CO, which according to (3) amounts to $$ {\partial \rho_{CO}^t \over \partial t } = f_t \ {\partial \rho_{H_2O}^a \over \partial t } \eqno(16)$$ The equations of continuity for H$_2$O and CO, involving time derivatives of both vapour and ice densities, may be written $$ {\partial \rho_{w}^v \over \partial t } + {\partial \rho_{H_2O}^i \over \partial t } + \nabla \cdot {\bf J}_w = 0 \eqno(17)$$ $$ {\partial \rho_{CO}^v \over \partial t } + {\partial \rho_{CO}^i \over \partial t } + {\partial \rho_{CO}^t \over \partial t } + \nabla \cdot {\bf J}_{CO} = 0 \eqno(18)$$ These equations are supplemented by the relevant equations for the rate of change of vapour or ice density, depending on the circumstances. For the case of H$_2$O, in the interior of the nucleus we have: $\rho_w^v=\rho_w^s(T)$ and in the boundary layer $\partial\rho_{H_2O}^i/\partial t$ is replaced by $q_w$ as computed from (8). For CO we have either $\rho_{CO}^v=\rho_{CO}^s(T)$ ($\delta_{CO}=1$) or $\rho_{CO}^i=0$ ($\delta_{CO}=0$). In the energy equation the following heat sources and sinks are considered: \noindent $\bullet$ the divergence of the heat flow due to the bulk thermal conductivity ($K$) of the solid matrix: $\nabla \cdot (K \nabla T)$; \noindent $\bullet$ the exchange of sensible heat between the gas and the solid matrix, assuming the two phases to be in local thermodynamic equilibrium: $-(C_w^v {\bf J}_w + C_{CO}^v {\bf J}_{CO}) \cdot \nabla T$, where $C_w^v$ and $C_{CO}^v$ are the specific heats of the two vapours; \noindent $\bullet$ the heat involved in the evaporation or condensation of ice: $H_w \partial \rho_{H_2O}^i / \partial t + H_{CO} \partial \rho_{CO}^i / \partial t$, where $H_w$ and $H_{CO}$ are the latent heats of sublimation of the two ices; \noindent $\bullet$ the latent heat of crystallization of amorphous ice $H_{cr}$ corrected for the energy expense ${{\cal{H}}}_{CO}(T)$ of the release of trapped CO: $\ -[H_{cr} - f_t {\cal{H}}_{CO}(T)] \ \partial \rho_{H_2O}^a / \partial t$. Thus the equation can be written $$ \eqalignno{ C \rho {\partial T \over \partial t } & = \nabla \cdot (K \nabla T) - (C_w^v {\bf J}_w + C_{CO}^v {\bf J}_{CO}) \cdot \nabla T \cr & \qquad + H_w {\partial \rho_{H_2O}^i \over \partial t } + H_{CO} {\partial \rho_{CO}^i \over \partial t } \cr & \qquad \quad- \bigl[ H_{cr} - f_t {\cal{H}}_{CO}(T)\bigr] {\partial \rho_{H_2O}^a \over \partial t } & (19)\cr}$$ In the left-hand member the overall specific heat per unit volume, $C \rho$, is given as a simple sum, $\sum_\alpha C_\alpha \rho_\alpha$, over all components ($\alpha$) of the material (solids and vapours). Thus, since $\rho = \sum_\alpha \rho_\alpha$ is the total density, $C$ is a mass-weighted average of the specific heats of the different components. For the computation of the thermal conductivity of the solid matrix, we make use of Greenberg's (1982) model for cometary nuclei, i.e., an agglomeration of $\mu$m-size grains, composed of a silicate and organic refractory core and a water ice mantle. The heat conductivity for a spherical grain with such a structure ($K_{cm}$) has been computed by Haruyama {\it et al.} (1993) by the following formula, which has the same form as the average dielectric function given by the Maxwell-Garnet theory: $$K_{cm} = K_w \left[ 1 \ + \ {3 g (K_d - K_w) \over K_d + 2 K_w - g (K_d - K_w)} \right] \eqno(20)$$ Here $g$ is a parameter related to the initial dust/ice ratio $R_{di}^0$ by $$ g = \left[ 1 + {\rho_{d,c} \over R_{di}^0 \rho_{w,c}^i}\right] ^{-1} \eqno(21)$$ Note that $K_{cm}$ is thus dominated by $K_w$, being of the same order of magnitude as the latter, independent of $K_d$, for any reasonable value of $R_{di}^0$. $K_w$ is computed as a mass-weighted average of the amorphous and crystalline bulk thermal conductivities. The thermal conductivity of the solid matrix is then computed as $$ K = \psi(p) \ {K_{cm} \ (\rho_d + \rho_{w}^i) + K_{CO}^i \ \rho_{CO}^i \over \rho_d + \rho_{w}^i + \rho_{CO}^i} \eqno(22)$$ including the correction for porosity by the factor $$\psi = 1 - p^\gamma \eqno(23)$$ (Smoluchowski 1981; Prialnik 1992). Different values of $\gamma$ are discussed in Sect. 2.7, as well as the temperature dependence of the bulk conductivity. We compute ${\cal{H}}_{CO}(T)$, following Prialnik (1992), as $$ {\cal{H}}_{CO} = \int_0^T \left(C_{CO}^v - C_{H_2O}^a\right) dT \eqno(24)$$ Values for the specific heats of vapours and ices, in particular that of amorphous ice, and results for $f_t {\cal{H}}_{CO}$ will be discussed in Sect. 2.7. In view of the above-mentioned approximation of saturated vapour densities, let us write down the explicit forms of the evolutionary equations taking the approximation into account. Considering the interior of the nucleus, i.e., excluding the surface boundary layer, we have $$ {\partial \rho_{H_2O}^a \over \partial t } = -\lambda(T) \ \rho_{H_2O}^a \eqno(25a)$$ $$ \left(C \rho + H_w {d \rho_w^s \over dT} + \delta_{CO} H_{CO} {d \rho_{CO}^s \over dT }\right) {\partial T \over \partial t } = \nabla \cdot (K \nabla T) $$ $$\eqalignno {\qquad & + \Bigl[F_w C_w^v \nabla \left(\rho_w^s \sqrt{T} \right)+ F_{CO} C_{CO}^v \nabla \left( \rho_{CO}^v \sqrt{T}\right)\Bigr] \cdot \nabla T \cr & + F_w H_w \nabla^2 \left(\rho_w^s \sqrt{T}\right) + \delta_{CO} F_{CO} H_{CO} \nabla^2 \left( \rho_{CO}^s \sqrt{T}\right) \cr & - \bigl[ H_{cr} - f_t {\cal{H}}_{CO}(T)\bigr] {\partial \rho_{H_2O}^a \over \partial t } & (25b) \cr } $$ $$\rho_w^v = \rho_w^s(T) \eqno(25c)$$ $${\partial \rho_{H_2O}^i \over \partial t } = F_w \nabla^2 \left( \rho_w^s \sqrt{T} \right) - {d \rho_w^s \over dT} {\partial T \over \partial t} \eqno(25d)$$ and for the CO densities, in addition to eq. (3) for $\rho_{CO}^t$, $$ \delta_{CO} = 1 \ : \cases{\rho_{CO}^v = \rho_{CO}^s & \hskip 15 true mm (26{\it a})\cr {\partial \rho_{CO}^i \over \partial t } = F_{CO} \nabla^2 \left( \rho_{CO}^s \sqrt{T} \right) \cr \qquad \qquad - {d \rho_{CO}^s \over dT} {\partial T \over \partial t} - {\partial \rho_{CO}^t \over \partial t } & \hskip 15 true mm (26{\it b}) \cr}$$ or $$\delta_{CO} = 0 \ : \cases{ {\partial \rho_{CO}^v \over \partial t } = F_{CO} \nabla^2 \left( \rho_{CO}^v \sqrt{T} \right) - {\partial \rho_{CO}^t \over \partial t } & \hskip 4 true mm (27{\it a}) \cr \rho_{CO}^i = 0 & \hskip 4 true mm (27{\it b}) \cr}$$ In (25$b$) the terms in the left-hand member involving the latent heats of sublimation are negligible compared to the $C \rho$ term. In (26$b$) the last term is included for completeness, but in fact the temperatures required for significant crystallization usually do not allow the existence of CO ice. The computations are started with $\delta_{CO} =0$, and this mode lasts until $\rho_{CO}^v$ reaches $\rho_{CO}^s$, when we instantaneously skip to $\delta_{CO}=1$. At this moment $\rho_{CO}^i$ starts increasing from zero to positive values, and the mode lasts until the ice density decreases to zero again. Whenever this occurs, $\rho_{CO}^v$ starts dropping below $\rho_{CO}^s$, and we instantaneously skip to $\delta_{CO}=0$. \titleb{Surface boundary layer} \medskip The only exception to the above-mentioned rule that the presence of ice implies saturated vapour is bound to occur next to the surface, where the gas is flowing out into the coma, i.e., an extreme low-pressure region. Adjacent to this very efficient sink there is a layer with unsaturated vapour, where the ice is sublimating at a high rate, thus feeding the outflow. We can estimate the thickness of this boundary layer ($\Delta$) to be of the order of several pressure scale heights, i.e. several mean free paths ($L$), which are of the order of the pore diameter $d_0$ according to (11). In the typical situation to be considered, this means that the boundary layer is extremely thin, i.e., $\Delta \la 100 \,\mu$m. Our model includes two regimes to be defined below. In one of these the surface temperature is so low that the sublimation rate cannot have any significant influence on the structure of the nucleus. Consequently no boundary layer is considered. In the other regime there is substantial H$_2$O outflow leading to a thin but finite boundary layer. As we shall see, the surface temperature is then high enough that it would not be reasonable to assume the presence of either CO ice or amorphous H$_2$O. Thus the boundary layer is characterized by $\delta_{CO} = 0 $ and $\rho_{H_2O}^a = 0 $. Let us denote the two regimes as the inactive and active ones, respectively. We have seen that the sublimation time scale $\tau_{sub}$ is extremely short. The relevant time scale to set up the structure of the boundary layer should be the one of gas diffusion through the layer ($\tau_d$). This can be estimated from the typical number of collisions with pore walls and the typical time elapsed between successive collisions. For our standard values, we find that $\tau_d \lta 10^{-3}$ s, so it is extremely short compared with the time scales characterizing the orbital motion and the processes occurring in the interior of the nucleus. As previously proposed by Prialnik (1992), we therefore consider a steady-state structure for the boundary layer. The steady state should naturally hold in a frame co-moving with the surface. For a typical erosional velocity at 1 AU of $v_e \sim 10^{-6}$ m/s, we find that the time scale for erosion to proceed through the whole boundary layer is much longer than $\tau_d$. Consequently we can use a simplifying approximation valid only in the boundary layer, considering the steady state to hold in a fixed frame as well. The continuity equations (17) and (18) for the vapours then yield $$ - q_w + \nabla \cdot {\bf J}_w = 0 \eqno(28)$$ $$ \nabla \cdot {\bf J}_{CO} = 0 \eqno(29)$$ replacing $\partial \rho_{H_2O}^i / \partial t$ by $-q_w$ as found from (8). According to (29) the CO vapour flux is constant throughout the boundary layer, so we only have to solve for the variations of water vapour flux and temperature. The heat balance equation becomes $$ \nabla \cdot (K \nabla T) - \nabla T \cdot (C_w^v {\bf J}_w + C_{CO}^v {\bf J}_{CO}) - H_w q_w = 0 \eqno(30)$$ Due to the very small thickness of the boundary layer, the solution of (28) and (30) can be found in plane-parallel geometry. Let us thus formulate the equations in terms of the depth $z$ below the surface and introduce the scalar variables $J_w = L \ \partial Q_w / \partial z$ and $G = K \ \partial T / \partial z$, denoting the outward mass flux of H$_2$O vapour and the outward conductive heat flux, respectively. Similarly, $J_{CO}$ denotes the outward flux of CO vapour. We get the following system of first-order ordinary differential equations for $T, G, Q_w$ and $J_w$: $$ {\partial T \over \partial z} = {G \over K} \eqno(31a)$$ $$ {\partial G \over \partial z} = - {G \over K} \bigl(C_w^v J_w + C_{CO}^v J_{CO}\bigr) + H_w S \bigl(Q_w^s - Q_w\bigr) \eqno(31b)$$ $$ {\partial Q_w \over \partial z} = {J_w \over L} \eqno(31c)$$ $$ {\partial J_w \over \partial z} = - S \bigl(Q_w^s - Q_w\bigr) \eqno(31d)$$ In these equations, $K$, $Q_w^s$ and $H_w$ are given functions of $T$, while $C_w^v$, $C_{CO}^v$, $S$ and $L$ are constants. \titleb{ Boundary conditions} In our spherically symmetric model for the interior the heat and gas fluxes, and thus the gradients of temperature and vapour densities, must vanish at the center. Therefore we have $${\partial T\over\partial r}\Bigr\vert_{r=0}=0\eqno(32)$$ $${\partial\rho_w^v\over\partial r}\Bigr\vert_{r=0}=0\eqno(33)$$ $${\partial\rho_{CO}^v\over\partial r}\Bigr\vert_{r=0}=0\eqno(34)$$ (33) is listed for completeness only, since it follows immediately from (25$c$) and (32). The corresponding conditions at the surface take somewhat different forms depending on whether the nucleus is in the active or inactive regime, i.e., whether the surface boundary layer is present or absent. A common boundary condition is used only for the CO flux, since this is conserved throughout the boundary layer. The form adopted is $$\rho_{CO}^v\bigr\vert_{r=R}=\delta_{CO}\bigr\vert_{r=R} \cdot\rho_{CO}^s(T_S)\eqno(35)$$ where $T_S$ is the surface temperature. As we have seen, in the active regime we always have $\delta_{CO}\bigr\vert_{r=R}=0$. Since the thickness of the boundary layer is negligible on the scale of the interior density gradient, it is immaterial whether we identify this level with the surface or the interface with the interior, which we shall call ``the boundary''. We adopt a pragmatic criterion for shifting between the active and inactive regimes, based on the fact that the concept of a boundary layer becomes ill-defined at too low surface temperatures. The critical temperature is $\simeq 150$--160 K. The inactive regime is characterized by very slow sublimation even from the surface, and whether in such a situation we treat the H$_2$O outgassing accurately by allowing for subsurface sublimation or we simply neglect the latter is not a matter of great importance. One further approximation used in the current model is a crude treatment of the return flux of H$_2$O molecules from the coma. In the inactive regime the assumption of free molecular outflow is a good one, implying no return flux, but going to the active regime we immediately shift to a hydrodynamic outflow with a constant ratio ($\alpha$) between the return flux and the outward H$_2$O flux at the surface. This means that the H$_2$O production rate per unit surface area is given by $${\cal{F}}_{H_2O} = \cases{Q_w^s(T_S) &; inactive\ reg.\cr (1-\alpha)\bigl[Q_w^s(T_S)+J_{w,S}\bigr] &; active\ reg. \cr} \eqno(36)$$ where $Q_w^s(T_S)$ represents the evaporation flux from the surface and $J_{w,S}$ is the surface value of the outward H$_2$O vapour flux in the boundary layer. Accounting for the heating effect of the return flux, we introduce another approximation, i.e., that all of this flux is absorbed at the surface and none of it penetrates into the boundary layer. Thus the surface energy balance equation becomes $${\cal I} = \cases{(1-A_{IR})\sigma T_S^4+K{\partial T\over\partial r}\bigr\vert_{r=R}+H_wQ_w^s(T_S) \cr \qquad \qquad +\delta_{CO}\bigr\vert_{r=R} H_{CO}Q_{CO}^s(T_S) & \hskip 3 true mm (37{\it a})\cr (1-A_{IR})\sigma T_S^4-G_S+(1-\alpha)H_wQ_w^s(T_S) \cr \qquad \qquad -\alpha H_wJ_{w,S} & \hskip 3 true mm (37{\it b})\cr}$$ for the inactive and active regimes, respectively. We assume a uniform average insolation for a spherical surface with visual Bond albedo $A_v$ at heliocentric distance $r_h$ $${\cal I}={1\over4}(1-A_v){\cal C}_\odot r_h^{-2}\eqno(38)$$ where ${\cal C}_\odot$ is the solar constant. The thermal infrared Bond albedo is $A_{IR}$, and $\sigma$ is the Stefan-Boltzmann constant. Equation (37a) accounts for the possibility of CO ice existing at the surface, even though this situation does not occur under the circumstances investigated in this paper. The surface boundary conditions for the H$_2$O vapour density in the two regimes are $$ \rho_w^v\bigr\vert_{r=R}=\rho_w^s(T_S)\qquad ;\quad {\rm inactive\ regime} \eqno(39a)$$ $$ J_{w,S}=Q_{w,S}\qquad\; ;\quad {\rm active\ regime} \eqno(39b)$$ In the active regime we need additional conditions at the boundary. The first of these is a temperature condition, expressing the requirement that the temperature and its gradient vary continuously across the boundary. Denoting by ``$\vert_{B^-}$'' the limiting value of any quantity when approaching the boundary from the interior and by ``$\vert_{B^+}$'' the limiting value from the side of the boundary layer, we have $$T\vert_{B^+}=T\vert_{B^-}\eqno(40a)$$ $$G\vert_{B^+}=-K{\partial T\over\partial r}\bigr\vert_{B^-} \eqno(40b)$$ Secondly, from the continuity of H$_2$O vapour density and temperature, we get $$Q_w\vert_{B^+}=Q_w^s(T\vert_{B^-})\eqno(41)$$ which, in combination with (31$d$) and (40$a$), yields $${\partial J_w\over\partial z}\bigr\vert_{B^+}=0\eqno(42)$$ We have an additional relation between $J_w\vert_{B^+}$ and $G\vert_{B^+}$, derived from (10), $$J_w\vert_{B^+}={L\over K}\ {dQ_w^s(T)\over dT}\ G\vert_{B^+} \eqno(43)$$ The ensemble of equations (32)--(43), which obviously are not independent, allows to specify the thickness of the boundary layer and the exact variation of temperature and vapour density as well as to couple this self-consistently to the interior temperature and vapour density profiles, as explained in Sect. 2.6. \titleb{Numerical scheme} \medskip \titlec{For the interior} \noindent In the spherically symmetric model described above, the cometary nucleus is characterized by its radius ($R$). The densities and temperature are functions of time and the distance $r$ to the centre. The processes modelled in this paper (i.e., crystallization, and H$_2$O and CO sublimation-recondensation) mainly occur in the outer part of the nucleus. Thus, the introduction of an exponential distribution of depth levels is justified (longer steps as one goes deeper into the nucleus), so we consider a new spatial variable ($x$), given by $$ x = \exp \left[\beta \left({r \over R} - 1\right)\right] \eqno(44)$$ where $\beta$ is a free scaling parameter to be fixed by, e.g., a desired depth step next to the surface. $x$ varies in the range: $\exp(-\beta) 61.55$ K \cr} \eqno(58)$$ The specific heats of vapours, as derived from the ideal gas law, are: $C_w^v = 1528 \ {\rm J kg}^{-1} {\rm K}^{-1}$ for H$_2$O, and $C_{CO}^v = 738 \ {\rm J kg}^{-1} {\rm K}^{-1}$ for CO. The latent heat of sublimation for H$_2$O is $$H_w = 2.888 \times 10^6 - 1116 \ T \ \ {\rm J kg}^{-1} \eqno(59)$$ (Delsemme and Miller 1971). A constant value for the latent heat of CO sublimation ($H_{CO} = 2.93 \times 10^5 \ {\rm J kg}^{-1}$) has been adopted, following Clayton and Giauque (1932). The latent heat of crystallization is $H_{cr} = 9 \times 10^4 \ {\rm J kg}^{-1}$ (Ghormley 1968). Using eqs. (24) and (57) with the above values for $C_{CO}^v$ and $f_t$, we get $f_t {\cal{H}}_{CO} \la 6 \times 10^3 \ {\rm J kg}^{-1}$ for $T \la 140$ K, showing that the release of CO consumes only a very small fraction of the crystallization heat (Prialnik 1992). The saturated vapour densities are given by $$\rho_x^s = A_x \cdot 10^{B_x(T)}/T \eqno(60)$$ with $A_w = 0.289 \ {\rm kg\ m}^{-3}{\rm K}$ and $A_{CO} = 4.49 \ {\rm kg\ m}^{-3}{\rm K}$, while the exponents are $$\eqalignno{ B_w & = - 2446 / T + 8.23 \ ^{10}\log(T) - 0.0167 \ T \cr & \qquad + 1.205 \times 10^{-5} \ T^2 - 6.757 &(61)\cr}$$ (Washburn 1928) and $$B_{CO} = \cases{- 418.2 / T + 4.127 \ ^{10}\log(T) - 0.0262 \ T + 1.474 \cr \qquad \qquad {\rm if\ } T \le 61.55 {\rm \ K} \cr - 425.1 / T + 7.823 \ ^{10}\log(T) - 0.00760 \ T \cr \qquad \qquad {\rm if \ } T > 61.55 {\rm\ K} \hskip 3 true cm (62) \cr}$$ (Clayton and Giauque 1932). For the nuclear surface characteristics, we use a low visual albedo of 0.04, as indicated by photometric and radiometric observations (see e.g. Jewitt 1991). For the IR albedo we use the same value, since typically this is also considered to be low (Delsemme 1985; Campins {\it et al.} 1987). The feed-back coefficient $\alpha$ is given the value of 0.25, which is typical of the hydrodynamic outflow of the collisional gas in an active comet near Earth's orbit (Crifo 1987). Our standard case for the orbital parameters is a perihelion distance of $q = 1.5$ AU and an aphelion distance of $Q = 6$ AU, typical for the observed sample of Jupiter family comets. The initial radius of the nucleus is 3 km, which is probably typical for many of these comets. The runs start at aphelion with an isothermal, low temperature nucleus ($T = 30$ K). Values between 0.2 and 0.4 m are used for the depth step at the surface. For the initial radius, the scale factor $\beta$ thus ranges from 30 to 15, respectively. \titlea {Results} \titleb{Comparison with previous models} The model just described differs in several respects from those previously published for the crystallization and chemical differentiation of cometary nuclei (Espinasse {\it et al.} 1989; Prialnik and Bar-Nun 1990; Espinasse {\it et al.} 1991; Prialnik 1992). Most important is probably the fact that we treat the H$_2$O and CO vapours as saturated in the interior of the nucleus wherever the corresponding ice is present. By contrast, the other authors have dealt with a formally unsaturated vapour, explicitly introducing eq. (8) as one of the basic equations governing the interior structure. Additional differences are that none of the other investigators included the erosional velocity of the surface in their equations, and that no hydrodynamic return flux from the coma was included into any previous model. Finally, while we assume Knudsen flow under all circumstances, Espinasse's model included a transition to hydrodynamic Poiseuille flow when the Knudsen number became too small. In view of such differences between the physical assumptions involved, we deem it essential to make comparisons with the previous models in order to show the effects of the differences, thus allowing uniform judgement of our results with the earlier ones within a common framework. There are important differences in the choice of physical parameters as well, so to make the comparisons meaningful we had to taylor our comparison models to fit the parameters chosen by the other authors as closely as possible. Let us start by a comparison with Prialnik's (1992) model. In this case we choose the orbit of P/Halley, a pore diameter of 20 $\mu$m, Klinger's amorphous ice conductivity of eq. (55$a$) and an infrared albedo of 0.5. The dust/ice ratio is 1/3 and the porosity is 0.63. The radius of the nucleus is 2.5 km. This model has been run for one orbital revolution as in Prialnik's paper. While the general behaviour is similar, involving a major crystallization spurt at $r_h \simeq 7$ AU pre-perihelion and several less dramatic post-perihelion spurts continuing to much larger distances, there are important differences as well. The fact that our H$_2$O vapour is saturated means that sublimation in the interior occurs only in response to changes of temperature and the diffusion of the saturated vapour, as expressed by eq. (25$d$). This limits the influence of the water sublimation term in eq. (19) and leads us to conclude that this term is not of major importance for the interior temperature profile -- an illustration is provided by Fig. 5 below. On the other hand, in Prialnik's case the modelling in terms of unsaturated H$_2$O vapour in principle allows for a possibility of much higher sublimation rates, should the deficit in vapour density be important enough. That such is actually the case was clearly demonstrated by Prialnik (1992), who stressed as a main result the flat temperature profile ensuing over a wide range of depths following the passage of the crystallization front. In her Fig. 1 this was compared with variant profiles obtained at the same heliocentric distance ($\sim 5$ AU pre-perihelion) by artificially cutting off different terms in the heat equation. \begfig 8 cm \figure{1}{P/Halley comparison model. Two temperature profiles at $r_h \simeq 5$ AU pre-perihelion. The full-drawn curve corresponds to an instant of active crystallization and the dotted curve is registered in a quiescent state 50 days later. The total radius of the comet is 2500 m. } \endfig We present in Fig. 1 two temperature profiles resulting from our model at $r_h \simeq 5$ AU pre-perihelion. The full-drawn curve corresponds to an instant of active crystallization and the dotted curve is registered shortly afterward in a quiescent state (see Sect. 3.2). We note a great similarity between our full-drawn curve and curve {\sl 'C'} of (Prialnik 1992, Fig. 1) which was obtained by cutting off the sublimation term in the heat balance equation. This indicates that the heat transfer rate due to H$_2$O sublimation in interior layers in Prialnik's stand-off model (curve {\sl 'A'} in her diagram) is much higher than in our model. As a consequence we do not reproduce her flat temperature profile, and the near-surface temperature gradient is negative ($T$ decreasing downward) rather than positive in most situations. Further discussion of the general features seen in Fig. 1 will be given in Sect. 3.2. A very important difference between our model and that of Prialnik (1992) is that she does not allow CO ice to condense below the crystallization front. The reason is not clear to us, since her equations do not include any obstacle to CO ice formation. It is true that an interior temperature of 30 K is on the high side as far as CO condensation in the solar nebula is concerned, but the vapour pressures resulting from the release of 10~\% CO during crystallization are far higher than the CO saturation pressure at 30 K. This causes a difference between our results and those of Prialnik (1992) in that we get crystallization of a layer of only 57 m thickness, whereas in her case this continued to 230 m. The necessity in our model to sublimate an important quantity of CO ice for the crystallization front to propagate downward causes a significant slowing down of the process, whereas in Prialnik's model this obstacle was not present. For the comparison with Espinasse {\it et al.} (1989, 1991) we ~ concentrate on ~ the ~ Jupiter ~ family ~ comet ~~~ P/Churyumov-Gerasimenko ($q=1.29$ AU and $Q=5.72$ AU), and to adapt the physical parameters to those of their model we treat a pure ice nucleus ($R_{di}^0=0$) with radius $R=1$ km, porosity $p=0.8$ and 5~\% initial trapped CO. We use Klinger's amorphous ice conductivity, and the models start with a uniform temperature of 10 K. Apart from the circumstance that our numerical procedure is entirely different from that applied by Espinasse, a fact whose consequences cannot be judged without going into a much more detailed comparison, we identify several differences of modelling assumptions in addition to those described above, which are specific to the comparison with Espinasse's model. The most important of these is probably that in her case the energetics of the release of trapped CO is treated like that of sublimation, so ${\cal H}_{CO}$ is replaced by $H_{CO}$, the latent heat of sublimation. Whereas in our case the release of CO consumes a very small fraction of the crystallization energy $H_{cr}$, this fraction is thus quite important in her models. Another significant difference is that the Espinasse models do not include any surface boundary layer. \begfig 7 cm \figure{2}{The evolution over the first 16 yr of P/Churyumov-Gerasimenko comparison model. The upper plot corresponds to the H$_2$O flux (full-drawn curve) and CO flux (dotted curve). In the lower plot we find, from top to bottom: total radius of the nucleus, radius where the amorphous ice is half the original value (dotted curve), radius of the outer border of the CO ice layer, radius of the inner border of the CO ice layer. The lower $x$-axes in both plots are expressed in years, while the upper ones in number of revolutions starting at aphelion. } \endfig The temperature profiles obtained show a broad similarity: temperatures a few meters below the surface in the crystalline layer range from 120 to 140 K in Espinasse's case, while we get a flatter profile near 140 K. The same peak structure in the temperature profile during a crystallization spurt (shown in Fig. 1 for a P/Halley model) is also found using both models. In the region where CO ice is recondensing we both get temperatures of $\sim 60$ K. Comparable amounts of CO ice are recondensed, reaching during the first revolution around 10~\% of the water ice density in the region below the crystallization front. Regarding the long-term propagation of this front, however, we find important differences. Fig. 2 shows the results of our comparison model over the first 16 years. As one might expect from a pure ice model with relatively little trapped CO, a furious rate of crystallization ensues. Within $\sim 10$ years more than half the mass of the nucleus has undergone the phase transition, and the CO outgassing rate is constantly very high -- about 10~\% of the perihelion H$_2$O outgassing rate. A similar runaway behaviour is seen in the results by Espinasse {\it et al.} though much less dramatic, and the most likely reason is that they consume a much larger fraction of the crystallization energy for CO release, as explained above. After 15 years they get crystallization of a layer of only 60 m thickness whereas in our case the depth of the front is 400 m. This is reflected in our CO flux being almost an order of magnitude higher than in their case. The difference between the temperature profiles in the crystalline layer appears to be another consequence of the same phenomenon, since our flatter profile reflects the fact that the crystallization zone is kept at a significantly higher temperature. Even the near-aphelion surface temperature is affected by this difference, at least in the beginning: we get 138 K whereas Espinasse {\it et al.} get 130 K. Due to the surface boundary layer, we also get a slightly lower near-perihelion surface temperature: 195 K for their model and 189 K for our model. \titleb{The standard model} The standard model was run for 500 years, a period of time comparable to $\sim 10$~\% of the typical dynamical lifetime of observable Jupiter family comets before ejection into orbits with larger $q$ (Lindgren 1992). Our description of these results will be divided into two parts: first we analyse a general picture of the depth structure of cometary nuclei typical of any stage of the evolution (a ``typical'' snapshot of the interior); then we consider how the structure evolves with time. The thermochemical characteristics of the cometary interior from the surface down to the center can be described as follows: \noindent$\bullet$ Starting at the surface, we first find the boundary layer, with a typical thickness of 30--40 $\mu$m (this layer is only solved for when the comet is close to perihelion and the surface temperature is higher than 150--160 K). Figure 3 shows three profiles for $T$, $\partial T/\partial z$, $\log Q_w$ and $\log J_w$ at three different orbital positions near the perihelion passage. The temperature gradient exhibits a steep decrease in absolute value over a short distance (a few $\mu$m) just below the surface, and it is by far not large enough to produce any considerable temperature contrast between the surface and the boundary. The gradient is smaller after than before perihelion, even though the surface temperature is higher, since the subsurface layers have been heated during perihelion passage. The flow of water vapour changes direction at $z \simeq 15\,\mu$m, from an outward flux in the outermost layer to an inward flux that proceeds via the boundary and far down toward the center. The magnitude of the latter is negligible compared to the flux that is escaping through the surface (typically, $-10^{-10}$ kg$\,$ m$^{-2}\,$s$^{-1}$ for the inward flux and $+10^{-5}$ kg$\,$m$^{-2}\,$s$^{-1}$ for the outward flux near perihelion). \begfigwid 16cm \figure{3}{Three boundary layer profiles of the standard model for $T$, $\partial T/\partial z$, $\log Q_w$ and $\log J_w$ at three different orbital positions near the perihelion passage (relevant information for each position is given on top of the plots). Note that, since $J_w$ changes its sign within the boundary layer, we plot $\log (|J_W|)$. The small thick vertical bars correspond to the depths where the change of sign occur. } \endfig \noindent$\bullet$ As illustrated by Fig. 4, below the boundary layer we first find a layer of crystalline water ice with no CO ice. The water vapour flows inward, while the CO vapour flows outward from the zone of crystallization. The H$_2$O flux is many orders of magnitude smaller than the CO flux, however. Near the surface the temperature profile reflects the orbital variation of the surface temperature, so that in the plotted case (near aphelion) there is an increase with depth, while near perihelion there is a rapid decline from a surface temperature reaching $\simeq 186$ K. From $\sim 10$ m below the surface down to the vicinity of the crystallization zone, the temperature decreases slowly independent of orbital position and stays, typically, in the range 120--130 K. \noindent$\bullet$ We then find a zone where crystallization occurs and large amounts of CO are being released and condensed or resublimated. Figs. 4 $a$ and $b$ illustrate how the structure of this zone changes between two modes: the quiescent one where crystallization is very slow and the active one with rapid crystallization, respectively. The zone has a typical thickness $\la 10$ m, within which the amorphous fraction goes essentially from 0 to 1 over a range of $\sim 1$ m (upper right panels). It acts as a source for CO vapour that diffuses away towards the surface and the interior. In the quiescent mode the crystallization front pushes a moderate amount of CO vapour ahead of it, increasing the density to $\simeq 0.2$ kg$\,$m$^{-3}$, whereas in the active mode the front is rushing downward and the peak reaches $\simeq 1$ kg$\,$m$^{-3}$ (lower left panels). In both cases this is accompanied by condensation of CO ice (lower right panels), which acts as a buffer for the large amount of released CO in the active mode (Fig. 4$b$). In addition, we note that the presence of CO ice envelopping the amorphous ice grains next to the crystallization front presents an obstacle to the propagation of the front: the sublimation of the CO ice consumes a considerable amount of energy. The latent heat released in the phase transition keeps the temperature profile flat on the crystalline side in the quiescent mode and causes a slight temperature peak in the active mode (upper left panels). On the amorphous side, the extensive CO condensation raises the temperature somewhat, thus accounting for the steep profile of saturated vapour density. \begfigwid 21 cm \figure{4}{Interior profiles of the standard model for $T$, $\rho_w^i$ (full line) and $\rho_{H_20}^a$ (dashed line), $\rho_{CO}^v$ and $\rho_{CO}^i$ during a quiescent (a) and active (b) phase of crystallization. In the water density plot, {\tt A} corresponds to the amorphous ice region and {\tt C} to the crystalline ice one. } \endfig \noindent$\bullet$ Below this zone the water ice is amorphous, and there is a region with condensed CO ice remaining from previous crystallization events, reflected in the wiggles on the ice density curves (lower right panels). The latent heat of these condensation events keeps the temperature at $\sim 50$--60 K over a considerable depth range. Note that the amount of CO ice in this layer is quite substantial, representing nearly 10~\% of the total ice density. At the lower boundary, the CO ice density drops to zero and the temperature reaches its central value of 30 K. There is a leakage of CO vapour across the boundary, and as a consequence there is CO ice formation and increase of the temperature. Hence the CO ice layer gradually spreads downward. In the inner ice-free zone the heat is again mainly conducted by the solid matrix and the CO vapour flows inward. To some extent this region may be considered as a low-temperature analogue of the near-surface layer with pure crystalline ice and inward-flowing water vapour. The heat transfer processes are illustrated by Fig. 5, where the values of the different terms in the heat balance equation (19) are plotted as functions of radius at almost exactly the same time as the snapshot in Fig. 4$a$. The full-drawn curve marks the rate of change of internal energy (left-hand member plotted as moved to the righ-hand side), which is balanced by the sum of all the other terms. Conduction via the solid matrix is seen to dominate among the latter throughout the crystalline ice layer. In the crystallization zone the latent heat release takes over as the dominant process, and below this the most important role is played by the latent heat of CO sublimation. \begfig 10 cm \figure{5}{The contribution of the different terms to the heat balance equation (19). The full-drawn line corresponds to minus the rate of change of internal energy (left-hand member moved to the righ-hand side), the dashed-dotted line to the conduction via the solid matrix, the long dashed line to release of crystallization latent heat, the short dashed line to the latent heat of CO sublimation, and the dotted line to the exchange of sensible heat between the CO vapour and the solid matrix (only noticeable at the crystallization front at $\simeq 2840$ m and at the inner boundary of the CO ice layer at $\simeq 2770$ m). The water sublimation and the water vapour exchange of sensible heat terms are not plotted since they always are many order of magnitudes smaller than the main contributors in each region of the plot. } \endfig Let us now consider the dynamics of the processes involved. The crystallization starts to speed up when the amorphous ice reaches a temperature $\simeq 115$--120 K, since it becomes self-sustaining: the phase transition releases enough heat to produce a local rise of temperature (due to the low thermal conductivity of amorphous ice, very little heat is conducted toward the interior). The time scale of crystallization is thus reduced and further phase transition occurs until, locally, all the ice becomes crystalline. This happens mostly in a quasi-continuous manner, and the time scale of local crystallization is $\sim 10^{\rm d}$, but the picture is from time to time perturbed by sudden crystallization ``spurts'' (cf. Prialnik and Bar-Nun 1987), when layers of several meters of amorphous ice become crystalline on a time scale of a few days. Associated with these events, we find CO bursts where the production rate increases by about half an order of magnitude in the same period. \begfigwid 11 cm \figure{6}{The evolution during $\simeq 20$ years of, from bottom to top, the rates of change of the integrated amounts (in 10$^6$ kg day$^{-1}$) of CO ice, CO vapour, trapped CO gas, and finally the amount of escaping CO gas per day. For the CO vapour and CO ice plots a line at 0 rate of change is drawn. The lower $x$-axes in all plots are expressed in years, while the upper ones in number of revolutions starting at aphelion. } \endfig Figure 6 shows the evolution during $\simeq 20$ years of, from bottom to top, the rates of change of the integrated amounts (in kg per day) of CO ice, CO vapour, trapped CO gas, and finally the amount of escaping CO gas per day. Starting with the crystallization spurt at $t = 87$ yr, which was illustrated in Fig. 4, a large quantity of trapped CO is released (the extremum is far out of the limits of the plot), producing a sudden increase of the CO gas and ice densities. The CO gas density reaches values up to $\sim 1$ kg$\,$m$^{-3}$ just below the crystallization zone (Fig. 4). The main crystallization spurt produces a perturbation of the diffusive equilibrium structure and is thus followed by secondary, less dramatic peaks as well as by several recondensation and sublimation events. These subside on the typical gas diffusion time scale of several years, as estimated using eq. (28) for a depth $\Delta \sim 100$ m. The CO outgassing flux starts to increase as soon as the large amount of CO vapour procuced at the crystallization zone reaches the surface and affects the density gradient there. The decrease of the outgassing flux after its maximum is smoother than the variations of the other plotted quantities, since considerable smearing ensues from diffusion through the crystalline layer. The following phase is quiescent. For several years, the scenario is that of slow, gradual crystallization whose rate varies smoothly with time, inducing slow variations in the amounts of CO vapour and ice -- an increase in crystallization rate eventually leads to an increased amount of CO vapour and a sudden increase of the condensation rate. This is the time when the temporary runaway starts, apparently caused by the slow heating of the transition layer reaching a critical temperature threshold. \begfigwid 9 cm \figure{7}{The standard model evolution over 500 yr of the water and CO flux, and the structural boundaries of the nucleus interior. The explanations are the same as in Fig. 2. } \endfig The long-term evolution of the model is illustrated by Fig. 7. During the first revolution, already at a large heliocentric distance on the inbound branch, a large flux of CO outgassing arises due to the crystallization of a layer several meters thick. Several crystallization spurts and CO bursts occur during this perihelion passage. The surface temperature rises to $T_S = 185$ K and the water flux to ${\cal{F}}_{H_2O} = 2.8 \times 10^{-5}$ kg$\,$m$^{-2}\,$s$^{-1}$ at perihelion, values that remain almost constant during all subsequent revolutions. After the first orbit, the crystallization front is already 33 m below the surface, but afterwards it slows down considerably. During the third revolution the CO outgassing fluxes during the quiescent phase and the burst are $6 \times 10^{-7}$ and $3 \times 10^{-6}$ kg$\,$m$^{-2}\,$s$^{-1}$, respectively, and from then to the end of the computation both fluxes drop by one order of magnitude. During the first century, the bursts with their associated crystallization spurts (not clearly visible on the scale of these plots) show a quasi-periodic appearance, but the period is longer than the orbital one so that the bursts exhibit a phase drift with respect to the orbital motion. This corresponds to the increasing diffusion time for the perihelion heat pulse at the surface to reach the amorphous ice, which is buried at increasing depth. When the time lag reaches one orbital period, the occurrence of the bursts takes on a chaotic appearance. The major ones always occur near perihelion, but an increasing fraction of the perihelia are skipped, so the bursts become rarer as time proceeds. Moreover, they gradually evolve from the sharp spikes seen in the beginning to a much more subdued appearance, as the diffusion time across the crystalline ice layer increases. \begfig 6 cm \figure{8}{The evolutions of the maximum (short-dashed line) and mean (full-drawn line) water flux over a revolution and maximum (dotted line) and mean (long-dashed line) CO flux over a revolution. } \endfig The CO/H$_2$O outgassing ratio varies enormously along the orbit, as noted earlier by Espinasse {\it et al.} (1991). The values at perihelion, after a few centuries, range between 0.1 and 1~\%, depending on whether or not a burst is occurring. However, the overall balance between the outgassing and condensation of CO leads to a mean CO outgassing amounting to $\sim 20$~\% of that of H$_2$O during the first two centuries, as shown by Fig. 8. At later stages in the evolution the CO outgassing rate, even as a mean over one orbit, becomes quite erratic due to the erratically occurring near-perihelion bursts. Thus, during some revolutions the amount of CO outgassed exceeds 15~\% of the H$_2$O production, but typically the ratio is $\sim 5$~\%. In Sect. 2.4 we defined $v_e$ as the erosional velocity of the surface, and we now define the crystallization velocity $v_{cr}$ as the time derivative of the mean radius of the crystallization front and the CO ice velocity $v_i$ as the time derivative of the radius marking the inner border of the CO ice layer. Figure 9 shows the evolutions of $v_e$, $v_{cr}$ and $v_i$ with time. Note the rapid decrease of $v_{cr}$ and $v_i$ in absolute values after the first fifty years. Both velocities continue to show a nearly linear decrease with time, however with a trend to flattening out; extrapolating this result leads us to estimate that complete crystallization requires $\sim 10\,000$ years in an orbit similar to that of our standard model. After several centuries we have $v_{cr} \simeq v_i$, showing that the thickness of the CO ice layer becomes nearly constant. We can estimate the effect of explicitly including the erosional velocity into the evolutionary equations by comparing the two right-hand members of eq. (46$b$). For $v_e \sim 0.25 m yr^{-1}$ and $\beta = 15$, the factor multiplying $\partial /\partial x$ in (46$b$) is $\sim 10^{-8}$. Considering the sharp depth gradients observed in the crystallization zone, we conclude that the contribution of the erosional velocity, though small, is not negligeable and it must be explicitly included. \begfig 6 cm \figure{9}{The evolutions of the erosional velocity $v_e$ (dashed line), crystallization velocity $v_{cr}$ (full-drawn line) and the CO ice velocity $v_i$ (dotted line) (see text). } \endfig \titleb{ Variant models} All the variant models described below are run for 150 years in the same orbit as the standard model. To facilitate comparisons, we replot in Fig. 10 the first 150 years of the standard model. \begfigwid 9 cm \figure{10}{First 150 years of the standard model evolution. The explanations are the same as in Fig. 2. } \endfig \titlec{Variation of the dust/ice ratio} \noindent As discussed above, several characteristics of the dust component are very uncertain, e.g., the ratio between silicates and organic refractories, the thermal properties of these materials and mixtures of them, the dust to ice ratio, etc. Of these, the dust to ice ratio in the comet interior is the least uncertain if it is constrained by the two usual assumptions of homogeneity and solar system abundances. Rocky and metallic alone can acount for only about 20~\% of the total mass. Based on comet P/Halley mass spectra and on the organic refractory component of interstellar dust, amounting to about 20~\%, and on additional carbonaceous component of 5~\%, one achieve a dust to ice ratio of the order of 0.82 (Greenberg 1982), consistent with our nominal choice of 1. On the other hand, the Giotto-DIDSY analysis of the P/Halley dust (McDonnell {\it et al.} 1987) gave clear indications that the total ratio of dust to gas in the material leaving the nucleus was in excess of unity and perhaps as large as 2--3; and also, recent analysis of comet dust trails has indicated values of the dust to gas ratios as high as 3 (Sykes and Walker 1992). If we interpret this as evidence for a very high dust/ice ratio in the nuclear material, we are permitted to adopt $R_{di}^0=2$ for our first variant model. We note however, that this involve abandoning the solar system abundances constraint. Neither the porosity, nor the total density is different from the standard model, but the partial densities become: $\rho_{H_2O}^a = 143 $ kg m$^{-3}$, $\rho_{CO}^t = 14 $ kg m$^{-3}$, and $\rho_d = 314 $ kg m$^{-3}$. \begfigwid 9 cm \figure{11}{The evolution of the variant model with dust to ice ratio $R_{di}^0=2$. The explanations are the same as in Fig. 2. } \endfig Fig. 11 shows the evolution of outgassing fluxes and structural boundaries for this model. The most striking difference with respect to the standard model is the slower progression of the crystallization front toward the interior. Just as in the standard model, it is accompanied very closely by a layer of condensed CO ice, but this layer is now much thinner. The reason for the slower crystallization appears to be that a larger percentage of the material is now bound up in the form of inert dust so the amount of latent heat $H_{cr}$ produced per unit mass is smaller. Nonetheless, in the beginning there are occasions when a runaway occurs, as seen from the CO bursts. This means that a critical temperature is again reached, although this should be somewhat higher than in the case of the standard model. The bursts look very pronounced, but their peak height is no larger than in the standard model. By contrast, the crystallization rate and CO production rate during the quiescent phases preceding the bursts now drop practically to zero. Again, the runaways disappear and the bursts take on a much more subdued appearance after some time, when the crystallization zone is buried deep enough. However, we see no evidence for the erratic occurrence of these later bursts that was observed for the standard model. They occur very regularly near each perihelion passage, and the CO outgassing rate is a factor 3 lower than in the standard model during both the bursts and the interludes. An interesting feature is that the phase drift of the early bursts with respect to the orbital motion is seen again but is now quicker than in the standard model. We recall our explanation for this phase drift: as the crystallization zone is buried deeper and deeper below the surface, it takes an increasing time for it to be reached by the perihelion heat pulse. But the amplitude of the temperature oscillation decreases exponentially as a function of depth, so in a range closer to the surface the phase lag in order to reach a certain temperature is a more sensitive function of depth than it is much deeper down. \titlec{Variation of the porosity} \noindent In this set of variant models the composition of the cometary nuclei as well as the physical characteristics of the materials are kept unchanged; we only vary the structure of the solid matrix. In particular, several recent analyses of nongravitational effects on short-period comets (Rickman 1989; Rickman {\it et al.} 1992) have advocated densities even lower than that of our standard model, with a preference around 300 kg m$^{-3}$. This is also consistent with the nucleus density derived by Greenberg and Hage (1990) based on the infrared emission by the observed mass (size) spectrum of comet P/Halley dust. Thus we designed a model for this purpose, increasing the porosity to $p=0.8$. The densities become: $\rho_{H_2O}^a = 122 $ kg m$^{-3}$, $\rho_{CO}^t = 12 $ kg m$^{-3}$, $\rho_d = 134 $ kg m$^{-3}$; yielding a total density $\rho_{nuc} = 269 $ kg m$^{-3}$. The specific surface increases by 23~\% with respect to the standard model and the mean free path of the molecules by 65~\%. The main effect is a decrease in the thermal conductivity of the whole solid matrix. The heat gained by the surface layer during the perihelion passage now meets a higher resistance to being conducted inward. As shown in Fig. 12, the advance of the crystallization front towards the center is only slightly slower than in the standard model, whereas due to the higher porosity the rate of surface erosion is higher. The CO bursts reach the same heights as in the standard model, but the production rate drops much deeper during the interludes. Their phase drift with respect to the orbital motion is again a striking feature; in this model the reduction in conductivity yields a smaller thermal skin depth and thus a phase lag that varies more rapidly with depth below the surface. This is in agreement with the more rapid phase drift seen in Fig. 12 as compared with Fig. 10. The average CO outgassing rate per revolution decreases slightly with respect to the standard model. \begfigwid 9 cm \figure{12}{The evolution of the variant model with a porosity $p=0.8$. The explanations are the same as in Fig. 2. } \endfig \begfigwid 9 cm \figure{13}{The evolution of the variant model with an amorphous ice conductivity $K_w^a$ according to Klinger (1980). The explanations are the same as in Fig. 2. } \endfig Another variant model involved a decrease of the exponent $\gamma$ governing the influence of porosity on thermal conductivity from the standard value of 2/3 to 0.1. For $p=0.65$ this means a reduction of conductivity similar to the one accomplished by increasing $p$ to 0.8 keeping $\gamma$ at its standard value. On the other hand, in the variant model with $\gamma=0.1$ the specific surface and mean free path keep their standard values, and the erosional rate of motion of the surface is also the same as in the standard model. These differences notwithstanding, we find that the $\gamma=0.1$ model gives a crystallization and CO outgassing behaviour very similar to that of the $p=0.8$ model, showing that the thermal conductivity is the most important parameter in this context. The last variant model of this group corresponds to a decrease of the porosity in order to explore the whole range of possible P/Halley densities found by Rickman (1989). A porosity $p=0.5$ was thus assumed, corresponding to a total density of $\rho_{nuc} = 689$ kg m$^{-3}$. In this case the results are very similar to those of the standard model. \titlec{Variation of the amorphous ice conductivity} \noindent In view of the very wide range of values suggested for the amorphous ice conductivity and our choice of a compromise for the standard model, we also consider two variant models using a much higher conductivity according to Klinger (1980) and a much lower one according to Kouchi {\it et al.} (1992) as explained in Sect. 2.7. The low conductivity model yields results very similar to the standard ones, indicating that our standard value can be taken to represent low-conductivity behaviour in general, but the results for the high conductivity model are more interesting. They are shown in Fig. 13, exhibiting several new features as compared to the previous models. The crystallization spurts and CO bursts are much more pronounced, and they occur erratically right from the beginning -- often but not always in pairs. Much more heat is now conducted into the amorphous ice, raising its temperature high enough to prevent the CO ice from approaching the crystallization zone. In the typical situation there are many meters of amorphous ice without any CO ice just below the transition level. This together with the fact that now the latent heat can be used more efficiently to propagate the front downward, makes the crystallization spurts penetrate by 10--20 meters until they are stopped by the CO-ice cold trap. During the very first revolutions, the CO bursts even reach values comparable to the perihelion H$_2$O production rate. Subsequently the amplitude of the CO peaks decreases and the mean production rate reaches its standard value. The long-term average crystallization rate of this model is similar to that of the standard model, but the CO ice layer becomes thicker. \titleb{Capture models} It is well known that the orbits of Jupiter family comets are usually very unstable due to severe perturbations experienced at close encounters with Jupiter. In particular, the perihelion distance sometimes jumps by several AU on such occasions. Thus, for example, some of the comets under observation were discovered shortly after being captured from orbits with much larger values of $q$. Already within several centuries from the present time, most Jupiter family comets have orbits that differ markedly from their present ones. For statistics and discussion of these features, see Carusi and Valsecchi (1987), Tancredi and Rickman (1992). Thus we decided to have a first look at how the crystallization behaviour may be affected by such orbital changes. Two models were thus considered. In each of them the comet starts by performing 20 revolutions in the pre-capture orbit. One model has $q_{ini}=3.5$ AU and $Q_{ini}=8$ AU, and the other has $q_{ini}=5.5$ AU and $Q_{ini}=9$ AU during this phase. After the 20:th perihelion passage the orbit is suddenly changed into $q=1.5$ AU and $Q=6$ AU, as in the standard model, at the time the comet reaches a heliocentric distance $r_h=6$ AU post-perihelion. After this capture the comet is followed for 10 additional revolutions. All physical parameters characterizing the nucleus are the same as in the standard model. \begfigwid 12 cm \figure{14}{The evolution of a capture model with $q_{ini}=3.5$ AU and $Q_{ini}=8$ AU. After 20 revolutions in this orbit, the comet jumps to an orbit with $q=1.5$ AU and $Q=6$ AU, and performs 10 more revolutions. The upper plot corresponds to the H$_2$0 flux (full-drawn curve) and CO flux (dotted curve). In the middle plot we find, from top to bottom: total radius of the nucleus, radius where the amorphous ice is half the original value (dotted curve), radius of the uppest CO layer, radius of the lowest CO layer. The lower plot corresponds to the evolution of the surface temperature. The lower $x$-axes in all plots are expressed in years, while the upper ones in number of revolutions starting at aphelion. } \endfig Fig. 14 shows the results for the $q_{ini}=3.5$ AU model. We note the initial occurrence of crystallization spurts and CO bursts reaching the same amplitude as seen in the previous low-$q$ models. After five revolutions the bursts disappear, but the crystallization rate is remarkably high, i.e., $\simeq 1$ meter per year, indicating that a 3-km radius nucleus in this orbit might crystallize completely within several thousand years, unless the front slows down later on. The rate of CO outgassing decreases very slowly, and over the centuries plotted it exceeds the perihelion H$_2$O outgassing rate by a factor 100. Does this mean that the activity of Jupiter family comets with $q \ga 3$ AU is CO-driven rather than H$_2$O-driven? A fair answer to that question requires that another, more favourable insolation geometry be considered for estimating the H$_2$O flux. Assuming that a large fraction of the nuclear surface area passes near the subsolar point during the rotation of the nucleus, thermal modelling results for such conditions prevailing on pure H$_2$O nuclei (Froeschl\'e and Rickman 1986) can be applied, and one then finds that the H$_2$O outgassing rate near perihelion at $q=3.5$ AU might be comparable to, or even somewhat higher than, our CO production rate. We thus have an indication that large-$q$ Jupiter family comets might have an important part of their activity spread out along the orbits due to CO outgassing from the deep interiors of the nuclei. Upon capture, the depth of the crystallization zone is already 300 m, and the subsequent evolution cannot exhibit any major CO bursts. In fact, the evolutionary pattern looks very similar to that of the standard model at a similar stage: subdued, minor CO bursts occur erratically near some perihelion passages. We conclude that the comet basically behaves in accordance with the burial depth of the crystallization zone for any given set of physical parameters, independent of which previous orbital evolution has led to this state. \begfigwid 12 cm \figure{15}{The evolution of a capture model with $q_{ini}=5.5$ AU and $Q_{ini}=9$ AU. After 20 revolutions in this orbit, the comet jumps to an orbit with $q=1.5$ AU and $Q=6$ AU, and performs 10 more revolutions. The explanations are the same as in Fig. 14. } \endfig The second capture model is illustrated by Fig. 15. Initially this behaviour is similar to the previous one, exhibiting three major crystallization spurts with impressive, associated CO bursts. After this the crystallization proceeds at a very gentle rate, compared with the previous model -- whereas in that case the surface temperature near perihelion was typically in the range where rapid crystallization ensues, the temperature now stays significantly lower. As a consequence, the CO outgassing rate settles on a very low level which, nonetheless, is more than five orders of magnitude higher than that of H$_2$O production. Again, however, one should note that the H$_2$O production rate for a nucleus with active spots near the subsolar latitude would be dramatically different from that of our model. The gentle oscillations correspond to slight crystallization events triggered by heat flowing periodically from the surface. The burial depth of the crystallization zone stays at $\simeq 50$ m, comparable to the case of the standard model after only two revolutions. Subsequent capture then produces a major crystallization spurt, and for several revolutions the comet outgasses substantially more CO than during later phases. Crystallization proceeds as in the standard model and would probably slow down, if the computation were continued further in time. \titlea{ Discussion and Conclusions} Our model differs in important regards from those previously published, as seen in Sect. 3.1, where we also found correspondingly, significant differences of results when treating similar cases. We have improved on the earlier models by accounting for the state of near saturation attained by the vapour inside the nucleus and by explicitly including the erosional velocity of the surface. In addition, we have improved on Prialnik's model by allowing for condensation of CO ice below the crystallization front, and we have improved on Espinasse's model by dealing correctly with the energetics of gas release and by applying a numerical technique that adequately copes with the PDE's at hand. The shortcomings of our model become apparent when considering the very large CO vapour densities obtained near the crystallization zone at times of active phase transition (Fig. 4b). We then get a typical Knudsen number of $\sim 0.2$, indicating that the gas flow is somewhat outside the Knudsen regime. We are still far from a classical hydrodynamic flow, however, so some slight improvement would be gained by adopting an interpolation scheme (Espinasse {\it et al.} 1991). More important is the fact that our assumption of mechanical stability of the porous matrix is likely to break down, since the expected tensile strength of cometary material (Tauber and K\"uhrt 1987; Sekanina 1983) is not large enough to withstand the very large pressure gradient of the CO vapour diffusing from the crystallization zone toward the surface. As a result, material yielding is to be expected, whereby the initially very narrow interstices of the fluffy grain aggregate are replaced by wider flow channels. Hence the escape of the vapour could possibly occur on a much shorter time scale than we have modelled. We will return to this point below, when discussing future developments. As far as physical parameters are concerned, our basic improvements on earlier models are: 1) the allowance for an important or even dominant non-volatile component of the solid matrix; 2) the representation of this matrix as an aggregate of micron-sized core-mantle grains, implying an improved formula for the resulting conductivity of the ice-dust mixture and a correspondingly small pore diameter, in accordance with ideas of cometary formation based on the chemistry of the presolar cloud (Greenberg and Hage 1990); and 3) the adoption of a very low thermal conductivity of the amorphous ice mantles, based on Kouchi {\it et al.}'s (1992) results. Several other parameters should still be considered as highly uncertain and remain to be explored. The thermal conductivity of CO ice is such a parameter, since it condenses as an outer mantle enshrouding the initial grains and thus the bulk conductivity of the CO ice region should be dominated by $K_{CO}^i$ just like that of other regions is dominated by $K_w$. Moreover, in this paper we adopt $\gamma = 2/3$ for the dependence of conductivity on porosity, but our picture of the matrix as a loose aggregate of grains with relatively small contact areas would be more consistent with a much smaller $\gamma$ -- e.g., $\sim 0.1$ as in our variant model. One might also, preferably, consider the high end of our porosity interval in view of current indications on cometary densities (Greenberg and Hage 1990, Rickman {\it et al.} 1992), and together with a small $\gamma$ this implies a preference for models with very low conductivities in both the crystalline crust and the amorphous core. When considering such improvements, however, we have even more close to take into account the influence of pressure-induced yielding in the crystalline crust on the thermal conductivity. The chemistry of our model should only be viewed as an example with respect to considering CO alone as the very volatile trapped specie. It is reasonable to consider CO as a first choice, since the in-situ analysis of the P/Halley coma indicated a production rate of this molecule from the nucleus amounting to several percent of that of H$_2$O (Eberhardt {\it et al.} 1987). A quantitative fit to this observation is outside the scope of this paper, since P/Halley is moving in an orbit quite different from those of Jupiter family comets, and it may have suffered erosion down to $\sim 1/2$ of the initial radius (Hughes 1985; Rickman 1989). We note, however, that CO$_2$ appeared to have a comparable production rate from the P/Halley nucleus (Krankowsky {\it et al.} 1986; Combes {\it et al.} 1986), and hence should be another first-rank candidate for a trapped volatile. So far the only indication about the relative abundances of CO and CO$_2$ in Jupiter family comets comes from HST observations of P/Hartley 2 (Weaver {\it et al.} 1993), which provided a detection of CO$_2$ but not of CO. Eventually, the modelling of Jupiter family comets should be linked to that of Oort cloud comets in order to provide more stringent limits to the initial temperature. Our choice of 30 K should be regarded as a first example only, and it should be noted that higher temperatures at some early stage would reduce the amount of CO ice condensed below the crystallization front and thus speed up the phase transition. The actual value is probably set by the conductivity of the amorphous ice (Haruyama {\it et al.} 1993) and the relative contribution by silicates to the nuclear material, which governs the rate of radiogenic heat production per unit mass. Concluding, on the behaviour of Jupiter family comets from the present results, we first note the broad similarity between our standard and variant models. Within the framework of our assumption of an initially CO-rich amorphous ice and an important dust component, the complete crystallization of a sizeable nucleus with an initial radius of several km is estimated to take $\sim 10^4$ years, a time scale at least as long as a typical dynamical visit into the observable Jupiter family (Lindgren 1992) and longer than the expected active lifetime (Tancredi 1993). This means that Jupiter family comets with our assumed properties should still retain their CO, although in most cases buried deep below the nuclear surface. The observation of CO outgassing rates from Jupiter family comets by UV spectroscopy is an interesting challenge for the future, according to our results, since at least the youthful objects (i.e., recently captured ones) should exhibit near-perihelion rates $\sim 10^{24} - 10^{25}$ molecules s$^{-1}$ for each km$^2$ of unmantled surface area (the effect of a dust mantle on the CO outgassing rate remains to be investigated). In fact, even near aphelion the unmantled CO outgassing flux should be no less than $10^{18}$ molecules m$^{-2}$ s$^{-1}$, which is enough to entrail grains of dust and/or H$_2$O ice of radii up to 10 $\mu$m from a 3-km radius nucleus (Rickman {\it et al.} 1990). Therefore, Jupiter family comets that have active regions on their nuclei insolated near aphelion could exhibit faint comae due to the scattering of sunlight off such grains, thus casting doubt on nuclear photometry that does not have the spatial or temporal resolution required to separate the bare nuclei from these comae. From our pre-capture results with $q = 5.5$ AU we tentatively assert that a permanent coma of P/Schwassmann-Wachmann 1, such as detected by Jewitt (1991), is consistent with the CO outgassing predicted by our model. Moreover, uv spectroscopy of quasi-Hilda type comets (e.g., P/Gehrels 3, P/Smirnova-Chernykh or P/Helin-Roman-Crockett) might reveal CO outgassing rates that compete with, or even dominate, the H$_2$O production rates. We do not obtain any CO bursts in the $q = 5.5$ AU pre-capture model that can be compared with the outbursts of P/Schwassmann-Wachmann 1. It remains to be investigated whether for any reasonable spin properties of the nucleus our model might be consistent with the scenario envisaged by Jewitt (1991), treating the local insolation patterns for separate regions on the nucleus. In general, at least for small-$q$ orbits, the spin modulation of the insolation is not expected to have any influence on the crystallization or CO outgassing behaviour, since the front should nearly always be situated far below the surface layer affected by spin-modulated heat flow. When confronting the model results with observations of real comets, however, one should take into account: 1) that real nuclei of Jupiter family comets are usually largely covered by dust mantles, which have a significant influence on the pattern of surface heat flow and erosion; 2) that the nuclei are probably non-spherical with significant surface topography, perhaps partly as a result of the uneven erosion due to uneven dust mantling. Furthermore, if there is lateral variation of the crystallization behaviour due to such inhomogeneities, lateral heat and vapour flows will ensue with as yet unexplored consequences. It must therefore be kept in mind that, in addition to the above-mentioned improvements of our model that appear necessary to reach physical consistency (i.e., treatment of the mechanical response of the material to the vapour pressure gradients building up in the interior), a truly realistic model would have to take dust mantling and surface topography into account. In view of this, our present model cannot claim to account for quantitative details like instantaneous CO/H$_2$O production ratios measured at specific points in the orbits of observed comets, or the amplitudes of outbursts detected by means of photometry or imaging of individual comets. The above conclusions, general as they may be, are the only ones that we find robust enough to emphasize at present, while the remainder of our results includes many interesting, though still tentative, indications. After developing the model in the directions mentioned, it will be worthwhile to look into the quantitative production rates of different gases for comparison with observations and thus to find constraints on the initial composition of the nuclei. At that stage it will also be desirable to undertake the ambitious program of a follow-up of Jupiter family comets during their entire lifetimes in ever-changing orbits, thus to find how crystallization proceeds throughout the nuclei and what stages of activity the comets might hence experience in the long run until their possible, final demise. \begref{References} \ref Bar-Nun A., Herman G., Laufer D., Rappaport M.L., 1985, Icarus 63, 317 \ref Campins H., A'Hearn M., McFadden L., 1987, ApJ 316, 847 \ref Carusi A., Valsecchi G.B., 1987, in: Interplanetary Matter, Z. Ceplecha, P. Pecina (eds), Publ. Astron. Inst. Czechosl. Acad. Sci. No. 67, p. 21 \ref Carusi A., Kres\'ak L., Perozzi E., Valsecchi, G.B., 1987, in: Interplanetary Matter,Z. Ceplecha, P. Pecina (eds.), Publ. Astron. Inst. Czechosl. Acad. Sci. No. 67, p. 29 \ref Clayton J., Giauque W., 1932, J. Am. Chem. Soc. 54, 2610 \ref Combes M., Moroz V.I., Crifo J.F., Lamarre J.M., Charra J., Sanko N.F., Soufflot A., Bibring J.P., Cazes S., Coron N., Crovisier J., Emerich G., Encrenaz T., Gispert T., Grigoryev A., Guyot G., Krasnopolsky V., Nikolsky Y., Rocard F, 1986, Nature 321, 266 \ref Crifo J., 1987, A\&A 187, 438 \ref Delsemme A.H., 1985, in: Ices in the Solar System, J. Klinger, D. Benest, A. Dollfus, R. Smoluchowski (eds), NATO ASI C:156, p. 367 \ref Delsemme A.H., Miller D.C., 1971, Planet. Space Sci. 19, 1229 \ref Eberhardt P., Krankowsky D., Schulte W., Dolder U., L\"ammerzahl P., Berthelier J.J., Woweries J., Stubbemann U., Hodges R.R., Hoffman J.H., Illiano J.M., 1987, A\&A 187, 481 \ref Espinasse S., 1989, Modelisation du comportement thermique et de la diff\'erenciation chimique des noyaux de com\`etes, Doctoral thesis, Univ. of Grenoble \ref Espinasse S., Klinger J., Ritz C., Schmitt B., 1989, in: Physics and Mechanics of Cometary Materials, J. Hunt and T.D. Guyenne (eds.), ESA SP-302, p. 185 \ref Espinasse S., Klinger J., Ritz C., Schmitt B. 1991. Icarus 92, 350 \ref Froeschl\'e C., Rickman H., 1986, A\&A 170, 145 \ref Ghormley J.A., 1968, J. Chem. Phys. 48, 503 \ref Greenberg J.M., 1982, in: Comets, L. Wilkening (ed.), Univ. Arizona Press, Tucson, p. 131 \ref Greenberg. J.M., Hage J., 1990, ApJ 361, 260 \ref Haruyama J., Yamamoto T., Mizutani H., Greenberg J.M., 1993, Subm. to J. Geophys. Res. \ref Hughes D., 1985, MNRAS 213, 103 \ref Jewitt D., 1991, in: Comets in the Post-Halley Era, R.L. Newburn Jr., M. Neugebauer,J. Rahe (eds), Kluwer Acad. Publ., Dordrecht/Boston/London, p. 19 \ref Klinger J., 1980, Science 209, 271 \ref Kochan H., Feuerbacher B., Jo\'o F., Klinger J., Seboldt W., Bischoff A., D\"uren H., St\"offler D., Spohn T., Fechtig H., Gr\"un E., Kohl H., Krankowsky D., Roessler K., Thiel K., Schwehm G., Weishaupt U., 1989, Adv. Space Res. 9, 113 \ref Kouchi A., Greenberg J.M., Yamamoto T., Mukai T., 1992, ApJ, 388, L73 \ref Krankowsky D., L\"ammerzahl P., Herrwerth I., Woweries J., Eberhardt P., Dolder U., Herrmann U., Schulte W., Berthelier J.J., Illiano J.M., Hodge R.R., Hoffman J.H., 1986, Nature 321, 326 \ref Lindgren M., 1991, in: Asteroids, Comets, Meteors 1991, A. Harris, E. Bowell (eds), Lunar Planet. Inst., Houston, p. 371 \ref McDonnell J.A.M., Alexander W.M., Burton W.M., Bussoletti E., Evans G.C., Evans S.T., Firth J.G., Grard R.L., Green S.F., Gr\"un E., Hanner M.S., Hughes D.W., Igenbergs E., Kissel J., Kuczera H., Lindblad B.A., Langevin Y., Mandeville J.-C., Nappo S., Pankiewicz G.S.A., Perry C.H., Schwehm G.H., Sekanina Z., Stevenson T.J., Turner R.F., Weishaupt U., Wallis M.K., Zarnecki J.C., 1987, A\&A 187, 719 \ref Mekler Y., Prialnik D., Podolak M., 1990, ApJ, 356, 682 \ref Mendis. D., Brin G., 1977, The Moon 17, 359 \ref Patashnick H., Rupprecht G., Schuerman, D.W., 1974, Nature 259, 313 \ref Press W., Flannery B., Teukolsky S., Vetterling W., 1987, in: Numerical Recipes, Cambridge University Press, Cambridge. \ref Prialnik D., 1992, ApJ 388, 196 \ref Prialnik D., Bar-Nun A., 1987, ApJ 313, 893 \ref Prialnik D., Bar-Nun A., 1988, Icarus 74, 272 \ref Prialnik D., Bar-Nun A., 1990, ApJ 363, 274 \ref Rickman H., 1986, in: Comet Nucleus Sample Return Mission, O. Melita (ed), ESA SP-249, Noordwijk, The Netherlands, p. 195 \ref Rickman H., 1991, in: Comets in the Post-Halley Era, R.L. Newburn Jr., M. Neugebauer, J. Rahe (eds). Kluwer Acad. Publ., Dordrecht/Boston/London, p. 733 \ref Rickman H., Fern\'andez J., Gustafson B., 1990, A\&A 237, 524 \ref Rickman H., Festou M.C., Tancredi G., Kam\'el L., 1991, in: Asteroids, Comets, Meteors 1991, A. Harris, E. Bowell (eds), Lunar Planet. Inst., Houston, p. 509 \ref Rickman H., Tancredi, G., 1993, in: Workshop on the Activity of Distant Comets, W.F. Huebner, H.U. Keller, D. Jewitt, J. Klimger and R. West (eds.), Southwest Research Institute report, San Antonio, Texas, p. 182 \ref Schmitt B., Espinasse S., Grim R., Greenberg J.M., Klinger J., 1989, in: Physics and Mechanics ~ of ~ Cometary Materials, J. Hunt and T.D. Guyenne (eds), ESA SP-302, p. 65 \ref Schmitt B., Espinasse S., Klinger J., 1991, in: Abstracts of the Asteroids, Comets and Meteors 1991, Flagstaff, Arizona \ref Sekanina Z., 1983, AJ 88, 1382 \ref Smoluchowski R., 1981, ApJ 244, L31 \ref Smoluchowski R., 1982, J. Geophys. Res. 87, Supp. A, 422 \ref Sykes M., Walker R. 1992. in: Asteroids, Comets, Meteors 1991, A. Harris, E. Bowell (eds), Lunar Planet. Inst., Houston, p. 587 \ref Tancredi G., 1993, Subm. to Planetary and Space Science \ref Tancredi G., Rickman H., 1991, in: Proc. of IAU Symp. 152: Chaos, Resonance and Collective Dynamical Phenomena in the Solar System, S. Ferraz-Mello (ed.), p. 269 \ref Tauber F., K\"uhrt E., 1987, Icarus 73, 499 \ref Washburn E., 1928, in: International Critical Tables, vol. III, p. 210 \ref Weaver H., Feldman P., McPhate J., A'Hearn M., Arpigny C., Smith T., 1993, in: Abstracts of Asteroids, Comets and Meteors 1993, Belgirate, Italy. \ref West R., Hainaut O., Smette A., 1991, A \& A, 246, L77 \ref Yabushita S., Wada K., 1988, Earth, Moon and Planets 40, 303 \endref \bye % finish your text with this command