Thermal modelling of water activity on comet 67P/Churyumov-Gerasimenko with global dust mantle and plural dust-to-ice ratio

Abstract

We perform a thermo-physical analysis on water activity of comet 67P/Churyumov-Gerasimenko (67P). The sublimation of water is assumed to occur from beneath a global, desiccated dust mantle over the irregular-shaped nucleus. The concept of two thermal models, the recipe of model formulation and the strategy of application to comet 67P are described. For an accurate and efficient evaluation of energy input by insolation and self-heating over the nucleus, a Landscape data base is devised based on polyhedral shape models of the nucleus. We apply the thermal models to investigate the impact of certain parameters of nucleus properties on water production. It is found that the measured water production of 67P can be overall attributed to sublimation of water ice with a mass abundance of a few to 10 per cent beneath a uniform dust mantle of several millimetres to one centimetre in thickness. Insofar as 67P is concerned, we argue against the necessity to invoke assumptions on localized water activity, or on the distinction of active/dormant surface areas.

solid state: refractory, solid state: volatile, methods: numerical, comets: general, comets: individual (67P/Churyumov-Gerasimenko)

1 INTRODUCTION

Comet nuclei are frozen mixtures of volatiles and refractories that crumble and burst in the heat of solar irradiation, producing prominent comae of escaping gases and floating dust fragments, when they encounter the Sun (Whipple 1950; Keller 1989). The most abundant volatile in comet nuclei is water, followed by other species, such as CO₂, CO, with higher volatility and, thus, preserved in deeper, more insulated layers (Festou, Keller & Weaver 2004). While the detailed mechanisms remain to be resolved, the principal characterization of cometary activity as dust ejection driven by volatile outgassing seems definitive. That is, the nucleus structure may be weakened or even disrupted from the interior, at least in part, by the sublimation of volatile ices via pressurization. The dynamics of fragments is then determined by the drag forces of the outflowing gases, gravity of the nucleus and, depending on the sizes and reflectance properties of the ejecta, solar radiation pressure (Huebner et al. 2006).

It is known that comet nuclei are dark objects in general, e.g. reflecting a few per cent of sunlight in the visible range. Exposure of pure water ice over nucleus surface is limited. Instead, water ice must be concealed beneath a prevalent, desiccated dust mantle that may significantly influence the cometary activity and evolution (Prialnik & Bar-Nun 1988). Since the pioneering space enterprise to 1P/Halley, when the spacecraft Giotto, Vega 1 and 2, caught the first-ever clear views of a comet, all nuclei¹ resolved thus far are irregular (non-spherical) in shape, with four being distinctly elongated. It has become increasingly clear that the mass of comet nuclei is dominated by dust with lesser amounts of volatiles (Keller 1989; Küppers et al. 2005).

There exists an entrenched belief that cometary water activity arises only from localized, isolated areas, presumably in correlation with the scanty exposure of water ice on the nucleus surface. This conception is somewhat reinforced by outcomes of thermal modelling of cometary water activity, assuming a bare icy surface that almost invariably overestimated the water production compared with observations. The model results are sometimes scaled down by a small factor, say a few per cent, that is interpreted to represent a macroscopic fraction of the ‘active’ nucleus surface. The necessary consequence of this intuition, namely the existence of significant ‘dormant’ or ‘inactive’ areas, is invalidated by observations of 67P/Churyumov-Gerasimenko (67P), collected by Rosetta spacecraft during its two-(and-plus-)year rendezvous with the comet. As a typical Jupiter-family member, the irregular-shaped nucleus of 67P has a dust-to-(water-)ice ratio of at least 6 and a volatile-free surface (Sierks et al. 2015; Capaccioni et al. 2015; Rotundi et al. 2015; Filacchione et al. 2016). Yet, the distribution of water activity over the nucleus is global, and the existence of inactive areas is elusive (Fougere et al. 2016b,a; Fulle et al. 2016a; Hu et al. 2017).

The formation of the dust mantle and its influence on cometary activity have been subjects of extensive theoretical and experimental research over decades (see e.g. Prialnik & Bar-Nun 1988; Grün et al. 1991; Kömle et al. 1992; Skorov & Rickman 1995; Gundlach, Skorov & Blum 2011, for a glimpse of the literature). The dust mantle restricts the diffusion of vapour, thus reducing the sublimation flux of water ice underneath. In addition, the icy interior is partly insulated by the overlying dust mantle, which regulates the variation of ice temperature at which sublimation occurs. These effects have been modelled by Keller et al. (2015) for estimating the water activity of 67P. On the other hand, the water production from a fully icy nucleus, as is sometimes assumed in thermal models, always indicates an upper bound, whereas the production from a ice–dust mixture ought to be lower. Therefore, the existence of a prevalent dust mantle and a significant dust-to-ice ratio, e.g. greater than 6 (Rotundi et al. 2015; Fulle et al. 2016c), will lower the water production from the nucleus.

The present study is aimed to assess the impact of a global dust mantle of uniform thickness and a plural dust-to-ice ratio on water activity of 67P. Towards this end, we introduce the formulation of two basic thermal models and the strategy of model parametrization for estimating water production of the comet (Section 2). We will explain in detail the method for efficient evaluation of the energy input over the nucleus of 67P, where irregular topography gives rise to complex effects of shadowing and nucleus thermal irradiation (Section 3). We compare some benchmark solutions of the two models to illustrate the influences of model parameters, such as dust mantle thickness and water ice abundance (Sections 4 and 5). The two thermal models will be applied to interpret the observed or in situ measured evolution of the water production of 67P, where results are compared with those of other thermal models for 67P (Section 6).

2 MODELING WATER SUBLIMATION FROM DUSTY-ICE OR DUST-MANTLED NUCLEUS

2.1 Bare, icy nucleus surface

In the simplest case of nucleus as a uniform mixture of ice and dust, as a dusty ‘snowball’ (Whipple 1950),² the outgassing of water is governed by the energy balance at the surface, namely (Fanale & Salvail 1984),

\begin{equation} Q_{(0)} = \left. -\kappa \frac{\partial T}{\partial x} \right|_{0^+} + \sigma \varepsilon T^4_{(0)} + \ell Z_{(0)}\quad ({\rm W}\, {\rm m}^{-2}), \end{equation}

(1)

where Q₍₀₎ is the energy input, i.e. the absorption of irradiation, on the nucleus surface. The first term on the right-hand side of the equation describes the energy flux conducted inwards from the nucleus surface, depending on the thermal conductivity of the nucleus, κ (W m⁻¹ K⁻¹), and the temperature gradient, |$\left(\partial T/\partial x\right)|_{0^+}=\left(\partial T/\partial x\right)|_{x=0^{+}}$|⁠. Depending on the surface temperature, T₍₀₎ = T_x = 0 (K), the nucleus cools off via thermal radiation, |$\sigma \varepsilon T^4_{(0)}$|⁠, where σ (W m⁻² K⁻⁴) is the Stefan–Boltzmann constant and ε ≈ 1 is the emissivity of the surface. The last term on the right-hand side accounts for the energy consumption by water sublimation, where Z₍₀₎ = Z(T₍₀₎) (kg m⁻² s⁻¹) is the mass flux of sublimation and ℓ (J kg⁻¹) is the latent heat of water ice.

The sublimation flux can be evaluated as

\begin{equation} Z_{(0)} = f Z^{\rm H-K}(T_{(0)}), \end{equation}

(2)

where

\begin{equation} Z^{\rm H-K} = \alpha P_V \sqrt{\frac{\hat{m}}{2{\pi} k_{\rm B} T}}, \end{equation}

(3)

is the Hertz–Knudsen formula, predicting the sublimation flux from the surface of pure solid ice into vacuum (Kossacki et al. 1999; Gundlach et al. 2011). P_V is the temperature-dependent saturation vapour pressure of water, given by

\begin{equation} P_V = a \exp \left({-\frac{b}{T}}\right)\quad ({\rm Pa}), \end{equation}

(4)

with constants a and b that can be located in such references as Fanale & Salvail (1984) and Gundlach et al. (2011). k_B in equation (3) is the Boltzmann constant, and |$\hat{m}$| is the mass of a water molecule. The factor α is the sublimation coefficient, accounting for the impinging and, thereby, sticking of gas molecules back on to the icy surface that reduces the net sublimation flux. α is evaluated as follows (Gundlach et al. 2011):

\begin{equation} \alpha = c_0 + \frac{c_1}{1+\exp (c_2-c_3/T)}, \end{equation}

(5)

where c_i for i = 0, 1, 2, 3 are constants.

Icy area fraction.

 Because the Hertz–Kundsen formula of equation (3) applies to pure-ice surface, the resultant flux must be in excess of the actual in the case of an icy nucleus with impurities of dust. For this reason, a factor of ‘icy area fraction’ is introduced in equation (2) to approximate this reduction. To illustrate this, one may consider a unit area of the nucleus surface or a horizontal cross-section of the nucleus subsurface at a certain depth. f = 1 corresponds to a fully icy surface in which case sublimation may occur everywhere. In a dust–ice mixture, a certain portion of the area is occupied by dust where no sublimation may occur. As expounded by Crifo (1997), this factor can be expressed as follows:

