© ESO, 2015

1. Introduction

The Jupiter family comet 67P/Churyumov-Gerasimenko (hereafter 67P) was discovered in 1969. Its current orbit is characterized by an eccentricity e = 0.64102, an inclination i = 7.0405°, a perihelion distance q = 1.2432 AU, and a semi-major axis of 3.4630 AU resulting in a period P = 6.44 years. In 1959 67P had a close encounter with Jupiter which changed its perihelion distance from 2.74 AU to its present value (see Krolikowska 2003). Triggered by a launch delay of the European Space Agency Rosetta mission comet 67P was chosen as a suitable replacement for its original target comet 46P/Wirtanen.

Launched in 2004 the ESA Rosetta spacecraft arrived at comet 67P in August 2014 entering the rendezvous phase. The scientific optical camera system OSIRIS comprises the Narrow Angle Camera (NAC) and the Wide Angle Camera (WAC; Keller et al. 2007). Analysis of ground-based and OSIRIS photometric observations revealed that the spin period of 67P had decreased by 0.36 h to 12.4043 ± 0.0007 compared to the value derived during its last apparition (Mottola et al. 2014). Once the cometary nucleus was resolved in July 2014, its spin axis orientation could be determined as RA = 69° and Dec = 64°, resulting in an obliquity of 52°. When Rosetta arrived at the comet only the northern hemisphere and a near-equatorial part of the southern hemisphere were illuminated. A digital terrain model based on almost 1000 OSIRIS/NAC images has a resolution between 5 and 10 m on the visible surface (Preusker et al. 2015; Sierks et al. 2015). The nucleus of 67P is characterized by a distinct bi-lobate appearance where the bigger lobe is connected to the smaller by a relatively small neck resulting in a huge cavity (Fig. 1). The lobes do not contact near their long axes as earlier surmised for comets 1P/Halley (Keller et al. 1986), 19P/Borelly (Oberst et al. 2004; Britt et al. 2004), 103P/Hartley 2 (Thomas et al. 2013), and 8P/Tuttle (Harmon et al. 2010). The north pole is located on the bigger lobe near the neck.

Fig. 1 Shape and surface of comet 67P as seen from the Sun at three rotational longitudes (30°, 150°, 270°) at northern solstice when the sub-solar latitude is + 52°, at the equinoxes, and at southern solstice. The smooth parts of the surface indicate yet unobserved territory on the south side of the shape.

The volume of the body and its surface are not yet known to high accuracy because the southern hemisphere is not yet observable. The volume derived from the current shape model is 21.4 km³ and combined with the mass 1.0 × 10¹³ kg yields a density of 470 [kg m^-3] (Sierks et al. 2015). This low density requires a high porosity of 70 to 80% considering that the density of an ice-dust mixture should be around 2000 [kg m^-3]. This suggests a highly porous interior of agglomerates confirming earlier indirect determinations (A’Hearn et al. 2005; Davidsson & Gutiérrez 2005). About 70% of the surface was accessible when constructing the shape model. Its total extent can only be estimated at this time to be 4.82 × 10⁷ m².

The surface of the neck area (called Hapi, see Thomas et al. 2015b; and El-Maarry et al. 2015, for region designations) is rather smooth and is slightly blue compared to the reddish overall spectrum (Capaccioni et al. 2015; Fornasier et al. 2015). Here, dust features were observed by OSIRIS (Vincent et al. 2015) when 67P was still at a heliocentric distance of 3.7 AU. Other, more strongly illuminated areas did not yet show signs of activity at the time. However, the mm/submm wavelength instrument MIRO observed water (Gulkis et al. 2015). Somewhat later H₂O, CO₂, and CO were detected in situ by the ROSINA mass spectrometer (Altwegg et al. 2015; Hässig et al. 2015). In spite of the ubiquitous water molecules in the coma, water ice has not been observed on the cometary nucleus surface. The VIRTIS spectrometer sets an upper limit of 1% for exposed water ice surface coverage (Capaccioni et al. 2015). However, De Sanctis et al. (2015) found water ice near light/shadow boundaries in VIRTIS spectra and interpret this as a re-condensation of water molecules. It was suggested that the overall energy input of Hapi is strongly enhanced by IR radiation emitted from illuminated surface elements of the Seth and Hathor regions of the big and small lobes (Sierks et al. 2015). Indeed the view factors of many surface facets are strongly restricted, as expected within such a concavity. The overall insolation of the neck (area between the lobes) is, however, low because this area is shadowed by one or the other lobe for most of the spin period when the comet approaches perihelion. The complex shape of the northern hemisphere requires careful modeling based on a high resolution shape model to assess the energy input and sublimation.

Fig. 2 Observational parameters during the present orbit of 67P. The left figure shows heliocentric distance of the comet and solar latitude vs. time. The right figure demonstrates the asymmetry of the insolation as a function of heliocentric distance.

The large obliquity (52°) and orientation of the spin axis causes strong seasonal effects. The orientation of the spin axis of 67P is such that southern solstice occurs only 34 days after perihelion at a heliocentric distance of r_h = 1.24 AU. Consequently the northern solstice is reached after aphelion, at r_h = 5.68 AU, last time in August 2013. The equinoxes are passed at r_h = 1.7 AU, 82 days before perihelion (May 2015), and at r_h = 2.6 AU, 7.5 months after perihelion (Fig. 2). This geometry leads to a dichotomy: the northern hemisphere is relatively weakly illuminated for a long time, when the comet is far from the Sun, whereas the southern hemisphere is strongly insolated during the short perihelion passage. During the respective winter seasons the poles are completely dark.

The high resolution shape model of the nucleus makes it possible to describe the diurnal insolation conditions in detail for areas such as the neck, crater-like formations, pits, and cliffs. Sublimation of volatiles requires insolation and essentially depends on the surface temperature. The measured thermal inertia of cometary surfaces is small (Groussin et al. 2007). In the case of 67P it was found by the MIRO team to be between 10 and 50 [J m^-2 K^-1 s^−1/2] (Gulkis et al. 2015). This low value requires a very low thermal conductivity in the range of 10^-4 to 10^-3 [W m^-1 K^-1]. Then the skin depth of the diurnal heat wave is only a few millimeters. Under these circumstances the characteristic time scale of heat diffusion is short. The sublimation rate then depends on the instantaneously available energy that reaches the water ice.

In Sect. 2 we describe the model calculations of the insolation in detail, followed in Sect. 3 by the water sublimation models, water ice on the surface or under a thin layer of porous dust. These conditions justify steady state and do not require taking thermal memory effects into account. The effects of the insolation and erosion models are discussed in Sect. 4. We compare the model results with the Rosetta observations in Sect. 5. The asymmetry of the seasonal insolation influences the shape and morphology of the nucleus (Sect. 6). Finally we discuss the validity of the simple steady state sublimation models and the assumption of a uniform surface by applying a model that includes thermal memory effects but can only be calculated for much worse spatial resolution (Sect. 7) and we summarize our results in Sect. 8.

2. Insolation model

