A&A 374, 326-336 (2001)
DOI: 10.1051/0004-6361:20010742

Effects of irregular shape and topography in thermophysical models of heterogeneous cometary nuclei

P. J. Gutiérrez - J. L. Ortiz - R. Rodrigo - J. J. López-Moreno

Instituto de Astrofísica de Andalucía-CSIC, Apartado 3004, 18080 Granada, Spain

Received 9 February 2001 / Accepted 21 May 2001

Abstract
Several improvements in the thermophysical model by Gutiérrez et al. (2000) have been included in a new code to specifically deal with fully irregular cometary nuclei. Also, the new code allows for the inclusion of regions with different ice to dust ratios, regions of different albedos and regions of different emissivity. The new model has been applied to groups of irregular bodies characterized by 3 statistical parameters, the so-called Gaussian random shapes. In simulations, these bodies rotate steadily around their maximum inertia moment axes. The results of the runs show that the main conclusions of Gutiérrez et al. (2000) still hold, and some new features are observed:
1) In general, very irregular objects have higher water production rates than spheres of the same radius for most of the orbital period. The fact that an irregular object has a larger area than the sphere cannot explain the differences in water production. The main differences appear to be a consequence of its topographic features. Also, topography can diminish the pre- and post-perihelion asymmetries in the lightcurves.
Concerning the results for plausible albedo and icy fraction area distributions, 2) the mean water production of a comet with an albedo distribution on the surface is equal to the water production of a homogeneous comet with an albedo equal to the mean albedo of the distribution. The same result is obtained for icy fraction area distributions. 3) Close to perihelion, objects with icy fraction area distributions have nearly the same productions as fully water ice objects. 4) The largest diurnal oscillations in the synthetic lightcurves result from the irregular shape, whereas albedo and icy fraction area inhomogeneities induce oscillations of only a few percent.

Key words: comets: general - solar system: general


1 Introduction


  \begin{figure} \par\begin{tabular}{ c c c} {\scriptsize OBJECT 1 } & {\script... ...sizebox{!} {4cm}{\includegraphics{1128f1c.eps}} } \end{tabular}\par\end{figure} Figure 1: Objects used.
Open with DEXTER

At present, there are a number of reasons to suspect that most of the comet nuclei are highly irregular. Clear evidence exists for some comets, such as the close up views of comet P/Halley taken by the Giotto spacecraft (Keller et al. 1986) and the images by the Vega spacecraft (Sagdeev et al. 1986), as well as the variability in light curves of some bare comet nuclei (e.g. Meech et al. 1997; Meech et al. 2000; Lamy et al. 2000). However, little effort has been made to develop thermophysical models for irregular nuclei. The important effect of an aspherical nucleus on the circumnuclear coma has been stressed by the pioneering work of Crifo & Rodionov (1997b), who showed remarkable coma asymmetries arising from a particular aspherical object. Other works have stressed the importance of negative relief (craters) in detailed thermophysical nucleus models (Colwell 1997), but it is only recently that a few authors have addressed the issue of modeling fully irregular comet nuclei of arbitrary shape (e.g. Enzian et al. 1999a; Gutiérrez et al. 1999; Gutiérrez et al. 2000).

Gutiérrez et al. (2000) (henceforth Paper I) presented a model which provides the outgassing curves and surface temperatures for two homogeneous irregular bodies which had some symmetries and were smooth enough so that some simplifying assumptions could be made. For any combination of orbital parameters, rotation axis orientation, spin period, and physical properties of the nucleus (albedo maps, emissivity, and thermodynamical properties) the heat diffusion equation in the surface normal direction was solved, with the energy balance at the surface (which included shadowing effects), as a boundary condition. This was done for all the cells in which the objects were divided and for a complete orbital revolution around the Sun. This allows the calculation of outgassing and temperatures for all the cells, as a function of heliocentric distance.

An improved model that deals with any arbitrary fully irregular shape is presented here. Furthermore, the model allows for the inclusion of albedo, emissivity and icy fraction maps; in other words, the surface does not have to be homogeneous. From the dynamical point of view, the bodies are allowed to spin around their maximum inertia moment axes. The model has been applied to several random Gaussian shapes which were used by some investigators in the past to model asteroidal shapes (e.g. Lagerros 1997; Muinonen 1998; Muinonen & Lagerros 1998).

We must emphasize that the model is not intended to fully describe the behavior of a particular comet nucleus, but to advance in several aspects which have to do with shape and to show the fundamental differences with respect to spherical nuclei. Thus, a wealth of phenomena have been oversimplified, such as dust mantling (e.g., Brin & Mendis 1979; Fanale & Salvail 1984; Rickman & Fernández 1986), gas flow through pores (e.g., Mekler et al. 1990; Prialnik 1992; Benkhoff & Boice 1997), sublimation of several volatiles (e.g., Fanale & Salvail 1987, 1990; Espinase et al. 1991, 1993), quasi 3-D heat and gas diffusion (Enzian et al. 1997).

2 The comet model


  \begin{figure} \begin{tabular}{ccc} {\scriptsize$\sigma =$0.2, $\Gamma =$ 15... ...{0.8cm}{!}{\includegraphics{1128f2d.ps}}} } \\ \end{tabular} \par\end{figure} Figure 2: Some albedo distributions used in this study.
Open with DEXTER

2.1 Physical processes

