A two-dimensional model of crevasses formed by cometary activity

Planetary and Space Science 87 (2013) 46–65

Contents lists available at ScienceDirect

Planetary and Space Science journal homepage: www.elsevier.com/locate/pss

A two-dimensional model of crevasses formed by cometary activity Barbara L. Krebl, Norbert I. Kömle n Space Research Institute, Austrian Academy of Sciences, Schmiedlstrasse 6, A-8042 Graz, Austria

art ic l e i nf o

a b s t r a c t

Article history: Received 5 February 2013 Received in revised form 30 July 2013 Accepted 12 August 2013 Available online 30 August 2013

Current images clearly show the surface of cometary nuclei to be highly structured as well as the active regions to be distributed unevenly and limited to certain areas, while the remaining surface stays dormant. The dominant part of the emitted volatiles seems to come from crevasses and ice pockets embedded in more or less non-volatile material, describing the scenario discussed in this paper and implemented by a numerical model. We consider a crack in the cometary surface, initially filled with 95% amorphous water ice containing inclusions of more volatile species (e.g. CO2, CO) and 5% crystalline water ice, sometimes covered by a dust matrix. The thermal evolution of such a system at low temperature conditions in response to diurnal irradiation by the Sun has been studied. It is found that there is little interaction between the progress of the crystallisation front and the erosion process due to sublimation at the surface. Moreover, our results indicate that in most circumstances the gas flux due to surface sublimation is orders of magnitude higher than the additional gas release caused by crystallisation of the amorphous ice. In view of these results the often heard suggestion that the crystallisation process of the originally amorphous ice could be a major reason for the frequently observed cometary brightness outbursts appears at least debatable. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Cometary surface Crevasses Ice crystallisation Sublimation Dust mantle evolution Numerical simulations