The shape of 67P is approximated by a polyhedron made of triangular facets. This model consists of about 108 000 facets derived by decimating a high resolution model based on images acquired during the SHAP4 orbital phase of Rosetta. This SHAP4 model has been derived by means of stereo-photogrammetric techniques (Preusker et al. 2015). The insolation of a facet is derived based on its orientation with respect to the Sun. For this purpose, the orbit around the Sun is discretized in steps of approximately 2.5° in true anomaly. The spin is discretized in steps of 2% in rotational phase to follow the rate of change of illumination. For each time step we check mutual occlusion of the facets – because the body is concave. In addition to the direct solar insolation, the energy balance on the comet is affected by sunlight reflected by mutually facing facets and by their thermal IR radiation, which result in a local flux enhancement in the presence of concavities (Davidsson & Rickman 2014). Our model includes these effects.They are modeled by computing the mutual view factors F_ij (the fraction of power radiated by facet i illuminating the surface of facet j) for all facets of the shape model (Fig. 3). Efficient computation of the view factors is achieved by implementing a hierarchical algorithm that significantly reduces the number of necessary checks, and is based on horizon-based occlusion detection, coupled with a variant of the painter’s algorithm (Hearn & Baker 1994). As the view factors depend only on the shape model, and not on the illumination/viewing geometry, they are pre-computed and stored, and accessed at each time step. The thermal emission of a facet depends on its surface temperature and hence on the amount of direct insolation and irradiation from neighboring facets. The surface temperature can be solved iteratively. As a first step the equilibrium temperature is computed by ignoring self-illumination. In the following steps this effect is added based on the temperature derived from the previous iteration. This algorithm converges very rapidly, such that only three iterations are needed to keep the convergence error below a few percentage points.

Fig. 3 Viewing factors of the model facets. A view factor of 0 (black) implies that a facet is not exposed to irradiation from other facets. A view factor of 1 implies that the field of view of a facet is completely filled by other facets of the body.

3. Sublimation models

The energy absorbed by the cometary surface calculated in the previous section can be partly used to sublime volatiles on or under the surface. The amount of absorbed energy is only a rough indication of the sublimation capability of a facet because other processes also consume energy. The temperature of the surface can be calculated by taking into account

• the balance of irradiation, and

• losses due to re-radiation, sublimation, and by conduction either vertically or laterally. We concentrate on water ice sublimation to provide examples of surface erosion due to activity.

Water ice on surface. The simplest assumption is that exposed water ice sublimates directly from the surface. The ice must be mixed with dust, considering the low albedo of cometary nuclei. An estimate of exposed ice detectable in the IR was determined for comet Tempel 1 (Sunshine et al. 2006) and an upper limit of 1% for 67P is provided by VIRTIS (Capaccioni et al. 2015).

For a flat surface, if heat conduction is neglected (justified by the low heat conduction coefficient observed; Gulkis et al. 2015) energy conservation leads to (1)where I is the solar irradiation, A_v = 0.01 the bolometric Bond albedo, ε = 0.9 is the emissivity, σ the Stefan-Boltzmann constant, T the surface temperature in K, Z the sublimation rate in [kg m^-2 s^-1], and L the latent sublimation heat in [J kg^-1]. We assume L constant, neglecting the slight temperature dependence. The sublimation rate is given by the Hertz-Knudsen formula Z(T) = 2P(T)/(πv_th), where the water vapor pressure is P = 3.56 × 10¹²exp(−6141.667 /T) [kg m^-1 s^-2] (Fanale & Salvail 1984) and the thermal velocity [m s^-1] with R the gas constant and μ the molar mass of water.

Facets of a nonconvex nucleus shape can be shaded but also receive energy from (nearby) facets that limit their unobstructed (2π) field of view (see Sect. 2) and hence the left-hand side of Eq. (1) needs to be rewritten. The total absorbed power in [J m^-2 s^-1] is (2)The additional terms R_vis(t) and R_IR(t) are the reflected visible and IR contributions, respectively, and E(t) is the emitted temperature-dependent IR input from the facets within the field of view of the reference facet (Rozitis & Green 2011). Nonlinear Eq. (2) has to be solved for the temperature T. The IR emitted and reflected components of the mutual facet illumination (self-illumination) are computed iteratively. Then the sublimation rate Z can be calculated and integrated over the time step. So far, only small bright spots have been found to be widespread (Sierks et al. 2015) on the surface of 67P, but they have not been clearly identified as water ice or rich in water ice. Hence all the water ice is covered by a dust layer. This layer has to be thin enough to still pass the solar irradiation to the water ice and also fragile enough to be removed by the vapor pressure of sublimating water ice so that activity is not quenched by a growing dust layer.

Two-layer dust-ice model. Recent models of cometary nucleus formation due to gravitational collapse (Johansen & Klahr 2011; Skorov & Blum 2012) but also the hierarchical agglomeration model (Weidenschilling 1997) lead to an extremely porous mixture of large ice and dust agglomerates of distinct sizes. An extremely low tensile strength of the cometary nucleus material is achieved by a hierarchical distribution of agglomerates (agglomerates of agglomerates). The surface is covered by a dark thin ice-free layer which forms when the ice starts to sublimate. After losing its ice content, the silicate dust layer slightly collapses, but its morphology stays similar to the initial morphology of the ice-dust mixture. Below this refractory layer the pristine ice-dust medium exists. The inert layer is removed by the increasing gas pressure of the subliming ice. By applying models for the heat transfer and permeability the activity of comets can be simulated (Blum et al. 2014).