The comet model is based on the following assumptions and hypotheses. The nucleus is a non-homogeneous irregular body composed of crystalline water ice and dust aggregates. The sublimation occurs only at the surface and thermodynamic equilibrium between the sublimated gas and the surface ice is assumed. The local sublimation rate, Z, is calculated by using the so-called Hertz-Knudsen sublimation rate, $Z_{\rm HK}$. With regard to the vapor pressure over the surface, the approximation for the Clausius-Clapeyron equation given by Fanale & Salvail (1984) is considered valid. This expression is

\begin{displaymath}% P_{\rm v} = 3.56 \times 10^{12} \cdot \textrm{e}^{{-6141.667}/{T_{{\rm s}}}} \end{displaymath} (1)

where $T_{\rm s}$ is the surface temperature. Currently, the model does not account for recondensation, collisions or photochemistry of the sublimated gas. Collisions play an important role in the inner coma gas flux and in the local production rate as well because they control the back-scattered flux toward the surface (e.g.: Crifo & Rodionov 1999). This must be specially considered with irregular bodies, where the concavity can lead to higher local pressure and temperature. In order to evaluate the effect of returning molecules to the nucleus on the gas production, a complex hydrodynamic coma model must be coupled to the nucleus model. Thus, the model presented here provides an approximation to the production by irregular bodies. Under these assumptions, the total gas production rate, Q, is
 
$\displaystyle % Q=\sum_{i=1}^{n_{\rm c}} s_{i}\,f_i\,Z_{{\rm HK},i} = \sum_{i=1... ...i\,\frac{P_{\rm v},i}{\sqrt{2\pi mkT_{{\rm s},i}}} \; \textrm{molec.~s$^{-1}$ }$     (2)

where the summation is over all the cells ($n_{\rm c}$) in which the surface is divided because of computational requirements, $P_{\rm v,i}$ and $T_{{\rm s},i}$ are the vapor pressure and surface temperature respectively, m the molecular mass of the gas and k is Boltzman's constant. si and fi are the area and the icy fraction area of cell i, respectively. The latter is a correction for the fact that the net sublimation rate from a dusty ice surface is lower than the one from a pure ice surface. Considering a dust to gas ratio pattern on the surface, $\Re (\theta,\phi)$ (where $\theta$ and $\phi$ are the cometocentric latitude and longitude respectively), and if the water ice and dust densities are $\rho_{\rm H_2O}$ and $\rho_{\rm d}$respectively, the icy fraction area, $f(\theta,\phi)$ (Crifo & Rodionov 1997a), on the surface is calculated as

\begin{displaymath}% f(\theta,\phi)=\frac {1}{1+\frac{\rho_{\rm H_2O}}{\rho_{\rm d}} \, \Re (\theta,\phi)}\cdot \end{displaymath} (3)

Currently, the fraction $f(\theta,\phi)$ is maintained constant during the orbital evolution because erosion and dust mantle development processes are not dealt with in the model.
The surface temperature for cell i is calculated according to the energy balance equation at the surface for that cell, i.e.:
 
$\displaystyle % \frac{S_{0}(1-A_{{\rm v},i})\textrm{cos} (\theta_{i})C_{i}}{r_{... ...\rm s},i}) + \kappa (T) \frac {\partial T}{\partial z} \Bigg \vert _{{\rm s},i}$     (4)

where
S0 is the solar constant;
$A_{{\rm v},i}$ the Bond albedo of cell i;
$\theta_i$ the Sun-zenith angular distance as seen from cell i;
Ci=1 if cell i is illuminated;
Ci=0 if cell i is not illuminated (at night or in shadow);
$\epsilon_i$ the IR emissivity of cell i;
$T_{{\rm s},i}$ the surface temperature of cell i;
$\sigma$ Stefan-Boltzman's constant;
H the sublimation latent heat;
$Z_{{\rm s},i}(T_{{\rm s},i})$ Hertz-Knudsen sublimated gas flux;
$\kappa (T)$ the thermal conductivity;
Fi,j the viewing factor;
fi the icy fraction area
and z the normal direction to the surface.
The first term of the lefthand side of Eq. (4) is the solar energy input. This term does not take into account the radiative transfer properties (attenuation, multiple scattering, emission) of the inner coma. This approximation is based on the results from Salo (1988). His calculations seem to suggest that, in general, the total energy input to the nucleus is only weakly dependent on the opacity of the coma, the radial distribution of the dust or the details of the extinction processes. Again, a complete description would need a radiative and gas kinetic coma model, especially in the regions of hight concavity.

As an improvement of the model described in Paper I, the energy balance equation (Eq. (4)) accounts for the absorbed radiation of scattered and reradiated energy from the comet surface itself. This additional heating source is the second term of the lefthand side of Eq. (4) where the summation is over all the visible surface cells from cell i. In this term, the factor Fi,j (defined as the fraction of radiative energy from one surface cell directly hitting another cell, Lagerros 1997) can be approximated by

 \begin{displaymath}% F_{ij}=v_{ij} \cdot \frac{\textrm{cos} (\alpha_{i}) \cdot \textrm{cos} (\alpha_{j}) }{\pi d^2} s_j \end{displaymath} (5)

where $\alpha_{i}$ is the angle between the normal to cell i and the line connecting this cell with cell j, d is the distance between the two cells, and si is the area of cell i. The factor vi,j will be 1 if cell j is visible from cell i. Expression (5) is approximately correct under the assumptions of small planar cells and Lambert's law of diffuse reflection.

Lagerros (1997) introduced the self heating parameter, $\chi$ as a measure of the importance of this energy source term. This parameter is defined as the ratio between the power directly striking the surface itself to the total power emitted from the surface. Thus, if the incident energy flux on the cell i is Gi and the emitted flux from that cell is Ji, the self heating parameter is