1. Introduction Already the first close fly-by missions at a comet nucleus (the Giotto and VEGA missions to P/Halley in the year 1986) revealed that gas and dust emission from the irradiated surface of a cometary nucleus is a highly complex phenomenon. Typically, the emissions are not evenly distributed over the irradiated surface, but rather are restricted to relatively small parts of it (a few percent of the total surface area). Later on, missions to other comets (like Deep Space 1 to comet Borelly in 2001, Stardust to comet Wild-2 in 2004, Deep Impact to comet Tempel-1 in 2005 and EPOXI to comet Hartley-2 in 2010) have revealed similar results (Boice et al., 2002; Tsou et al., 2004; A'Hearn et al., 2005, 2011). Moreover, dust emissions from active regions were observed to take the shape of strongly focussed “dust jets” with little tendency to expand laterally across the surface. All these observations point towards the conclusion that cometary activity is not a pure surface phenomenon and that many of the gas and dust emissions come from the interior of the

n

Corresponding author. Tel.: þ 43 316 4120651; fax: þ 43 316 4120690. E-mail address: [email protected] (N.I. Kömle).

0032-0633/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.pss.2013.08.006

comet nucleus and are only triggered by the solar radiation absorbed near the surface. Moreover, the strong lateral heterogeneity suggests that volatiles stored in the interior flow through pre-existing pathways and are expelled from the nucleus through crevasses and holes. These findings call for a new generation of thermal models for comet nuclei. Most of the currently published models describing the energy and mass balance of the comet nucleus are only onedimensional in the sense that they describe a fast spinning spherical nucleus on an excentric orbit around the Sun. Examples are the models published by Espinasse et al. (1991) and by Prialnik and BarNun (1990). Such global models frequently assume implicitly that gas and dust emission takes place homogeneously and changes only with time, depending on the distance of the comet from the Sun. Variations caused by the rotation of the nucleus are usually not considered, nor are spatial inhomogeneities taken into account. At this point we want to refer to the numerous analyses and calculations done for the latest comet mission Rosetta, which is scheduled to rendezvous with its target comet 67P/ChuryumovGerasimenko, on May 22, 2014. Overall portraits of 67P can be found in de Sanctis et al. (2006) and Lamy et al. (2007), studies on productivity and activity are the topics of Rosenberg and Prialnik (2010) and de Almeida et al. (2009). Dust mantel formation is modelled in Rosenberg and Prialnik (2009), while the modelling of

B.L. Krebl, N.I. Kömle / Planetary and Space Science 87 (2013) 46–65

orientation and structure of 67P is studied in Kossacki and Szutowicz (2006). Davidsson and Gutierrez estimated 67P's properties from non-gravitational force modelling in Davidsson and Gutierrez (2005). Further work to be referred to in this connection are de Sanctis et al. (2010) and Kelley et al. (2009). Examples for studies on the Rosetta mission itself are, e.g. Biele and Herfort (2012) and Ulamec et al. (2006).

47

In the present work we try to include spatial inhomogeneities by developing a two-dimensional model of the comet nucleus surface, taking into account the presence of cracks and holes in the surface. We study the evaporation of ices residing inside such cracks and the possible interaction of processes associated with the warm-up of these ices, like amorphous/crystalline phase change, ice sublimation, and the build-up of a porous dust mantle covering the ices.

Fig. 1. (a) Whipple's icy conglomerate model, (b) fractal aggregate model, (c) primordial rubble pile model, (d) icy-glue model (adapted from Weissman et al. (2004), page 340 in the book Comets II).

Fig. 2. Crevasses and holes on a cometary surface, as illustrated in Houpis (1990) (p. 185).

48

B.L. Krebl, N.I. Kömle / Planetary and Space Science 87 (2013) 46–65

The following pages contain a short general introduction into the topic as well as a more detailed theory section, describing the mathematical formalism on which the involved physics is based. These sections are followed by the implementation of the numerical model and the results it yielded. The text concludes with a discussion of the presented results and a listing of conclusions drawn from them. 1.1. Crevasses as source regions of cometary activity The current view on the general build-up of cometary nuclei is based on the assumption that comets have formed during the accretion phase in the outer regions of the solar system. Their basic building blocks are planetesimals in the 10–100 m range, which may have condensed and accreted first, when the temperature in the region was still to high to allow the condensation of gases. Later on, the cometary nuclei, which are typically bodies of several up to several tens of kilometres size, accreted by slow collisions of these planetesimals, leaving holes and cracks along the boundaries. These cavities were finally filled by condensing water vapour and gases of higher volatility, as CO2 and CO. Various types of such “rubble pile” models were suggested in the literature. Some of them are shown in Fig. 1. Here we mention in particular the so-called “icy glue” model proposed by Gombosi and Houpis (1986). This model predicts the presence of weak regions along planetesimal boundaries on the cometary surface, where the icy constituents should be concentrated and which therefore should be the preferred regions of “cometary activity”. Fig. 2 shows a schematic sketch of such a crevasse region on a cometary surface, with all features that might occur. Although it is known that cometary nuclei are quite diverse, there is a broad agreement that they contain three types of chemical species in a solid form, namely:

 Rocky components composed of silicates, which resemble





aggregates of interplanetary dust particles and the layers constituting, e.g. asteroids. These components are rather inert when they are irradiated by the Sun, except in extreme cases when the comet comes very close to the Sun or crashes into the solar atmosphere. They do not exhibit significant phase changes like melting or evaporation when the comet passes through its perihelion. Organic components composed of hydrocarbons, which sometimes cover the mineral grains. They may have a consistence similar to tar or paraffines and should behave “semi-volatile” under intense solar radiation. This means that they could partly melt or evaporate, and thus they could act as a “glue” between particles and therefore contribute to an increase of cohesiveness. Icy components composed of water ice, but also more volatile frozen gases, like CO2, CO, etc.

Although the existence of these species in comets has been largely confirmed by various space missions, it turns out more and more that their distribution inside and near the surface of the nucleus is quite different than originally thought. It is far from being homogeneous, at least near the surface. It appears that only small parts of the cometary surface become “active” in the sense that they start to emit gases and particles when heated up by the Sun as the comet approaches its perihelion. Rather, gas and dust emissions are restricted to a small percentage of the illuminated surface, the much larger part remains inert and does not show any measurable activity. The activity is concentrated along restricted areas, which may delineate surface cracks or holes.

1.2. Dust mantle build-up Formation and growth of a dust mantle may vary widely on a real comet. It depends on the initial dust content of the ice and on the topological structure of the ice/dust agglomerates, for which there are in principle two possible configurations: (i) in the form of a “dirty iceball”, where ice forms the basic structure, with single dust particles embedded, or (ii) in the form of an “icy dirtball”, where the topology is the other way round. Aside of this, the formation and growth of a dust mantle mainly depends on the local gravity acceleration and on the particle size distribution. In our actual model we assume two limiting cases: (i) A No Dust Model (NDM), which implies that any enclosed dust particle is dragged away by the gas stream immediately after it has been liberated from the surrounding ice. In this case no dry porous dust mantle will form on the surface and the sublimed water vapour is always free to flow into the low pressure coma. (ii) A Dust Model (DM), where the comet nucleus surface remains at its original level despite the action of ice sublimation. In this case both the water vapour formed at the sublimation front and the gases liberated in deeper layers (for example due to the crystallisation process of amorphous ice) cannot freely escape, but are forced to diffuse through the porous dust mantle before they can finally escape into the coma. Moreover, as we will see below, such a dust mantle forms a thermal barrier for the penetrating heat wave, since its bulk thermal conductivity is expected to be much lower than that of the underlying ice/dust mixture.

2. Theoretical basis In this section the most important differential equations and formulae for the material properties used to describe heat and mass transfer in the cometary surface material are given. This includes source and loss terms due to phase changes (sublimation) and the release of trapped gases associated with the crystallisation of amorphous ice. The ice filling the crevasse is assumed to be porous and porosity is given as a prescribed parameter, which modifies other material parameters like heat capacity and thermal conductivity by including a volume fraction factor. Generally, as illustrated in more detail in the next section, the computational domain consists of a free flow region (forming a part of the inner coma), and the domain inside the crack filled by the porous ice. The gas flow in the porous domain is described by the Brinkman Equations. In the free domain the flow regime is assumed to be turbulent and is calculated by the k–ϵ model (with k describing the turbulent kinetic energy and ϵ describing the dissipation rate of the turbulence), which is applied to a Reynolds-averaged Navier–Stokes (RANS) equation. This method is a standard approach for describing turbulent gas flows in various regimes (see e.g. Kuznetsov and Sheremet, 2010). Sublimation is implemented by calculating the surface pressure developing due to the evaporating gases, using the saturation pressure, the sublimation rate and the thermal velocity. The equations to implement the process of crystallisation are those derived by Schmitt et al. (1989), Ghormley (1968), and Klinger (1980, 1981). Table 12 in the Appendix contains a descriptive and quantitative listing of variables and constants used in the model.

B.L. Krebl, N.I. Kömle / Planetary and Space Science 87 (2013) 46–65

The text describing the heat transfer, Brinkman and turbulent flow summarises the most important points regarding the s implementation at hand as given in the COMSOL Multiphysics 1 theoretical documentary of the physics modules and the subsequent sources cited therein.

2.1. Heat transfer All forms of heat transport (conduction, convection and radiation) are present in different degrees. Heat transfer by radiation is only simulated by using the blackbody emission term modified by an IR-emissivity factor in the heat flux boundary condition. The following general heat transfer equation incorporates solid state heat transfer via the solid parts of the porous matrices (separate ice and dust matrices in our case) and heat transfer via the present gases in the pores: ðρC p Þeq


 ∂T þ ρC p u  ∇T ¼ ∇  αeq ∇T þ Q HT : ∂t

ð1Þ

Here ρ denotes the density of the solid, Cp is the heat capacity by mass at constant pressure, T is the temperature, u is the velocity of the fluid, and α is the thermal conductivity. QHT is a collective term representing all heat sources. The index eq stands for “equivalent” and means that these parameters were modified to include porosity. As inside the porous ice the convective term containing u is usually negligible compared to the conductive and the heat source and loss terms, the heat transfer equation can be simplified by setting u ¼ 0. Equivalent to the porosity of the material, the volume fraction θ is used: θ ¼ 1ψ , where ψ is the porosity. The volumetric heat capacity at constant pressure (ρC p ) and the thermal conductivity are then modified to account for porosity: ðρC p Þeq ¼ θs ρs C ps þ θf ρC p ;

2.2. Fluid flow 2.2.1. Brinkman equations The Brinkman equations describe the velocity and pressure conditions for a flow of certain viscosity through a porous medium, characterised by its porosity and permeability. The first equation handles mass conservation with Q as the mass source term
 ∂ðψρÞ þ∇ ρu ¼ Q ; ∂t

and

αeq ¼ θs αs þ θf αf :

ð4Þ

The indices s and f stand for “solid” and “fluid”. Surface boundary condition: There is only one heat flux to implement as a boundary condition, namely the heat flux at the surface. It contains the Stefan–Boltzmann law (blackbody radiation), the energy loss due to sublimation and the incoming energy due to irradiation from the Sun. The complete energy balance at the nucleus' surface can then be written as q ¼ ð1AÞ 

qSun d

2

 cos ζ ps ðT Þ  H sub 

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi M H2 O em  sT 4 : 2π  R  T

ð5Þ

q is the heat flux, A is the albedo, qSun is the solar constant, d is the heliocentric distance, cos ζ is the incident angle, ps(T) is the saturation pressure at temperature T, Hsub is the latent heat of sublimation, M H2 O is the molar mass of water, R is the universal gas constant, em is the emissivity of the surface and s is the Stefan– Boltzmann constant. All boundaries other than the nucleus surface are thermally isolated. 1

Finite Element software which was used for the numerical implementation. For more detailed information refer to the webpage www.comsol.com.

ð7Þ

and is related to (6) via its stress and force terms. ∇P is the pressure gradient, η the viscosity of the fluid, I the stress tensor, κ the permeability of the porous material and F sums up volume and gravity forces. Because incompressible flow is assumed, the first Brinkman equation simplifies to the stationary form

ρ∇  u ¼ Q :

ð8Þ

Also the inertial term – ðu  ∇Þu=ψ on the left-hand side of Eq. (7) is neglected. 2.2.2. Turbulent flow To calculate turbulent flow in the free domain, the Reynoldsaveraged Navier–Stokes (RANS) model for incompressible, Newtonian flow is used. The turbulent character of the flow is described by the k–ϵ model. The essential equations needed are

 the RANS-equation: h i
 ∂U þ ρU  ∇U þ ∇  ρu′  u′ ¼ ∇P þ ∇η ∇U þ ð∇UÞT þ F; ∂t

ρ

where ð3Þ

ð6Þ

while the second one describes a momentum balance   ρ ∂u u þ ðu  ∇ Þ ¼ ∇P ψ ψ ∂t      2  1 η η ∇u þ ð∇uÞT  ηð∇  uÞI  þ Q u þ F þ∇  3 ψ κ

ð2Þ

θs þ θf ¼ 1;

49



ð9Þ

where U is the sum of averaged and fluctuating velocity parts and u′ is only the fluctuating part of the velocity, and the Reynolds stress tensor, which is the relation between the RANS and the k–ϵ model h i 2 ð10Þ ρðu′  u′Þ ρk ¼ ηT ∇U þ ð∇UÞT : 3

The latter introduces two new quantities: ηT , called the turbulent or eddy viscosity and k, the turbulent kinetic energy. Together with the dissipation rate of the turbulence energy ϵ, one obtains the basic equations of the k–ϵ turbulence model. It treats the fluctuating part of the velocity as a product of cumulation and dissipation of turbulent kinetic energy. Together they compose a function for the turbulent viscosity 2

k

ηT ¼ ρC η : ϵ

ð11Þ

C η is one of the five empirical model constants, given in Table 1 and explained in further detail, e.g. in Kuznetsov and Sheremet Table 1 Table of empirical constants for the k–ϵ turbulence model. Constant

Value

Constants

Value

Cη C ϵ1 C ϵ2

0.09 1.44 1.92

sk sϵ

1.0 1.3

50

B.L. Krebl, N.I. Kömle / Planetary and Space Science 87 (2013) 46–65

(2010). With the relations given above one obtains the following differential equations for the variables k and ϵ:   ∂k η ρ þ ρu  ∇k ¼ ∇ η þ T ∇k þ P k ρϵ ð12Þ ∂t sk

the crystalline ice present) and the function

  ∂ϵ η ϵ ϵ2 ρ þ ρu  ∇ϵ ¼ ∇ η þ T ∇ϵ þ C ϵ1 P k C ϵ2 ρ ∂t sϵ k k