\begin{equation} f = \left(1+\frac{\rho _{\rm i}}{\rho _{\rm d}}\mu \right)^{-1}, \end{equation}

(6)

with ρ_i, ρ_d (kg m⁻³) indicating the densities of water ice and dust, respectively. μ is the dust-to-ice ratio of the nucleus in mass. Hence, outgassing from a dusty nucleus surface with higher μ is weaker than that from an icy nucleus. Note that f in this context measures the microscopic areal abundance of water ice at given depth, where the dust is assumed to be in thermal equilibrium with ice. This parameter is adopted in other works, such as Gutiérrez et al. (2001), Davidsson & Skorov (2002) and Groussin & Lamy (2003).

For the sake of clarity in the subsequent discussion, we refer to the above formulation of the problem where water sublimation occurs from an icy nucleus surface as the ‘dusty-ice’ model.

2.2 Mantled nucleus surface

The validity of the dusty ice model rests on the prevalent exposure of water ice on the nucleus surface. Comet 67P is known to be overall desiccated at the surface with scarce exposure of water ice (Capaccioni et al. 2015; Filacchione et al. 2016). Hence, a more realistic characterization of the nucleus is a dust–ice mixture overlain by a dry dust mantle. An alternative formulation to the dusty ice model is needed. This model is termed ‘dust mantle model’, where the sublimation of water ice is assumed to occur underneath dry mantle. The energy balance at the dry surface is formed by omitting the energy consumption by water sublimation from the budget on the right-hand side of equation (1) (Kührt & Keller 1994),

\begin{equation} Q_{(0)} = \left. -\kappa \frac{\partial T}{\partial x} \right|_{0^+} + \sigma \varepsilon T^4_{(0)}. \end{equation}

(7)

The expense should be instead explicated at an additional boundary of the ice front at x = X, i.e. (Kührt & Keller 1994),

\begin{equation} {\left. -\kappa _{\rm d} \frac{\partial T}{\partial x} \right|_{X^-} = \left. -\kappa _{\rm i} \frac{\partial T}{\partial x} \right|_{X^+} + \ell Z_{(X)},} \end{equation}

(8)

where κ_d and κ_i denote, respectively, the thermal conductivities of the dust mantle and the icy interior. Note that the mantle corresponds to x ∈ [0, X). The sublimation flux beneath a dust mantle is given by

\begin{equation} Z_{(X)} = \Psi f Z^{\rm H-K}(T_{(X)}), \end{equation}

(9)

where T_(X) denotes the temperature at the ice front. Ψ accounts for the permeability of the mantle to gas flow that may reduce outgassing with respect to a bare, icy surface. We assume that the dust mantle comprises uniform spherical dust aggregates with diameter d_D and has random packing (Skorov & Blum 2012). As determined by Gundlach et al. (2011) via diffusion experiments, the reduction or ‘quenching’ factor is given by

\begin{equation} \Psi = \frac{1}{1+p\,H}, \end{equation}

(10)

where p is a constant, and where H = X/d_D measures the ‘height’ of the mantle in diameters of constituent particles.

The icy area fraction, f, in equation (9) is as defined by equation (6). Here, f > 0 holds only for the icy interior; f = 0 applies in the dry dust mantle, otherwise.

2.3 Heat transport in nucleus

Neglecting ice sublimation and gas diffusion that arise from the porosity of the nucleus, the energy transport in the interior is described by the 1D heat equation,

\begin{equation} c \rho \frac{\partial T}{\partial t} = \frac{\partial }{\partial x}\left(\kappa \frac{\partial T}{\partial x}\right), \end{equation}

(11)

where c (J K⁻¹) and ρ are the specific heat capacity and the mass density of the material, respectively. The thermal inertia of the material is specified via,

\begin{equation} \mathbb {I} = \sqrt{\kappa c \rho }\quad ({\rm W}\,{\rm K}^{-1}{\rm m}^{-2}{\rm s}^{1/2}). \end{equation}

(12)

Depending on the time-scale or periodicity of the energy input, the heat flux becomes negligible from a certain depth, i.e.

\begin{equation} {\left.-\kappa \frac{\partial T}{\partial x} \right|_\mathbb {X}= 0, } \end{equation}

(13)

or equivalently, |$T_{(x)} = T_{(\mathbb {X})}$| for |$x \ge \mathbb {X}$|⁠, such that the nucleus becomes nearly isothermal from depth |$\mathbb {X}$|⁠.

The adiabatic condition of equation (13), together with either equation (1) or the set of equations (7) and (8) distinguishing the structures of the nucleus surface layers, constitute the boundary conditions for solving equation (11). The dust-mantle model requires one more boundary condition at the bottom of the mantle or ice front.

2.4 Note on literature

The dusty ice model has been frequently adopted in studies of cometary activity and evolution, though numerical treatment varies greatly with different works. References can be found in Smoluchowski (1981), Weissman & Kieffer (1981), Froeschle, Klinger & Rickman (1983), Kührt (1984), Kührt (1999), among numerous others.

It is conclusive that exposure of water ice over cometary nuclei is limited (Filacchione et al. 2016). The effect of this dust mantle on heat transport and gas diffusion through the nucleus as well as the conditions of dust activity has been studied in great detail. On this matter, the reader is referred to Mendis & Brin (1977), Brin & Mendis (1979), Brin (1980), Fanale & Salvail (1984), Kömle & Steiner (1992), Kömle et al. (1992), Kührt & Keller (1994), Skorov & Rickman (1995), Davidsson & Skorov (2002), Davidsson & Gutiérrez (2005), Kossacki & Szutowicz (2008), Skorov et al. (2011), Gundlach et al. (2011), Gundlach & Blum (2012) and Keller et al. (2015), if only for just a partial review of the immense literature.

What some might consider to be a drastic simplification with both the dusty-ice and dust-mantle models is the implicit assumption that the sublimation occurs from a surface of pure solid ice, in contradiction with cometary nuclei being porous objects in reality. Porosity of the nucleus, and in particular, of the icy interior gives rise to complex behaviour of volatile transport inside the nucleus. Most notably, the gas diffusion does not occur unilaterally; while a certain fraction of the gas will seep upwards through the dust mantle and escape, some gas will also likely diffuse inward and re-condense upon cooling. This phenomenon requires and, indeed, was first revealed by substantially more sophisticated thermal modelling (see e.g. Prialnik & Bar-Nun 1988; Spohn & Benkhoff 1990; Mekler, Prialnik & Podolak 1990; Prialnik & Mekler 1991; Yabushita 1995; Enzian, Cabot & Klinger 1997; De Sanctis et al. 1999; Capria et al. 2000; Gortsas et al. 2011; González, Gutiérrez & Lara 2014). These models account explicitly for not only heat transport but also mass transfer that is governed by the sublimation and re-condensation of ices in the nucleus interior. The application of these more comprehensive models is beyond the scope of this work, but clearly merits further investigations in the near future.

3 APPLICATION OF THERMAL MODEL TO COMET 67P

3.1 Evaluation of energy input

Insolation.

 The primary source of energy input to the nucleus is insolation, given by

\begin{equation} {\rm Q}_{\odot } = \left(\frac{r_{\odot }}{1\,{\rm au}}\right)^{-2} {\rm C}_{\odot } (1-\mathcal {A})\, \delta (\alpha _{\odot },\varphi _{\odot }) \sin \varphi _{\odot } , \end{equation}

(14)

where C_⊙ = 1361 W m⁻² is the solar constant; |$\mathcal {A}$| is the bolometric bond albedo of the nucleus surface; r_⊙ is the heliocentric distance of the nucleus; α_⊙, φ_⊙ (degree) denote, respectively the azimuth and altitude of the Sun measured with respect to the local topocentric coordinate system at the given location. Variable and, in particular, concave topography may give rise to obstruction of sunlight, or shadowing, even if φ_⊙ > 0. Hence, it is defined that δ = 1 if the local surface is illuminated, and δ = 0 otherwise. The strategy for the determination of δ is discussed in Section 3.2.

Self-heating.

 On the other hand, topographic variations tend to impede cooling of the nucleus via enhanced absorption of thermal radiation from the surroundings, a phenomenon known as ‘self-heating’ (Colwell & Jakosky 1987; Colwell et al. 1990; Lagerros 1997; Gutiérrez et al. 2001; Ivanova & Shulman 2006; Davidsson & Rickman 2014). The effect is pronounced in the rugged landform, such as valleys, coves, cliffs, compensating for a deficiency of insolation due to shadowing. We assume that the nucleus radiates isotropically. Suppose a surface area of A_F (m²) around a point F on the nucleus subtends a small solid angle, dΩ ≈ A_Fcos θ/d², from another surface point, F^΄, where θ is the emission angle of the thermal radiation at F towards F^΄ and d is the distance in between. It can be shown that the thermal irradiation at F^΄ by A_F with surface temperature T₀ is given by (Davidsson & Rickman 2014),