\begin{displaymath}% \chi=\frac {\sum_i G_i s_i}{\sum_i J_i s_i}=\frac {\sum_i \sum_{j\not=i}F_{ij} s_i}{\sum_i s_i} \end{displaymath} (6)

where the summation is over all the surface cells and si is the area of cell i. It can be seen that this parameter depends on geometrical factors only. Thus, this parameter can be considered as a measure of the concavity of an irregular object. Also, since the radiative and dynamical effects of the coma on the surface are more important as the concavity increases, this parameter could be a rough measure of the importance of these effects.

The righthand side of Eq. (4) includes the terms of thermal reradiation (first), sublimation (second) and thermal diffusion to/from the interior (third). Horizontal heat diffusion and the flow of gas through pores are not considered. In order to estimate the energy that diffuses to or from the interior of the nucleus, the unidimensional heat diffusion equation must be solved, i.e.:

 \begin{displaymath}% \rho \cdot C(T) \frac{\partial T}{\partial t} = \frac {\par... ...} \Bigg [ \kappa(T) \cdot \frac{\partial T}{\partial z} \Bigg] \end{displaymath} (7)

where $\rho$ is the comet density, C(T) the specific heat capacity, $\kappa (T)$ the thermal conductivity, T the temperature and z the depth from the surface. Equation (7) must be integrated with two boundary conditions. The first boundary condition, located at the surface, is Eq. (4) and the second one is placed taking into account the thermal skin depth ($\delta$). $\delta$ is the depth at which a sinusoidal temperature variation is reduced by a factor "e''. Thus, the heat diffusion equation is integrated down to a depth equal to several times the thermal skin depth for a heating cycle of years. At that depth, the thermal gradient can be assumed to be 0. Therefore,

\begin{displaymath}% \frac{\partial T}{\partial z} \bigg \vert _{n \, \delta} = 0 \end{displaymath} (8)

is the second boundary condition needed to solve Eq. (7). With regard to the specific heat and thermal conductivity of the dust-ice mixture, the effective values used are

 \begin{displaymath}% C=\frac {\tilde{\rho}_{\rm H_2O}} {\rho} \, C_{\rm H_2O } % + \frac {\tilde{\rho}_{\rm d }} {\rho} \, C_{{\rm d}} \end{displaymath} (9)