as determined from measurements by Schmitt et al. (1989). Using the inverse of the factor 1.05  1013 in Eq. (19) yields the time necessary to complete crystallisation at temperature T:

with P k ¼ ηT



i 2 2 ∇u : ∇u þ ð∇uÞT  ð∇uÞ2  ρk∇u 3 3 h

ð13Þ

ð14Þ

2.3. Water ice sublimation The conditions prevailing at the surface of a cometary nucleus dictate that due to the low surface pressure water is not expected to exist in a liquid form. Rather we may assume that – if it has condensed in a low temperature environment as generally expected – it originally existed in the form of amorphous ice and crystallised as the temperature raised above  120 K for the first time. When the temperature approaches 200 K, the crystallised ice starts to evaporate. To describe this process mathematically, the so-called Clausius–ClapeyronEquation can be applied, which describes the saturation vapour pressure (ps evolving in a closed vapour space which is in contact with ice walls of temperature T): ð15Þ 12

with a ¼ 3:56  10 Pa and b ¼6141.667 K (see Fanale and Salvail, 1986). For the case of free sublimation into a vacuum this formula can be used to calculate the water vapour emission rate Z [kg m  2 s  1] from an icy surface with temperature T: rffiffiffiffiffiffiffiffiffiffiffiffi M Z ðT Þ ¼ ps  : ð16Þ 2π RT To relate the actual gas pressure at the surface to the gas emission rate due to sublimation one has to consider the mass balance at the surface together with the ideal gas law, as given, e.g. in Espinasse et al. (1991). For a porous crystalline ice, where the gases are allowed to diffuse through open pores driven by the pressure gradient, one arrives at the following expression for the surface pressure: pv ∂p ¼ ψ 2=3 D surf þZ ðT Þ ð17Þ kB T ∂x Hereby the velocity v is usually taken to be of the order of the sound speed, which is a mere function of temperature. D is the diffusion coefficient.

Ea ; kB T

ð19Þ

ð20Þ

where Ea is the activation energy for the onset of the phase change. The constant ðEa =kB Þ has a value of 5370 K. The heat and mass source terms depend on the crystallisation rate, which is a function of temperature. They are given by Eqs. (21) and (22).2 Q HC ¼ 1:05  1013  exp

5370 ρs Hcry ; T

ð21Þ

with Hcry as the latent heat of crystallisation. The used mass source is given by Eq. (22) and was taken from Prialnik (1991). Q MC ¼ X gas ρ

∂X a ; ∂t

ð22Þ

with Xgas being the amount of trapped gas in the amorphous ice matrix. The threshold temperature for the onset of crystallisation is around T ¼ 120 K. Above a temperature of 166 K it can be assumed that the crystallisation process is completed. Compared to sublimation, the phase change of the ice from the amorphous into the more common crystalline form shows two particular features, which need to be taken into account properly in the numerical implementation: (i) It is an exothermal process, i.e. heat is released during the crystallisation process, while during sublimation heat is consumed and (ii) the process is irreversible, i.e. in every mass element it takes place only once, no matter weather or not the ice is cooled down again after the transformation process has been completed. The peculiar consequence is that the crystallisation front can always only move into the comet and never propagate in upward direction, even when the comet nucleus cools down again to temperatures lower than the crystallisation temperature, when it moves away from the Sun on the outbound branch of the orbit. Furthermore, the thermal conductivity of the ices is calculated differently. For amorphous ice: C ices  vsound  Λ  ϱA : 4

ð23Þ

Cices is the specific heat capacity of the ices, vsound is the sound velocity, Λ is the mean free path and ϱA is the density of the amorphous ice. For crystalline ice, the thermal conductivity is calculated as follows: aαC T

ð24Þ

with aαC as a constant of 576 W m  1 K  2. The specific heat capacity, which is taken to be the same for amorphous and crystalline ice and is calculated by

2.4. Amorphous ice crystallisation The rate of crystallisation is given in the form of an exponential decay law ∂X a ¼ λðT ÞX a ; ∂t

T

t c ðT Þ ¼ 9:54  1014 exp

Eqs. (12) and (14) are describing the transport and production of k, the turbulent kinetic energy, while Eq. (13) handles dissipation, also by a transport equation, which has the same form as (12). Pk itself is the production term. The k–ϵ turbulence model assumes the Reynolds number to be high enough so that turbulence develops and equilibrium in boundary layers holds (∂k=∂t ¼ ϵ).

ps ¼ aeb=T

Ea

λðT Þ ¼ 1:05  1013  exp kB ;

C ices ¼ aC  T þ bC :

ð25Þ

ð18Þ

with X a being the amorphous ice mass fraction (values from 0–1, corresponding to the percentage of amorphous ice with respect to

2 The subscripts HC and MC stand for Heat source of Crystallisation and Mass source of Crystallisation, respectively.

B.L. Krebl, N.I. Kömle / Planetary and Space Science 87 (2013) 46–65

6

51

y {m}

5 3

ANTLE UST M

D

0,002 m

TALLIN CRYS

0,5 m

2

E ICE

4

0,5 m

1 0 y S ICE

PHOU

9m

AMOR

-1 -2 -3 -4 -5 -6

x

z

-7 -8 -9

x {m} 0

1

2

3

4

5

Fig. 3. (a) 3D-sketch of the model geometry for the dust models. (b) 2D-cut across the crevasse geometry showing in detail the domains and boundaries used in the COMSOL Multiphysics model. This figure shows the geometry for the Dust models. In the No Dust models the dust mantle is part of the turbulent domain. The individual domains and boundaries referred to in Tables 2–5 are indicated.

3. Numerical model s

The model, built and calculated with COMSOL Multiphysics 4.2., represents a symmetric rift in 2D (y is the symmetry axis) on a comet's surface, where the computational domain represents onehalf of the total cross section of a real crevasse, with x being the lateral and y the vertical coordinate (y¼ 0 corresponds to the comet nucleus' surface and  y is the depth coordinate counted from the surface). It is filled with amorphous and crystalline ice, embedded in non-volatile material and either covered by a dust mantle (cohesive dust matrix) of 2 mm initial thickness (DM – Dust Model) or exhibiting a free sublimating surface (NDM – No Dust Model). The model geometry and its dimensions are shown in Fig. 3 (a and b). The sublimation/dust front is not at even ground with the nucleus' surface, which had to be done in order to avoid mesh distortions. Since we assume symmetry across the x ¼0 plane, only the x 40 domain of the full crevasse needs to be computed. As the comet moves along its orbit, one side is illuminated by the Sun and thus heated, while the other is not, depending on the comets rotational period. Starting at approximately 120 K, the rift's amorphous ice layer crystallises and a certain amount of trapped