The ice-dust mixture is characterized by the filling factor (1 − porosity) and by the permeability of the refractory dust layer. The individual agglomerates have a volume filling factor of ϕ_a. The volume filling factor of the ice-free dust agglomerate packing structure is ϕ_p and the volume filling factor of the entire dust layer is then given by ϕ_aϕ_p. We assume (following Zsom et al. 2010) that ϕ_a = ϕ_p = 0.4. Thus we get a volume filling factor of the ice free surface layer of ~0.16 (which is between the values corresponding to the two extreme scenarios considered by Skorov & Blum (2012). The thickness of the ice-free dust-agglomerate layer and the agglomerate size are free parameters. The limited permeability of a porous refractive layer greatly reduces the flow of gas for all regimes of gas diffusion (from Knudsen diffusion to a continuous flow). This is compensated by an increase in pressure according to the energy available. A collisionless gas flow (Knudsen regime) can be assumed for the mass transport through the porous dust layer, and we use a modified Clausing formula (Skorov et al. 2011) and an empirical experimental formula given by Gundlach et al. (2011). Thus the effective gas production under a porous refractive layer can be evaluated for a given layer thickness, grain size, and resulting surface temperature of the ice.

The amount of energy available for sublimation is determined by the heat conductivity through the dust layer. Generally the conductivity through the solid medium, the conductivity via thermal emission inside the pores, and conductivity via gas diffusion are considered. The last component is neglected based on the low gas density for the considered conditions. A detailed description of the first two mechanisms is presented by Gundlach & Blum (2012). We simulate the effective conductivity considering the hierarchic structure of the medium (porous layer consisting of porous aggregates).

A stationary approach to the energy balance neglecting thermal history effects can be applied if the characteristic time scales for gas diffusion and heat conduction are smaller than the characteristic time of irradiation changes and the integration step. In this case one has only to consider energy balances on the top and on the bottom of the refractive porous layer in order to evaluate the resulting gas flow and the resulting gas pressure above the sublimating layer and one can neglect the history. Thus Eq. (2) can be generalized to the following nonlinear algebraic equation with the unknown temperatures T_s and T_i determined by energy conservation at the surface and at the sublimation front below the dust layer, (3)where (4)Here Ψ_perm is the permeability as a function of the agglomerate radius r_A and thickness of the dust layer l, is the effective heat conductivity of the dust layer, and controls the flux into the core of the nucleus. The equilibrium orbital temperature T_core and orbital skin depth l_core are evaluated following McKay et al. (1986). This system of equations can be solved for the most general case when the effective conductivity is a nonlinear function of temperature (due to radiation inside the pores; Gundlach & Blum 2012) and the sublimation energy of ice is a function of its temperature T_i (Cowan & A’Hearn 1979).

The stationary case is only valid as long as the characteristic time of heat diffusion τ_hd through the dust layer is short compared to the time step (15 min) used to calculate the received energy input onto a facet (see Sect. 2). The effective heat conductivity of the dust layer, , for the parameters r_A and l used in our calculations does not exceed 2 × 10^-3 [W m^-1 K^-1] even at perihelion for 67P. This is in good agreement with the observed low thermal inertia between 10 and 50 [J m^-2 K^-1 s^−1/2]. The characteristic time τ_hd is then 5 min or less when l< 1 mm, smaller than the integration time step. The heat conductivity of the bulk dust-ice mixture is about 50% higher than .

Fig. 4 Surface temperature as function of the irradiance for a gray body and models A to C. Model A represents the lower limit for the temperature of an active surface. The temperatures of model C come close to that of an inactive gray body.

Fig. 5 Water production rates vs. time for models A to C with (black lines) and without (red lines) self-illumination (SI).

Fig. 6 Erosion and insolation integrated over one orbit. The left figure shows the erosion (model A). The maximum water erosion on the (unseen) south side is 4.8 [m]. The right figure shows the insolation on the same color scale with the maximum of 7.54 × 109 [J m-2].

Model parameters. In addition to the case of water ice on the surface (model A) we calculate the effects of a refractive porous dust layer. We use (model B) a thin layer of small aggregates with a radius r_A = 5μm and a thickness of l = 50μm (5 aggregate layers) and (model C) a thick layer (limited by τ_hd to be short) with agglomerates of 50 μm and l = 1 mm (10 aggregate layers). The temperatures T_s and T_i and the water production rate are tabulated as a function of the absorbed energy over the range of heliocentric distances of comet 67P. This range is extended for values to cover additional energy input from neighboring illuminated facets and on the lower end to cover weak absorption due to large inclination angles.

Model A with dark water ice on the surface sets the upper limit for the sublimation rate but is unrealistic because water ice on the surface of 67P (and of other comets) has not been detected or only in minor (<1%) amounts (Capaccioni et al. 2015). Model B uses a very thin layer of dust above the ice/dust mixture of the nucleus, about the minimum to make the ice invisible to spectroscopy. Near perihelion the influence of this dust layer on the sublimation rate is minimal. The dust layer of model C comes close to the diurnal thermal skin depth of about 5 mm for [W m^-1 K^-1]. The upper limit of the surface temperature of a gray body is reduced by the heat needed to sublime volatiles (water ice) near the surface. A lower limit is given by model A with ice directly on the dirty surface (Fig. 4). For model B the surface temperature is slightly higher than that of model A. Model C approaches the gray case representing a relatively inert dust surface as observed during the beginning of the Rosetta rendezvous phase with comet 67P. It could be used to evaluate the maximum effect of self-illumination and to simulate inert surfaces. Near perihelion the thin dust layer of model B hardly restricts the sublimation, which almost matches the production rate of model A. The peak production rate of model C is about a factor of 4 to 5 lower than that of dirty ice. The water production integrated over the whole nucleus is slightly enhanced at larger heliocentric distances if self-illumination is taken into account (Fig. 5). The enhancement coincides with northern summer. The small increase in surface temperature enhances the albeit minute sublimation rate at large heliocentric distances. The dirty ice is most sensitive to the additional power input by self-illumination during the aphelion passage.

Erosion rate. The local erosion rates of cometary material at a specific surface element (facet) can be estimated by Δx/ Δt = Z(T) /ρ where Z(T) is measured in [kg m^-2 s^-1] and ρ is the density of water ice in the nucleus. Δt is the time step of integration over which Z(T) stays constant. ρ can be found by correcting the density of the nucleus (ρ_N) by the dust to ice ratio, (γ) as ρ = ρ_N/ (1 + γ).

4. Discussion of insolation, sublimation, and erosion

The southern hemisphere is strongly heated during and near the perihelion passage of 67P over 10 months from May 2015 to March 2016 when the sub-solar latitude is negative because southern solstice occurs 34 days after perihelion (Fig. 2). However, this relatively short intense insolation is counterbalanced by the long time when the comet is outside 2.6 AU. Sublimation depends not only on the amount (integral over time) of insolation, but particularly at large heliocentric distances also on the intensity of the insolation. Figure 6 compares the total erosion and insolation of the northern hemisphere, both integrated over the whole orbit. Both representations use the same color scale. The maxima are of course not visible because they lie on the southern hemisphere. The comparison shows that only about 25% of the maximum erosion is reached on the northern hemisphere while the total insolation reaches 80%. This clearly demonstrates the nonlinearity of the sublimation with temperature.

Fig. 7 Production rate for model A. Dashed lines indicate the times of perihelion and maximum production.

Therefore it is not surprising that about 75% of the total water production, integrated over a whole orbit, is produced during southern summer from equinox pre to equinox post perihelion (Table 1). This value varies only very slightly from model to model. The asymmetry of the sub-solar latitude variation with respect to perihelion leads to an asymmetry of the water production rate. Its maximum is reached 6 days after perihelion (Fig. 7). As a consequence the water production is slightly higher post perihelion at the same heliocentric distance. This trend reverses around 2 AU and at large heliocentric distances the production after aphelion is slightly higher than before. This trend is independent of the sublimation model. The asymmetry is weak considering the complex shape of 67P and the strong obliquity of its rotation axis. Overall regional variations of the insolation, accumulated over one orbit, do not exceed a factor of two. As expected, model A with water ice on the surface results in the strongest production rate over the whole heliocentric distance range. Even at aphelion the water production is still about 2 [kg s^-1] (Fig. 8), which may mainly derive from neglecting the heat conduction (see also Sect. 7). The more realistic model B with a thin layer of dust on top and heat conduction shows an increasing reduction if compared to model A at heliocentric distances beyond 2.5 to 3 AU. Near perihelion the water production of model B almost reaches the rates of model A. The “resistance” of the thin dust layer to the heat penetration is so small that all energy is used for sublimation as in model A. At large heliocentric distances the thin dust layer quenches the production almost as strongly as in model C where a thicker dust layer is considered. In the case of model C the production rate of water is subdued by about a factor of 4 to 5 at heliocentric distances below about 2.5 AU and runs essentially parallel to model A. There is no upturn visible as for model B. It is interesting to note that the enhancement of the sublimation due to self-illumination for model C is in the same range as for models A and B although the considerably higher surface temperature (Fig. 4) leads to enhanced heating by the IR radiation of the neighboring facets. The total sublimation over the whole orbit increases by 3% for models A and B and 5% for model C. Overall, self-illumination smoothens the contrast between highly and less illuminated morphologies (Fig. 9) and hence leads to more uniform sublimation.

Fig. 8 Water production rates of models A to C versus heliocentric distance.

The spatial resolution of our model calculations using the underlying SHAP4 model is good enough to describe the illumination within the big pits found in the Seth region (Fig. 3). Here one finds of course the strongest self-illumination but on an overall low insolation level. Inside such a pit with a diameter of 220 m and a depth of 185 m (Vincent et al. 2015) all facets have high mutual view factors. If higher resolution shape models will become available the computational requirements to proceed with the present approach will be too high to be useful. Focus will be given to spatially restricted but detailed comparisons of observations and models. Our method can easily be applied to cutouts of the surface as long as the total number of facets is on the order of 10⁵. We will provide an example in Sect. 5 (Fig. 14). Usually we integrate the absorbed insolation of each facet over the whole orbit, over part of the orbit, or over one spin of the comet. If the results of each integration step are saved we provide the diurnal variation of the insolation or the temperature of the surface element. Here again the assumption is that the temperature corresponds to the instantaneous insolation and does not depend on any temperature history.

The total production of water calculated by integrating over one orbit and the 4.82 × 10⁷ m² cometary surface is 6.5 × 10¹⁰ kg, 4.9 × 10¹⁰ kg, and 1.5 × 10¹⁰ kg for models A, B, and C, respectively (see Table 1). We assume the surface to be uniform and have to keep in mind that the southern side of 67P is not yet illuminated and hence we do not know its morphology or even topography. During the early activity phase of 67P in July to September 2014 a dust to gas mass ratio γ = 4 was derived from observations in the coma (Rotundi et al. 2015). Then the mean erosion of surface material per orbit, averaged over the whole surface, reaches values of 14.5, 11.0, and 3.4 m for the three models, respectively. The corresponding mass loss is a fraction 3.3 × 10^-2, 2.5 × 10^-2, and 7.5 × 10^-3 of the nucleus mass of 1.0 × 10¹³ kg (Sierks et al. 2015), respectively, on the order of 1% per orbit. These numbers refer to a fully active surface. Comparison with observations suggests that the water production reaches only about 6% (see Sect. 5).

Usually we refer to the pure ice erosion rate and leave it up to the reader to apply the appropriate, possibly local dust to ice ratio. One should keep in mind that our models provide erosion rates under specific model assumptions. These will not be fulfilled on many localities due to the morphology of the surface such as coverage by thick dust (refractory) layers. In particular model A provides an overestimated upper limit of the erosion, better the potential for erosion based on the actual insolation. The average and maximum erosion rates are listed in Table 1. They vary between maxima of 5 to 1 m. The comet could lose up to about 20 m per orbit at specific locations (for γ = 4). The peak water erosion of the comet (model A) in August 2014 was about 1 mm per spin (Fig. 10). Until the end of 2014 the comet should have lost up to 1.5 m (for 300 spins and γ = 4) on well insolated spots, comprising a small percent of the northern areas. Changes of the morphology should become detectable for the science cameras OSIRIS. An admixture of super volatiles such as CO₂ or CO will increase the erosion.

The dichotomy of erosion between the northern and southern side of 67P is very distinctive. The results for model B are shown in Fig. 11. The south side is eroded nearly uniformly by close to 4 m per orbit (possibly because of the lack of known topography) whereas the north facing regions erode much less rapidly but with variations by a factor of 2 due to the rough topography. The equatorial areas that are inclined towards the south are well eroded and display a smoother topography than the north.

Fig. 9 Effect of self-illumination. The left image shows the energy absorbed from direct sunlight integrated over the whole orbit. For the right image the energy reflected and radiated (thermal) from other facets in the field of view of each facet is included.

Fig. 10 Water erosion (model A) integrated over one cometary spin on 9 August 2014. The northern hemisphere with the neck area is shown. The erosion is enhanced along the transition from the valley bottom (Hapi) to the cliff wall (Hathor).

The water production rate in August 2014 was 1.2 [kg s^-1] (Gulkis et al. 2015), when the heliocentric distance of 67P was 3.6 AU. This is 2% of the 61 [kg s^-1] that our calculation yields for model B. During the last apparition of comet 67P a few days after perihelion the maximum water production was determined by the SWAN instrument onboard SOHO (Bertaux et al. 2014) to 1.73 × 10² [kg s^-1]. This is 6% of 2.7 × 10³ [kg s^-1], the maximum of model B. In other words only 6% of the southern hemisphere will show the activity predicted by model B. The increase in the water production from 3.6 AU to perihelion in our models is less steep (even for model B, see Fig. 8) than indicated by the observations. A perfect match cannot be expected because the topography of the southern hemisphere is not yet known. The surface activity is not necessarily homogeneously distributed (we discuss this in Sect. 7). Within the validity of model B one could conclude that the active percentage on the northern hemisphere is only a third of the value on the south. On the northern hemisphere most horizontally-oriented areas are covered by what looks like dust (Thomas et al. 2015b). Here the activity is quenched due to the very low thermal conductivity. Model A and C provide even higher production rates at 3.6 AU relative to their perihelion activity, than Model B – their slopes with heliocentric distance are less steep. This is to be matched by an even smaller percentage of active area on the northern relative to the southern hemisphere. The more intense insolation during perihelion passage can maintain activity on larger parts of the surface.

The in situ instruments on board Rosetta may be interested in the modeled diurnal variation of the water sublimation rates. We calculated samples for dates in 2014 and 2015. The diurnal water production rate (Fig. 12) varies slightly at the beginning of activity when the Sun is still at high northern latitudes and shines continuously on the north polar region. With time the Sun’s latitude decreases and the insolated parts of 67P vary more and more with the cross section of its elongated shape (see Fig. 11). The variation of the water density observed by ROSINA (Hässig et al. 2015) in August 2014 reflects the visibility of illuminated portions of the nucleus as seen from the spacecraft.

5. Comparison with observations

The neck (Hapi, Hathor, Seth). The overall shape of 67P deviates strongly from an essentially convex ellipsoid, being dominated by the large concavity, the neck, where the two lobes join. This shape is somewhat reminiscent of comet 1P/Halley that appears bi-lobate in heavily processed images taken during the Vega flyby missions (Sagdeev et al. 1986). The neck valley sets itself morphologically apart from the rest of the northern hemisphere by being flat and smooth, seemingly dust covered. It covers a gravitational low area that makes accumulation of refractory material (dust) probable (Pajola et al. 2015). It girds half of the comet on a length of more than 2 km separating the small from the big lobe. Its width varies around 500 m. On the side of the small lobe it is flanked by an almost perpendicular (relative to the local gravity field) wall, Hathor, of almost 1 km height. On the other side, the Seth region of the big lobe rises steeply with terraces. It is obvious that direct insolation of the shielded valley is less than on the convex tops of the lobes. Why then was dust activity of 67P detected here by OSIRIS in early August 2014 (around a heliocentric distance of 3.7 AU)?

The north pole of 67P lies close to the Hapi (neck) region (see Fig. 16). Therefore Hapi and the whole north facing surface are only illuminated during northern summer when the comet is at heliocentric distances >1.6 AU pre-perihelion and >2.6 AU post-perihelion (see Fig. 2).Our calculations show that the total energy input in the neck area is strongly enhanced by self-illumination (reflected light and IR radiation) from the walls of the valley (Fig. 3). In August and September 2014 (r_h< 3.6 AU) the erosion is relatively strong where Hapi and the high wall of Hathor join. Only a few minor, isolated places show higher erosion rates (see Fig. 10). Figure 14 shows a detailed comparison of the model output and an image taken by the OSIRIS NAC camera. The erosion is about 30% higher at the yellow band along the Hathor wall than in the center of the valley. The surface temperature along this band of strong sublimation is locally enhanced. A detailed comparison of the modeled surface temperatures and the VIRTIS observations is currently in progress based on model C which comes closest to the temperature variation of an inert surface (see Fig. 4). This relative enhancement seems sensitive to the specific insolation geometry and weakens and may disappear when the sub-solar latitude decreases.

Fig. 11 Three different views of the shape showing total (integrated over one orbit) water erosion (model B). Left: northern areas and neck, middle: equatorial view, right: southern, yet unobserved side.

Fig. 12 Daily variation of water production (model A) vs. rotational phase. The dates are 9 August 2014, when jets started to be observed at the neck area, 14 November 2014 (landing of Philae), 31 January 2015, 2 April 2015, 9 May 2015 (equinox), and 9 August 2015 (perihelion).

Fig. 13 Dust activity of 67P. The left image is a combination of two WAC images taken on 10 January 2015. The dominating fan-like dust feature above the neck region is strongly enhanced by integration along the line of sight, parallel to the fan plane from the active neck valley. This effect is demonstrated by the small images on the right. The lower image was taken with a similar line of sight as the left image, whereas the nucleus is slightly rotated in the upper view taken half an hour earlier. The enhancement by the integration effect does not appear.

Fig. 14 A comparison of the neck area (Hapi with Hathor and Seth). The model shows erosion (model A) integrated over one spin on 12 September 2014. The strongest erosion is on a cliff in the Ma’at area (see Fig. 15) The erosion along Hapi is concentrated at the foot of the Hathor wall. The enhancement is visible as a yellow band. The surface temperature (based on model C) is as expected also high in that region, depending on the local time of day.

Fig. 15 Activity in the Ma’at area. In the left figure (Vincent et al. 2015) the pits 1 and 2 were identified as exhibiting dust activity. Indeed, our model for 12 September 2014 (right figure) reveals small areas of strong erosion in both pits. Below Ma’at 1 there is a cliff that shows the strongest erosion of the entire illuminated surface (compare with Fig. 10).

When the comet approaches equinox pre-perihelion and the insolation increases, then Hapi is mainly illuminated with low incident angles while the Sun shines parallel to the valley. Over most of the rotation the area is shadowed either by the small or the big lobe. Integrated over the whole orbit Hapi is one of the least insolated areas of 67P’s surface (see Fig. 9). It helps little that self-illumination within the cavity is strong (see Fig. 3). The comparison of the total insolation with and without self-illumination shows the strongest relative effect on the large cliff of the Hathor area. The self-illumination including IR radiation roughly increases the energy input by about 50% (compare Fig. 9), but even then the received energy reaches only about 70% of the peak energy on the north facing areas. Earlier calculations based on a less detailed shape model produced qualitatively similar results but overestimated the energy input in the Hathor and Hapi regions (Sierks et al. 2015).

Dust activity. At the beginning of the Rosetta rendezvous dust features were first detected by the OSIRIS cameras when the line of sight was integrating along the neck valley. This does not mean that sublimation or even particularly strong sublimation occurred only there. MIRO (Gulkis et al. 2015) and the in situ mass spectrometer ROSINA (Hässig et al. 2015) detected water long before the first dust activity was seen. Analysis of the measured water densities shows that the water production originates from the insolated surfaces (Bieler et al. 2015) and that the rates are roughly proportional to the strength of the insolation. This suggests that the activity of the observed nucleus surface is rather homogeneous at least on a large scale. This conclusion is also suggested by the observed dust activity in December 2014 and January 2015. All illuminated parts of the (northern) areas are active and produce dust features. The structures within the near surface dust probably reflect the morphological variation of the surface with crater-like features and cliffs. The comparison of the small images of Fig. 13 demonstrates the effect of integration along the neck valley where the dust features resemble a fan.

Local activity spots. The resolution of the SHAP4 model is good enough to single out “hot spots” with locally strongly enhanced erosion. One striking activity spot is located in the Ma’at area on the small lobe clearly showing up on the diurnal erosion calculations for 9 August and 12 September 2014 (Figs. 11 and 15). Vincent et al. (2015) discussed dust activity in connection with pits and investigated some specific locations on the surface of 67P. Small jets seemingly arise from the Ma’at area. We calculated the erosion for the date (12 September 2014) of their image applying an enhanced (less aggressively decimated) shape model. We find a spot of strong activity on a cliff in front of their Ma’at pit 1 (Fig. 15). One should keep in mind that the model calculations show erosion where insolation is ample assuming a homogeneous surface. Our model calculations combined with a high resolution shape provide a robust tool to localize active spots, even tiny ones within craters or pits. This combination will indicate where to look for changes of morphological details caused by the erosion of surface features.

6. Erosion and morphology

The intense southern summer around perihelion leads to the strong dichotomy that ? of all eroded material is lost during this season. The least erosion is found at the directly north facing parts of Ma’at on the small lobe and Ash on the big lobe. The erosion by water sublimation (model B) is here typically 10% (and at Hapi even less) of the maximum reached on the south side of the nucleus (Fig. 11). Erosion increases for the areas of Ash and Ma’at that slope towards the south. A gradient is visible on Seth where the part close to the pole is least illuminated because it is shaded by the small and big lobes during the diurnal spin. Seen from the north, almost all of the visible surface experiences little erosion by water sublimation (blue and deep blue color in Fig. 11). Only the periphery on the small lobe receives somewhat more insolation. Generally the edge of this area is so steep that surface areas illuminated when the Sun approaches the equator are not visible from the north. In particular for the big lobe one has the impression of looking down on a slightly overhanging “roof”. The transition from the “roof” to the well eroding equatorial areas is abrupt and forms a distinct ridge particularly at the transition from Seth to Babi/Aten but also to the other side to Anubis (Fig. 16). Here the “roof” of Seth even juts out as demonstrated by the increased erosion along the underside of the overhang, a phenomenon that has been interpreted as an indication of multiple layers (Massironi et al. 2015). The view from the south is similar. The preliminary shape model indicates that the south facing part of the nucleus shape is rather flat. Equatorial terrains bulge only slightly out. The exception is the big relatively flat Imhotep plain.

Ring-like or horseshoe-like indentations with diameters on the order of several 100 m with steep cliff-like walls cover the Seth surface. The indentations are flat (not bowl-shaped) and horizontal relative to local gravity whereas the cliffs are vertical. They are arranged like steps leading upslope in the Seth region. Here, essentially all horizontal areas are covered by what seems to be air fallen dust. Most probably the dust is liberated during the perihelion passages on the south and settles on the then completely dark and inert northern regions. The present dust layer seems thick enough (on the order of meters; Thomas et al. 2015a) to quench any activity from the original nucleus surfaces. On the side of Seth that slopes southward towards Anubis the indentations open quite often towards the southeast direction, towards the direction from which the strongest insolation comes. The cliffs facing in this direction show by far the strongest erosion potential (Fig. 16) of the Seth area, reaching about 60% of the maximum in the south. Deep pit-like holes are frequent (Thomas et al. 2015b; Vincent et al. 2015). The morphology of the north facing parts of Ma’at resembles that of Seth. The predominant sublimation (erosion) of the cliffs will lead to surface flattening due to sublimation-driven slope retreat and this flattening could potentially be the process responsible for planation of the equatorial and southern surfaces, for example, Imhotep (see below).

The Seth area is reminiscent of the extremely rough surface of comet 81P/Wild 2 (Brownlee et al. 2004). 81P’s perihelion distance recently changed from 4.9 to 1.5 AU by a close Jupiter encounter in 1974 (Sekanina & Yeomans 1985). The surface of comet 81P is covered with flat-floored craters very similar to those found on Seth (see Fig. 7 in Brownlee et al. 2004). This comet’s morphology was discussed by Basilevsky & Keller (2006). The slope distributions of Seth and of 81P’s surface seem both quite extreme, dominated by slopes near zero and near vertical. This similarity could suggest that both surfaces are relatively pristine if compared to stronger eroded surfaces and comets such as 9P/Tempel 1 (A’Hearn et al. 2005). This, however, would imply that the orientation of the nucleus of 67P has been unfavorable to good insolation of the (now) northern areas for many orbits. It would take at least 10 perihelion passages with inverted south-north insolation to smooth the present surface of Seth. Several cliffs are higher than 100 m.

Imhotep. We do not yet know the morphology of the south side of the nucleus because it is still unilluminated. The large flat surface almost at the end of the large lobe, Imhotep, extends across the equator and is therefore illuminated during most of the orbit of 67P. The angle between the rotation axis and the plain of Imhotep is about 50° similar to the sub-solar latitude at northern solstice. During that season Imhotep sees barely any light. Nevertheless, the erosion potential of this 10 km² large area reaches about 70 to 80% of the maximum found on the southern regions. Again the cliffs in the north of Imhotep undergo stronger erosion, reaching up to 90%. This is a good example of how the current insolation tends to flatten the surface. Figure 17 shows the daily (integrated over one cometary spin) erosion at Imhotep at three different seasons. On the left Imhotep awakens from polar night when the Sun starts to touch the highest mountains. The center model depicts the strong erosion around summer solstice right after perihelion. Watch the scales of the color bars. The maximum here is 180 times larger than on the left side. The large cliff at the north (indicated by arrows) shows the strongest erosion. This is not anymore the case when the Sun stays above the equator and illuminates the surface uniformly albeit with much lower intensity because the distance to the Sun is already 2.6 AU. Then the steep slope of the cliff is hit by the rays with a large incident angle at local noon, leading to less absorbed energy per surface area. This interplay of the varying insolation drives the cliffs across the Imhotep plain from the south to the north and finally leads to planation of the surface.

On the opposite side of the nucleus is Hatmehit, also a flat surface on the equator, though rather circular with a diameter of about 800 m. Its orientation to the rotation axis and hence the irradiation situation are very similar to that of the Imhotep plain. At the present stage where we mostly see the cometary surfaces facing north, the Imhotep region is interesting and peculiar for several reasons. This region, studied in detail by Auger et al. (2015), is located on the equator and will soon be strongly illuminated, i.e. one can expect that it is similar to strongly insolated southern regions that are still in polar night. It comprises a large (3000 m × 3500 m) flat horizontal plain that is slightly inclined towards the south. Its northern rim is limited by a low mountain range. The region shows a large diversity of terrains (smooth or rocky), specific areas (accumulation basins and bright spots) and remarkable morphological features (linear features, roundish features, and boulders (Pajola et al. 2015). The relatively strong insolation (compared to the presently visible surfaces) integrated over the orbit and the corresponding erosion potential lead to morphological diversification. This makes Imhotep a role model for the unknown southern territories. Imhotep is a good candidate for being a region that has been shaped by a combination of several processes driven by relatively strong insolation including collapse, air falls, degradation, thermal stress and most important activity (Auger et al. 2015). The most striking evidence of recent activity processes is likely the presence of roundish features, located in the gravitational low of the region, depleted in boulders and in fine material deposits contrary to the rest of the visible nucleus. Although their origin is uncertain, these roundish features are likely the remnants of recent strong activity. Another evidence of recent activity is the presence of fields of bright boulders, bluer than the surrounding terrains and composed of dust and water ice according to Pommerol et al. (2015). Although probably not important in terms of water production rate, these bright boulders formed recently, most likely at the previous perihelion passage, suggesting that activity indeed occurs today on Imhotep.

Constrained activity of a homogeneous nucleus. Constrained areas of notably high activity cannot be maintained over many orbital revolutions of the comet. After 10 orbits the active spots will lie in depressions of 150 to 200 m depth (model A). Sunlight would reach the active bottoms of the pits only during shorter and shorter intervals and hence activity would become marginal. Constrained activity is digging its own grave. A way out of this conundrum is for the depressions (pits) to expand laterally by activity of their cliff walls. Activity of the cliffs can be more easily maintained because there it is not quenched by remnant or “air fallen” dust. In this way, activity can erode areas that cannot maintain activity on horizontal plains. Consequently, activity will be found predominantly at cliffs of evolved cometary nuclei. This is in good agreement with prior observations of cometary nuclei such as 19P/Borrelly (Britt et al. 2004), 81P/Wild 2 (Sekanina 2004), and 9P/Tempel 1 (Farnham et al. 2007, 2013), for which a correlation of activity with steep slopes was observed. OSIRIS observations of dust features above the northern surfaces strongly suggest a correlation with the ubiquitous cliffs (Vincent et al. 2015). The total surface represented by cliffs and their irradiation are sufficient to support the currently observed activity as we discussed earlier.

Consequently, this scenario, where only a small part of the surface is active, does not mean that most of the comet is inert and cannot support activity due to lack of volatiles. A thin, accumulated, surface layer of refractory material is enough to prevent direct activity. Activity on an evolved surface requires that the active spots can wander over most of the surface. The nucleus has to be homogeneous (supporting activity) otherwise activity will decay very fast. The surface is eroded, layer by layer, when the cliffs wander across the plains. An example of an early phase of this erosion process is Seth, a later stage can be observed on Imhotep.

Fig. 16 Upper image: the nucleus as seen in the equatorial plane. Lower image: the nucleus as seen from the north almost parallel to the rotation axis pointing out of the plane towards the observer (near the red bullet). This view shows also the neck, which is very little eroded by water sublimation.

Fig. 17 Erosion per spin (model B) on the large Imhotep plain when the sub-solar latitude is 48° (left image, May 2014), at southern solstice (center image, September 2015), and at equinox (right image, March 2016). The view is in the equatorial plane. The rotation axis is inclined by about 50° relative to the Imhotep plain. At southern solstice the cliff (indicated) is insolated and eroded much stronger than at equinox, when the Sun is high but farther away.

7. Global versus localized activity

So far we have assumed that the cometary surface is homogeneous with respect to activity and erosion. This seems to be reasonable − at least at large scale − as suggested by the overall dust activity and the water molecule distribution in the inner coma. This does not imply that activity is based on water sublimation only. The ROSINA in situ observations show that the densities of super volatiles (CO, CO₂) relative to the dominating water vary by factors with time and/or position of the Rosetta spacecraft (Hässig et al. 2015). Compositional variability was also observed during the flybys of comets 9P/Tempel 2 (Feaga et al. 2007) and 103P/Hartley 2 (Protopapa et al. 2014). However, even at the present low activity phase the water sublimation models show an overall activity that is about 16 (model A) to 4 (model C) times higher than observed. There are three possibilities for reducing the overall sublimation: a) the surface is uniformly covered by an even thicker dust layer than assumed for model C, so the activity is further quenched; b) only parts of the surface are active, varying on a relatively small scale like the fields on a chessboard. The portion of the currently insolated surface parts covered by cliffs is about 10 to 15%. Hence the observed activity could come from the cliffs only; and c) surface activity may vary (to some degree) on a large scale. An indication of this situation could be derived from the fact that models A to C predict a shift of the maximum production by only 6 days after perihelion rather than the observed about 20 days. In addition, the observed maximum is more peaked than predicted. We therefore briefly discuss a model developed earlier that assumes dirty water (like model A) distributed in spots on the surface.