where $\tilde{\rho}_{\rm H_2O}$ and $\tilde{\rho}_{\rm d}$ are the mass of water ice and dust, respectively, per unit volume in the nucleus comet, and $C_{\rm H_2O }$ and $C_{\rm d}$ the specific heat of water ice and dust, respectively. On the other hand,

 \begin{displaymath}% \kappa = h \left (\frac {\tilde{\rho}_{\rm H_2O}} {\rho_{\r... ...tilde{\rho}_{\rm d}} {\rho_{\rm d}} \, \kappa_{\rm d} \right ) \end{displaymath} (10)

where $\kappa_{\rm H_2O}$ and $\kappa_{\rm d}$ are the heat conductivity of water ice and dust, respectively, and $\rho_{\rm H_2O}$ and $\rho_{\rm d}$ their densities. The parameter h, the so-called Hertz factor, is introduced because the intrinsic thermal conductivity of a porous medium is smaller than that of the corresponding compact material. This parameter is a function of the density, size of contact areas between grains and porosity.

2.2 The irregular shape

With regard to the nucleus shape, the current model works with fully irregular shapes. In this study, the so-called Gaussian random shapes (Peltonieni et al. 1989; Muinonen 1996a) have been used to define the nucleus surface.

2.2.1 The Gaussian random shape

The main idea of these irregular shapes is to describe the radii as Nrandom variables (where N is the number of grid points) $\bf {r}=(r_1,...,r_N)$ with corresponding spherical coordinates $\Omega = (\theta_1,\phi_1;...;\theta_N,\phi_N)$, obeying the lognormal multivariate statistics. The distribution is totally characterized by the means, <ri>, the variances, $\sigma_i^2<r_i>^2$, and the covariance matrix, $\Sigma_{\rm s}$. The elements of this matrix are defined in terms of the correlation function, $C_{\rm r}$. We have assumed, as Lagerros (1997) and Muinonen and Lagerros (1998), that the radii have equal means, <r>, equal variances, $\sigma^2<r>^2$, and that the correlation function as a function of the angle between two directions, $\gamma$, is

\begin{displaymath}% C_{\rm r}(\gamma,\Gamma)=\frac{(1+\sigma^2)^{C_{\rm s}(\gamma,\Gamma)}}{\sigma^2}-\frac{1}{\sigma^2} \end{displaymath} (11)

where

\begin{displaymath}% C_{\rm s}(\gamma,\Gamma)=(1-\omega) P_2(\textrm {cos}(\gamm... ...rm {sin}^2(\gamma / 2)}{\textrm {sin}^2 (\Gamma / 2)} \right ) \end{displaymath} (12)

where $\omega$ is a so-called "weight" and $\Gamma$ is a correlation angle. $P_2(\textrm {cos}(\gamma))$ is the Legendre polynomial of second degree. Muinonen & Lagerros (1998) show that this autocorrelation function can be fitted to the autocorrelation function of real asteroidal shapes. It can be shown that the N radii can be obtained from

\begin{displaymath}% r(\theta,\phi)=\frac {<r>\textrm{exp}(s(\theta,\phi))}{\sqrt{1+\sigma^2}} \end{displaymath} (13)

where
 
$\displaystyle % s(\theta,\phi)=\sum_{l=0}^\infty \sum_{m=0}^l P_l^m(\textrm{cos}(\theta)) (a_{lm} {\rm cos}(m\phi) + b_{lm} {\rm sin}(m\phi))$     (14)

provided that the coefficients alm and blm of a spherical harmonics expansion for $s(\theta,\phi)$are independent Gaussian random variables with zero means and equal variances, $\beta_{l,m}$,

 \begin{displaymath}% \beta_{lm}=(2-\delta_{m0}) \frac {(l-m)!} {(l+m)!} \, c_l \cdot \textrm{ln}(1+\sigma^2). \end{displaymath} (15)

In Eq. (14) Plm are the associated Legendre functions and in Eq. (15) cl are the coefficients of a Legendre expansion for $C_{\rm s}(\gamma)$. A more complete description of the statistic of Gaussian random shapes is done in Muinonen (1996a) and Muinonen (1996b) where the author shows how these shapes can be generated. Moreover, a discussion about the features of the autocorrelation function is given in these references. An important feature of these irregular shapes is that they are fully characterized by the four parameters $\sigma$, $\Gamma$, $\omega$ and <r>.

2.3 Albedo and icy fraction area patterns

For the albedo and icy fraction distributions, we have adopted similar distributions as for the radii. Thus, a normal multivariate statistic describes the albedo distribution on the surface.

\begin{displaymath}% A(\theta,\phi)=\frac {<A>}{\sqrt{1+\sigma_{A}}} {\rm exp} (s(\theta,\phi)) \end{displaymath} (16)

where <A> is the mean albedo, $\sigma_A^2 <A>^2$ the variance and s verifies the same properties as before. As for the radii, these distributions are characterized by <A>, $\sigma_A$ and $\Gamma$ (for the albedo and icy fraction area, $\omega=0$). For the icy fraction area distribution, we have used the distributions:

\begin{displaymath}% f(\theta,\phi) = 1- A(\theta,\phi). \end{displaymath} (17)

One of aims of this study is to describe in qualitative terms the effect of albedo and icy fraction area patterns on the surface.
 

 
Table 1: Orbit and rotation parameters. The orbital elements correspond to the orbit of a typical Jupiter family comet. The spin axis orientation is given following the notation of Sekanina (1981): I, obliquity of the orbit plane to the equator, and $\Phi $, the angle of the subsolar meridian at perihelion from the ascending node.
Semimajor axis   3.5 UA
Eccentricity   0.6
Comet rotation period   24 h
Rotation axis orientation Case 1 $I=90\hbox{$^\circ$ }$ $\Phi= 0\hbox{$^\circ$ }$
  Case 2 $I=0\hbox {$^\circ $ }$ $\Phi = 90\hbox {$^\circ $ }$
  Case 3 $I=45\hbox {$^\circ $ }$ $\Phi = 60\hbox {$^\circ $ }$



   
Table 2: Characteristics of the objects considered.
  OBJECT 1 OBJECT 2 OBJECT 3
RADIUS: (m) 1000 1000 1000
VOLUME: (m3) $4.56 \times 10^9$ $5.09 \times 10^9$ $4.98 \times 10^9$
VOL/VOL SPHERE: 1.09 1.22 1.20
AREA $({\rm m}^2)$: $1.52 \times 10^7$ $1.63 \times 10^7$ $1.58 \times 10^7$
AREA/AREA SPHE: 1.21 1.30 1.26
SELF HEATING: $1.96 \times 10^{-2}$ $3.31 \times 10^{-2}$ $8.0 \times 10^{-3}$

3 Parameters and model calculations

In the previous model (Paper I), the input of the code was the nucleus shape and the parameters defining the rotation and the orbit of the comet. In the new code, the albedo and the icy fraction area patterns must be specified as well. Table 1 shows the orbital and rotation parameters used in the calculations. To define the spin axis orientation, two angles must be specified. Thus, following the notation of Sekanina (1981), the rotation axis orientation is selected introducing the angle $\Phi $ of the subsolar meridian at perihelion from the ascending node, and the obliquity, I, of the orbit plane to the equator.
 

 
Table 3: Physical parameters.
PARAMETER   VALUE
Bond albedo $A(\theta,\phi)$ 0.1, 0.05 or distrib.(Fig 2)
Emissivity $\epsilon(\theta,\phi)$ 1 - A
Icy fraction area $f(\theta,\phi)$ 1.0, 0.9, 0.95 or distrib.
Comet bulk density $\rho$ 500 Kg/m$^{\rm 3}$
Ice water density $\rho_{\rm H_2O}$ 917 Kg/m$^{\rm 3}$
Dust density $\rho_{\rm d}$ 3500 Kg/m$^{\rm 3}$
Sublimation latent heat L 48600 J/mol
Water Ice heat capacity $C_{\rm H_2O }$ $0.09+ 0.00749 \cdot T \; {\rm J g^{-1}\ K^{-1}}$
Dust heat capacity $C_{\rm d}$ 1.2 $ {\rm J g^{-1}\ K^{-1}}$
Water ice conductivity $\kappa_{\rm H_2O}$ 567/T ${\rm W m^{-1}\ K^{-1}}$
Dust conductivity $\kappa_{\rm d}$ 4.2 ${\rm W m^{-1}\ K^{-1}}$
Hertz parameter h 0.18


With regard to the shape, the Gaussian random shapes used in this study as well as comet nucleus surface and their parameters are shown in Fig. 1. These surfaces are divided in 2320 triangular cells (30 latitude bins and 40 longitude bins). Table 2 summarizes some of the characteristics of the objects of the Fig. 1. Figure 2 also shows some of the albedo distributions studied and their parameters.

The physical parameters used in the model are summarized in Table 3. The comet bulk density used is close to the densities calculated by Rickman (1986, 1989) for P/Halley and by Solem (1994) and Asphaug & Benz (1994) for Shoemaker-Levy 9. Concerning the water ice thermal parameters, $C_{\rm H_2O }$ and $\kappa_{\rm H_2O}$, the values given by Klinger (1981) have been adopted. For the dust, we have used the values for $C_{\rm d}$ and $\kappa_{\rm d}$ given by Ellsworth and Schubert (1983).

Concerning the Hertz parameter, h, Mendis & Brin (1977) assigned typical values for h of 10-3 to 10-5. Nevertheless, theoretical calculations to fit experimental results by Seiferlin et al. (1996) (considering gas flow throw pores also) give values for h that range from 0.001 to 0.01. This interval for h is the usual one used in some thermophysical models (e.g. Enzian et al. 1999b; Kossacki et al. 1999). In Paper I we used a larger value for h, equal to 0.1, in order to compare our results with the results from Colwell (1997). He adopted that value based on the discussion of Smoluchowsky (1982), which is previous to the experimental studies, about the thermal conductivity of porous water ice. In order to test the new code, we have maintained this high value for the Hertz factor. Thus, if thermal conductivity is given by Eq. (10), a reduction factor of 0.1 for a porous ice free dust nucleus requires a Hertz factor of 0.18 for the densities given in Table 3.

 

 
Table 4: Numerical parameters.
PARAMETER VALUE
Number of surface cells 2320, triangular
Time steps ${\approx}500$ per rotation period
Depth of temperature 5 $\delta$
calculations  
Number of depth bins in 64
heat diffusion equation  


3.1 Model calculations

The code starts at aphelion with an initial nuclear temperature of 70 K. At every time step, the new orbital position and the solar incidence angle for each cell are calculated. After that, the algorithm checks for night time or projected shadows (See Paper I) and the surface temperature is calculated for every cell. The correct calculations of the surface temperature for each time step require an iterative loop between the solutions of the non-linear system of the surface energy balance equations (because the surface temperature of a cell depends on the surface temperature of other cells) and the solutions of the heat diffusion equations. Since this calculation scheme is very time consuming, a different method has been adopted. For a cell i, the temperature calculated in the previous time step is used to calculate the term of reradiation from the other cells hitting cell i (second term on the lefthand side of Eq. (4)). If the time step is very small, this permits us to decouple the non-linear system of surface energy balance equations with only a small error. After that, the energy balance equation at the surface is solved for every cell by means of the Newton-Raphson method. Using the value obtained for the surface temperature, the heat diffusion equation is solved. This equation is integrated using Kirchoff's transformation and a Crank-Nicholson scheme.

With the value of the surface temperature obtained with this scheme and the time step adopted, the energy balance equation at the surface, Eq. (4), is very close to zero, always less than 1% of the solar energy input.

4 Results


  \begin{figure} \par\resizebox{8.8cm}{!}{\includegraphics{1128f3.ps}}\end{figure} Figure 3: Ratio of the water production rate of irregular objects 1, 2, and 3 to the water production rate of an equal area sphere. The spin axis orientation is $I=45\hbox {$^\circ $ }$, $\Phi = 60\hbox {$^\circ $ }$. Solid line: object 1, dash-dot-dot-dot line: object 2, dash-dot line: object 3. It should be noted that the data in all the figures are shown every 2000 time steps ( ${\approx }4$ rotation periods). For this reason the plots appear to show a longer period than the rotational period. All the curves are smoothed (averaging in a box of 10 points) to enhance the mean trend.
Open with DEXTER

4.1 Production of irregular objects


  \begin{figure} \par\resizebox{8.8cm}{!}{\includegraphics{1128f4.ps}}\end{figure} Figure 4: The same as Fig. 3 with a different spin axis orientation. In this plot, $I=0\hbox {$^\circ $ }$, $\Phi = 90\hbox {$^\circ $ }$. Solid line: object 1, dash-dot-dot-dot line: object 2, dash-dot line: object 3. All the curves are smoothed (averaging in a box of 10 points) to show the mean trend. The real short term oscillations are around 50% the mean production.
Open with DEXTER


  \begin{figure} \par\begin{tabular}{ c c c c} {\scriptsize OBJECT } & {\scripts... ...resizebox{1cm}{!}{\includegraphics{1128f5g.ps}}} } \\ \end{tabular}\end{figure} Figure 5: Sinusoidal projections of the surface temperatures of the sphere and of the object 3 for three orbital positions are shown. The spin axis obliquity is $45^\circ $.