gas (in this case CO23) is released. Crystallisation is an exothermal process, which means that heat is generated in the respective regions. According to the Clausius–Clapeyron Equation (see formulas (15) and (16)) the resulting crystalline ice layer is then expected to sublimate around approximately 200 K. As sublimation is an endothermic process, it consumes energy in order to sustain the phase change. Note that they only take place at the fronts of their respective domains. This section will explain how the fluid flow, heat transfer and phase changes were modelled. For the numerical implementation a moving mesh with automatic remeshing (when grid distortions become too large) is used. By employing this feature it can be insured that during the computations the moving fronts in the domains are always well resolved and a reasonable mesh quality is maintained at any time. For a more detailed description of the model building process see Krebl (2012).

3 Besides CO, CO2 is one of the usual gas constituents known to exist in comets. Also the material parameters of CO2 are better known than those of CO.

52

B.L. Krebl, N.I. Kömle / Planetary and Space Science 87 (2013) 46–65

Table 2 Domain details.

Table 5 Boundary conditions for the turbulent flow in the coma.

Domain

Description

Dimensions (x [m]  y [m])

Used equations

1 2 (CD) 3 (DD) 4

Amorphous ice Crystalline ice Dust mantle Turbulent flow

0.5  9.0 0.5  0.5 0.5  0.002 5.0  6.0

Eqs. (9), (12)–(14)

1, 2, 3

Brinkman flow

0.5  9.502

Eqs. (1), (6), (7)

Table 3 Boundary conditions for the heat transfer through the porous media. Boundary condition (see Fig. 3) for boundary numbers

NDM boundaries

DM boundaries

Thermal isolation – heat flux ¼0 W/m2 Heat flux on the surface: Eq. (5) Heat sink of sublimation (only DMs): Eq. (5), sublimation term Heat source of crystallisation: Eq. (21) Symmetry axis

2, 11 and 12 6

2, 11, 12 and 13 8

– 4 1, 3

6 4 1, 3 and 5

Table 4 Boundary conditions for the Brinkman flow through the porous matrices. Boundary condition (see Fig. 3) for boundary numbers

NDM boundaries

DM boundaries

Pressure: Eq. (17) Walls: set to slip condition Symmetry axis

6 4 and 12 1, 3

– 4, 12 and 13 1, 3 and 5

3.1. Heat transfer Heat transfer is calculated in domains 1 and 2 (NDMs) and also in domain 3 in the DMs at an initial temperature of 90 K. The fluid flow is taken to be isothermal in the free domains. Therefore no heat transfer is implemented in these domains. As mentioned before, heat transfer via fluids is neglected in our model. The boundary conditions used are listed in Table 3.

Boundary condition (see Fig. 3) for boundary numbers

NDM boundaries

Inlet pressure: Eq. (17) Inlet velocity:

– 8 – velocity field of adjacent Brinkman interface 10, 16 10, 16 13, 14, 15 and 17 14, 15 and 17 5, 7 and 9 7, 9

Outlet pressure¼ 0 Pa Walls: slip condition Symmetry axis

DM boundaries

6 –

Turbulent flow: Domains 3 and 4 in the NDMs and only 4 in DMs are representing a small section of the coma above the comet's surface. The initial pressure and velocity conditions match the ones for the Brinkman equations. The initial turbulent energy k and dissipation rate ε are set to their default values as implemens ted in COMSOL Multiphysics . 3.3. Sublimation and crystallisation The velocities of the sublimation and the crystallisation front are calculated from their respective sublimation/crystallisation rates. The crystallisation rate resulting from Eq. (18) is already given as a percentage, which then is applied to the current volume. This percentage is the currently crystallised amount of ice and is used to calculate all other crystallisation related quantities, such as the mass flow of CO2 and the heat source of crystallisation. The velocity of the crystallisation front can only be negative and thus the front can never move upwards. In this way the irreversibility of the crystallisation process is taken into account. The fluid flow is only implemented in domain 2 (no multi-gas flow by mixing of water vapour and CO2 ) with an initial velocity of 0 m/s and is caused by a boundary gas source, implemented at boundary 4 and calculated by Eq. (22). The sublimation rate is used to calculate the velocity of the sublimation front – the mass it yields is taken as the respective mass of ice evaporated and is then recalculated into a volume using the density of the ice. It is implemented by the pressure boundary condition on boundary 6. 4. Results The calculated results were plotted to depict

3.2. Fluid flow The fluid behaves like an ideal gas – all model parameters and external calculations use the ideal gas law to describe the gas properties. However, as we assume incompressible flow, an average constant density, calculated on the basis of the expected scales of the pressure and temperature conditions, was used. Recondensation of the free gases is not considered. The permeability of the ice and dust layer is an average value calculated on the basis of previous heat transfer calculations and is mainly dependent on the thickness of the actual layer. The permeability is of the order 10  13 m2 for the dust and 10  17 m2 for the ice layers in both model varieties. Brinkman flow: It is implemented in the same domains as the heat transfer is. The initial pressure is set to the saturation pressure at the current temperature and the initial velocity is set to 0 m/s. Boundary conditions are listed in Table 4. A “wall” describes every outer boundary of a domain, where no explicit boundary condition has been set.

    

the temperature development in and at the surface of the rift, the movement of the sublimation and crystallisation front, the volume change of the dust mantle, the development of the gas flow, caused by sublimation and crystallisation and the sublimation and crystallisation rates and heat sources.

As irradiation intensity the solar intensity at a solar distance of 1.29 AU (perihelion distance of comet ChuryumovGerasimenko), 5.74 AU (aphelion distance) and at 2.7 AU is taken. A distance of 2.7 AU corresponds to the solar distance of the comet at the planned landing time of the Rosetta Lander Philae. Due to the 12.7 h rotation period of the comet (see Lamy et al., 2007), the irradiation intensity is modulated by a sinusoidal variation at daytime and set to zero at night time. The rotation axis is assumed to be perpendicular to the orbital plane. Starting with constant initial conditions for temperature (T ¼90 K) and gas flow (p ¼ water vapour saturation pressure at temperature T, v ¼0 m/s), we calculate the

B.L. Krebl, N.I. Kömle / Planetary and Space Science 87 (2013) 46–65

Fig. 4. Development of the temperature at different depth levels and at the sublimation and crystallisation front at a heliocentric distance of 1.29 AU.

evolution of the system over 40 rotation periods. After approximately 10 periods a quasi-steady state solution should have established, which may already be close to the real situation. The plots are arranged into panels containing the development of certain quantities at one of the calculated heliocentric distances. The notation of the plot titles is as follows: D or ND for “Dust” or “No Dust” distance – 1.29, 2.7 or 5.74 AU – and the variable displayed. Tables 6–9 are attached at the end of their corresponding plot series to summarise the data displayed in the graphs. The data given in the tables display the development after the 10th period, because the plots indicate a settling of the system to a quasi-equilibrium state. Furthermore, because the fronts stopped moving at a heliocentric distance of 5.74 AU, the fluid flow, sublimation and crystallisation plots for that distance are excluded as they hold no useful data due to non-existing phase changes. The velocity of the sublimation and crystallisation front and their heat sources had to be based on a finite extension of the rift in z-direction, thus the phase change and gas flow plots have to be regarded as calculations for a 3D-rift with a 20 m extension in z-direction. Two examples for this necessity are the velocity of the sublimation front and the heat source of crystallisation.

53

Fig. 5. The same as Fig. 4 but for a heliocentric distance of 2.7 AU.