\begin{equation} {\rm d}Q_{\rm F} \approx \sigma \varepsilon T^4_{(0)} \cdot \frac{\cos \theta ^{\prime } \cos \theta }{{\pi} d^2}A_{\rm F} , \end{equation}

(15)

where θ^΄ is the incidence angle of irradiation at F^΄. The total energy flux of self-heating at F^΄ is then evaluated as

\begin{equation} {\rm Q}_\otimes = \int _{{\rm F} \in \otimes , {\rm F} \ne {\rm F}^{\prime }} {\rm d}Q_{\rm F}. \end{equation}

(16)

It follows that the total energy input on the left-hand side of equation (1) and (7) is

\begin{eqnarray} Q_{(0)} = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}{\rm Q}_{\odot }, \quad {\rm \quad without\, self{-}heating;}\nonumber\\ {\rm Q}_{\odot } + {\rm Q}_\otimes ,\quad {\rm with\, self{-}heating}. \end{array}\right.\\ \end{eqnarray}

(17)

There are other sources of energy input. The effects of multiscattering, such as reflection of sunlight and re-absorption of nucleus thermal radiation, are neglected in this work. Arguably, this simplification should be reasonable for low-albedo objects. We note that these complex phenomena have been taken into account in various thermal analyses (see e.g. Gutiérrez et al. 2001; Davidsson & Rickman 2014; Keller et al. 2015).

3.2 Landscape data base

The evaluation of Q_⊙ demands the knowledge of r_⊙, α_⊙ and φ_⊙ at the epochs of interest. r_⊙ can be obtained directly from the SPICE kernels for Rosetta (Acton 1996). α_⊙, φ_⊙ need to be defined and calculated locally at given position on the nucleus. The determination of Q_⊗ requires totalling the thermal radiation from all surface areas visible to the local observer. Therefore, the provision of such information on the ‘visibility’ of the surrounding topography can facilitate the determination of shadowing, as well.

In order to account adequately and efficiently for the impact of the irregular shape on the energy input of the nucleus, an assortment of records depicting the landscape of the nucleus is organized. We call the collection ‘Landscape’ data base for 67P nucleus. The concept has been implemented by Lagerros (1997) and, in the case of 67P, by Keller et al. (2015). The elements of the data base are enumerated and explained in this section.

3.2.1 Shape model supplementary

Shape models of 67P have been developed by Preusker et al. (2015) and Jorda et al. (2016), represented by a polyhedron whose surface consists of a number of triangular facets. Let us denote the vectors of the respective three vertices as |$\boldsymbol {v}_k$| with k = 1, 2, 3, ordered counterclockwise with respect to an external observer. The surface normal is defined by the cross product, |$\boldsymbol {n}$| = |$(\boldsymbol {v}_2-\boldsymbol {v}_1) \times (\boldsymbol {v}_3-\boldsymbol {v}_1)$|⁠, and always directed outwards.

It is advantageous to augment the original form of the shape model with some supplementary materials that will greatly facilitate the thermal analysis.

Centroid of facet.

The position vector for a given facet is defined as pointing to the centroid of the facet, i.e.

\begin{equation} \boldsymbol {r}_{\rm F} = \frac{1}{3} \left( \boldsymbol {v}_1 + \boldsymbol {v}_2 + \boldsymbol {v}_3 \right), \end{equation}

(18)

which uniquely locates the position of a facet, irrespective of its orientation in space. As will be discussed, |$\boldsymbol {r}_{\rm F}$| is the defined origin for a local coordinate system associated with each facet.

Area of facet.

The area of the facet is useful for assessing the thermal radiation from the nucleus. Depending on the topography or, more exactly, concavity, thermal radiation can be re-absorbed by the nucleus at other locations. The facet area is computed by the Heron’s formula,³

\begin{equation} A_{\rm F} = \sqrt{s(s-a)(s-b)(s-c)}, \end{equation}

(19)

where a, b and c denote the lengths of the edges, and where |$s = \frac{1}{2}(a + b + c)$|⁠.

3.2.2 Visibility

The Landscape data base provides the mutual visibility between two locations on the surface, indicated by the centroids of the respective encompassing facets of the shape model.

It is straightfroward to introduce the ‘mutual-visibility factor’, v_{k, l}, for a pair of facets k and l (Lagerros 1997),

\begin{eqnarray} v_{k,l} = v_{l,k} = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\, 1,& {\rm if \,facets}\,k\, {\rm and}\,l\, {\rm are\, mutually\, visible},\nonumber\\ \, 0,& {\rm otherwise}. \end{array}\right.\\ \end{eqnarray}

(20)

The information is visualized in Fig. 1. The fraction of the nucleus surface visible from a given location varies notably across the nucleus. Concavities, such as the ‘neck’ between the two lobes, result in large visible surface area (Figs 1a–c), in contrast to the convex topographies, e.g. on the far sides of the lobes, that are less obstructed by the surrounding landscape (Fig. 1d). Additionally, it is worth noting that the distribution of the visible facets may be scattered, which likely results from the irregularity in topography reflected even in a low-resolution shape model (Figs 1a and c).

Figure 1.

Visible surface portions from four respective locations on the nucleus, as inferred from the data base of Landscape: Visibility. The shape of the nucleus is represented by an SPC-shape model with one thousand facets (Jorda et al. 2016). The local facet is marked in red; the visible facets are in green. (a) The local facet is at the floor of a valley between the two lobes, with visible facets distributed across the walls. (b,c) The local facet is on the wall on the big lobe looking into the valley, but loses sight of facets on the other side almost entirely, where topography drops below the horizon. (d) The facet situated on the far side of the small lobe with nearly convex topography sees no other facets.

Some auxiliary information can be derived from the above analysis with little extra cost, as detailed below.

Distance between facets.

The distance between mutually visible facets is required for evaluating self-heating via equation (15). Therefore, for each pair of visible facets, k and l, the distance, |$d_{k,l}=\left| \boldsymbol {r}_{{\rm F}\, k} - \boldsymbol {r}_{{\rm F}\,l} \right|$|⁠, should be stored. The distances between invisible pairs are discarded.

Mutual view of facets.

Let θ_{k, l} denote the orientation of (the line-of-sight to) facet k with respect to the surface normal to local facet l (Fig. 2), such that

\begin{equation} \cos \theta _{k,l} = \frac{\hat{{\boldsymbol n}}_l \cdot \left( {{\boldsymbol r}}_{{\rm F}\,k} - {{\boldsymbol r}}_{{\rm F}\,l} \right) }{ d_{k,l} }, \end{equation}

(21)

where θ_{k, l} ≠ θ_{l, k} in general. A particularly useful quantity is

\begin{equation} w_{k,l} = v_{k,l} \frac{ \cos \theta _{k,l} \cos \theta _{l,k} }{\pi d^2_{k,l}}, \end{equation}

(22)

defined here as the ‘mutual-view factor’ considering the obvious symmetry, i.e. w_{k, l} = w_{l, k}. The preceding mutual-visibility factor, v_{k, l}, on the right-hand side suggests that only values for those visible pairs need to be stored. w_{k, l} differs subtly from the well-defined ‘view factor’ (see Lagerros 1997, and reference therein), and is related to the latter via F_{l, k} = w_{k, l} A_F k; the symmetry in this case is given by A_F lF_{l, k} = A_F kF_{k, l}. Note again that the use of w_{k, l} in equation (15) is only strictly valid if A_F k ≪ d². This condition may not be fulfilled for nearby, and in particular, adjacent facets. Violation of it may lead to erroneously large F_{l, k}. The remedy is to refine F_{l, k} via subdivision of facet k (Davidsson & Rickman 2014). This procedure is not applied in this work.

Figure 2.

Relative orientation of two facets. The incidence angle, θ, is measured between the line-of-sight towards the distant facet and the surface normal of the local facet. θ is complementary to the altitude φ. |$\boldsymbol {n}$| indicates the surface normal.

3.2.3 Horizon

Another element of the Landscape data base is the record of local horizons over the nucleus. The horizon on a convex object, such as a sphere, is defined conveniently by cos θ = 0 in every direction. The local horizon on 67P nucleus is highly irregular, making it difficult but necessary to delineate the skyline in order to account for shadowing effects (Gutiérrez et al. 2000). Trivially, determination of shadowing demands finding the intersection of a ray of sunlight emitted towards the given location with topography in the foreground. However, this approach requires the search for intersections of a line and a polyhedral shape model at all locations on the nucleus and at each epoch of analysis. This computation can be exceedingly demanding when high temporal and spatial resolutions are desired.