Open with DEXTER

Figures 3 and 4 show the water production rates of the irregular objects studied relative to the water production rates of equal area spheres for two spin axis orientations. In general, it is expected that the different total area and the different cross sections of the different bodies viewed from the Sun will lead to very distinct production rates. In Figs. 3 and 4 it can be seen that, indeed, the production rate is strongly dependent on the geometrical shape and on the spin axis orientation.

The results obtained when the spin axis obliquity is 45$^\circ$ are shown in Fig. 3. With this spin axis orientation, all the objects studied have cross sections close to 25% of their total areas (close to 26% for object 3). The irregular objects are more productive at large heliocentric distances than equal area spheres. If the curves of Fig. 3 are multiplied by the area factors shown in row 5 of Table 2, the total production relative to the production of a 1 km sphere is obtained. Hence, it is clear that these irregular bodies have a larger water production than a sphere of the same radius. Thus, the radius of the sphere that produces the same outgassing as the irregular objects must be higher than the mean radius of the objects.

Nevertheless, the differences in total water production from the different bodies cannot be explained by the differences in the total area. Table 2 shows that the three objects have similar total areas and Fig. 3 shows very different water productions at large heliocentric distances. The large differences could be a consequence of topographic features. The irregular bodies can have several regions whose surface normal directions are nearly parallel to the Sun direction. These regions will have similar temperatures and gas production to that of the subsolar point and, therefore, these regions will control the sublimation at large heliocentric distances.