The velocity of the sublimation front was calculated from the amount of sublimated ice, given by the sublimation rate. Since the sublimation rate is a rate per area rather than volume, a transformation by using the density of the ice had to be performed to be able to specify a velocity for the sublimation front in the negative y-direction. The heat source of crystallisation is given in W/m3 and depends on the amount of crystallised ice. In order to change the dimension to fit the requirements of a boundary heat source, the total amount of freed energy by crystallisation had to be divided by a defined surface area, needing of course a finite value for the extension in z-direction. 4.1. Temperature plots Figs. 4–10 are showing the development of the temperature with time at the surface and at the various fronts. They are also depicting the development of the temperature with depth at every third maximum starting with the 11th (first curve from right) and closing with the 38th (first curve from left) maximum of the cycle made up of 40 periods with a duration of 12.7 h each. 4.1.1. Rift without a dust cover See details of Figs. 4–6.

54

B.L. Krebl, N.I. Kömle / Planetary and Space Science 87 (2013) 46–65

Fig. 6. The same as Fig. 4 but for a heliocentric distance of 5.74 AU.

4.1.2. Rift with a dust cover See details of Figs. 7–10. 4.2. Front movement plots The change of the negative y-coordinate and the volume change of the dust mantle is shown in Figs. 11–15.

Fig. 7. Development of the temperature at the dust front (surface) for the calculated heliocentric distances 1.29, 2.7 and 5.74 AU.

4.3.2. Rift with a dust cover See details of Figs. 18 and 19.

4.4. Crystallisation and sublimation plots Figs. 20–23 are depicting the crystallisation rate and heat source for the calculated distances. The rate of sublimation and its heat source is presented in Figs. 24–27.

4.2.1. Rift without a dust cover See details of Figs. 11 and 12.

4.4.1. Crystallisation – rift without a dust cover See details of Figs. 20 and 21.

4.2.2. Rift with a dust cover See details of Figs. 13–15. 4.3. Gas flow plots

4.4.2. Crystallisation – rift with a dust cover See details of Figs. 22 and 23.

In this section we present the combined as well as the separate mass flow at the surface due to crystallisation and sublimation in Figs. 16–19.

4.4.3. Sublimation – rift without a dust cover See details of Figs. 24 and 25.

4.3.1. Rift without a dust cover See details of Figs. 16 and 17.

4.4.4. Sublimation – rift with a dust cover See details of Figs. 26 and 27.

B.L. Krebl, N.I. Kömle / Planetary and Space Science 87 (2013) 46–65

Fig. 8. Development of the temperature at different depth levels and at the sublimation and crystallisation front at a heliocentric distance of 1.29 AU.

55

Fig. 9. The same as Fig. 8 but for a heliocentric distance of 2.7 AU.

5. Discussion 4.5. Lowered thermal properties of the dust mantle – temperature and front movement The assumed values for the thermal conductivity and the heat capacity of the dust mantle of 0.1 W m  1 K  1 and 1300 J kg  1 K  1 respectively are valid for a ice/dust matrix as suggested by Rosenberg and Prialnik (2010) and Merk and Prialnik (2006). However, the thermal conductivity of a loose, dust-only aggregate in a vacuum is 1–2 orders lower than the one mentioned above (see references Koemle et al., 1996; Seiferlin et al., 1996). The thermal conductivity was lowered from a value of 0.1 Wm  1K  1 to 0.01 Wm  1K  1 and the heat capacity was reduced by approximately 40% to a value of 770 J kg  1 K  1 as used by Merk and Prialnik (2006). Figs. 28–30 are showing the temperature development in the crevasse and at its various fronts. The development of the dust mantle thickness and the front movements is presented in Figs. 31–33. The relevant parameters for this case are listed in Tables 10 and 11. 4.5.1. Temperature plots See details of Figs. 28–30. 4.5.2. Front Movement plots See details of Figs. 31–33.

Based on the results presented in Section 4, we now discuss the main findings from the simulations. 5.1. Thermal behaviour 5.1.1. No dust models With an initial temperature of 90 K, only the models calculated at a heliocentric distance of 1.29 AU develop a temperature pattern at the surface that stays roughly constant from period to period, as can be seen in Fig. 4. At 2.7 AU heliocentric distance, activity due to water sublimation is negligible as the temperature still rises from period to period and has not yet reached the sublimation point of 200 K at the end of the whole cycle. For a distance of 5.74 AU, the temperature at the surface only rises to 114 K and to 112 K at the crystallisation front over the 40 periods. Therefore it is to be expected that all activity is suppressed at this distance, as the temperatures at the various fronts reach neither the sublimation nor the crystallisation point. Looking at Figs. 11, 12, 14 and 15 this assumption seems to be confirmed. However, crystallisation continues for a distance of 1.29 and 2.7 AU – the temperature settles around 140 K in each case – of course the smaller the distance the sooner this point is reached. While at a distance of 2.7 AU the sublimation front movement decreases by more than 2 orders in comparison with the movement at 1.29 AU, the distance covered by the crystallisation front is only  1 m less at

56

B.L. Krebl, N.I. Kömle / Planetary and Space Science 87 (2013) 46–65

2.7 AU. As mentioned in Section 2.4, the thermal conductivity of amorphous and crystalline ice is calculated differently. Due to this circumstance the graphs depicting the temperature development with depth are also showing the boundaries and the growth/shrink of the two layers (crystalline and amorphous ice). In the crystalline domain the temperature decreases more smoothly with depth, while the temperature gradient in the amorphous domain is steeper. Figs. 11 and 12 are showing the front movement for several distances. While the sublimation front movement allows for recognising the periods and their maxima and minima, the movement of the crystallisation front is a smooth line.

5.1.2. Dust models Although the temperature at the surface is high compared to the temperature at the sublimation front and the dust mantles

Fig. 10. The same as Fig. 8 but for a heliocentric distance of 5.74 AU.

initial thickness is small (2 mm), it has a quenching effect on sublimation that sets in quite soon. These circumstances are controlled by the thermal properties of the dust mantle, such as its thermal conductivity (0.1 W/m K) and heat capacity (1300 J/ kg K). In Figs. 8 and 9 a vertical line around y ¼0 m can be recognised, representing the dust mantle and its steep temperature gradient, which indicates the influence of its thermal properties. As the dust mantle grows the temperature at the sublimation front decreases over the 40 periods. Since there is no mechanism built in to decrease the dust mantles thickness in the current model, it will grow and quench sublimation further until its total suppression. The behaviour at a distance of 2.7 AU is similar to that of the model without a dust mantle at the same distance. Temperatures

Fig. 11. Movement of the sublimation front at heliocentric distances of 1.29, 2.7 and 5.74 AU.

Table 6 Temperature ranges of the thermal plots between the 10th and 40th maximum. Model_distance

ND_129

ND_27

ND_574

D_129

D_27

D_574

T T T T

163–200 141–144 90–200

140–167 135–142 90–167

100–114 101–112 90–114

300–320 165–187 141–144 90–330

165–175 140–170 135–142 90–175

104–114 98–112 98–110 90–112

dust (K) sub. (K) cry. (K) (y) (K)

B.L. Krebl, N.I. Kömle / Planetary and Space Science 87 (2013) 46–65

57

Fig. 12. Same as Fig. 11 but for the movement of the crystallisation front.

Fig. 14. Movement of the sublimation front at heliocentric distances of 1.29, 2.7 and 5.74 AU.

Fig. 13. Volume change of the dust mantle calculated for the heliocentric distances 1.29 and 2.7 AU. The plot for a distance of 5.74 AU is constant at 2 mm.