A more cost-effective strategy is to trade memory for computational speed, i.e. by pre-storing the variation of the local horizon as a function of azimuth (Lagerros 1997). These data specify altitudinal intervals (i.e. along meridian) at each azimuth in which no obstruction occurs so that illumination is possible. Thus, the repetition of laborious search for intersections of a line and a polyhedron is eliminated at the expense of memory space for data storage.

For each facet of the shape model, we define a local horizontal coordinate system, whose origin coincides with the centroid of the facet, |$\boldsymbol {r}_{\rm F}$|⁠. The z-axis is aligned with the local surface normal and points along the zenith; the direction of the local meridian, i.e. the x-axis, is arbitrary within the horizontal plane (see the Appendix). The azimuth and altitude of an object in the corresponding spherical coordinate system are given by |$\alpha = \arctan (y/x)$| and |$\varphi =\arcsin (z/r)$|⁠, respectively, where r² = x² + y² + z².

A survey is performed over the topocentric hemisphere on each facet. The hemisphere is discretized and binned in azimuth and altitude, such as, α_n and φ_m for non-negative integers n and m. The visual obstruction of topography is indicated by any intersections of a ray from the facet with the shape model in the direction specified by α_n and φ_m.

The above analysis results in a spherical point grid that samples the visual (non-)obstruction by topography in all directions for every facet, namely

\begin{eqnarray} \delta (\alpha _n,\varphi _m) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\, 1,\quad {\rm if\, open\, sky\, is\, along\, line\, of\, sight},\nonumber\\ \, 0,\quad {\rm if\, topography\, is\, along\, line\, of\, sight}. \end{array}\right. \\ \end{eqnarray}

(23)

The local surface within the facet can be illuminated by the Sun from (α_n, φ_m) if δ = 1. These data are stored as intervals of altitude, i.e. (φ_a, φ_b), such that δ(α_n, φ_a < φ < φ_b) = 1 for all α_n at a given facet. Note that there might be multiple intervals of (φ_a, φ_b) at a certain α_n.

The horizons at four locations indicated in Fig. 1 are depicted in Fig. 3. The panoramic view is projected over the upper hemisphere, discretized into a spherical point grid at 5-deg interval in altitude and 4-deg interval in azimuth.⁴ The black dots indicate nucleus surface along the line of sight, whereas the unfilled areas indicate the open sky. The complex topography distorts the local horizon from a ‘horizontal’ line. In the cases of Figs 3(a)–(c), the view of the facets is obstructed by the two lobes and confined along the valley in between. In Fig. 3(b), the facet is mostly ceiled by the wall over the small lobe, with a view similar to the outlook from a cave. Open view is not uncommon, e.g. on the far side of the small lobe where global-scale topographic variation is absent, as in Fig. 3(d). It is evident that the shadowing effect needs to be taken into account for at least a significant portion of the nucleus surface.

Figure 3.

Shape of local horizon in the local horizontal coordinate system of azimuth and altitude, as inferred from the Landscape: Horizon data base. (a) The horizon is distinctly rugged; the open and obstructed views alternate more than once along the same meridian, e.g. at around 30 deg. (b) Cave view with visual obstruction over zenith. (c) Another case of multiple alternations of open and obstructed views along the meridian around 180 deg. (d) Open view on a convex topography. The location of the observer is within the red triangle on the shape model in the left-hand panel. The x-, y-, and z-axes of the local horizontal coordinate system for the indicated facets (red triangles) are in blue, green and red, respectively.

3.3 Numerical procedure

The 67P nucleus is subject to periodic heating of insolation and self-heating due to both rotation and orbital motion of the nucleus around the Sun. The periodicity of the energy input is expressible by Q₍₀₎(t) = Q₍₀₎(t ± t_P), with t_P being the period of variation of insolation. It is then imperative to solve for periodic variations of temperatures and outgassing flux of the nucleus, e.g. T(t) = T(t ± t_P), according to the physical time-scale of the problem at hand.

The thermal skin depth of the nucleus can be defined by Huebner et al. (2006),

\begin{equation} X_\mathcal {S} = \sqrt{ \frac{t_{\rm P} \kappa }{{\pi} \rho c} }, \end{equation}

(24)

which approximates the depth at which the temperature variation decays by a factor of e⁻¹ with respect to the surface temperature in response to periodic heating. With a rotation period of 12.4 h (Preusker et al. 2015), the diurnal variation of water activity is considered to be mostly affected by the top few centimetres of the subsurface. On the other hand, the temperatures and heat transport through the first few metres of the subsurface need to be resolved in order to describe the orbital behaviour of water activity on 67P over about 6.5 yr.

3.3.1 Crank–Nicolson method

We solve equation (11) via the Crank–Nicolson method. The nucleus subsurface down to the isothermal depth is discretized as x_ȷ = x₁ + (ȷ − 1)Δx, where ȷ = 1, 2, ⋅⋅⋅, ȷ_max index the discrete layers and Δx denotes depth increment. The numerical treatment for the dust-mantle model differs from that for the dusty-ice model, in which the temperature profiles in the mantle and in the icy interior are solved in separate sets for the former. For instance, the energy balance at the lower boundary for the solution of mantle coincides with the upper boundary condition for interior (equation 8). We may define the dust mantle as occupying the top |$\mathcal {J}$| layers.⁵ Let a set of Crank–Nicolson solution be denoted as {u_j, j = 1, 2, ⋅⋅⋅, j_max}. In the dust-mantle model, we obtain two solutions, i.e. {u_j, j = ȷ}_mantle for temperatures in the mantle and |$\left\lbrace u_j,\,j=\jmath -\mathcal {J}\right\rbrace _{\rm interior}$| for those in the interior. The two sets are concatenated to yield a full temperature profile as {T_ȷ}. On the other hand, temperatures in a uniform nucleus are given by one set, i.e. {u_j, j = ȷ}, in the dusty-ice model. If the adiabatic condition is imposed at the bottom boundary as given by equation (13), it must be ensured that the maximum depth of the numerical solution always exceeds the thermal skin depth.

The temperatures at the boundaries, i.e. T₍₀₎ at the nucleus surface and, if applicable, T_(X) at the bottom of the mantle, are not output directly by the Crank–Nicolson solution; instead, they are solved iteratively via the Newton–Raphson method once the Crank–Nicolson solutions are completed at each time-step.

3.3.2 Local and global iterations of solution

To initialize the solution, we assume that the nucleus is isothermal and let |$T_\jmath (t_0) = T_{(\mathbb {X})}$| at the starting epoch, t₀. As noted, we seek a solution where the nucleus temperatures are equilibrated and repeat with the periodic variation of energy input. The temperatures and, depending on the depth of the ice front, the sublimation fluxes are propagated or solved at every step in time, Δt, until a full period, t_P, is reached. The deviation between the temperature profiles one period apart is evaluated as,

\begin{equation} D_T = \sqrt{ \frac{1}{\jmath _{\rm max}} \displaystyle \sum ^{\jmath _{\rm max}}_{\jmath =1} \left[T_\jmath (t_0) - T_\jmath (t_0+t_{\rm P}) \right]^2 }. \end{equation}

(25)

The solution needs to be iterated, for example, re-initialized by |$T^{(i+1)}_\jmath (t_0) = T^{(i)}_\jmath (t_0+t_\mathrm{P})$|⁠, where i counts the number of iterations performed. The solution is considered convergent once D_T < ε is fulfilled for some small threshold value of ε of choice. The amplitude of the temperature variation at depth x_ȷ is given by |ΔT_ȷ| = max [T_ȷ(t)] − min [T_ȷ(t)] for t ∈ [t₀, t₀ + t_P]. The ‘numerical’ thermal skin depth is then found as the first layer below the surface where |ΔT_ȷ| ≈ 0.

The iterative procedure described above applies to the model solution with insolation as sole energy input. Hence, it can be performed for each facet independently. We call this procedure ‘local iteration’. In case self-heating is included, the thermal irradiation depends on the instantaneous surface temperature across the nucleus. Consequently, another ‘global iteration’ is necessitated following the local iterations, where Q_⊗ in equation (17) is evaluated via equations (16) and (15) and updated in each iteration. The criterion of convergence for the final solution is the same as equation (25) and not elaborated here.

4 MODEL RESULTS AND INTERPRETATION

4.1 Model parametrization and output

We present and discuss some basic simulations of nucleus temperatures and water outgassing flux of 67P via the dusty-ice and dust-mantle thermal models. The model parameters are summarized in Table 1. We first solve for the orbital temperature fluctuations at various locations of the nucleus, assumed to be initially isothermal at T = 15 K at aphelion on the present orbit of 67P. Orbital solutions indicate that temperature varies marginally around 100 K at the depth of about 1 m. Hence, we adopt |$\mathbb {X}=1\ {\rm m}$| as the isothermal depth for the diurnal solution. Because |$\mathbb {X}$| far exceeds the diurnal skin depth (as approximated by equation 24), the diurnal solution is insensitive to the choice of |$T_{(\mathbb {X})}$| for the temperature at the bottom boundary.