When the spin axis is in the orbital plane ( $I=90^\circ$), we obtain similar results. Irregular objects have a larger water production than equal area spheres, and, therefore, than a sphere with the same radius. These results confirm those of Paper I, obtained with a simplified model and "less irregular'' bodies.

Nevertheless, some exceptions to the above picture do occur for some spin axis orientations. With the spin axis perpendicular to the orbital plane (Fig. 4) the subsolar point circles the cometary equator and the cross section of the three irregular objects viewed from the Sun oscillate around 21%, 22% and 19% of their total area, respectively. As the relative cross sections are smaller than the relative cross section of the sphere (25%), the objects studied are less productive than equal area spheres. Among the three bodies, only object 2 (the most irregular one) has a larger total production than a 1 km sphere. This can be verified multiplying the ratio for the object 2 shown in Fig. 4 (about 0.9) by the ratio of the total area of object 2 to the area of the sphere shown in Table 2, (1.30). Thus, and because all the objects studied have the same mean radius (1 km), object 2 produces more outgassing than a sphere with the same radius. Looking at Table 2, it is clear that the difference in total area cannot explain totally the large water production of object 2. This body has a total area similar to object 3, but the former has a much larger total production than the latter. The differences in water production arise from topographic features.

In summary, from Fig. 4 it can be stated that, although in general irregular objects are more productive than spheres of the same radius, and, also, than equal area spheres, some irregular objects for certain spin axis orientations can have a lower water production than spheres of the same radius (despite the fact that their total area is much larger than that of the sphere).

As an example of surface temperatures patterns, Fig. 5 is shown. This figure depicts sinusoidal projections of the surface temperatures of the sphere and of the object 3 for three orbital positions when the spin axis is 45$^\circ$. Concerning the asymmetry of the lightcurves around perihelion, it is interesting to discuss the case of obliquity 45$^\circ$. The curves of Fig. 3 have asymmetry parameters, E, (as defined by Festou et al. 1990) given in Table 5.

 

 
Table 5: Asymmetry parameters of the lightcurves shown in Fig. 3.
  Sphere Obj. 1 Obj. 2 Obj. 3
Orbit 0.022 0.055 0.55 0.094
$r_{\rm h} < 2.5$ AU 0.042 0.089 0.084 0.139
$r_{\rm h} > 2.5$ AU -0.020 -0.034 -0.033 -0.045


With this spin axis orientation, the cross sections of the irregular objects change a similar amount close to perihelion. Thus, it appears that the differences in the asymmetry parameter arise from the topography. The presence of small regions with a production rate close to that of the subsolar point will maintain the production rate of objects 1 and 2 more symmetric or "constant'' around perihelion than the production rate of object 3, which has a much smoother surface than the other irregular objects.

Concerning the short term variability, we can observe that changes in the cross section viewed from the Sun as the body rotates give rise to very important oscillations in the lightcurve. For the objects studied these can be up to 50% of the mean production. This large amplitude is a particular result. Different surfaces will lead to different amplitudes.

4.1.1 Importance of self heating in irregular objects

  \begin{figure} \par\resizebox{8.8cm}{!}{\includegraphics{1128f6.ps}}\end{figure} Figure 6: Ratio of the water production rate of object 2 considering self heating to the water production rate of the same object without self heating. The spin axis orientation is $I=0\hbox {$^\circ $ }$, $\Phi = 90\hbox {$^\circ $ }$.