at the sublimation front have not yet reached the sublimation point – the front movement plots show an increase of the dust mantle by only half a millimetre. Heating will continue until the sublimation point is reached – after that a similar behaviour to the model at a distance of 1.29 AU is expected and the temperature will fall. Crystallisation seems to be insensitive to the presence of the dust mantle. Front movement remains roughly the same in the models with and without a dust cover at the same solar distances. The values indicate that a suppression of sublimation at the surface in fact furthers crystallisation in deeper layers if the temperature is given time to rise to the several phase change temperature points (transformation from amorphous to crystalline ice around 140 K). This asynchrony between sublimation and crystallisation is already evident when investigating the temperature development in the models without a dust cover. The temperature increase at the crystallisation front is caused by the sublimation process itself. As long as the sublimation front is fed enough energy to sustain sublimation, the temperature will be roughly constant around 200 K and less energy is transported towards the crystallisation front. As soon as this “energy barrier” is dissolved due to diminishing sublimation, more energy gets transported towards the crystallisation front. Of course this also depends on the distance separating the fronts, which initially is 0.5 m in all models.

58

B.L. Krebl, N.I. Kömle / Planetary and Space Science 87 (2013) 46–65

5.2. Gas flow 5.2.1. No dust models The development of the singular mass flux due to crystallisation/sublimation and the total mass flux displayed in Figs. 16 and 17 for distances 1.29 and 2.7 AU matches the thermal behaviour of the models. Without a dust cover sublimation is free to develop and continues if not slowed down by recondensation or shadowing by the rift walls. While the mass flux due to sublimation decreases by 3 orders when the distance increases from 1.29 to 2.7 AU, the mass flux due to crystallisation roughly stays the same when developing towards a quasi-equilibrium state within the last

Fig. 15. Same as Fig. 14 but for the movement of the crystallisation front.

10–20 periods. As already discussed in the thermal section, at a distance of 2.7 AU, the temperature at the sublimation front will increase until there is enough energy for sublimation – the corresponding mass flux behaves in the same way. The total mass flux for each distance is the sum of the mass flux due to sublimation and due to crystallisation. Thus, the shape of the graph is determined by the dominant mass flux at the surface, which also depicts the composition of the gas leaving the comet at a certain distance. The mass flux due to sublimation prevails at both heliocentric distances in the NDmodels.

Fig. 16. Total and singular mass flow at the surface caused by sublimation and crystallisation at 1.29 AU.

Table 7 The data describe the total vertical loss/gain of the respective domain. The minus sign indicates that the front can only move downwards (irreversibility of crystallisation and no recondensation of sublimated vapours). Model_distance

ND_129

ND_27

ND_574

D_129

D_27

D_574

Sublimation (m) Crystallisation (m) Dust mantle thickness (m)

 0.23  1.8 –

 0.001  0.87 –

– – –

 0.038  1.7 þ 0.38

 0.0005  0.92 þ 0.0005

– – –

B.L. Krebl, N.I. Kömle / Planetary and Space Science 87 (2013) 46–65

59

5.2.2. Dust models As in the ND-models, the mass fluxes behave according to their temperature development graphs. For a distance of 1.29 AU the mass flux due to sublimation rapidly decreases with increasing thickness of the dust mantle. Table 8 shows a decrease in mass flux due to sublimation by 9 orders if a dust mantle is present at a heliocentric distance of 1.29 AU. When increasing the distance to 2.7 AU, this causes the mass flux to drop by another 2 orders. Like in the NDmodels, the sublimation point is not reached at the end of the 40 periods at a distance of 2.7 AU and thus will grow until it does. The temperature development and front movement plots already showed that the ongoings on the surface matter little to the process of crystallisation if a phase change temperature point is reached (140 K

in this case). The mass fluxes due to crystallisation are roughly the same for a distance of 1.29 and 2.7 AU and when reaching a quasiequilibrium state, which is the case in the last 10–15 periods. If a dust mantle is present, sublimation is strongly influenced and suppressed. Therefore the plotted total mass flux at the surface is dominated by the CO2 mass flux pattern due to crystallisation. Note that, in the case of the DMs, the total mass flux is the sum of the mass fluxes due to sublimation at the surface and the mass flux due to crystallisation at the sublimation front. Since a mixed gas flow is not implemented, the CO2-flux is only calculated in domain 2 and has its outlet at the sublimation front. However, a comparison of the mass flux due to sublimation at the sublimation front and at the surface has shown the effect on the flow pattern to be minimal.

Fig. 17. Same as Fig. 16 but for heliocentric distance of 2.7 AU.

Fig. 18. Total and singular mass flow at the surface caused by sublimation and crystallisation at 1.29.

Table 8 Singular and total mass flux at the surface due to crystallisation and sublimation from the 10th to 40th maximum. Model_distance

ND_129 2

Mass flow sub (kg/m s) Mass flow cry (kg/m2 s) Total mass flow (kg/m2 s)

ND_27 7

(22–24)  10 (7–13)  10  13 (22–24)  10  7

D_129  10

(15–95)  10 (5–11)  10  13 (15–95)  10  10

D_27  16

(0.5–2.5)  10 (6–11)  10  13 (6–11)  10  13

(12–80)  10  18 (5–11)  10  13 (5–11)  10  13

60

B.L. Krebl, N.I. Kömle / Planetary and Space Science 87 (2013) 46–65

Fig. 21. Same as Fig. 20 but for a heliocentric distance of 2.7 AU. The rate is given in volumetric percent crystallised per second based on the respective remaining amount of amorphous ice.

Fig. 19. Same as Fig. 18 but for a heliocentric distance of 2.7 AU.

Fig. 22. Crystallisation rate and heat source at a heliocentric distance of 1.29 AU. The rate is given in volumetric percent crystallised per second based on the respective remaining amount of amorphous ice.

5.3. Phase changes

Fig. 20. Crystallisation rate and heat source at a heliocentric distance of 1.29 AU.

5.3.1. Crystallisation Naturally the onset of crystallisation is shifted to a later time if the heliocentric distance is increased. Still, the rate of crystallisation and its heat source for a distance of 2.7 AU do not differ much from the rate and heat source at 1.29 AU when looking at

B.L. Krebl, N.I. Kömle / Planetary and Space Science 87 (2013) 46–65

61

Fig. 25. Same as Fig. 24 but for a heliocentric distance of 2.7 AU. Fig. 23. Same as Fig. 22 but for a heliocentric distance of 2.7 AU.

Fig. 26. Sublimation rate and heat source at a heliocentric distance of 1.29 AU. Fig. 24. Sublimation rate and heat source at a heliocentric distance of 1.29 AU.

the last 10–15 periods. This is true for models with and without a dust cover. Crystallisation rates are varying between 7  10  6 and 10  10  6%/s and heat sources are varying between 0.008 and 0.01 W/m2 in the last third of the whole cycle. Since the crystallisation is a exothermal and irreversible process, one would expect it to be self-sustaining – a so-called runaway crystallisation. As it turns out, the energy amounts produced by the heat source are too small to affect the temperature in the vicinity. This conclusion was confirmed by additional model calculations described in Krebl (2012), which demonstrated explicitly that a switch-off of the boundary heat

source of crystallisation has no visible effect on the calculated results. 5.3.2. Sublimation At a distance of 1.29 AU and with no dust mantle present, sublimation is an ongoing process settling at a rate of approximately 19  10  5 kg/m2 s and a heat sink of 500–550 W/m2. Sublimation related values increase for a distance of 2.7 AU over the whole 40 periods – again presumably to the point of the onset of sublimation. In the presence of a dust mantle, sublimation will be quenched after a certain time, depending on the thermal properties of the dust mantle.

62

B.L. Krebl, N.I. Kömle / Planetary and Space Science 87 (2013) 46–65

The rate of sublimation levels out at about 10–15  10  6 kg/m2 s in the last 10–15 periods and has a declining tendency as the thickness of the dust mantle increases. Naturally the same is true for the heat sink of sublimation, which is decreased by a factor of  6 over the 40 periods ending up at  30 W=m2 . The trend of the heat source is also declining. The results of the Dust Models for a distance of 2.7 AU are similar to those of the No Dust Models, as the temperature at the respective front has not yet risen to a point where sublimation becomes important.