Table 1.

Parameters of thermal models for diurnal solutions.

Parameters	Symbols	Values
Step in depth	Δx	1	mm
Step in time	Δt	tP/1200
Rotation period	tP	44650	s
Bond albedo	|$\mathcal {A}$|	0.05
Emissivity	ε	1
Heat conductivity	κ	2 × 10−3	W m−1 K−1
Specific heat capacity	c	1000	J kg−1 K−1
Density	ρ	500	kg m−3
Diameter of dust particlea	dD	1	mm
Latent heat of water ice	ℓ	2.3 × 106	J kg−1
Saturation vapour pressureb	a	3.23 × 1012	Pa
	b	6134.6	K
Sublimation coefficientc	c0	0.146
	c1	0.854
	c2	57.78
	c3	11580	K
Quenching factord	p	0.14
Bottom depthe	|$\mathbb {X}$|	1	m
Interior temperaturee	|$T_{(\mathbb {X})}$|	100	K

Notes.^aSee Section 5.3 in Hu et al. (2017) for a justification for adopting this value.

^bEquation (4) (Gundlach et al. 2011).

^cEquation (5) (Gundlach et al. 2011).

^dEquation (10) (Gundlach et al. 2011).

^eObtained by the orbital solution.

The temperature profiles at two locations, hereafter designated as A and B, on 67P nucleus at UTC 00:00:00 on 2014 June 15, are shown in Fig. 4. The shape model of 1000 facets is used (see Fig. 1 or 3). Unless specifically noted otherwise, self-heating is neglected hereafter.⁶ For the sake of brevity, the epoch is referred to as t₀. With a low thermal inertia of ∼30 W K⁻¹ m⁻² s^1/2 in the nucleus subsurface (Gulkis et al. 2015), only temperatures in the topmost centimetres of the nucleus affect diurnal variation of the water activity. The coincidence of the profiles at t₀ and after one rotation at t₀ + t_P suggests that convergence of the numerical solution is achieved. Location A in the neck region (Fig. 4a) was illuminated at t₀, as reflected by steep negative temperature gradient in the top centimetres below the surface. The profile starts to flatten at the depth of around 1.5 cm marking a pivot point on the curve. Correspondingly, the temperature profile at B in shadow at t₀ shows a crest at roughly the same depth (Fig. 4b). It can be inferred that diurnal variations of temperatures nearly vanish below this depth of insulation; a steady heat flux is directed towards the nucleus interior underneath. The diurnal insulation depth of a mantled nucleus is greater than that of a bare nucleus. The presence of the dust mantle results in a notable increase of temperature in the interior.

Figure 4.

Profiles of nucleus temperatures at two locations, designated as A and B, on 67P as derived via the thermal models. The profiles for the nucleus with a 5 mm dust mantle at t₀, UTC 00:00:00 on 2014 June 15, are indicated by the solid red curves; the dashed black curves, nearly coincident with the red curves, indicate the temperatures after one rotation period, t₀ + t_P. The blue curves indicate temperatures in a uniform nucleus at t₀, derived via the dusty ice model for comparison. The locations are marked by a red triangle on the shape model in the respective panels. The pattern of instantaneous insolation according to equation (14) is in grey-scale.

4.2 Thermal skin depth

The surface temperature of the dust-mantled nucleus is higher by about 10 K than that for an icy surface from the dayside, but lower from the nightside (Fig. 4). As shown in Fig. 5, the surface temperature on the dust-mantled nucleus fluctuates more strongly than that on the icy nucleus surface. During the daytime, more energy is used to warm up the nucleus when the sublimation rate is reduced by the dust mantle; conversely, at night, the ice front beneath the dust mantle remains warmer and, thus, maintains stronger sublimation than a bare icy surface, in which case less energy from the interior can escape from the surface. This also explains why diurnal heat waves penetrate to larger depths in the mantled nucleus (Fig. 4).

Figure 5.

Diurnal variation of temperatures at different depths. The results in (a) and (b) are derived via the dusty-ice and dust-mantle models, respectively, for location A as indicated in Fig. 4(a). The results in (c) and (d) are for location B as indicated in Fig. 4(b). The mantle thickness is 5 mm.

The temperature variation attenuates with depth. According to the dust mantle model, for instance, the temperature at the ice front of x = 5 mm varies by roughly 1/2 of the range at the surface. At the depth of 10 mm, the temperature range attenuates to 1/4. Therefore, the diurnal skin depth that denotes the decay by e⁻¹ lies between 5 and 10 mm. This conclusion holds for the bare icy nucleus, as well. The numerical result is consistent with the analytic approximation of about 8 mm by equation (24). Note that the diurnal skin depth is always smaller than the insulation depth (about 1.5 cm, Fig. 4).

4.3 Energy consumption

The conservation of energy can be expressed as follows:

\begin{equation} Q_{(0)} = q_{\varepsilon \,(0)} + q_{Z} + \dot{U}, \end{equation}

(26)

where Q₍₀₎ is the energy input (recalling equation 17) and |$q_{\varepsilon \,(0)}=\varepsilon \sigma T^4_{(0)}$| is introduced as a symbol for thermal radiation. q_Z = ℓZ stands for energy flux consumed by sublimation of water ice. U is, loosely speaking, the internal energy of the system, which varies with temperature, i.e. ΔU = cρ · ΔT (J m⁻³) for a unit volume of material. Over a unit area at a certain location on the nucleus, the change of internal energy in the subsurface is expressed by

\begin{equation} \Delta U = \int _0^{\mathbb {X}} c\rho \Delta T {\rm d}x\quad ({\rm J}\,{\rm m}^{-2}), \end{equation}

(27)

where the integration is carried out formally from the surface down to the isothermal depth, |$\mathbb {X}$|⁠, where the heat flux vanishes, according to equation (13). The density and specific heat capacity are assumed to be constant as a simplification in this work. The rate of change of the internal energy is evaluated numerically by summing up contributions from discrete layers, i.e.

\begin{equation} \dot{U} = \displaystyle \sum _{\jmath =1}^{\jmath _{\rm max}} c \rho \dot{T}_\jmath \Delta x_\jmath , \end{equation}

(28)

where |$\dot{T}_\jmath$| denotes the rate of temperature variation of the ȷth layer with thickness Δx_ȷ. In the case of uniform discretization, Δx_ȷ is a constant.

While tantamount to equation (1) or equation (7) indicating the surface energy balance for the respective thermal models, equation (26) explicates the conservation of energy and should not be considered a mere repetition; it can also be used to check numerical consistency.

The diurnal variations of energy fluxes in the nucleus at the locations A and B on 2014 June 15 are shown in Fig. 6. A distinction between the results of two models lies with the energy consumption of sublimation on the dayside. The sublimation from an icy nucleus surface consumes up to about 15 W m⁻² (light blue curves in Figs 6a and c), whereas the consumption is negligible on the nucleus with a dust mantle (Figs 6b and d). Accordingly, the mantled nucleus cools off via stronger thermal radiation compared with a bare nucleus on the dayside. The energy consumption in a mantled and bare nuclei is similar during the night (where the red curves trace zero in Fig. 6).

Figure 6.

Energy fluxes inside the nucleus over one comet rotation on 2014 June 15. Results are for locations A and B on the nucleus as indicated in Fig. 4. Energy input of insolation is indicated by the red curve. Thermal radiation is indicated by the dark blue curve. The energy consumed by sublimation corresponds to the light blue curve. The change rate of nucleus internal energy is given by the black curve.

4.3.1 Global energy consumption and variation with heliocentric distance

The total outgassing flux of water from the nucleus is obtained by integrating the contribution from all facets of the shape model, such as,

\begin{equation} \Sigma Z = \displaystyle \iint _A Z \cdot \mathrm{d}A = \displaystyle \sum _k Z_k \cdot A_{{\rm F}\,k}, \end{equation}

(29)

where Z_k indicates the sublimation flux from the kth facet of the shape model with surface area A_F k. The notation of ‘Σ’ applies to other quantities, e.g. Σq_Z = ∑_kq_Z k · A_F k, denotes the total expense of energy for water ice sublimation from the nucleus.

Fig. 7 illustrates the global energy input and consumption of the nucleus on 2014 June 15, 2014 December 15 and 2015 June 15, respectively, when 67P was at the corresponding heliocentric distances of about 3.9, 2.8 and 1.4 au. The 1000-facet shape model is used. Consistent with the results for single locations (Fig. 6), the energy dissipation is dominated by thermal radiation at 3.9 au, regardless of the existence of the dust mantle (dark blue curves in Fig. 7a). At 2.8 au, the water outgassing becomes a significant source of energy loss from a bare icy nucleus (light blue curve in the left-hand panel of Fig. 7b). Close to perihelion at 1.4 au, the bare nucleus cools off predominantly via water sublimation (light blue curve in the left-hand panel of Fig. 7c). Water outgassing accounts for little loss of energy from the dust-mantled nucleus that dissipates most of the heat via thermal radiation at all heliocentric distances (right panels of Figs 7a–c). In all cases, the dayside temperature, and hence, thermal radiation from the mantled nucleus is higher than from an icy surface.