Open with DEXTER


  \begin{figure} \par\begin{tabular}{cc} \rotatebox{90}{\resizebox{6cm}{!}{\inclu... ...esizebox{1cm}{!}{\includegraphics{1128f7b.ps}}} & \end{tabular}\par\end{figure} Figure 7: Surface temperature of object 2 considering self heating minus surface temperature of the same object without self heating at perihelion. The spin axis orientation is $I=0\hbox {$^\circ $ }$, $\Phi = 90\hbox {$^\circ $ }$. In this figure, black areas represent a difference smaller than 1 K and white areas larger than 5 K.
Open with DEXTER

Figures 6 and 7 show an example of the effect of self heating. In Fig. 6, the ratio between the water production rate of the object 2 considering self heating and the water production rate of the same object without self heating is represented. The spin axis orientation considered was $I=0\hbox {$^\circ $ }$, $\Phi = 90\hbox {$^\circ $ }$. Considering the whole surface, the total heat input into the nucleus because of self heating is around 3% of the solar energy input. Nevertheless, locally, in surface regions of hight concavity, the effect of self heating can reach 20%, and even higher, of the solar energy input. These numerical results correspond to the particular cases studied. Different objects, with different self-heating parameters, will have different numerical values. Figure 6 shows that this local effect can increase water production rates by a 20%. Close to perihelion, the relative effect of reradiation is diminished because nearly all the direct solar input energy goes into sublimation. In Fig. 7 the difference between the surface temperature of the object 2 considering self heating and the surface temperature of the object 2 without self heating at perihelion is shown. In this figure, black areas represent a difference smaller than 1 K and white areas represent differences larger than 5 K. Since the water production rate scales exponentially with the temperature, small differences in surface temperature can lead to significant differences in water production rate. It can be seen that, at aphelion, the importance of self heating is a little bit higher than at perihelion. The reason is the relative importance of the reradiation term in the energy balance equation at the surface.

4.2 Effects of albedo and emissivity distributions


  \begin{figure} \par\resizebox{8.8cm}{!}{\includegraphics{1128f8.ps}}\end{figure} Figure 8: Ratio of the water production by a sphere with a distribution of albedo on the surface (dashed line: <A> = 0.05, $\sigma^2~=~0.6$, $\Gamma =20^\circ $, thin solid line: <A> = 0.1, $\sigma^2~=~0.2$, $\Gamma =20^\circ $) to the water production rate of a sphere with a uniform albedo. Thick solid lines are fits of the amplitude of the oscillations following Eq. (18). The spin axis obliquity is $0^\circ $.
Open with DEXTER


  \begin{figure} \par\resizebox{8.8cm}{!}{\includegraphics{1128f9.ps}}\end{figure} Figure 9: The same as Fig. 8 for the irregular object 2.
Open with DEXTER


  \begin{figure} \resizebox{8.8cm}{!}{\includegraphics{1128f10.ps}}\par\end{figure} Figure 10: Ratios of the water production rate of a comet with a uniform icy fraction area, f, (solid lines: f=0.90, dashed lines: f=0.95) to the water production rate from the same comet with no dust on the surface (f=1). The non-oscillating curves correspond to the sphere and the oscillating ones to the irregular object 1. The spin axis is perpendicular to the orbital plane.
Open with DEXTER


  \begin{figure} \resizebox{8.8cm}{!}{\includegraphics{1128f11.ps}}\par\end{figure} Figure 11: Ratio of water production rate by a spherical comet nucleus with an icy fraction area distribution characterized by $<f>\ = 0.9$, $\sigma = 0.2$ and $\Gamma = 15$ to the water production rate of a sphere with a uniform icy fraction area equal to 0.9. The spin axis is perpendicular to the orbital plane.
Open with DEXTER

Figures 8 and 9 are examples of some effects of the albedo distributions on the total production rate. All the spin axis orientations and all the objects show the same general behavior. The first result that can be drawn from these figures is that the mean production of the heterogeneous comet used in this study is equal to the production of a uniform comet with an albedo equal to the mean of the albedo distribution. This behavior is a consequence of the linearity of the production rate with the albedo for the range of albedos considered in the distributions (A is always less than 0.4). Furthermore, it can be also seen that the albedo distributions used produce small oscillations in the gas production curve (always smaller than 10% of the mean production). Concerning the amplitude, the larger the variance of the distribution, the larger the oscillations; as it was expected. For all the spin axis and all the objects studied, the general trend of the variation of amplitude of the oscillations with heliocentric distance is the one shown in Figs. 8 and 9. First, where the sublimation is not important from an energetic point of view ( $r_{\rm h} >\ \approx 2.5~{\rm AU}$) , the oscillations are larger and have an increasing amplitude toward perihelion. As sublimation acquires more importance when the comet approaches the Sun ( $r_{\rm h} \approx\ < 2.5~{\rm AU}$), the amplitude of the oscillations decreases. In general, for the albedo distributions studied the amplitude of the oscillations will depend on the rate of change of the mean production, $Q_{{\rm esf},<A>}$, with heliocentric distance. For the spin axis shown in Fig. 8, the amplitude of the variations can be suitably fitted to a linear expression of the mean production rate of change with heliocentric distance, i.e.:

 \begin{displaymath}% {\rm Envelope} = m\frac{1}{Q_{{\rm esf},<A>}} \, \frac{\partial{Q_{{\rm esf},<A>}}}{\partial{r_h}}+b \end{displaymath} (18)

where m and b are parameters that can be fitted. Thick solid lines in Fig. (8) (which follow Eq. (18)) show the variation of the amplitude of the oscillations with heliocentric distance for the case considered (sphere spinning around an axis perpendicular to the orbital plane).

In all the runs shown, the emissivity pattern on the surface was assumed to be equal to 1 $-A(\theta,\phi)$ where A is the albedo. Several additional runs have been done to study the effect of an emissivity different to 1-A. In those runs (which are not shown here and assuming a constant emissivity of 0.9) the main differences with respect to the 1-A emissivity case appear at large heliocentric distances, where the reradiation term is the most important one in the energy balance equation at the surface. There are no differences close to perihelion. The general effect, close to aphelion, is to increase the amplitude of oscillations which can reach 15% of the mean production for some albedo distributions.

4.3 Effects of icy fraction area distributions

The main effects of the icy fraction area distributions, f, are shown in Figs. 10 and 11. In Fig. 10, ratios between water production rate of a comet with a uniform icy fraction area (equal to 0.9 (solid lines) and to 0.95 (dashed lines)) and the water production rate from the same comet with f=1 are shown. The spin axis is perpendicular to the orbital plane. The non oscillating curves correspond to the sphere and the oscillating ones to the irregular object 1. This figure shows that the water production rate from the dusty comet of this study is larger than the one from the pure ice comet for most of the orbital period. Close to perihelion, it can be seen that water production rate does not depend on the local icy fraction area. This result is observed for all the objects, all the spin axis orientations and all the icy fraction area distributions studied.

It might be expected that far away from perihelion, and if thermal diffusion into the nucleus is not taken into account, the presence of dust on the surface can lead to a smaller water production rate than the one from a pure water ice nucleus. The obvious reason would be that the dusty nucleus has a smaller sublimation area and that close to aphelion the equilibrium temperature is controlled by the reradiation term in Eq. (4). This would be true even for a dusty nucleus with very similar thermal properties to the ones of a pure water icy nucleus. Nevertheless, if the thermal diffusion is taken into account and thermal properties of the dusty nucleus are represented by Eqs. (9) and (10), the results are quite different. Figure 10 shows that the water production rate from our particular dusty comet is larger than the one from a pure water ice comet. The reason is that thermal conductivity of a dusty nucleus is smaller than the conductivity of a pure water ice nucleus. This means that the former diffuses into the interior a smaller amount of energy than the latter. Thus, the dusty nucleus reaches a higher surface temperature than the pure water icy nucleus.

The differences between both water production rates start to diminish when the sublimation term becomes the dominant one from an energetic point of view. Close to perihelion, where the effect of thermal diffusion is negligible and the equilibrium temperature is given, essentially, by sublimation, comets produce similar amounts of gas regardless the mean icy fraction area. That means that the dusty nucleus reaches an equilibrium temperature which increases as thermal conductivity decreases (or, in this case, the local icy fraction area is smaller). This increase in temperature nearly compensates the reduction of sublimating area because of the presence of dust on the surface.

The oscillations shown in Fig. 10 can be explained as follows. We have observed that the water production rate oscillations and the cross section oscillations due to rotation of an irregular body are out of phase. This lag depends on the value of thermal conductivity. Thus, water production curves from comets with different icy fraction area are out of phase because thermal conductivity depends on the icy fraction area.

Figure 11 shows the ratio between the water production rate by a spherical comet nucleus with an icy fraction area distribution characterized by $<f>\ = 0.9$, $\sigma = 0.2$ and $\Gamma = 15$ and water production rate of a sphere with a uniform icy fraction area equal to 0.9. Similar plots are obtained for different icy fraction area distributions regardless the comet shape or the spin axis orientation. It can be seen that icy fraction area distributions produce very similar effects to those derived from the albedo distributions. Mean water production rate by a comet depends only on the mean value of the icy fraction. The inhomogeneities of the ice distribution studied may only introduce oscillations smaller than 8% of the mean production rate.

5 Summary and conclusions

With the model presented here, surface temperatures and water production rates for a fully irregular cometary nucleus with a heterogeneous surface can be calculated for any combination of orbit and rotation parameters. Since the code can work with realistic cometary shapes, a new tool is now available to analyze the data obtained by future cometary missions. The model has been applied to several simulated cometary shapes with different albedo and local icy fraction area patterns on the surface. Some results from the simulations have been presented.

With regard to the effects of the nucleus shape on water production rates, in general, the water production of a very irregular comet is larger than that of a sphere with the same radius, and also than equal area spheres for most of the orbital period. However, some irregular objects, for certain spin axis orientations, have lower water production than equal spheres because they show a relatively low cross section as viewed from the Sun. Also, it can be concluded that the differences in total area cannot explain the differences in total water production. Irregular objects with similar total areas have very different water productions. Therefore, topographic features appear to be an important cause of large differences in water production of irregular objects. Very irregular objects can have more regions whose surface normal directions are nearly parallel to the Sun direction than spherical comets. These regions will have similar temperatures and gas production to that of the subsolar point.

Irregular shape can lead to large seasonal asymmetries in lightcurves. Elongated objects rotating around a principal axis which is inclined toward the orbital plane will change their cross sections as viewed from the Sun close to perihelion. As a general trend, the larger and more asymmetric the change in cross section around perihelion, the larger the seasonal asymmetry in the lightcurve. On the other hand, topography acts in the contrary way. Topography tends to diminish this pre- and post-perihelion asymmetry. The presence of small regions with production rates close to the production rate of the subsolar point maintain the total gas production of very irregular objects more symmetric around perihelion.

These conclusions are in agreement with the main conclusions obtained in Paper I, and with the conclusions obtained by Colwell (1997) for a spherical nucleus with craters. Nevertheless, these results could be particular cases and cannot be totally generalized. Further studies about the effects of irregular shape in water production are necessary.

Concerning the effects of a heterogeneous surface, the albedo and local icy fraction area distributions considered in this study produce similar effects, regardless of the nucleus shape. These effects appear at large heliocentric distances. When sublimation becomes the most important term from an energetic point of view, the effects of the albedo distributions or of the mean local icy fraction area diminish. As expected, the water production rate of a comet with an albedo or an icy local fraction distribution on the surface is equal to the water production of a homogeneous comet with the mean albedo or the mean icy fraction area. From the simulations of a comet with dust on the surface, we have found that a dusty comet (with an icy fraction area equal to 0.9) produces similar amounts of gas than a pure water icy comet at perihelion and even more at large heliocentric distances. The presence of dust on the comet diminishes the thermal conductivity of the pure water ice leading to higher surface temperatures than in a pure water ice comet. This result depends indeed on the parameterized expression for thermal conductivity. Thus, different results can be obtained from different expressions for thermal conductivity in terms of the local icy fraction area.

It can be observed that geometrical shape, albedo distributions and local icy fraction area patterns on the surface give rise to short term variability in comet lightcurves. From our simulations, the changes in cross section because of rotation result in diurnal oscillations which are larger than those caused by albedo or icy fraction area inhomogeneities. The amplitude of the former can reach up to 50% of the mean production. Oscillations induced by the albedo or local icy fraction area distributions on the surface used in this study are always lower than 10% of the mean production. Different values can be obtained from different distributions.

References

 

Copyright ESO 2001