We consider three alternative models, called the homogeneous (D) and spotted (E and F) models. These models are less realistic than models A−C in the sense that a much coarser shape model is used (SHAP4 reduced to 5000 facets), self-illumination is ignored, and shadowing is included in an approximate manner. However, they are more realistic than models A to C in the sense that the effect of thermal inertia is considered explicitly (models D, E, and F), and allow for a heterogeneous distribution of ice (activity) on the nucleus (models E and F).

In models D and E the one-dimensional heat conduction equation is solved, (5)with upper and lower boundary conditions (6)Since the specific heat capacity c, mass density ϱ, and heat conductivity κ here are considered temperature-independent, the surface temperature and sublimation rate do not depend on the individual values of c, ϱ, and κ but only on the thermal inertia . We here apply Γ = 30 and 100 [J m^-2 K^-1 s^−1/2] for models E and F, respectively. The upper boundary condition balances direct solar radiation (solar constant S_⊙ = 1367 J m^-2 s^-1, Bond albedo A = 0.032¹, incidence angle i, and heliocentric distance r_h), thermal emission (emissivity ε = 1), sublimation rate of ice (sublimation rate Z, latent heat of sublimation L_ice), and heat conduction (depth x, temperature T). The vanishing temperature gradient used as lower boundary condition is applied at a depth of ~10 skin depths.