Figure 7.

Diurnal variations of energy fluxes integrated over the nucleus surface on 2014 June 15 (a), 2014 December 15 (b) and 2015 June 15 (c). In each case, results in the left-hand panel are for the bare icy nucleus, and those on the right-hand panel are for the mantled nucleus. The mantle thickness is 5 mm. The red curves indicate global energy input of insolation. The black curves indicate variation of internal energy of the nucleus. The dark blue curves correspond to thermal radiation. The light blue curves are energy fluxes consumed by water sublimation.

5 INFLUENCE OF MANTLE THICKNESS, DUST-TO-ICE RATIO AND SELF-HEATING ON GLOBAL WATER PRODUCTION

5.1 Dust mantle thickness

The thickness of a dust mantle significantly influences the water activity over the nucleus. Simply put, the dust mantle reduces the sublimation flux with respect to free sublimation of ice on the surface; in addition, it overlies the ice front and moderates the temperature of sublimation.

Fig. 8 illustrates the total water outgassing flux of the nucleus covered by a uniform dust mantle of varying thicknesses over one cometary rotation. A fully icy (sub)surface is assumed (f = 1). The 1000-facet shape model is used. Water sublimation from below a dust mantle of 2 mm is less than half of that from an icy nucleus surface; an even thicker mantle of 5 mm restricts the outgassing flux to less than 20 per cent of the surface flux. The range of diurnal variation attenuates with increasing mantle thickness. For instance, the outgassing flux varies by less than 5 kg s⁻¹ under a 5-mm-thick mantle, in contrast to the variation of about 20 kg s⁻¹ in the case of ice sublimation from the surface.

Figure 8.

Modeled water outgassing flux of 67P over one nucleus rotation on 2014 June 15. The dark blue curve indicates the water production over an icy nucleus without a dust mantle. The medium and light blue curves correspond to the dust mantle thickness of 2 and 5 mm, respectively.

More insight is gained into the specific role of the dust mantle in inhibiting the sublimation flux by comparing the quenching factor given by equation (10), which arises solely from the resistance of the porous mantle to gas diffusion, with the numerical results. For instance, the sublimation flux below a 2-mm dust mantle is reduced to 78 per cent of that from a bare nucleus at the same temperatures; the quenching factor is 59 per cent for 5-mm dust mantle. The stronger decrease of sublimation fluxes with mantle thickness than predicted by equation (10), as shown in Fig. 8, suggests that the reduction of sublimation results from the lower temperature below the dust mantle at which sublimation occurs relative to the surface temperature of a bare icy nucleus.

It would be misleading to conclude that the presence of the dust mantle always results in lower temperature at the ice front than that on a bare icy nucleus (on the day side). In addition to heat conductivity, the modulation of ice temperature depends critically on the permeability of the dust mantle (Kömle et al. 1992): An impermeable mantle, one that consists of finer, micron-sized dust grains, for example, would quench more effectively water sublimation underneath (see equation 10) and, thus, could induce ice temperature higher than that on a bare nucleus. This phenomenon, known as ‘pressure cooking effect’, has been substantiated both experimentally and numerically (Kömle & Steiner 1992; Kömle et al. 1992). In our simulations, the dust mantle comprises millimetre-sized particles and the pressure cooking effect does not play a notable role. Instead, the low heat conductivity of the dust mantle lowers the temperature at the ice front compared with that on an icy nucleus surface (Kührt & Keller 1994).

The maximum water production rate below the dust mantle of 2-mm thickness lags behind the peak production of the icy nucleus by about 20 min (Fig. 8), whereas the lag in the case of a 5-mm-thick mantle increases to over one and half hours. This delay reflects the time-scale of thermal conduction to the ice front, which can be approximated by a simple formula (Huebner et al. 2006):

\begin{equation} \tau \sim (x\mathbb {I}/\kappa )^2\quad ({\rm s}). \end{equation}

(30)

Therefore, the outgassing flux from beneath a thicker dust mantle would peak even later.

5.2 Dust-to-ice ratio

The icy area fraction of the nucleus is governed by the dust-to-ice ratio (equation 6). The presence of non-volatile components decreases the portion of the (sub-)surface area where sublimation may occur, leading to reduced sublimation flux with respect to a fully icy nucleus (sub-)surface.

However, icy area fraction may not be a direct measure of the moderation of outgassing in the presence of non-volatile impurities. That is, the sublimation flux from a nucleus (sub-)surface with icy fraction, f, is not in proportion to that from a pure-ice nucleus (with f = 1). In Fig. 9, the total outgassing flux of 67P nucleus (approximated by the 1000-facet shape model) with respective icy area fractions of 100 per cent, 10 per cent, and 1 per cent are derived via the dusty ice thermal model. f = 0.1 results in a reduction to about 20 per cent of the production from a fully icy surface. The outgassing flux in the case of f = 0.01 is around 3 per cent of that from a fully icy surface.

Figure 9.

Diurnal variation of total outgassing flux from 67P influenced by the icy area fraction of the nucleus. Three values are considered, i.e. 100 per cent, 10 per cent and 1 per cent. The dusty ice thermal model is applied for analysis.

The under-reduction of the outgassing with respect to f is due to the compensation by increasing temperature of ice sublimation. Referring to the energy balance given by equation (1) (or equation 8 for the dust mantle model), one sees that a reduction by f would result in a decrease of energy consumption by ice sublimation. Hence, more heat will be transported into the interiors, thus raising the temperatures at ice front. As a result of such negative feedback, sublimation will be partly compensated. While warming of the nucleus will counter the reduction, it does not overcompensate sublimation flux such that the net result is still a decrease in production.

We note that only in the case where the nucleus is locally active, i.e. with substantial dormant areas, is it justified to directly (but still arbitrarily) scale down the total outgassing flux with the icy fraction of f, which is, in effect, the macroscopic measure of the active fraction of the total nucleus surface. In this study, the nucleus is assumed to be globally active, in which case f is a microscopic quantity regulating water outgassing from a homogeneous surface area.

5.3 Self-heating

The nucleus of 67P is among the most irregular-shaped Solar system objects ever observed. The globally bi-lobed shape and locally abrupt topography need to be carefully considered for evaluation of energy input of insolation and thermal irradiation of the nucleus (equation 17) (Sierks et al. 2015). To demonstrate the impact of self-heating, we apply a higher-resolution shape model consisting of 15 000 facets (Jorda et al. 2016), such that the size of an individual facet is negligible with respect to the dimension of the nucleus. The valley or ‘neck’ region on the nucleus is overall poorly illuminated due to shadowing of the lobes; however, it is where the effect of self-heating is most probable. Self-heating accounts for up to 20 per cent of additional energy input relative to the maximum insolation over the concavity (Fig. 10a). The maximum increase coincides with abrupt topographic variations, e.g. tracing (accentuating) the edges of scarps and recesses.

Figure 10.

Effect of self-heating over 67P nucleus at t₀ (2014 June 15). (a) Increase of energy input. (b) Increase of surface temperature. (c) Enhancement of sublimation flux. The results are derived via the dusty ice model. In each case (row), the left-hand panel shows the result without self-heating; the middle panel shows the result with self-heating; the right-hand panel shows the net increase due to self-heating. Note that a narrower colour-scale needs to be used to exhibit the pattern of net change, always smaller in magnitude than the total in the right-hand panel in all cases.

The valley can be warmed by nearly 40 K over the shadowed areas, absorbing thermal radiation from the illuminated areas across the concavity (Fig. 10b). The pattern of temperature increase is anticorrelated with distribution of the temperature itself and, naturally, pattern of instant illumination (comparing left-hand and right-hand panels in Fig. 10b). Hence, the impact of self-heating is more pronounced over the shadowed areas, especially those adjacent and exposed to the well-illuminated areas and, thus, prone to thermal irradiation.

The enhancement of water production shows a somewhat different pattern (Fig. 10c). The maximum increase still occurs within the concavity, but over the illuminated area (Fig. 10c, right-hand panel). We infer that this contrast arises from the high temperature for substantial water sublimation (e.g. above 180 K) and the exponential dependence of water sublimation on temperature (equation 4). While self-heating could induce a temperature increase by as much as 40 K within the shadow, the warming contributes little to the water production. On the contrary, even a small increase of temperature by a few K could significantly enhance sublimation flux over illuminated and warm surface areas.