5.4. Lowered thermal properties of the dust mantle For both heliocentric distances, 1.29 and 2.7 AU, the temperature at the sublimation front is  5–10 K lower if the ability

Fig. 27. Same as Fig. 26 but for a heliocentric distance of 2.7 AU.

Fig. 28. Same as Fig. 7 but for a lowered heat capacity and thermal conductivity of the dust mantle.

Fig. 29. Same as Fig. 8 but for a lowered heat capacity and thermal conductivity of the dust mantle.

Table 9 Rates of the phase changes and their heatsources/sinks. Values from 10th to 40th maximum. Model_distance Cry. rate (%/s) Cry. heat source (W/m2) Sub. rate (kg/m2 s) Sub. heat sink (W/m2)

ND_129

ND_27 6

(10–16)  10 0.03–0.012 (17–19)  10  5  500 to  550

D_129 6

(6.5–13)  10 0.024–0.008 (10–80)  10  8  0.3 to  2.2

D_27 6

(8–14)  10 0.025–0.01 (30–12)  10  6  90 to  30

(7–13)  10  6 0.024–0.008 (10–75)  10  8  0.3 to  2.2

B.L. Krebl, N.I. Kömle / Planetary and Space Science 87 (2013) 46–65

63

Fig. 31. Same as Fig. 13 but for a lowered heat capacity and thermal conductivity of the dust mantle.

Fig. 30. Same as Fig. 29 but for a heliocentric distance of 2.7 AU.

Table 10 Same as Table 6 but for lowered thermal properties of the dust mantle. Model_distance

D_129

D_27

T T T T

320–340 185–180 118–116 90–340

165–180 147–162  118 90–180

dust (K) sub. (K) cry. (K) (y) (K)

of the dust mantle to store and conduct heat is lowered. As a result the quantity of ice sublimated and crystallised is reduced by  50% – see Figs. 28, 32 and 33. The difference between temperatures at the crystallisation front is 2–5 times higher when compared with the conditions at the sublimation front in this and the original model. This can be explained by the step wise crystallisation of the amorphous ice.

6. Conclusions and summary Summarising our computational results, the following main conclusions can be drawn:

 Activity due to sublimation: At a distance of 1.29 AU the comet is active, even if the surface is covered by dust.

Fig. 32. Same as Fig. 14 but for a lowered heat capacity and thermal conductivity of the dust mantle.

However, activity will be suppressed after several more periods if the dust mantle is not thinned out by the dragforce of the vapours. At a distance of 2.7 AU, sublimation is negligible but has to be taken into account if the rift is heated further. Counted from period 10 the total increase would be 10 K over 30 periods, which means the temperature at the surface would increase with a rate of 0.3 K/ period. If the sublimation point is reached depends on the

64

B.L. Krebl, N.I. Kömle / Planetary and Space Science 87 (2013) 46–65

amount of matter transformed by phase changes, since the temperatures of the ice domains are lower too. Acknowledgment The authors would like to thank the reviewer Prof. Jacek Leliwa-Kopystynski from Poland for his very careful review work, which contributed significantly to improve contents and presentation of this paper.

Appendix A. Model parameters The following Table 12 provides a list of the most relevant parameters and variables used in our model. Table 12 Parameters used in the model – constants and averaged quantities. Parameter Value

Description

ϱs

Density of ices (kg m  3)a

ϱ Fig. 33. Same as Fig. 15 but for a lowered heat capacity and thermal conductivity of the dust mantle.

Table 11 Same as Table 7 but for lowered thermal properties of the dust mantle. Model_distance

D_129

D_27

Sublimation (m) Crystallisation (m) Dust mantle thickness (m)

 0.018  3.2 þ 0.018

 0.0002  2.4 þ 0.0002



 





orbit of the comet and if it moves towards its perihelion or towards its aphelion, respectively. Activity due to crystallisation: Crystallisation related quantities remain roughly unchanged in all the models once a quasi-equilibrium state is reached and the temperature settles around 140 K. Crystallisation starts at 120 K and needs much less energy than sublimation. At a distance of 2.7 AU crystallisation is still ongoing and the flow pattern of the expelled CO 2 dominates over the flow pattern caused by sublimation, if a dust mantle is present and quenches sublimation. Asynchronous behaviour: Crystallisation is more likely to progress if there is no sublimation taking place on the surface (energy loss due to the heat sink caused by sublimation). Inactivity: If the activity is only based on water ice as a volatile and on freed gases upon crystallisation (as assumed in these model calculations), the comet should be completely inactive at a distance of 5.74 AU. Positive total heat flux at the surface: The overall temperature in the rift is increasing for all three heliocentric distances, which means that the total energy loss of the rift at the surface is smaller than the respective gain. The positive heat flux from the illumination of the surface surpasses the negative ones, which are the latent heat sink caused by sublimation and the radiative heat loss (IR-radiation) from the surface. Influence of thermal properties of the dust mantle: In general a lower thermal conductivity of the dust mantle reduces the

ψ

θs

θf cos ζ em A s d Cp Cd Cices aC bC αs

vsound Λ aαC α M Hcry Hsub χa Xgas kB R Na QHC and QMC Na ηT Pk a

940 (crystalline ice) 910 (amorphous ice) 3250 (dust) 1.2041  10  9 (H2O vapour) 2.814  10  7 (CO2) 0.6 (crystalline ice) 0.3 (amorphous ice) 0.8 (dust) 0.4 (crystalline ice) 0.7 (amorphous ice) 0.2 (dust) – 1

Dust density (kg m  3)b Density of fluids (kg m  3)c

Porosity

Volume fraction of solid material

Same as porosity of a solid Cosinus of local solar zenith angle (vertical illumination) 0.98 (ice), 0.90 (dust) Emissivity of surfaced 0.04 Comets albedob 5.670373  10  8 Stefan Boltzmann constant (W m  2 K  4) 1.29, 2.7 and 5.74 Heliocentric distance (AU) 1850 (H2O vapour), Spec. heat cap., P const. (J kg  1 K  1)d 850 (CO2) 1300 Spec. heat cap. dust (J kg  1 K  1)b See Eq. (25) Spec. heat cap. ices (J kg  1 K  1) 7.49 Const. 1, spec. heat cap. (J kg  1 K  2)b 90 Const. 2, spec. heat cap. (J kg  1 K  1)b 0.1 Thermal cond. of solid (dust) (W m  1 K  1)b Eqs. (23) and (24) (ices) 2500 Velocity of sound (m s  1)e 5 Mean free path (Å)e 567 Const., thermal cond. cry. ice (W m  1 K  2)b 0.025 (H2O vapour) Thermal cond. of fluid (W m  1 K  1)d 0.01465 (CO2) 0.018 (H2O vapour) Molar mass of water vapour (kg mol  1) 0.044 (CO2) 9  104 Latent heat of crystallisation (J kg  1)b 2.8345  106 Latent heat of sublimation (J kg  1)b Amorphous ice fraction 0.05 Trapped gas mass fraction  23 1.3806  10 Boltzmann constant (J/K) 8.314 Universal Gas constant (J K  1 mol  1) 23 6.02214199  10 Avogadro Constant [mol  1] Heat (J m  2) and mass source (kg m  3) of crystallisation 6.02214199  1023 Avogadro Constant (mol  1) Eq. (11) Turbulent viscosity (Pa s  1) Eq. (14) Production term for k (kg m  1 s  2)

http://www.lsbu.ac.uk/water/. Rosenberg and Prialnik (2010). c Calculated from ideal gas law with the models temperature profile and averaged over the periods. d http://www.engineeringtoolbox.com/. e Espinasse et al. (1991). b