The water production rate is calculated for each nucleus facet as a function of time t, accounting for i – variation during daytime and gradual cooling during nighttime, but ignoring transient shadowing. Values are stored at every 10° of nucleus rotation. For these specific nucleus rotational phases, the shadowing conditions are calculated for each facet using ray casting. In case a facet is formally on the day side, but is experiencing transient shadowing, its gas production rate for that specific rotational phase is set to zero. A rotationally averaged production rate is then calculated for each facet, and this is done for the interval 150 days pre-perihelion to 150 days post-perihelion, with a 10 day temporal resolution. We note that models D−F (like model A) formally are valid only if there is exposed water ice on the surface, while all comet nuclei observed so far are covered by dust. However, according to model B, the addition of a dust layer just thick enough to suppress spectral absorption features of the ice, does not have a strong influence on the water production rate near perihelion. Thus, the lack of water absorption features on 67P does not necessarily invalidate Eqs. (1) and (2).

Figure 18 shows the measured water production rate of 67P according to Hanner et al. (1985), Crovisier et al. (2002), Feldman et al. (2004), Schleicher (2006), Ootsubo et al. (2012), Bertaux et al. (2014). It also shows model D as a dashed curve, obtained by allowing all facets of the shape model to be active, but reducing the production rate to 6.9% of its actual calculated value to minimize the residual with respect to observations. For model A the reduction is 6.5% (see Sect. 4), hence both models are well in agreement. Consideration of the thermal history due to the inertia (models D−F) and a somewhat higher albedo value reduce sublimation slightly farther away from the Sun. The additional increase from 30 (E) to 100 [J m^-2 K^-1 s^−1/2] (F) applying the same model leads to a reduction within the overall uncertainties of the model(s) (see Fig. 18). Model D can be seen as a nucleus where 7% of each facet is covered by water ice while 93% is covered by inert dust, and the nucleus is homogeneous in the sense that all facets have the same ice coverage. As can be seen in Fig. 18, the homogeneous models D and A reproduce the near-perihelion rates, but yield pre- and post-perihelion rates that are too high and with too shallow slopes. From this comparison, it is clear that the comet nucleus does not produce water vapor in constant proportion to its daily averaged geometric cross section, as seen from the Sun.