We adopt the dusty ice thermal model to investigate the enhancement of total water production from 67P nucleus due to self-heating. The results are for three full nucleus rotations when the comet was inbound at 3.9, 2.8 and 1.4 au from the Sun, respectively. While self-heating accounts for an enhancement by nearly 50 per cent at 3.9 au (Fig. 11a), it yields a moderate contribution of 10 per cent at 2.8 au (Fig. 11b). Its role steadily diminishes as the comet approaches perihelion. At 1.4 au, the effect of self-heating is imperceptible, where the total water production of 67P is dominated by intense insolation close to perihelion (Fig. 11c). Overall, self-heating enhances but does not dominate total water production of 67P from a heliocentric distance of 4 au inwards. This conclusion applies to the mantled nucleus, as well.

Figure 11.

Enhancement of total water production from 67P over three full nucleus rotations due to self-heating, when the comet was inbound at 3.88 (a), 2.77 (b) and 1.43 au (c) from the Sun, respectively. The results are derived via the dusty ice model.

6 GLOBAL WATER PRODUCTION OF 67P AROUND PERIHELION

We apply the dust mantle thermal model to estimate the global water production rate of 67P as function of heliocentric distance. The shape model of 15 000 facets is adopted for the analysis (Jorda et al. 2016). Self-heating is neglected. Specifically, we average the water production over one nucleus rotation, such as

\begin{equation} n_Z(t_0) = \frac{1}{\hat{m}} \left( \frac{1}{t_{\rm P}} \displaystyle \int ^{t_0+t_{\rm P}}_{t_0} \Sigma Z {\rm d}t \right) \quad ({\rm molecule} \cdot {\rm s}^{-1}), \end{equation}

(31)

where |$\hat{m}$| is the mass of a water molecule (as in equation 3).

Hereafter, water production rate is uniquely defined as a time-averaged quantity and given in number of molecules per second. Thus, production rate is distinguished from outgassing flux, with the latter quantifying the instantaneous sublimation flux of water ice from the entire nucleus.

6.1 Evolution of water production rate of 67p

In Fig. 12, we show the water production rates from nucleus overlain by a dust mantle of X = 5 and 10 mm in thickness. Note that the choices of mantle thickness are not random; rather, the values are chosen to bound the estimate of 6 mm by Shi et al. (2016) derived from the observed duration of the dust activity that continued after sunset, probably sustained by thermal lag in the shallow subsurface where sublimation occurred. Furthermore, we adopt two values, f = 0.1 and 0.01, for the icy area fraction of the nucleus interior in order to characterize the evolution of the water production around perihelion. The f = 0.1 roughly accommodates the dust-to-water ratio of about 6, as reported by Rotundi et al. (2015) and Fulle et al. (2016c), and a somewhat larger value of 8.5 for the bulk nucleus, as reported by Fulle et al. (2016b). A low icy area fraction, f = 0.01, may not be surprising. It is plausible that the expansive dust deposits over the northern hemi-nucleus of 67P, formed by the deposition of ejecta from the south, contained a few per cent of water ice in mass (Hu et al. 2017).

Figure 12.

Total water production of 67P as function of time since perihelion in 2015, or heliocentric distance, estimated with the dust mantle thermal model for different mantle thickness and icy area fractions. The results are plotted in comparison with various measurements as presented in Fougere et al. (2016b). The estimates of water production based on DSMC (Direct Simulation Monte Carlo) are included (Fougere et al. 2016a). The first two rows of figure legend read ‘mantle thickness – icy fraction area’. The reader is referred to Fougere et al. (2016b), where the description of the measurement data and sources thereof are provided.

The model results are compared with the in situ measurements of water production rates by Rosetta instruments as well as estimates based on ground-based observations, as collected in Fougere et al. (2016a,b). We note that the two combinations of X = 5 mm with f = 0.01 and X = 10 mm with f = 0.1 could largely reproduce the evolution of the water production from 67P as measured, with notable underestimation only within about 50 d around perihelion (Fig. 12). It should be stressed that the data interpretation is non-unique and that the model assumption of constant mantle thickness is ideal and likely questionable, i.e. the ice front may vary on both diurnal and seasonal time-scales (De Sanctis et al. 2015; Fornasier et al. 2016). Nevertheless, the solutions of the thermal models clearly indicate the distinct roles of a prevalent dust mantle and limited ice abundance in regulating the water production of the nucleus.

6.2 Comparison with other thermal models

It is instructive to compare the model results presented above with some previous investigations devoted to characterizing the activity of 67P (see e.g. Davidsson & Gutiérrez 2005; Fulle et al. 2010; Keller et al. 2015, among others).

Keller et al. (2015) applied six thermal models to reproduce some historical measurements of water production of 67P around perihelion. Specifically, models A, B and C employ the ‘stationary’ approach, under the assumption that local thermal equilibrium is established instantly with varying insolation. The three models differ in the assumed depth of the ice front. Model A is a dusty ice model where the energy balance at the surface is similar to equation (1). Models B and C address the influence of a thin dust mantle; the mantle in model B is 50 μm in thickness and composed of dust particles of 10 μm in diameter (i.e. H = 5 in equation 10); a thicker dust mantle of 1 mm is assumed in model C that consists of 100 μm particles (for H = 10).

Three other models, D, E and F, are not restricted by the assumption of instant thermal equilibrium and account for the effect of thermal inertia, or thermal lag, of the nucleus subsurface in response to changing insolation. All three models fall in the class of dusty ice thermal models for bare icy nucleus surface. It is noted that model D applies to a homogeneous, or globally active, nucleus, whereas models E and F simulate water activity arising sporadically and randomly from the nucleus.

The dusty ice model in the present study coincides with model D in Keller et al. (2015) for a fully icy nucleus with f = 1 (in equation 2). However, a direct comparison between our dust mantle model and models B and C of theirs is difficult, to say the least. The thickness of the dust mantle investigated here exceeds the maximum of 1 mm considered by Keller et al. (2015).

The curves of water production derived via the dusty-ice and dust-mantle models presented in this work are superposed on results in Keller et al. (2015) (Fig. 13). We note that our dusty-ice model produces nearly identical results to those of model D (orange curve in Fig. 13), although estimates from both models have to be reduced to match measurements.⁷ It is worth noting that, even with this reduction, these estimates still exceed the measurements until well within 100 d before perihelion. On the other hand, our dust mantle thermal model, adopting X = 10 mm with f = 0.1 and X = 5 mm with f = 0.01, respectively (same as in Fig. 12), both yield better approximation of overall trend of increasing water production. It is (re-)emphasized that these results are not arbitrarily scaled for matching measurements. The results suggest that it is, at least, not necessary to invoke the assumption of strong inhomogeneity in terms of the distribution of active/dormant areas over the nucleus.

Figure 13.

Water production of 67P derived via dust mantle thermal model with varying mantle thicknesses and icy area fractions, in comparison with model results reported in Fig. 18 in Keller et al. (2015). The estimates of water production reported by Fougere et al. (2016a) are appended and marked by black dots. The orange curve is derived by the dusty ice thermal model with icy area fraction of 100 per cent, where results are scaled down by a factor of 6 per cent. This curve nearly overrides the dashed blue curve for model D in Keller et al. (2015) factored by 7 per cent. The thick solid curve in bright green is produced with the dust mantle model with mantle thickness of 10 mm and icy area fraction of 10 per cent. The thick solid curve in dark cyan corresponds to the mantle thickness of 5 mm and the icy area fraction of 1 per cent.

Two thermal models, developed respectively by De Sanctis, Capria & Coradini (2006) and Gortsas et al. (2011), were adopted in an earlier effort to study water production of 67P by Fulle et al. (2010). The distinction of these more sophisticated modelling approaches is briefly described in Section 2.4. We compare our results based on the dust mantle model for different parameterizations with those presented by Fulle et al. (2010) in Fig. 14.

Figure 14.

Water production of 67P derived via dust mantle thermal model with varying mantle thicknesses and fixed icy area fraction of 10 per cent, in comparison with model results reported in Fulle et al. (2010). The legend in the upper-right corner reads ‘mantle thickness – icy area fraction’. The colored curves are superposed on the original monochromatic figure of fig. 20 in Fulle et al. (2010). The reader is referred to the cited work for explanation of the data points and references thereof. The dashed dark curve is based on the thermal model of De Sanctis et al. (2006); the dotted dark curves are modelled by Gortsas et al. (2011) for different parameters. The red curve for the mantle thickness of 5 mm and the icy area fraction of 10 per cent coincides with the upper dotted curve.

It is somewhat interesting and, indeed, reassuring to note the coincidence of the curve for X = 5 mm and f = 0.1 (red curve in Fig. 14) with one of the curves produced by the model of Gortsas et al. (2011). It seems that another curve based on X = 10 mm is also in general agreement with the modelled water production within the heliocentric distance of 3 au in Fulle et al. (2010). The agreement suggests that the difference in the model formulations does not result in substantial, intractable discrepancies in simulations, and that the simplified treatment of the physical mechanisms in the present thermal models, e.g. neglecting mass diffusion of gas within the porous nucleus interior, is reasonable for analysing the seasonal evolution of the water production on 67P.