B.L. Krebl, N.I. Kömle / Planetary and Space Science 87 (2013) 46–65

References A'Hearn, M.F., Belton, M., Delamere, W., Feaga, L., Hampton, D., Kissel, J., Klaasen, K., McFadden, L., Meech, K., Melosh, H., Schultz, P., Sunshine, J., Thomas, P., Veverka, J., Wellnitz, D., Yeomans, D., Besse, S., Bodewits, D., Bowling, T., Carcich, B., Collins, S., Farnham, T., Groussin, O., Hermalyn, B., Kelley, M., Kelley, M., Li, J.-Y., Lindler, D., Lisse, C., McLaughlin, S., Merlin, F., Protopapa, S., Richardson, J., Williams, J., 2011. EPOXI at comet Hartley-2. Science 332, 1396–1400. A'Hearn, M.F., Belton, M., Delamere, W., Kissel, J., Klaasen, K., McFadden, L., Meech, K., Melosh, H., Schultz, P., Sunshine, J., Thomas, P., Veverka, J., Yeomans, D., Baca, M., Busko, I., Crockett, C., Collins, S., Desnoyer, M., Eberhardy, C., Ernst, C., Farnham, T., Feaga, L., Groussin, O., Richardson, J., Wellnitz, D., White, R., 2005. Deep impact: excavating comet Tempel-1. Science 31, 258–264. Biele, J., Herfort, U., 2012. Rosetta Project Comet Reference Models – comet 67P/ Churyumov-Gerasimenko. Boice, D.C., Soderblom, L., Britt, D., Brown, R., Sandel, B., Yelle, R., Buratti, B., Hicks, M., Nelson, R., Rayman, M., Oberst, J., Thomas, N., 2002. The Deep Space 1 encounter with comet 19P/Borelly. Earth, Moon and Planets 89, 301–324. Davidsson, B.J., Gutierrez, P.J., 2005. Nucleus properties of comet 67P/ChuryumovGerasimenko estimated from non-gravitational force modeling. Icarus 176 (2), 453–477. de Almeida, A., Sanzovo, D.T., Sanzovo, G., Boczko, R., Torres, R.M., 2009. Comparative study of productivity of the Rosetta target comet 67P/Churyumov-Gerasimenko. Advances in Space Research 43, 1993–2000. de Sanctis, M., Capria, M., Coradini, A., 2006. 67P/Churyumov-Gerasimenko nucleus model: portrayal of the Rosetta target. Advances in Space Research 38 (9), 1906–1910. de Sanctis, M., Lasue, J., Capria, M., Magni, G., Turrini, D., Coradini, A., 2010. Shape and obliquity effects on the thermal evolution of the Rosetta target 67P/Churyumov-Gerasimenko cometary nucleus. Icarus 207 (1), 341–358. Espinasse, S., Klinger, J., Ritz, C., Schmitt, B., 1991. Modeling of the thermal behavior and of the chemical differentiation of cometary nuclei. Icarus 92 (2), 350–365. Fanale, F.P., Salvail, J.R., 1986. A model of cometary gas and dust production and nongravitational forces with application to P/Halley. Icarus 66 (1), 154–164. Ghormley, J., 1968. Enthalpy change and heat-capacity changes in the transformation from high-surface-area amorphous ice to stable hexagonal ice. Journal of Chemical Physics 48, 503–508. Gombosi, T.I., Houpis, H., 1986. An icy-glue model of cometary nuclei. Nature 324, 43–44. Houpis, H., 1990. Models of cometary nuclei. Comet Halley Investigations Results Interpretations 2, 173–188. Kelley, M.S., Wooden, D.H., Tubiana, C., Boehnhardt, H., Woodward, C.E., Harker, D.E., 2009. Spitzer observations of the comet 67P/Churyumov-Gerasimenko at 5.5– 4.3 AU from the Sun. The Astronomical Journal 137, 4633–4642.

65

Klinger, J., 1980. Influence of a phase transition of ice on the heat and mass balance of comets. Science 209, 271–272. Klinger, J., 1981. Some consequences of a phase transition of ice on the heat balance of comet nuclei. Icarus 47, 320–324. Koemle, N., Kargl, G., Thiel, K., Seiferlin, K., 1996. Thermal properties of cometary ices and sublimation residua including organics. Planetary and Space Science 44, 675–689. Kossacki, K.J., Szutowicz, S., 2006. 67P/Churyumov-Gerasimenko: modeling of orientation and structure. Planetary and Space Science 54, 15–27. Krebl, B.L., 2012. Numerical Model of the Gas Flow and Temperature Development in an Ice Filled Crevasse on the Surface of a Comet's Nucleus. Diploma Thesis, Located in the library of the Karl-Franzens-Universitaet Graz. Kuznetsov, G., Sheremet, M., 2010. Numerical simulation of natural convection in a rectangular enclosure having finite thickness walls. International Journal of Heat and Mass Transfer 53, 163–177. Lamy, P.L., Toth, I., Davidsson, B.J.R., Groussin, O., Gutierrez, P., Jorda, L., Kaasalainen, M., Lowry, S.C., 2007. A portrait of the nucleus of comet 67P/ChuryumovGerasimenko. Space Science Reviews 128, 23–66. Merk, R., Prialnik, D., 2006. Combined modeling of thermal evolution and accretion of trans-neptunian objects – occurrence of high temperatures and liquid water. Icarus 183 (2), 283–295. Prialnik, D., 1991. A model of gas flow through comet nuclei. In: Koemle, N.I., Bauer, S.J., Spohn, T. (Eds.), Theoretical Modelling of Comet Simulation Experiments. Verlag der Oesterreichischen Akademie der Wissenschaften, Wien, pp. 1–10. Prialnik, D., Bar-Nun, A., 1990. Gas release in comet nuclei. Astrophysical Journal 363, 274–282. Rosenberg, E.D., Prialnik, D., 2009. Fully 3-dimensional calculations of dust mantle formation for a model of comet 67P/Churyumov-Gerasimenko. Icarus 201 (2), 740–749. Rosenberg, E.D., Prialnik, D., 2010. The effect of internal inhomogeneity on the activity of comet nuclei – application to comet 67P/Churyumov-Gerasimenko. Icarus 209, 753–765. Schmitt, B., Espinasse, S., Grim, R., Greenberg, J., Klinger, J., 1989. Laboratory studies of cometary ice analogues. Physics and Mechanics of Cometary Materials – ESA Special Publication 302, 65–69. Seiferlin, K., Koemle, N., Kargl, G., Spohn, T., 1996. Line heat source measurements of the thermal conductivity of porous H2O ice, CO2 ice and mineral powders under space conditions. Planetary and Space Science 44, 691–704. Tsou, P., Brownlee, D., Anderson, J., Bhaskaran, S., Cheuvront, A., Clark, B., Duxbury, T., Economou, T., Green, S., Hanner, M., Horz, F., Kissel, J., McDonnell, J., Newburn, R., Ryan, R., Sandford, S., Sekanina, Z., Tuzzolino, A., Vellinga, J., Zolensky, M., 2004. Stardust encounters comet 81P/Wild-2. Journal of Geophysical Research 109, E12S01. Ulamec, S., Espinasse, S., Feuerbacher, B., Hilchenbach, M., Moura, D., Rosenbauer, H., Scheuerle, H., Willnecker, R., 2006. Rosetta Lander – Philae: implications of an alternative mission. Acta Astronautica 58 (8), 435–441. Weissman, P.R., Asphaug, E., Lowry, S.C., 2004. In: Festan, M.C., Keller, H.U., Weaver, H.A. (Eds.), Structure and Density of Cometary Nuclei. Comets II. Arizona University Press, Tucson, pp. 337–357.