The spatial distribution of water in the inner coma of Comet 9P/Tempel 1: Comparison between models and observations

(1)

The spatial distribution of water in the inner coma of Comet 9P/Tempel

1: Comparison between models and observations

S. Finklenburg

a,⇑

_{, N. Thomas}

a

_{, C.C. Su}

b

_{, J.-S. Wu}

b a

Physikalisches Institut, University of Bern, Sidlerstrasse 5, CH-3012 Bern, Switzerland

b

Department of Mechanical Engineering, National Chiao Tung University, 1001 Ta-Hsueh Road, Hsinchu 30010, Taiwan

a r t i c l e

i n f o

Article history:

Received 14 November 2013 Revised 6 March 2014 Accepted 24 March 2014 Available online 2 April 2014

Keywords: Comet Tempel 1 Comets, coma Comets, nucleus

a b s t r a c t

The near nucleus coma of Comet 9P/Tempel 1 has been simulated with the 3D Direct Simulation Monte Carlo (DSMC) code PDSC++_{(Su, C.-C. [2013]. Parallel Direct Simulation Monte Carlo (DSMC) Methods for}

Modeling Rarefied Gas Dynamics. PhD Thesis, National Chiao Tung University, Taiwan) and the derived column densities have been compared to observations of the water vapour distribution found by using infrared imaging spectrometer on the Deep Impact spacecraft (Feaga, L.M., A’Hearn, M.F., Sunshine, J.M., Groussin, O., Farnham, T.L. [2007]. Icarus 191(2), 134–145.http://dx.doi.org/10.1016/j.icarus.2007.04.038). Modelled total production rates are also compared to various observations made at the time of the Deep Impact encounter. Three different models were tested. For all models, the shape model constructed from the Deep Impact observations by Thomas et al. (Thomas, P.C., Veverka, J., Belton, M.J.S., Hidy, A., A’Hearn, M.F., Farnham, T.L., et al. [2007]. Icarus, 187(1), 4–15.http://dx.doi.org/10.1016/j.icarus.2006.12.013) was used. Outgassing depending only on the cosine of the solar insolation angle on each shape model facet is shown to provide an unsatisfactory model. Models constructed on the basis of active areas suggested by Kossacki and Szutowicz (Kossacki, K., Szutowicz, S. [2008]. Icarus, 195(2), 705–724.http://dx.doi.org/ 10.1016/j.icarus.2007.12.014) are shown to be superior. The Kossacki and Szutowicz model, however, also shows deficits which we have sought to improve upon. For the best model we investigate the properties of the outflow.

1. Introduction 1.1. General introduction

Comet 9P/Tempel 1 was discovered in 1867 and is a Jupiter-family comet with an orbital period of 5.5 years and a perihelion distance of 1.51 AU. It was encountered by the Deep Impact space-craft on 2005 July 4, one day before perihelion. The peak outgas-sing currently occurs before perihelion and the water production rate is then about 5–10 1027_{molecules per second (}_A’Hearn

and Combi, 2007). The thermal inertia was found to be very low (Groussin et al., 2007a, 2007b).Sunshine et al. (2006)discovered water ice patches on the surface but they were too small and too few to explain the total water production rate. These results imply gas emission from ice just below the visible surface or from a surface layer dominated by dark material intimately mixed with the ice.

Tempel 1 is the ﬁrst comet for which there is a detailed shape model (Thomas et al., 2007, 2013) as well as observations of the

gas distribution in the innermost coma (Feaga et al., 2007). This allows one to study the importance of the surface shape for the gas distribution and make detailed comparisons between model results and observations. Crifo et al. (2002)suggested that many of the observed features in the near nucleus coma of comets can be explained by an irregular shape. They suggested that, in the case of Comet 1P/Halley, the main structure of the coma could be sim-ulated by low order spherical harmonics of the shape. This con-tradicted the paradigm originating in the initial interpretation of the Halley Multicolour Camera data (e.g.Keller et al., 1987) which concluded that there were active and less active regions on the nucleus.

The shape inﬂuences the gas coma in two ways. On the one hand, the surface features cause different illumination conditions and therefore different gas production rates if one assumes a homoge-neous composition and outgassing as a simple function of insolation. Secondly, the topography itself causes structures in the gas distribu-tion as has been demonstrated by, for example, (Keller et al., 1994) where crater-like features ‘‘collimated’’ the outﬂow.

In this paper, we compare the measured gas distribution in the near nucleus coma of Comet Tempel 1 (Feaga et al., 2007) with

⇑ Corresponding author. Fax: +41 31 631 4405.

E-mail address:susanne.ﬁnklenburg@space.unibe.ch(S. Finklenburg).

Contents lists available atScienceDirect

Icarus

(2)

simulated gas distributions generated using the Direct Simulation Monte Carlo (DSMC) method. These distributions have been gener-ated in full 3-D using the shape model of Tempel 1 (Thomas et al., 2007). Firstly, we assume homogeneous surface properties, mean-ing that the sublimation rate is proportional to the Hertz–Knudsen flux at each surface element. In this way, we compare the observed gas distribution with a model which follows the Crifo hypothesis (Crifo et al., 2002) of outgassing proportional to local insolation. In a second model, we adopt an activity distribution suggested byKossacki and Szutowicz (2008). We will call this ‘‘the KS model’’. These authors used the activity of Tempel 1 over its orbit derived from ground-based observations and concluded that the data would be best matched if only certain areas were active. They cal-culated the production rate with the so-called ‘‘long tube’’ formula and suggested different thicknesses of the dust crust to explain the different production rates. The long tube formula consists of the Hertz–Knudsen flux scaled with a dust layer dependent constant. It is therefore in principle similar to the proportionality factor, C, used in the homogeneous model and the icy area fraction used byCrifo (1997). The only difference is, that we use here the values suggested by Kossacki and Szutowicz who found them by fitting the observed production rates over whole orbits from Comet 9P/Tempel 1 to their model. Thirdly, we tested an ad hoc model

modiﬁed from the KS model which was an attempt to ﬁt H2O

column densities as measured byFeaga et al. (2007). This will be called the ‘‘Estimated model’’. Finally we show results of a varia-tion of the Estimated model, where the producvaria-tion rate is the same as before, but the surface temperature is everywhere constant at 200 K. This was done to show the impact of the surface tempera-ture on the column densities and the radiation map.

In the next section, we shall discuss some of the inputs needed. Then the DSMC method and the speciﬁc implementation that we use will be described. We will also explain why DSMC is required here (rather than use of a ﬂuid approach) and how the high density regions can be treated. In Sections2.1–2.3we explain the three dif-ferent boundary conditions (models) which were simulated. The results are shown in Section3 and we draw our conclusions in Section4.

1.2. Source values, inputs, and constraints

The total production rate of Comet 9P/Tempel 1 at the time of the Deep Impact ﬂy-by (2005-07-04) was estimated to be 5 1027_H

2O molecule/s (Schleicher, 2007), 9.4 ± 0.69 1027

H2O molecule/s (DiSanti et al., 2007) and 4.6 1027H2O

mole-cule/s (Feaga et al., 2007) using independent methods. The total surface area is approximately 120 km2_{which is an output of the}

shape model produced by Thomas et al. (2007)using the Deep

Impact camera measurements. The smallest facet of the shape model has a size of 100 m2_{; the largest one is 10,000 m}2_{; the}

average facet size is around 4000 m2_{. The surface temperatures}

are between 272 K and 336 K on the sunlit hemisphere (Groussin et al., 2007a, 2007b). Groussin et al.’s temperature map has a resolution of about 120 m.

The effusion velocity of a gas which leaves the surface with a half-Maxwellian is given by

v

eff¼

ffiffiffiffiffiffi

2kT

pm

q

. Note that, as we consider here only a half-Maxwellian and not a full one, there is a mean velocity greater than zero. The half-Maxwellian is directed in a way, that all molecules leaving the surface have a velocity compo-nent pointing away from the surface. We assume this velocity dis-tribution function (VDF) in the three models. With an approximate temperature of 300 K and a molecular mass of 3 1026_{kg for}

H2O,

v

eff 300 m/s. If we assume that only half of the nucleus is

active, the local production rate would be between 8 1019 _and

1.6 1020molecules/m2/s for homogeneous dayside emission depending on the surface temperature. With the estimated

effusion velocity, the number density, n, at the surface is then between 3 1017 _{and 5 10}17_molecules/m3_{. Note that we used}

here the expression for a half-Maxwellian above the surface and neglect all backscattered molecules. The mean free path of a water molecule in the VHS (variable hard sphere;Bird, 1994) model is k¼pffiffiffi2

p

d2_ref1 T

Tref x0:5

1

n. For the VHS parameters of water, we

use dref= 6.15 1010m, Tref= 300 K, and

x

= 1.1 (Crifo, 1989).

x

is the temperature exponent in the viscosity power law. This results in mean free paths above the sunlit surface of 1.2–2.4 m. Note that we have used for this estimation the simpler VHS model instead of the VSS model which we use later. The parameters are chosen in a way that the results of the two models are similar. On the night side, as well as further away from the surface, the mean free paths are greater (up to more than 5 km on the night-side) because of the lower densities. If the gas production is more localised on smaller areas, the local densities can be higher and the mean free path therefore smaller.

The Deep Impact observations of the water distribution in the vicinity of the nucleus provide the principal constraint on our model. It should be noted that we use two aspects of the results inFeaga et al. (2007). Firstly, the image of the water distribution in the imme-diate vicinity of the nucleus is used (Feaga et al., 2007;Fig. 11, second line in this paper; Fig. 5 inFeaga et al., 2007) for qualitative compar-ison. In order to be able to compare Feaga’s observations, which is a radiance map, with our simulation result, which is a density, several computational steps and assumptions have to be made. First we integrate our densities over the line of sight of the observation conditions. This gives the column density maps which will be shown later. Then we assume that the water band at the comet is optically thin. It is not immediately obvious that this is sufﬁcient because estimations byFeaga et al. (2007)show that the strongest water line is most probably optically thick close to the nucleus. However, the observed water band comprises many lines so that the total emission comes from lines which are mostly optically thin. Furthermore, our model predicts higher rotational tempera-tures than assumed by Feaga et al. which should reduce the optical thickness. With this assumption, our radiance maps (Figs. 7, 9 and 11) give upper limits. We will argue below that the assumption of low optical thickness is probably sufﬁcient.

Clearly, Feaga et al.’s radiance map constrains the source distri-bution and the position of putative active regions. However, we also use the determination of water column densities at 7 km resp. 3.8 km distance from the nucleus which are given numerically in the Feaga et al. paper as a quantitative constraint.

We seek to test various hypotheses against these constraints using 3D-DSMC modelling.

1.3. DSMC

The DSMC method was described byBird (1994). The basic idea is that large numbers of real gas molecules are represented by a smaller number of simulation molecules. These simulation molecules move and collide with each other in a statistical way. All molecules within one collision cell at the time of the collision step can collide with each other. The collision pairs are chosen ran-domly. The macroscopic gas properties such as density, velocity and temperature are calculated by averaging the appropriate quantity over all simulation molecules in a sampling cell. The detailed formulae are found for example in Bird (1994). The rotational temperature is calculated by giving each simulation molecule an ‘‘internal temperature’’ and energy is redistributed partially between the internal and external degrees of freedom in each collision (Larsen–Borgnakke model (Borgnakke and Larsen, 1975)). The efﬁciency of this redistribution is parameterized by the rotational relaxation collision number which is given by the

(3)

Parker model (Parker, 1959). The molecular collisions are calcu-lated with the variable soft sphere (VSS) scheme (Koura and Matsumoto, 1991). The model parameters were found by ﬁtting the viscosity and the self-diffusivity description of the VSS model to the results from an extended corresponding states principle analysis (Xiang, 2001) of water. The details are described in

Finklenburg (2014). The temperature exponent is 0.923. The VSS parameter is 0.929 and the reference diameter at a reference tem-perature of 300 K is 5.072 Å.

It is clear that the separation of the microscopic behaviour of a gas in a translational step and a collision step as required by DSMC is only valid if the mean distance between two colliding molecules is smaller than the mean free path. With separation, we mean that the movement and collision of simulated particles in the DSMC simulation are uncoupled such that particles ﬁrst move and then collide with each other based on an instantaneous collisional probability. The ratio of the mean collisional distance (the mean distance between two colliding simulation molecules in the simu-lation) to the mean free path has to be smaller than 1.Bird (1994)

proposed a value of 0.3 as an appropriate limit. This constrains the size of the collision cells (the cells in which the collisions are calculated).

The DSMC program used in this paper is PDSC++ ₍_{Su, 2013}₎

which is based on the PDSC code under development by Wu’s group since 2000 (Wu and Lian, 2003; Wu and Tseng, 2005). The PSDC++code allows a parallel simulation of 2D and 3D ﬂows on unstructured grids. The unstructured grid is deﬁned with the help of the software GRIDGEN™ before the DSMC simulation is started. The parallel performance is optimised by the use of a domain re-decomposition method.

The domain re-decomposition method includes four steps: (1) The initial DSMC simulation to reach a steady state using very few particles (e.g., one particle per cell) based on the domain decomposition with approximately an equal number of cells in each domain; (2) re-decompose the domain based on the particle number distribution obtained in Step (1); (3) decrease the particle weighting to increase the total number of particles for a better statistical accuracy (e.g., lower statistical uncertainties), and (4) perform the ﬁnal DSMC simulation. We have found that for most DSMC simulations which look for steady-state solutions, this approach is very powerful and simple. For the study presented here, the approach is appropriate.

In addition, this version of PDSC++_{can employ hybrid grids in}

two- and three-dimensional cases including either triangles, quadrilaterals or mixed for the former, and either tetrahedrons, pyramids, prisms, hexahedrons or mixed for the latter.

The main difference between DSMC and fluid dynamics meth-ods such as Navier–Stokes or Euler Equations solvers is that DSMC does not have any constraints on the degree to which the gas might be in a non-equilibrium state. Over every outgassing surface there is a non-equilibrium layer called the Knudsen layer. In the 2D case (evaporation from an infinite, flat surface), it takes about 200 mean free paths to reach 99% of the equilibrium value (Cercignani, 2000). Cercignani approximates the velocity distribution function of the gas in the Knudsen layer by a combination of a drifting Maxwellian and two semi-Maxwellians, one pointing towards the surface and one away from it. We speak here from ‘‘reaching 99% of the equi-librium value’’ when the contribution from the drifting Maxwellian to the total velocity distribution function is 99%.

The case of evaporation from a sphere was studied byLukianov and Khanlarov (2000). They characterize the ﬂows by their global Knudsen number which is the ratio of the mean free path to the radius of the source. In our case, the ‘‘global’’ radius of Tempel 1 is about 3 km. With a mean free path of the order of meters, the global Knudsen number is of the order of 103_{. Locally the}

curva-ture can be greater and therefore this estimate is a lower limit.

Lukianov and Khanlarov (2000) found that, for global Knudsen numbers larger then 103_{, the non-equilibrium boundary layer}

(Knudsen layer) goes over into a zone where there is

non-equilibrium because of the expansion. Their criterion to deﬁne local thermal equilibrium is that the difference between the temperatures has to be less than 5%. (They use either the ‘‘temperatures’’ measured in radial and normal direction or the difference between translational and rotational temperature to determine non-equilibrium.)

We conclude that while there may be regions in the flow field where CFD (computational fluid dynamics) methods could be used, these regions are in fact small. The exact size and shape of these regions depends on the CFD method used. The Euler equations

assume local thermal equilibrium everywhere while the

Navier–Stokes equations allow small deviations. Also the degree of accuracy one wants to have influences the regions suitable for CFD methods. It is also difficult to couple the two solution methods and to define suitable boundary conditions for the fluid approach

although this was done by e.g. Burt and Boyd (2009),

Schwartzentruber and Boyd (2006), Wu et al. (2006)and others. There are Knudsen layer correction factors (e.g.Cercignani (1980)

for polyatomic gas) but the thickness is unclear and depends on the local production rate (as explained above, the Knudsen layer thickness is a function of the mean free path of the gas and this depends on the density and therefore on the production rate). Typically one assumes local thermal equilibrium at the sonic surface, the surface at which the Mach number becomes 1. This is the largest thickness for which one can calculate the Knudsen

layer with the approach suggested by Cercignani (1980) and

Ytrehus (1975). Local thermal equilibrium has to be assumed there in order to have well defined input conditions for the CFD part. The shape of this surface, which is typically used as a starting surface in CFD solvers, is only in special cases similar to the shape model (Crifo et al., 2003,Fig. 4). Given these difficulties we have chosen to use an alternative method for optimising the DSMC method for simulation of higher density gas flows.

1.4. Variable time-step (VTS) scheme

This procedure is basically the same as the one proposed by

Kannenberg and Boyd (Kannenberg and Boyd, 2000; Kannenberg,

1998). If the cell size is set to be proportional to the local mean free path, then it is inversely proportional to the number density, n, for the three-dimensional case. The volume of such a cell in three dimensions is therefore proportional to n3_{. As the number of}

par-ticles per cell is proportional to the cell volume times the number density, in this idealised case it is proportional to n2_{. Low density}

regions therefore have more particles than needed for a good DSMC simulation. This will lead to great computational waste, especially in regions far away from the comet and above the night side. To avoid this problem, a variable time-step (VTS) scheme (Wu et al., 2004) has been proposed to obtain a more uniform spatial distribution of simulated particles. The basic idea of the VTS scheme is to enforce the conservation (mass, momentum and energy) of moving simulated particles when crossing the interface between two neighbouring cells. If neither the number of particles nor their velocity should change this implies that the ratio of a weighting factor to the local time step has to be constant over the whole simulation region. As the number density is a priori unknown one assumes again that the cell size is proportional to the mean free path and therefore to n1_{. One assumes therefore}

that the number density is proportional to 1_ffiffiffi V

3

p, where V is the cell

volume. The cell size is thereby assumed to be proportional to the local mean free path, even though it is not strictly true everywhere in the computational domain. One has to remember that the grid is constructed before the simulation run and it is

(4)

therefore up to of the user to build a mesh with cells of a volume that depends on the local mean free path (which is unknown and requires estimating). The local weighting factor, Wi, should be

inversely proportional to the local number density in order to have a uniform distribution of simulation particles. Therefore the local weighting factors are set to be Wi¼ Wmin

ffiffiffiffiffiffiffi Vi Vmin 3 q . The timestep,Dt is then found by the condition that the ratioWi

Dti¼

Wmin

Dtminand therefore

Dti¼Dtmin ffiffiffiffiffiffiffi_V i Vmin 3 q

. One of the advantages in implementing the VTS scheme is to reduce both the number of simulated particles and the number of transient time steps to reach steady state, at which point the simulation particle property sampling normally starts in DSMC. This will result in appreciable time saving for the steady DSMC simulation. The VTS scheme does not ensure the same phys-ical time everywhere in the computational domain. Thus, it can only be used for obtaining the steady-state solution, which most DSMC simulations do. However, it is highly efﬁcient in reducing the total number of particles and number of steps towards the steady state.

In this paper, the VTS scheme is used for all simulations. 1.5. Transient adaptive sub-cell (TAS) method

As stated above, to have reliable results with the DSMC method, the distance between two colliding molecules (the mean collisional spacing, mcs) should be smaller than the mean free path (mfp) of the gas. In standard DSMC procedures, the two colliding molecules are chosen randomly from a cell. Their position inside the cell is ignored. This means that the cell size has to be small compared to the mean free path in order to fulfil the requirement. This leads, especially in a 3-D flow field, to an enormous number of cells. In our case, this would be more than needed for the required resolu-tion. A further minor issue is that the result files become large which makes them hard to handle. The way around these issues is to use ‘‘transient adaptive subcells’’ or TAS (Su et al., 2010).

The TAS method is similar to Bird’s subcell method employed on a structured grid, in which the number of subcells depends on the ratio of local mean free path to the local cell size. This has been further extended to the case of unstructured grids (Su et al., 2010). To the best knowledge of the authors, Ivanov’s group has not used this idea to improve the collision quality; instead, this group used their own so-called Majorant Collision Frequency Scheme (Venkattraman et al., 2012) to improve the collision quality. TAS divides every collision cell dynamically into many subcells based on the number of simulation particles and the collision partners are chosen to be in the same or in adjacent subcells. There are on average 2 simulation particles per subcell. This allows one to find nearest neighbours for the collisional step and reduces thereby the mcs. Note that the subcells are only used to find the nearest neighbours. The result files still have the coarser resolution of the sampling cells.

There might be some concern because normally many more than two simulation particles per cell are required. This require-ment is correct for the sampling cells, but we are talking here about the sub-cells. In PDSC++_{, we dynamically divide the sampling}

cell (or background cell) into different numbers of subcells based on local conditions in the sampling cell (number of particles). We divide the sampling cell into many sub-cells. We select randomly the first particle from the particles in the cell and then select the second by searching through the particles in the same subcell. If not, then neighbouring subcells are searched for the second parti-cle until it is found. In this way, we can guarantee the almost near-distance 2nd particle can be selected. Fig. 1shows a typical 2D example. The number of sub-cells in each coordinate direction is estimated as ffiffiffiffiffiffiffiffiffiffiffiffiffiffiNcell=P

p

. Here, Ncelland P is the number of particles

per background cell and the number of particles per sub-cell. Based on our experience, 2 particles per subcell are enough to guarantee a good quality of collision, even though the cell size is much larger than the local mean free path. Indeed, there are many more than 2 particles in the background cell in general. For more details one can refer toSu et al. (2010).

In these simulations, it turns out that the mcs/mfp ratio is, even with the transient adaptive subcells, larger than 1 in some regions. We tested the influence of this problem for our case by running one simulation with a coarser grid and without TAS (1,119,719 cells and max. 7,000,000 particles) and then the same case with a finer grid and TAS (6,191,856 cells and max. 10,000,000 particles). Note that ‘‘max.’’ gives the maximum number of particles which can be calculated. This number is not equal to the number of particles used in the simulation. The only information given about the real number of simulation particles is that this is smaller than this max-imum number. By doing this we wanted to test the two extremes. Above the region with the highest production rate on the dayside, the mcs/mfp is about 40 in the TAS case and more than 100 in the coarser grid case. The largest difference in the column densities in the orientation we use for comparison with Deep Impact observa-tions is about 8%. This is in a small region on the night side where both simulations have mcs/mfp values of less than 0.1. Therefore, the dominant differences occur in the rarefied regions which should not be caused by the grid. We suspect that the large differ-ences in the rarefied regions are caused by numerical scattering together with the fact, that the density in this region is very low so that small absolute differences produce high relative errors. The mean difference over the whole flow field is about 2% which therefore provides an estimate of the numerical error caused by the high mcs/mfp ratio. Therefore we conclude that the presented results allow us to compare with the measured data to an accuracy which is probably better than 10%. For the presented models, the fine grid with TAS was used throughout.

It should be noted that use of TAS does not have a major impact on computational speed (computation times are around 5% longer using TAS) which is also a signiﬁcant advantage given the long com-putation times necessary for 3D rareﬁed gas dynamics calculations.

2. The simulated models

For all models, we used the 3D shape model fromThomas et al. (2007). The coordinate system is deﬁned inThomas et al. (2007). The coordinate origin is at the centre of mass of the model assum-ing a homogeneous internal mass distribution. The spin vector, which deﬁnes the north–south axis, is along the calculated maxi-mum moment of inertia. Finally, the prime meridian is set at the centre of a 350 m crater. SPICE (Acton, 1997) gives the sub-solar point as the interception point between the comet surface and the comet–Sun line at the time 2005-07-04T05:39:18 to be 12.6°W (=347.4°E) and 13.7° latitude. This is the time of observa-tion 9000039 which was used for the high resoluobserva-tion temperature map (Groussin et al., 2007a, 2007b). This is used as the position of the sub-solar point in our modelling. The vector describing the viewing direction for the Feaga et al. observations is [0.358, 0.830, 0.427] in our cometocentric coordinate system.

The simulations were typically run with a maximum of 107

particles in around 6 106_{cells. The smallest timestep was}

1 103_{s. 3 10}4_{timesteps were simulated. Unfortunately, there}

is no strict criterion by which we could decide if steady state was reached. Beside the fact that 30 s is more time than needed by the average molecule to reach the end of the simulation region, we checked if the total number of simulation molecules was still changing with time. As there no systematic increase of this number in the last seconds was visible, we concluded that a steady state

(5)

was reached. We used TAS as well as VTS to optimise the compu-tation. Given that each simulation run requires 12 h of computa-tion on a 32 CPU parallel system, it should be clear that a comprehensive investigation of the possible parameter space is not computationally feasible and hence some ‘‘guesswork’’ has been used with respect to the source distributions to limit the number of runs while still retaining consistency with the con-straints. Besides the computation time, it should be noted that a surface distribution of sources of unknown extent cannot be described by a single parameter which can be varied to cover all possibilities. Therefore each case has to be developed by hand. This excludes an automatic parameter study.

2.1. Homogeneous surface model

In the homogeneous model, we assume that the gas ﬂux is everywhere proportional to the Hertz–Knudsen ﬂux, i.e.

Q ðTÞ ¼ C ffiffiffiffiffiffiffiffiffiffiffiffi m 2

p

kT r psatðTÞ ðkg=m2sÞ ð1Þ psatðTÞ ¼ i1e i2 T _ðPaÞ _ð2Þ With i1= (3.23 + 0.73/0.38) 1012Pa and i2= 6134.6 ± 17.0 K

(Gundlach et al., 2011). Q is the local water vapour ﬂux, m is the mass of a water molecule. k is the Boltzmann constant, T is the local temperature of the sublimating surface, psatis the local saturation

vapour pressure at the given temperature and C is a proportionality constant smaller than 1. C is determined by the total production rate which is set to 9.4 1027_{molecules/s (}_{DiSanti et al., 2007}_{) in}

this model.

The temperature is the equilibrium temperature under the assumption that the thermal conductivity into the interior is neg-ligible and only IR emission and the net sublimation play a role, i.e.

Sc R2h

ð1 AbÞ maxðcosð

a

Þ; 0Þ ¼

r

2T4þ QðTÞH ð3Þ

H ¼ k

mi2ðJ=kgÞ ð4Þ

Here, Sc= 1360 W/m2is the solar constant, Rh= 1.51 AU is the

helio-centric distance of the comet at the time of the encounter, Ab= 0.014 is the visible Bond albedo (Groussin et al., 2007a,

2007b),

a

is the solar zenith angle,

r

¼ 5:6704 108W m2_K4

Ste-fan–Boltzmann constant,

¼ 0:9 for consistency with Groussin et al. (2007a, 2007b) and H(T) is the enthalpy of sublimation (Gundlach et al., 2011).

The shape model is quite smooth so that there are no obvious regions where shadow could play a role other than on the regions

which face away from the Sun and they are included in Eq. (3). Also, there is no diffusion from the interior considered.

To calculate the model we used the solar zenith angles given by the SPICE system (Acton, 1997). For each surface element, its cen-tre was used as a reference point for which the solar zenith angle and the temperature was calculated. This value was then assumed to be the same for the whole element. The resulting boundary con-ditions are shown inFig. 2a (Temperature) andFig. 2b (Production rate). The temperatures are between 6 K and 316 K (measure-ments: ‘‘272–336 K on the sunlit hemisphere’’; Groussin et al., 2007a, 2007b). The minimum temperature has been set artiﬁcially to be consistent with no gas production. There is simply no gas production on the night side. The maximum temperature (316 K) arises because in matching the desired production rate, the propor-tionality constant, C, becomes less than 1, and the requirement for energy balance, implicit in Eq.(3), forces the temperature to this value.

Simulation molecules returning to the surface are absorbed.

2.2. The KS model

Here we simulate a nucleus with a dust layer height distribu-tion as suggested byKossacki and Szutowicz (2008). We refer to this as the ‘‘KS model’’ and the paper from which it is derived as ‘‘KS’’.

The temperature map was calculated by neglecting all heat transfer into the interior and assuming that all incoming energy is emitted. The emissivity,

, was calculated such that the equilib-rium temperature at the sub-solar point, which was the hottest region on the nucleus, is 336 K – the maximum measured (Groussin et al., 2007a, 2007b). This gives an emissivity of 0.82. On the night side the minimum temperature was arbitrarily set to 80 K. KS do not discuss the nightside temperature, but as most of it is covered by a thick dust layer in our model, the exact value of the nightside temperature does not play a big role. Even if it would be much higher, the 2 m thick dust layer (see below) would, on the one hand, not conduct enough heat to the subsurface and, on the other hand, the gas could hardly diffuse through it in the relevant time scales.

The relevant equation is

Sc R2h

ð1 AbÞ maxðcosð

a

Þ; 0Þ ¼

r

2T4 ð5Þ

Fig. 3compares our temperature model with the one published by Groussin et al. (2007a, 2007b). In order to facilitate the comparison, the temperature scales are the same. Note again that

(6)

the lowest temperature in the model has been set to 80 K (on the night side where there are no observations) and not 272 K which was the lowest measured temperature.

The temperature of the subliming layer just below the surface is calculated by the condition that all energy conducted to the inte-rior goes into sublimation. The thermal conductivity of the surface is constant and assumed to be

j

= 2 102_{W/m/K, as given in KS}

for a non-uniform dust mantle.

KS describe the ﬂow coming through a dust layer thicker than 0.5 cm with the long tube formula:

Q ðTÞ ¼8 3 rd

W

d

s

D

d ffiffiffiffiffiffiffiffiffiffiffiffi m 2

p

kT r psatðTÞ ðkg=m2sÞ ð6Þ

with rdis the effective pore radius (0.10 mm),Wdis the porosity of

the dust layer (0.73),

s

is the tortuosity (assigned a value ofpffiffiffi2) and

Ddis the thickness of the dust layer (KS). With, rd,Wd, and

s

held

constant, the production rate is a function of the assumed temper-ature and the dust layer thickness. As input for our simulation, only the gas density inside an imaginary subsurface reservoir and its temperature are needed as input conditions. To speculate about

the physical processes which could produce these conditions is out of the scope of this paper.

In the KS paper, the nucleus is represented by an ellipsoid and is divided into 32 sectors. According to their analysis, only 3 of these 32 sectors should be active. KS use 4 longitude sectors, 90° apart with the ﬁrst sector starting at 0°. There are 8 latitude sectors with centres at ±61°, ±30°, ±21° and ±12°. The 3 active sectors ‘‘are cen-tred at the latitudes 12°, 21° and 30° at any longitudes among four that are analysed in the model with the 32 sectors. In the selected 3 sectors, the dust layer is only 1 cm thick, while the remaining surface is covered by a thick isolating mantle.’’ As KS use the equatorial bulge as zero longitude the longitude values are shifted by about 40° with respect to the grid deﬁned by

Thomas et al. (2007).

We have chosen the longitudes from those proposed by KS in a way to ﬁt the observed gas distribution best because we have no rotational phase constraints. Also it is clear in KS’s paper that low activity at high latitudes does not change the total production rate signiﬁcantly at any position along the orbit. Such an activity is observed byFeaga et al. (2007)but is not present in the KS analysis. Therefore in the KS model presented here we have added a small

Fig. 2. Top row: Surface temperature in the ‘‘homogeneous surface’’ case. In the ﬁrst picture, the Sun is to the right and the shown axes give the orientation of the cometocentric coordinate system as deﬁned in ‘‘tempel1_coordinates.lbl’’. The z-axis is the rotation axis and the x-axis points toward the prime median. This orientation is the one in which the observations at the time 2005-07-04T05:39:18 UTC were made. The second and the third picture show the view on the rotation poles of the comet. Second row: Production rate in molecules/m2

/s in the ‘‘homogeneous surface’’ model. The orientation of the images is the same as in the ﬁrst row.

Fig. 3. (a) Temperature distribution in the KS model. The minimum temperature is 80 K. The emissivity is 0.82 – a value set by the condition that the maximum temperature should be 336 K at the sub-solar point. The albedo is 0.014 (Groussin et al. (2007a, 2007b)). Tempel 1 was at a heliocentric distance of 1.51 AU. (b) Temperature map of Groussin et al. (2007a, 2007b). The dark blue regions are the ones without measurements or where no temperature data were given. (For interpretation of the references to colour in this ﬁgure legend, the reader is referred to the web version of this article.)

(7)

source at high latitude, called the north polar source in the follow-ing, to attempt agreement without breaching the observational constraints.

The active sectors and their location as well as their dust layer thicknesses are given in Table 1. Sectors 2, 3, 4 are the ones described by KS. All other sectors (numbers 1, 5 and 6) were added to ﬁt the observed gas distribution. Note that sector 6, the north polar source, is smaller than a single cell and a smoothing function has been applied. The parameters inside each surface element are constant. The smoothing makes the boundary between the differ-ent sectors less sharp by changing the dust layer thickness over about 5 surface elements. The resulting dust thickness distribution is shown inFig. 4.

Outside the sectors listed inTable 1, the dust layer is 2 m thick and the surface therefore practically inactive. The heights of the dust layer at the sides of the sectors were smoothed to minimise structures appearing because of jumps in the boundary conditions at sector boundaries. The gas density directly above the surface also depends on the local temperature. Therefore the density maps are not simply proportional to the dust layer thickness. The density map can be seen inFig. 5.

In total 32.5% of the total surface is ‘‘active’’ and produces 96% of the total gas ﬂux. The total amount of gas produced is 7.65 1027

H2O molecule/s. This is in the middle of the range of gas production

rates found by various authors (see above andTable 4). Note that we deliberately used here the long-tube formula and the best ﬁt values from the KS paper to compare the results from this paper as directly as possible with Feaga et al.’s results.

Simulation molecules hitting the surface are absorbed. 2.3. Estimated model and its constant temperature variant

A case was calculated in which the dust layer thickness was chosen ‘‘by hand’’ in a way to reproduce the observed gas distribu-tion to a higher degree of ﬁdelity. Only water vapour was simu-lated. In this case, we focused on reproducing the lower limits of the column densities at 3.8 km and 7 km distance given inFeaga et al. (2007).

The distribution of the sectors on the nucleus in this ‘‘Estimated’’ model is similar to the KS model. But sector 6 was changed and some sectors are essentially ‘‘inactive’’ because their dust layer thickness is the same (0.4 m) as the one of the ‘‘inactive’’ regions. The returning molecules are treated in the same way as before. The difference is that the dust layer thicknesses have been changed and the temperature model is different. The position of the sectors and the dust layer thicknesses are given inTable 2.

The temperature was taken fromGroussin et al. (2007a, 2007b)

(Fig. 3b). Everywhere where there was no temperature

measurement, the minimum measured (272 K) was assumed. This gives much higher night side temperatures than the modelled

temperature maps used in the KS model. However, the long-tube formula we are using here has two parameters (the surface temperature and dust layer thickness) to control the production rate. Hence, the outgassing rate is controlled by a dust layer of varying thickness. In the Estimated model, the maximum thickness assumed is L = Lmax= 0.4 m so that, together with the higher

night-side temperatures, emission from the nightnight-side is non-negligible.

Fig. 6shows the resulting initial number densities and the local production rates.

The total production rate is here 4 1027_{molecules/s and}

therefore slightly below the observed production rates.

Finally, we simulated a situation in which the production rate is everywhere the same as in the Estimated model. Therefore also the total production rate is the same. The difference is, that in this variation, the initial temperature of the gas is everywhere constant at a value of 200 K. By comparing the Estimated model with its variant one can therefore estimate the importance of the initial temperature to the column densities and the radiation maps. 3. Results

Feaga et al. (2007)measured the distribution of water vapour, CO2and dust in the vicinity of the nucleus using the High

Resolu-tion Instrument Infrared Spectrometer (HRI-IR). In their paper, two observations are of immediate interest to us.

Firstly, radiance maps of the spectral lines of the different gases are shown. We concentrate here on the water vapour which is assumed to be the dominant gas component and not inﬂuenced by ﬂows of other species. The water vapour radiance map from

Feaga et al. (2007)is shown inFig. 11(second line) to provide com-parison with our best model. Assuming that the vapour is not opti-cally thick, the radiance is proportional to the column density along the line of sight. We assume an optically thin atmosphere and calculate the radiance resulting from a given number density with the formula (Feaga et al., 2007):

Radiance ¼Ng r2 hc k 1 4

p

with N being the number density in particles per square metre and g = 2.85 104_s1_{being the g-factor of water at 1 AU. r = 1.51 AU is}

the heliocentric distance at the time of the observation, h is the Planck constant, c is the speed of light, k ¼ 2:66 106_{m is the}

cen-tral wavelength of the emission band. The radiance is then given in (W/m2_/sr).

Secondly, Feaga et al. gave lower limits on the column density at 3.8 km resp. 7 km distance from the centre of the nucleus. These are listed inTable 3.

In order to compare the simulation results with Feaga’s et al. results we need to compute column densities. When computing the column density, we must account for the ﬁnite size of the sim-ulation region. Simply adding up the densities leads to artefacts at the edges. The problem is, that when the simulation region is a sphere and you add up all cells which lie behind each other in the line of sight, then there are more cells in the lines going through the centre of the sphere and only very few cells on the sides. In the extreme case, on the edges, there is only a single cell. In order to avoid these artefacts, a radial r2_{dependence of the gas}

density is assumed outside of the simulated 10 km sphere so that a simple integration of the column can be made. The densities along the line of sight are summed up and then displayed in a colour plot. The radiance is calculated from these column densities to be compared directly withFig. 11.

In addition, the column densities at 3.8 km resp. 7 km distance from the centre of the nucleus have been calculated to be compared with the lower limits measured byFeaga et al. (2007).

Table 1

Boundary conditions for the KS model. L is the thickness of the dust layer in the active regions and it is L = 102

m thick. The dust layer of the inactive regions is 2 m thick. Sector nr. Mid latitudea (min; max) Mid longitude (min; max) Thickness of dust layer 1 61° (90°; 45.5°) 0° (0°; 360°) 2L 2 30° (45.5°; 25.5°) 85° (40°; 130°) L 3 21° (25.5°; 16.2°) 0° (315°; 45°) L 4 12° (16.2°; 0°) 0° (315°; 45°) L 5 (45.5°; 70°) 270° (225°; 315°) L 6 (89°; 90°) 0° (0°; 360°) 10L

We assumed that the boundary between two adjacent sectors was in the middle: e.g. for 30° and 21° the boundary should be at (30° to 21°)/2 = 25.5°.

a

The given value for the mid-latitudes is the one given inKossacki and Szutowicz (2008).

(8)

Finally, the simulated production rates are compared to the observations (Table 4).

Fig. 7shows the column densities in the case of the homoge-neous model. There is only one region with high density and this is close to the subsolar point. There is little emission from the poles

which is contrary to the measurements. Fig. 8 compares the

measured column densities at 3.8 km and 7 km distance to the nucleus with the ones given as lower limits byFeaga et al. (2007). It is obvious that this model fails to match either the near-nucleus radiance map or the column density limits at 3.8 km and 7 km. It can be immediately concluded that an outgassing from a homoge-neous surface with a half-Maxwellian emission from each facet of the given shape model is not compatible with observations.

Figs. 9 and 10show the column density map, the radiance map and the column densities at 3.8 km and 7 km distance for the KS model. In this model, there is more emission in the direction of the poles and the column density map shows a number of similar-ities with the observations. Specifically, the highest column densi-ties are well represented in the model. The peak of the water column density in the model is no longer above the sub-solar point but accurately representing the more southerly activity. The added north polar source fits the observation reasonably well although a south polar source might be additionally needed. On the basis of the column density map alone one might conclude that the KS model (with slight modification) is actually a good fit to the data. However, there are significant differences in the column densities

Fig. 4. The dust layer thickness of the KS model seen from different perspectives. (a) Observation perspective, (b) North pole, and (c) South pole.

Fig. 5. First row: Initial number density of the KS model shown from different perspectives. Left: Observation perspective, middle: North pole, right: South pole. The scale is logarithmic and goes from 1 1017

to 1 1019

particles per m3

. Second row: Production rate of the KS model. The orientation is the same as in the ﬁrst row. The scale is logarithmic and goes from 5 1019

to 2.5 1021

particles per m2

per second.

Table 2

Boundary conditions for the Estimated model. L = 0.4 m dust layer thickness. The dust layer thickness of the inactive areas is also 0.4 m and therefore much smaller than in the KS model. Note that the sectors 2, 3, 5 and 6 are not different from the inactive area any more.

Sector nr. Mid latitude (min; max) Mid longitude (min; max) Thickness of dust layer 1 61° (90°; 45.5°) 0° (0°; 360°) 1.5L 2 30° (45.5°; 25.5°) 85° (40°; 130°) L 3 21° (25.5°; 16.2°) 0° (315°; 45°) L 4 12° (16.2°; 0°) 0° (315°; 45°) 0.1L 5 (45.5°; 70°) 270° (225°;3 15°) L 6 (70°; 90°) 0° (0°; 360°) L

(9)

at 3.8 km and 7 km. The North pole and night side emission are below the lower limits set by Feaga et al. and the angular position of emission peaks do not ﬁt.

Figs. 11, 12 and 14show the results of the Estimated model and

Fig. 11(second line) shows Feaga et al.’s radiance map for compari-son. Here, the column density map does show the high density region south of the subsolar point and the ﬁt is generally rather good. There are two points where the simulated radiation map differs sig-niﬁcantly from the observed one. One the one hand, there is a region close to the surface on the nightside where limited emission was observed. Further out from this point, the observed radiation is

higher again. In our simulation, we ﬁnd a higher radiation close by the surface and a steady decrease going outwards. The lack of radi-ation close to the nucleus may be caused by shadowing effects or inhomogeneous nightside emission not accounted for in the model. On the other hand, there is the south polar region (below the nucleus in Fig. 11). Here our simulation shows less radiation than Feaga et al.’s observations. In the KS model simulations (Fig. 9) we ﬁnd more radiation in this region, but even more than observed. We con-clude therefore, that both proposed models are not ideal in this region. A better model will have an intermediate production rate in this region. The column density at 3.8 km and 7 km distance from the centre of the nucleus (Fig. 14) is now much more similar to the measured values and everywhere close to or above the lower limit. Also the total production rate of 4 1027_{molecules per second is}

within 15% of the measured production rate.

We conclude that this model is superior to the other two mod-els of outgassing. It should be clear however that we have selected sources and source strengths in a somewhat ad hoc manner and that there is improvement possible. The problem is almost cer-tainly severely under-constrained and the solution we present is undoubtedly not unique. On the other hand, we can, at this point, show the properties of the ﬂow should this particular solution be valid.

One alternative solution is the variation with the constant ini-tial temperature shown in Fig. 13. In the column density map, the difference between the Estimated model and its variant is barely visible. In the radiance map, which is simply proportional to the column density, the differences are clearer but still small. Due to the lower initial temperature, the gas expands slower and the density is slightly higher. The gas velocity 10 km from the cen-tre of the source is in this case, about 700 m/s compared to the 900 m/s in the Estimated model. But the main features of the col-umn density and radiation maps are very similar in both Estimated models although the initial temperature distribution varies signiﬁcantly.

Fig. 15shows the speed of the gas along a line along the vector [1, 0, 0] (in cometocentric coordinates). Because we used an unstructured grid, it is highly unlikely to have points exactly on this vector. We show therefore the points in the vicinity of this

Fig. 6. First row: Initial number density in molecule/m3

for the Estimated model. The colour scale is logarithmic and the same as in the ﬁgure for the KS model (1017

molecules/m3

in dark blue to 1019

molecules/m3

in red). Second line: Production rates in molecules per square meter per second for the Estimated model. The colour scale is logarithmic and goes from 5 1019

molecules/m2

/s (dark blue) to 2.5 1021

molecules/m2

/s (red). (a) Observation perspective, (b) North pole, and (c) South pole. (For interpretation of the references to colour in this ﬁgure legend, the reader is referred to the web version of this article.)

Table 3

Table 3 fromFeaga et al. (2007)giving the lower limits for the water column densities in 7 km distance from the centre of the nucleus together with the estimated position of the region. The angle to the x-axis is in the projected view asFig. 11, second line. The angles are not given by Feaga et al. but simply estimated from the ﬁgures.

Region H2O (molecules m2) Distance from the

centre of the nucleus (km) Angle to x-axis (°) Ecliptic north >6.9 1020 3.8 60–135

Positive rot. pole >2.8 1020

7 15–60

Sunward >4.0 1020

7 0–15; 315–360

Ecliptic south >7.1 1020

3.8 240–315

Negative rot. pole >2.9 1020

7 195–240

Anti-sunward >2.0 1020

7 135–195

Table 4

Overview of the different production rate measurements and calculations at the time of the Deep Impact encounter.

Total production rate estimate (molecules/s) Feaga et al. 4.6 1027 Schleicher 5 1027 DiSanti et al. 9.40 ± 0.69 1027 Homogeneous model 9.4 1027 KS model 7.65 1027 Estimated model 4 1027

(10)

vector. (We calculate the scalar product of the position vector of the centre point of each cell scaled to have norm 1 and the direction vector which is also scaled to have norm 1. If this product is larger than 0.999, the point is considered to be along the given vector.)

It is obvious that the gas does not reach its ﬁnal velocity in the simulated 10 km domain. Also this ﬁnal velocity will be clearly more than the 500 m/s assumed in theFeaga et al. (2007)paper to calculate production rates. We note thatGombosi et al. (1986)

using a ﬂuid approach also found higher gas velocities than the 500 m/s assumed. We note that the column densities used by Feaga et al. to calculate the total production rate vary by a factor of two within their two ‘‘bins’’. Therefore there is some imprecision

in the given value. Comparing our velocities to values found by other groups and for other cases does not seem to be useful as the velocity depends on the initial gas temperature and this varies for the different thermophysical nucleus models and the different heliocentric distances for which the various models were calculated.

One can try to estimate the final gas velocity by approximating the comet by a sphere. Then the gas flow field is given by the formulas in Sone and Sugimoto (1993). For the free molecular (collisionless) flow approximation, the final velocity is simply u1;c¼ 2 ffiffiffiffiffiffiffi 2kT0 pm q which is about u1;c¼ 34:3 ffiffiffiffiffiT0 p in (m/s) for water

Fig. 7. Column densities in the case of the homogeneous model. The numbers on the axes are the distances from the centre of the nucleus in meters.

Fig. 8. Column densities at 3.8 km and 7 km distance from the centre of the nucleus in the case of the homogeneous model (given in molecules per m2

). The simulation results are shown as crosses (7 km) or diamonds (3.8 km); the solid line represents the column densities estimated byFeaga et al. (2007)for 7 km distance, the dashed lines are Feaga et al.’s data for 3.8 km distance.

Fig. 9. Column densities in the case of the KS model. The orientation is again with the Sun to the right and the ecliptic North pole to the top. The numbers on the axes are the distances from the centre of the nucleus in meters.

Fig. 10. Column densities at 3.8 km and 7 km distance to the nucleus for the KS model (given in molecules per m2

). The crosses represent the simulation results, the solid lines areFeaga et al.’s (2007)estimated column densities for the 7 km distance from the source. The diamonds and the dashed lines show the simulated values resp. theFeaga et al. (2007)estimates for the 3.8 km distance case.

(11)

vapour. T0is the initial gas temperature (K), k the Boltzmann

con-stant and m the mass of a single molecule in (kg).

In the ﬂuid limit and for three internal degrees of freedom, the formulae given inSone and Sugimoto (1993)have to be slightly

changed (see for example, Finklenburg, 2014). They give the

formulae for gases with no internal degrees of freedom. In the fluid limit, the Mach number of the flow goes to infinity at an infinite distance from the source. For this case, one finds that u1;f¼

ffiffiffiffiffiffiffiffiffi

28kT

3m

q

. Tis here the temperature at the sonic surface (K).

As we look at the ﬂuid limit, we neglect the Knudsen layer thick-ness and assume that the Knudsen layer can be described by its 1D model described inCercignani (2000). For a gas with 3 internal degrees of freedom, T

T0¼ 0:814 and one ﬁnds for water that u1;f¼ 59:2

ffiffiffiffiffi T0

p m=s.

In the Estimated model, the initial gas temperature in the region close to the [1, 0, 0] vector is about 320 K. Therefore the ﬁnal

velocities are u1,f= 614 m/s and u1,f= 1059 m/s. The velocity at

the end of the simulation region at this position is about 900 m/s but there is clearly still an acceleration going on.

In the 200 K surface temperature variant of the Estimated model, the initial gas temperature is 200 K and one ﬁnds the values u1,f= 485 m/s and u1,f= 837 m/s. At the end of the simulation

region, the velocity is at about 700 m/s and again there is still an acceleration noticeable.

We note also that the pure fluid description is probably not cor-rect as the gas rarefies with expansion. Therefore values below the fluid estimation are expected. The simulation results are therefore consistent with expectation.

Fig. 11. Top: Estimated model, radiance map (W m2

sr1

). The Sun is to the right and the ecliptic north to the top. The comet nucleus is marked with white dots, the black squared underneath it, are the regions without data. The distances are given in meters. Bottom: Radiance map from Fig. 5 inFeaga et al. (2007). The colour bar indicates the radiance values (W m2_sr1_{) along the line of sight. The ﬁeld of view}

is 43 km 10 km. The Sun is to the right. (For interpretation of the references to colour in this ﬁgure legend, the reader is referred to the web version of this article.)

Fig. 12. Estimated model, column densities (m2

). The Sun is to the right and the ecliptic north to the top. The comet nucleus is marked with white dots, the black squared underneath it, are the regions without data. The distances are given in meters.

Fig. 13. Estimated model with constant surface temperature of 200 K. The radiance map and the column densities are given with the same scales as used for the Estimated model inFigs. 11 and 12. Due to the lower initial temperature, the gas densities are slightly higher as can be best observed in the radiance map.

Fig. 14. The column densities at 3.8 km and 7 km distance from the centre of the nucleus for the Estimated model (given in molecules per m2

). The simulation results are marked with crosses (7 km) or diamonds (3.8 km), the lower limits fromFeaga et al. (2007)are shown as solid line (7 km) or dashed line (3.8 km).

(12)

Fig. 16 shows different gas parameters in cuts through the origin of the cometocentric coordinate system. The local number density is the ﬁrst parameter shown inFig. 16. The second row shows the mean free path of the gas. It is nearly everywhere larger than 1 m but, inside the simulation region, smaller than 100 m. That means that it is everywhere large compared to dust grains but small compared to the cometary size. This, in turn, implies that a free molecular aerodynamic approach in modelling the dust ﬂow would be acceptable (e.g.Bird, 1994).

The third parameter shown inFig. 16 supports this result. It shows the degree of non-equilibrium by presenting the ratio of the rotational temperature to the translational one. This value is everywhere close to, but larger than 1. That it is larger than 1 is expected because there are no mechanisms to cool the rotational degrees of freedom other than collisions. The translational temper-ature drops because of the expansion. Therefore the rotational degrees of freedom are always warmer or equal in energy with respect to the translational ones. That the simulation result is everywhere close to one, indicates that the gas is close to a local thermal equilibrium.

The total temperature is shown in the last row ofFig. 16.

Fig. 17shows the speed and the gas velocities in the radial direction and perpendicular to it as well as the ratio of these two velocities. By comparing the plots, it is immediately clear that the radial flow is dominant outwards from about 1 nucleus radius above the surface. But there are perpendicular velocities which are not negligible. The flow field <1 nucleus radius from the surface is extremely complex. The perpendicular velocities are especially high close to the surface. This is the region where dust is mostly accelerated as the gas is densest there. Therefore a complex spatial distribution of dust might be expected.

Going back toFig. 6, second row, one can compare the typical nightside to a typical dayside production rate of the Estimated model. Note that the colour bar is logarithmic. The highest dayside production rates are about 2.5 1021_molecules/m2_{/s. The}

night-side production rate is only around 1 1020_molecules/m2_{/s and}

hence considerably lower than the local dayside production but clearly non-negligible.

4. Conclusions and discussion

Three different models (and one variant) of the outgassing of Comet 9P/Tempel 1 were simulated and compared with Deep Impact measurements. The measurements used were a radiance map in an H2O emission line, an estimation of a lower limit of

the column densities at different points around the nucleus in the inner coma, and the total production rate.

It is clear that a nucleus, uniformly active in response to solar insolation, cannot explain the observed gas distribution around Tempel 1. Active regions are required which is in broad agreement with the conclusions ofFarnham et al. (2012)and further supports the concept of inhomogeneous surface properties of cometary nuclei.

A model suggested byKossacki and Szutowicz (2008)has also been studied. Their result appears to reproduce roughly the posi-tions of active regions on the surface but with a number of limita-tions. The simulated near-nucleus radiance map is similar to the Deep Impact radiance map ofFeaga et al. (2007)but the KS model requires additional polar sources to provide a better ﬁt. However, the additional constraint provided by Feaga et al.’s numerical values for the column densities 3.8 km and 7 km from the centre of the nucleus is poorly reproduced by this model in particular because of the low nightside emission.

A modiﬁcation of the KS model, the Estimated model, which uses a different temperature model and different dust layer thick-nesses has been constructed to reproduce the column density at 3.8 km and 7 km distance from the nucleus. The approach was here rather successful in reproducing both the main features of the radi-ance map and the inferred column densities.

Finally a small variation of the Estimated model was also simu-lated. In this variation, the temperature on the surface is every-where equal to 200 K. The column densities are slightly higher compared to the Estimated model, but the features and the distri-bution are not influenced significantly. We conclude that gas emis-sion from the nightside of the nucleus is necessary to match the column densities given by Feaga et al. The calculations also show (e.g. in the KS model) that lateral flow (i.e. roughly parallel to the surface) of gas, which is a known and expected phenomenon (e.g.Keller and Thomas, 1989) and also present in our simulations, is not sufficient to populate the inner coma above the nightside to the required density. Raising the surface temperature to allow H2O

outﬂow from the nightside surface in the model is required. The observations ofFarnham et al. (2012)suggest that dust emission from the nightside does indeed occur and that emission may occur from regions which have been in darkness for at least 4 h. Hints of nightside activity have been seen previously in ﬂy-by observations of other comets (e.g.Sagdeev et al., 1987). Here, in the absence of further constraints, we have used an approach which implies a hemispheric, homogeneous outgassing from the nightside. How-ever, it should be clear that this approach does not offer a unique solution to the observed distribution.

The kinetic energy of the gas when it has its final velocity is, at the source, given in the form of a temperature. We have used here the surface temperature distribution of Groussin et al. (2007a, 2007b) to set this temperature as a constraint. As pointed out above, we are in essence assuming that the gas is emitted from an ice surface at approximately 200 K but immediately accommo-dates to the observed surface temperature as it flows through the surface layer before leaving the comet. The relationship between effusion velocity and initial temperature has been given above and would suggest a decrease in the outflow velocity of around 20% if the temperature of the gas exiting the surface is close to the free sublimation temperature of water ice.

As shown in Fig. 15, the velocity of the gas reaches almost 900 m/s at the edge of the simulation field. This is considerably higher than assumed byFeaga et al. (2007), for example, who con-sider 500 m/s in their analysis. The model also shows that lateral flows near the surface can locally reach 300 m/s with ratios of lat-eral to radial outflow velocity approaching 1 in some areas. This has immediate implications for simulations of the near nucleus region. For example, one cannot simulate only part of the surface unless the simulation encompasses all emission/activity. There are two main reasons. First, lateral flows entering the simulation

Fig. 15. Speed of the gas along the model x-axis. The speed is given in m/s. The data are for the Estimated case.

(13)

Fig. 16. Flow ﬁeld parameters in cuts through the origin of the coordinate system. From left to right, one sees the yz-plane, the xz-plane and the xy-plane. Note that none of these plains is the observation plane, but the xz-plane comes closest to it. Top row: The number density: The colour bar is logarithmic and goes from 1015_#/m3_{(dark blue) to}

1018_#/m3_{(red). 2nd row: The mean free path: The colour bar is logarithmic and goes from 1 m (dark blue) to 100 m (red). 3rd row: The non-equilibrium parameter (T} rot/Ttr):

The colour bar is linear and goes from 0.9 (dark blue) to 1.2 (red). Bottom row: The total temperature: The colour bar is linear and goes from 50 K (dark blue) to 250 K (red). (For interpretation of the references to colour in this ﬁgure legend, the reader is referred to the web version of this article.)

(14)

Fig. 17. Velocity components in cuts through the origin of the coordinate system. From left to right, one sees the yz-plane, the xz-plane and the xy-plane. Top row: Total speed in m/s and stream lines. The colour bar is linear and goes from 200 m/s (dark blue) to 900 m/s (red). Second row: Radial velocity: The colour bar is linear and goes from 150 m/ s (dark blue) to 900 m/s (red). Third row: Velocity perpendicular to the radial direction: The colour bar is linear and goes from 20 m/s (dark blue) to 300 m/s (red). Bottom row: Velocity ratio: The colour bar is logarithmic and goes from 0.01 (dark blue) to 1 (red). (For interpretation of the references to colour in this ﬁgure legend, the reader is referred to the web version of this article.)