7 CONCLUSIONS FOR 67P: NEAR-CENTIMETRE THICK DUST MANTLE AND PLURAL DUST-TO-ICE RATIO

We presented two thermal models for characterizing the water activity on comet 67P. The dusty ice model applies to water sublimation from a bare, icy nucleus surface, while the dust mantle model addresses the impact of an insulating, gas-flow-resistant dust mantle on water outgassing. In order to apply the thermal models to 67P with an irregular-shaped nucleus, we constructed a Landscape data base based on the shape models, which facilitates and expedites the evaluation of energy input of insolation and self-heating over the nucleus.

We performed a variety of simulations to isolate and assess the influences of certain parameters of nucleus subsurface properties. The following conclusions were drawn:

• The thickness of dust mantle plays a vital role in modulating the sublimation flux underneath. The sublimation flux decreases with increasing dust mantle thickness. The mechanisms are twofold: first, the increasing thickness due to accumulation of more particle layers reduces the permeability of the mantle and, thereby, directly inhibits outgassing (Gundlach et al. 2011). The flow resistance of the mantle gives rise to the ‘pressure cooking effect’, where an increase of the ice temperature must occur with the quenching of outgassing by the mantle, that partly compensates for the sublimation flux (Kömle & Steiner 1992; Kömle et al. 1992). However, this effect is not pronounced for a dust mantle comprising coarse grains, for instance, those millimetres in size, as assumed in the present simulations. On the other hand, the insulation effect, or the low heat conductivity, of the dust mantle moderates the ice temperature underneath compared with that on an icy surface (on the dayside) that further reduces sublimation (Kührt & Keller 1994). In the case of 67P, the outgassing flux from beneath a dust mantle of 5 mm in thickness is 20 per cent of that from a bare nucleus at the heliocentric distance of 3.9 au.

• The dust-to-ice ratio influences the icy area fraction of the nucleus surface or subsurface from which sublimation occurs. As a result, ice abundance in the nucleus interior significantly affects water outgassing. In the case of a bare icy nucleus, an icy area fraction of 10 per cent restricts outgassing flux to about 20 per cent of that from a pure-ice surface.

• The irregular-shape of the nucleus gives rise to non-negligible self-heating effect. This phenomenon is more pronounced at larger heliocentric distances. However, self-heating probably does not dominate total water production that is still directly driven by insolation.

• We parametrized and applied the dust mantle model to approximate the evolution of observed and measured water production rates of 67P around the 2015 perihelion. Conclusions from previous studies tend to interpret the overestimation of water production by thermal models as indicating localized activity, or the existence of significant dormant surface areas over the nucleus. Our results show that this overproduction of thermal models could simply reflect a lack of deliberation on the role of a prevalent dust mantle or that of low ice abundance limiting the total water production of the nucleus.

With due caution in mind that data interpretation is rarely unique, we propose that the water production of 67P was sourced from a nearly uniform dust mantle whose thickness is between several mm to about 1 cm and that the ice abundance is between a few and about 10 per cent on average. In the case of 67P, the assumption of a locally active nucleus surface may not be necessary at all for model interpretation of water activity.

Acknowledgements

We thank Dr. Norbert Kömle for his careful assessment of an earlier manuscript and providing constructive criticisms and detailed suggestions. We are particularly grateful for his advices that enabled us to improve the discussion about the influence of the dust mantle on cometary water activity.

XH acknowledges financial support of the International Max Planck Research School for Solar system Science at the MPS in Göttingen, Germany.

1

They are, in chronological order of explorations, 1P/Halley (Keller et al. 1986), 19P/Borrelly (Soderblom et al. 2002), 81P/Wild 2 (Brownlee et al. 2004), 9P/Tempel 1 (A’Hearn et al. 2005), 103P/Hartley 2 (A’Hearn et al. 2011) and 67P/Churyumov-Gerasimenko (Sierks et al. 2015).

2

Whipple himself did not suggest the use of equation (1) for energy balance in his dusty snowball model of cometary nucleus. In fact, he specifically envisaged the presence of a dust mantle due to de-volatization in the topmost layers of the nucleus.

3

The area of a facet can also be derived directly as |$|\boldsymbol {n}|/2$|⁠.

4

Higher resolution data base is used in practice, where altitude and azimuth are discretized at 1- and 2-deg intervals, respectively.

5

That is, |$X=x_\mathcal {J}+\frac{\Delta x}{2}$| in equation (8).

6

The impact of self-heating is discussed in Section 5.3. The effect is neglected in all other analyses in this work.

7

The dusty-ice model for f = 1 overestimates the total water production of 67P and is therefore scaled by a factor of 6 per cent. Model D in Keller et al. (2015) likewise needs to be downscaled by 6.9 per cent.

8

Defined in accordance with the Cheops reference frame (see Preusker et al. 2015).

APPENDIX A: LOCAL HORIZONTAL COORDINATE SYSTEM

Let us establish a right-handed Cartesian coordinate system centred at a given point on the nucleus surface with z-axis aligned with the surface normal (Fig. A1). The direction of the local prime meridian, i.e. the x-axis, is arbitrary within the horizontal plane. With a shape model, the origin is then chosen to be the centroid of the given facet, |$\boldsymbol {r}_{\rm F}$|⁠. The three coordinate axes are specified by three unit column vectors as

\begin{eqnarray} {\boldsymbol u}_z &=& \hat{{\boldsymbol n}} = \frac{ (\boldsymbol {v}_2-\boldsymbol {v}_1) \times (\boldsymbol {v}_3-\boldsymbol {v}_1) }{| (\boldsymbol {v}_2-\boldsymbol {v}_1) \times (\boldsymbol {v}_3-\boldsymbol {v}_1) |},\nonumber\\ {\boldsymbol u}_x &=& \frac{ {\boldsymbol v}_2 - {\boldsymbol v}_1 }{ \left|{\boldsymbol v}_2 - {\boldsymbol v}_1\right| } ,\nonumber\\ {\boldsymbol u}_y &=& \frac{ {\boldsymbol u}_z \times {\boldsymbol u}_x }{ \left| {\boldsymbol u}_z \times {\boldsymbol u}_x \right| }, \end{eqnarray}

(A1)

with |${\boldsymbol v}_1$| and |${\boldsymbol v}_2$| locating the first and second vertices of the facet, respectively. |$\hat{{\boldsymbol n}}$| denotes the unit normal vector to the facet. Note that all quantities are referred to the body-fixed (BF) coordinate system⁸ of 67P. It can be easily shown that the transformation of a vector, |${\boldsymbol r}$|⁠, from the body-fixed coordinate system of 67P to the local horizontal coordinate system is given by

\begin{equation} {\boldsymbol r}^{[{\rm H}]} = {\boldsymbol R}^{{\rm [H]}}_{[{\rm BF}]} \left( {\boldsymbol r} - {\boldsymbol r}_{\rm F} \right), \end{equation}

(A2)

where the 3 × 3 rotation matrix is formulated as

\begin{equation} {\boldsymbol R}^{{\rm [H]}}_{[{\rm BF}]} = {\left(\begin{array}{ccc}\, {\boldsymbol u}_x &\, {\boldsymbol u}_y &\, {\boldsymbol u}_z \end{array}\right)}^{\rm T}. \end{equation}

(A3)

The superscript ‘[H]’ is adopted to express that a vector is referred to the horizontal coordinate system.

Figure A1.

Local horizontal coordinate system defined for an individual facet.

For the purpose of projecting an object on the topocentric hemisphere, it is natural to introduce a spherical coordinate system are radius, r, altitude, φ, and azimuth, α, such as

\begin{eqnarray} r &=& \sqrt{\left(x^{[{\rm H}]}\right)^2 + \left(y^{[{\rm H}]}\right)^2 + \left(z^{[{\rm H}]}\right)^2},\nonumber\\ \alpha &=& \arctan \left(\frac{y^{[{\rm H}]}}{x^{[{\rm H}]}}\right),\nonumber\\ \varphi &=& \arcsin \left(\frac{z^{[{\rm H}]}}{r}\right). \end{eqnarray}

(A4)

For practical calculations, the two-argument function, atan2(x, y) ∈ ( − π, π], should be used to evaluate α, i.e.

\begin{eqnarray} \alpha = {\rm atan2} \left( x^{[{\rm H}]},y^{[{\rm H}]} \right) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\, \arctan \left(\frac{y^{[{\rm H}]}}{x^{[{\rm H}]}}\right), & {\rm if}\quad x > 0, \nonumber\\ \, \arctan \left(\frac{y^{[{\rm H}]}}{x^{[{\rm H}]}}\right) \pm {\pi} , & {\rm if}\quad x < 0 , \nonumber\\ \, {\rm sgn}\left(y^{[{\rm H}]}\right) \frac{{\pi} }{2}, & {\rm if}\quad x = 0,\\ \end{array}\right. \\ \end{eqnarray}

(A5)

with sgn being the sign function.