Fig. 18 Water production rates of comet 67P measured during earlier apparitions near perihelion. Models A and D for a homogeneous nucleus are scaled to the observed maximum. Model E is a best fit for spotted activity. Model F is calculated for the same distribution but with a thermal inertia of 100 [J m-2 K-1 s−1/2].

A possible way to change the pre- and post-perihelion slopes is to de-activate some facets completely, while increasing the availability of ice elsewhere. Model E is an attempt to construct such a heterogeneous or “spotted” nucleus. This is done by considering surface activity maps (Davidsson & Gutiérrez 2004). Here, an active spot is constructed by placing the tip of a cone at the nucleus center, selecting its symmetry axis orientation by using randomly generated azimuth and elevation angles, and allow all surface facets located within the cone to be active, while those being outside are considered inactive. For each surface activity map, a random number of spots are created (at most ten), with randomly selected opening angles (at most 15°). For each map, the corresponding water production rate is calculated. In the current simulation, a total of 10⁶ surface activity maps were considered, of which 29 maps were selected based on their small residuals with respect to observations. The average production rates of these 29 maps are shown in Fig. 18 as red squares, where error bars denote the standard deviation among the different maps. As can be seen, such spotted model nuclei are capable of reproducing the observations, i.e. an activity that starts off lower but increases more steeply pre-perihelion than for model D, peaks about 10 days after perihelion, then falls off more rapidly than model D post-perihelion.

In order to investigate what the 29 successful surface activity maps have in common, the probability of activation was calculated for each facet. This probability is shown as a color coding on the shape model in Fig. 19 for four different nucleus rotational phases separated by 90°, as seen from the Sun at perihelion. The region of activity is remarkably confined. One center of activity is found at the interface between the small lobe and the large lobe, on the opposite side of the Hapi region. The large lobe is essentially inactive, except for an extended region that is attached to the first one. Therefore, the success of these maps relies on the large active area being poorly illuminated, except for a ~100 day window centered on the perihelion passage. The low level of activity on the sides and bottom of the large lobe (defining the top of the large lobe as the side containing the Seth and Ash regions) are necessary in order not to provide too much water vapor before and after perihelion.

Fig. 19 Comet 67P represented by a low resolution shape model, as seen from the Sun at perihelion at four different rotational phases separated by 90°. The color indicates the probability that a facet is active.

We note that the region of activity coincides with the parts of the nucleus for which the Sun is permanently high in the sky during perihelion. It cannot be excluded that the nucleus is heterogeneous, and that the amount of water ice simply is higher in the active region than elsewhere. However, there may also be an alternative explanation. If the removal of quenching dust is an inefficient process, the surface layers may clog up easily when illumination conditions are relatively weak. This may be the reason why most of the nucleus appears to be relatively inactive, but in areas where the solar insolation is continuously strong throughout the perihelion passage, insolating dust layers may not have time to grow but are constantly being removed or thinned out. Such efficient removal of dust may keep a constant supply of water ice available at or near the surface within a restricted area, and cause a relatively peaked water production curve around perihelion.

Our prediction for the August 2015 perihelion passage is therefore that the gas production curve of 67P depends less on cross section variation, and more on nucleus inhomogeneity and/or the efficiency of dust removal. Fortunately, the chemically and physically complex process of strong comet activity will be studied in detail by Rosetta.

8. Summary and conclusions

We approximate the complex shape of 67P by a high resolution polyhedron (shape model) and model the insolation of each facet taking into account the self-illumination by other facets in its field of view. In a first step the reciprocal viewing factors of the facets are calculated and stored. Then the surface temperatures of the facets are calculated so that the thermal radiation next to the reflected visible and IR radiation can be iteratively considered. The calculations are fast enough to consider more than 10⁵ facets. The surface temperatures depend on the physical properties of the surface such as bolometric albedo, emissivity, and sublimation activity. We can neglect the thermal history as long as the characteristic time scales for heat conduction are smaller than the time steps. A comparison with a model that includes the effects of thermal history shows that our simplification is reasonable at least near perihelion when activity is observed. We consider three models for the sublimation of water ice, on the surface (dirty ice), under a thin layer of dust, and a thicker layer that quenches the water sublimation rates by a factor of 4 to 5 if compared to the dirty ice model.

During the ten months near perihelion, when the sub-solar latitude is south the models predict that the amount of water sublimation is three times higher than what is produced during the rest of the orbit. The assumption of a homogeneously active surface leads to production rates that are 4 to 16 times higher than what was observed during earlier perihelion passages of the comet. Consequently only a small part of the nucleus surface is active. We tested various distributions of active spots on the south side of the nucleus and could match the variation of the observed water production rates. If the model production rates are normalized at perihelion, the activity predicted for August 2014 is more than three times what was observed. This suggests that the fraction of active surface on the north side has to be further reduced. The low insolation level restrains activity and cannot activate as much surface as on the strongly insolated side. The gradient of the water production rate is steeper than expected from the models. Nevertheless, activity seems to be globally distributed. In December 2014/January 2015 dust features are observed above all well illuminated regions. The in situ measurements of water also show that the production rates follow the insolated surface. Therefore we conclude that active and inactive areas are distributed on small scales like on a chessboard and most of the currently observed area is inactive. Active spots make for only about 5 to 10% of the surface. This is roughly in agreement with the amount of surface that is occupied by cliffs that are distributed over the whole northern side. Our calculations show that the cliffs are more strongly insolated than the horizontal surfaces, which are mostly covered by refractory material. This activity pattern leads to planation of the surfaces. A good example for this effect is Imhotep, a flat plain bordered in the north by a mountain range that is driven north by the strong sunshine from the south during perihelion.

Just taking water sublimation into account the maximum erosion on the south side could reach one to five times (γ + 1) meters, with γ being the dust to ice ratio (around 4). Thus the activity from the south side will have quite a different quality. During the perihelion passage dust jets are produced that can be observed with telescopes from earth. The northern regions are in polar night during this season. They form a cold trap for water molecules and icy (dust) particles and blocks that do not escape from the nucleus. In addition the region around the pole and the neck between the lobes are in a gravitational low where particles can easily be collected. This suggests that the northern terrains that are less steep (relative to local gravity) than the angle of repose are covered by global “airfall” from the high activity in the south. This agrees with the current observations.

We expect the topography on the south side to be smoother than what is observed on the north. Imhotep could be a role model. Our test calculations also show that the activity is concentrated near the center of the south side. The variation of the water production rate near perihelion is not proportional to the illuminated cross section of the nucleus.

Acknowledgments

OSIRIS was built by a consortium of the Max-Planck-Institut für Sonnensystemforschung, in Göttingen, Germany, CISAS-University of Padova, Italy, the Laboratoire d’Astrophysique de Marseille, France, the Instituto de Astrofísica de Andalucia, CSIC, Granada, Spain, the Research and Scientific Support Department of the European Space Agency, Noordwijk, The Netherlands, the Instituto Nacional de Técnica Aeroespacial, Madrid, Spain, the Universidad Politećhnica de Madrid, Spain, the Department of Physics and Astronomy of Uppsala University, Sweden, and the Institut für Datentechnik und Kommunikationsnetze der Technischen Universität Braunschweig, Germany. The support of the national funding agencies of Germany (DLR), France (CNES), Italy (ASI), Spain (MEC), Sweden (SNSB; Grant No. 74/10:2), and the ESA Technical Directorate is gratefully acknowledged. H. Rickman was also supported by Grant No. 2011/01/B/ST9/05442 of the Polish National Science Center. We thank the ESA teams at ESAC, ESOC and ESTEC for their work in support of the Rosetta mission.

References

All Tables

All Figures

	Fig. 1 Shape and surface of comet 67P as seen from the Sun at three rotational longitudes (30°, 150°, 270°) at northern solstice when the sub-solar latitude is + 52°, at the equinoxes, and at southern solstice. The smooth parts of the surface indicate yet unobserved territory on the south side of the shape.
In the text

	Fig. 2 Observational parameters during the present orbit of 67P. The left figure shows heliocentric distance of the comet and solar latitude vs. time. The right figure demonstrates the asymmetry of the insolation as a function of heliocentric distance.
In the text

	Fig. 3 Viewing factors of the model facets. A view factor of 0 (black) implies that a facet is not exposed to irradiation from other facets. A view factor of 1 implies that the field of view of a facet is completely filled by other facets of the body.
In the text

	Fig. 4 Surface temperature as function of the irradiance for a gray body and models A to C. Model A represents the lower limit for the temperature of an active surface. The temperatures of model C come close to that of an inactive gray body.
In the text

	Fig. 5 Water production rates vs. time for models A to C with (black lines) and without (red lines) self-illumination (SI).
In the text

	Fig. 6 Erosion and insolation integrated over one orbit. The left figure shows the erosion (model A). The maximum water erosion on the (unseen) south side is 4.8 [m]. The right figure shows the insolation on the same color scale with the maximum of 7.54 × 109 [J m-2].
In the text

	Fig. 7 Production rate for model A. Dashed lines indicate the times of perihelion and maximum production.
In the text

	Fig. 8 Water production rates of models A to C versus heliocentric distance.
In the text

	Fig. 9 Effect of self-illumination. The left image shows the energy absorbed from direct sunlight integrated over the whole orbit. For the right image the energy reflected and radiated (thermal) from other facets in the field of view of each facet is included.
In the text

	Fig. 10 Water erosion (model A) integrated over one cometary spin on 9 August 2014. The northern hemisphere with the neck area is shown. The erosion is enhanced along the transition from the valley bottom (Hapi) to the cliff wall (Hathor).
In the text

	Fig. 11 Three different views of the shape showing total (integrated over one orbit) water erosion (model B). Left: northern areas and neck, middle: equatorial view, right: southern, yet unobserved side.
In the text

	Fig. 12 Daily variation of water production (model A) vs. rotational phase. The dates are 9 August 2014, when jets started to be observed at the neck area, 14 November 2014 (landing of Philae), 31 January 2015, 2 April 2015, 9 May 2015 (equinox), and 9 August 2015 (perihelion).
In the text

Fig. 13 Dust activity of 67P. The left image is a combination of two WAC images taken on 10 January 2015. The dominating fan-like dust feature above the neck region is strongly enhanced by integration along the line of sight, parallel to the fan plane from the active neck valley. This effect is demonstrated by the small images on the right. The lower image was taken with a similar line of sight as the left image, whereas the nucleus is slightly rotated in the upper view taken half an hour earlier. The enhancement by the integration effect does not appear.

In the text

Fig. 14 A comparison of the neck area (Hapi with Hathor and Seth). The model shows erosion (model A) integrated over one spin on 12 September 2014. The strongest erosion is on a cliff in the Ma’at area (see Fig. 15) The erosion along Hapi is concentrated at the foot of the Hathor wall. The enhancement is visible as a yellow band. The surface temperature (based on model C) is as expected also high in that region, depending on the local time of day.

In the text

	Fig. 15 Activity in the Ma’at area. In the left figure (Vincent et al. 2015) the pits 1 and 2 were identified as exhibiting dust activity. Indeed, our model for 12 September 2014 (right figure) reveals small areas of strong erosion in both pits. Below Ma’at 1 there is a cliff that shows the strongest erosion of the entire illuminated surface (compare with Fig. 10).
In the text

	Fig. 16 Upper image: the nucleus as seen in the equatorial plane. Lower image: the nucleus as seen from the north almost parallel to the rotation axis pointing out of the plane towards the observer (near the red bullet). This view shows also the neck, which is very little eroded by water sublimation.
In the text

Fig. 17 Erosion per spin (model B) on the large Imhotep plain when the sub-solar latitude is 48° (left image, May 2014), at southern solstice (center image, September 2015), and at equinox (right image, March 2016). The view is in the equatorial plane. The rotation axis is inclined by about 50° relative to the Imhotep plain. At southern solstice the cliff (indicated) is insolated and eroded much stronger than at equinox, when the Sun is high but farther away.

In the text

	Fig. 18 Water production rates of comet 67P measured during earlier apparitions near perihelion. Models A and D for a homogeneous nucleus are scaled to the observed maximum. Model E is a best fit for spotted activity. Model F is calculated for the same distribution but with a thermal inertia of 100 [J m-2 K-1 s−1/2].
In the text

	Fig. 19 Comet 67P represented by a low resolution shape model, as seen from the Sun at perihelion at four different rotational phases separated by 90°. The color indicates the probability that a facet is active.
In the text