© S. Benseguane et al. 2022

1 Introduction

Comets are among the least processed remnants of the early stages of our planetary system. The study of comets thus provides critical information to help us better understand the physical processes that lead to planet formation, and the early material that formed the protoplanets (Festou et al. 2004; Cochran et al. 2015). Jupiter-family comets (JFCs) are a subpopulation of comets, with short-period orbits dominated by the gravitational influence of Jupiter. They are thought to originate from the Kuiper Belt and the scattered disk (Brasser & Morbidelli 2013), where they got destabilized owing to interactions with Neptune. They evolved through the giant-planet region and toward the inner solar system where they are observed nowadays (Levison & Duncan 1997; Di Sisto et al. 2009; Nesvorný et al. 2017). Cometary activity starts beyond the orbits of Jupiter and Saturn for long-period comets (Meech et al. 2017; Jewitt et al. 2017; Hui et al. 2017, 2019; Yang et al. 2021; Farnham et al. 2021) and for Centaurs, which are the precursors of JFCs (Jewitt 2009; Lin et al. 2014; Epifani et al. 2017, 2018; Steckloff et al. 2020; de la Fuente Marcos et al. 2021). This implies that JFCs, which have been studied so far by space missions, mostly have evolved surfaces (e.g., Gkotsinas et al. 2022). In this framework, the European Space Agency (ESA)’s Rosetta mission aimed to study how a comet’s surface might be modified through cometary activity. Indeed, by understanding the physical processes that currently reshape the nucleus, we might reconstruct the properties that it would have had at the time of its formation.

Significant geological heterogeneity was observed at the surface of 67P/Churyumov-Gerasimenko (hereafter 67P). In addition to the presence of terraces, strata, fractures (Massironi et al. 2015), goose-bump features (Sierks et al. 2015), and wind-taillike features (El-Maarry et al. 2019), the observation of surface depressions, linked with cometary activity, thus offered the opportunity to look into the characteristics of the subsurface and the thermophysical processes actively modifying them (Sierks et al. 2015; Vincent et al. 2015). Two main types of depressions can be distinguished on the surface of 67P based on their dimensions, that is to say their diameter and their depth (mainly the latter): shallow depressions and deep depressions, as we detail below.

First, shallow depressions of only a few meters in depth are generally observed on smooth terrains, in many regions across the nucleus. These might be seasonal depressions shaped during perihelion passages, reported to be mainly driven by sublimation activity in the current orbits of the comet. For instance, Vincent et al. (2016) and El-Maarry et al. (2017) propose that the surface was reshaped via scarp retreat, and Pajola et al. (2016) infer that these shallow structures could indicate a future cliff collapse. Groussin et al. (2015) and Bouquety et al. (2021b) suggest that they could be seasonal structures shaped by progressive erosion, induced by activity sustained close to the perihelion approach. Bouquety et al. (2021a) named these depressions cometary thermokarst depressions due to their morphometrical analogy with thermokarstic lakes on Earth and scalloped terrain on Mars. These depressions will not be further studied in this work.

In this study, we are interested in the second kind of surface depressions, characterized by steep walls with depths of tens to hundreds of meters (Massironi et al. 2014; El-Maarry et al. 2015, 2019; Thomas et al. 2015b). These larger-scale structures include pits, as well as cliffs or alcoves (see Fig. 1 retrieved from Rosetta/OSIRIS’; NAC Narrow Angle Camera). They are mostly present on 67P’s northern hemisphere and are generally concentrated in some regions. For instance, the Maftet geological unit displays irregular-shaped pits of 1020 m deep and 100–150 m in diameter (Thomas et al. 2015a). The Seth region is dominated by multiple series of circular, flat-floored pits (Besse et al. 2015) and contains a pit chain similar to the one observed on Ma’at (Thomas et al. 2015a, see Fig. 1). It also contains cliffs that are tens to hundreds of meters high, which are also observed in Hathor (El-Maarry et al. 2015). Furthermore, Vincent et al. (2015) report the detection of cometary activity in the form of localized dust jets in some of these pits. Additionally, they note that active pits have a high depth-to-diameter ratio compared to inactive ones (Besse et al. 2015). Our study thus specifically considers the following pits: Seth_01, Seth_02, Seth_03, Seth_04, Seth_05, Seth_06, Ma’at_01, Ma’at_02, Ash_03, Ash_04, Ash_05, and Ash_06 (see Fig. 1 from Vincent et al. 2015). We also include additional pits, with similar geomorphological characteristics to the ones studied in (Vincent et al. 2015). We further study cliffs or alcoves, as these might be construed as deteriorated pits (Vincent et al. 2015).

As a result, our study focuses on features of at least a few tens of meters in depth, and a few hundreds of meters in diameter: the smallest depth and diameter are 35 and 130 m, respectively. In the rest of the paper, we indifferently use the term “pit” for the sake of simplicity.

Such pits have been observed on most comets directly studied by space missions, for example 19P/Borelly (seen by Deep Space 1, Soderblom et al. 2002), 81P/Wild 2 (seen by Stardust, Brownlee et al. 2004), 9P/Tempel 1 (seen by Deep Impact and Stardust-NExT, Belton et al. 2013), and 103P/Hartley 2 (seen by EPOXI, Syal et al. 2013). The mechanism at the origin of these structures is still a matter of debate (Brownlee et al. 2004; Belton & Melosh 2009; Belton et al. 2013; Thomas et al. 2013). Holsapple & Housen (2007) and Vincent et al. (2015) argue that impacts on cometary surfaces are expected to produce features with a morphology distinct from these observed pits, and thus they should be a signature of some process related to cometary activity rather than the result of collisions. Vincent et al. (2015), Kossacki & Czechowski (2018) and Leliwa-Kopystynski (2018) propose the formation of pits by sinkhole collapse due to subsurface cavities, either primitive or formed as a result of subsurface depletion of volatiles by ice sublimation. Massironi et al. (2014) argue that a sublimation process can lead to slope retreats and material ablation at the pit’s location, while Thomas et al. (2015b) argue that mechanisms such as ice sublimation or sinkhole collapse would not likely lead to the material structure giving rise to the quasi-circular aspect of pits.

Mousis et al. (2015) explored the possibility such structures forming due to phase transitions (i.e., sublimation, amorphous water ice crystallization, and clathrate destabilization). They showed that the time required to produce features of the spatial scale observed by Rosetta on the surface of 67P is long, on the order of a thousand years or more. Guilbert-Lepoutre et al. (2016) further showed that it is very unlikely that pits form with the current illumination conditions, as no known mechanism could carve the surface to form pits with a depth of ~200 m, and a diameter ranging from 100 to 300 m over short timescales. Moreover, since 67P’s previous orbits had perihelion distances farther away from the Sun (Maquet 2015), it is unlikely that quasi-circular pits were formed by the progressive effect of one phase transition. Because the distribution law of the pits’ size frequency is similar at the surface of 67P, 9P/Tempel 1, and 81P/Wild 2, Ip et al. (2016) suggested that they might have been formed with the same mechanism operating on many JFCs. Comparing the orbital history of those comets, they inferred that such processes might have carved pits before these comets entered the inner solar system.

With these arguments in mind, we want to understand how the progressive modification of 67P’s surface due to cometary activity might have affected the characteristics of these depressions. In particular, we are interested in understanding whether signatures of the formation mechanism at the origin of pits can still be found. Our goal is thus to quantify the amount of erosion sustained by pits at the surface of 67P, under the current illumination conditions that periodic, daily, and seasonal cycles entail. In Sect. 2 we present the surface energy model and the thermophysical evolution model used to address this issue. We present the results of the influence of several key parameters in our study in Sect. 3, and the results of the thermal simulations in Sect. 4. Finally, we discuss our results in Sect. 5.

Fig. 1 Image from OSIRIS/NAC of a part of the Seth region on which we illustrate the type of depressions we study: pits and alcoves (half circular-pits, Vincent et al. 2015).

2 Methods

Our study follows the work of Mousis et al. (2015) and Guilbert-Lepoutre et al. (2016), who constrained the thermal evolution of 67P’s subsurface over the recent past. They studied the possible formation of pits in general, averaging the energy input across 67P’s surface. To go beyond these studies, our goal is to quantify how the energy received locally, at a small scale on the surface, translates into phase transitions. We aim to quantify the extent of these phase transitions over multiple perihelion passages. Eventually, we want to assess whether these can be the origin of the formation of pits, or their evolution into structures with the spatial scale observed by Rosetta. To do so, we used a Stereo-PhotoGrammetric (SPG) shape model of 67P’s nucleus (Preusker et al. 2017, see Sect. 2.2.1 for details). The spatial resolution was chosen to provide several facets for each geometric portion of a pit (i.e., the bottom, the cliffs, and the plateau surrounding it). The energy received by each facet was computed and used as the surface boundary condition of a thermal evolution model. Each step of this method is detailed in the following sections.

2.1 Thermophysical evolution model

2.1.1 Main equations

A thermophysical evolution model was applied to each facet of this SPG model. The following aspects needed to be taken into consideration when choosing our numerical scheme:

Each facet gets its own boundary condition at the surface. Therefore, a 1D thermal evolution model is the best option, further justified by the results of Macher et al. (2019), who found that temperature differences at the surface of 67P between a 1D and a 3D thermal simulation (i.e., accounting for lateral heat fluxes) amount to only ~0.1%.

Physical processes included in the model should be standard to thermal evolution models developed over the past few decades (Prialnik et al. 2004; Huebner et al. 2006): heat and gas diffusion, phase transitions for volatile species, drag of dust particles by the escaping vapor phase, and the formation of a dust mantle at the surface.

Finally, thermal evolution models necessarily rely on a number of free or poorly constrained thermophysical characteristics. Exploring the free-parameter-space ought to be rapid from a computational point of view, so as to provide an insight into the robustness of our results. Therefore, simple expressions of thermophysical characteristics such as the thermal conductivity for example are preferred to complex ones, as these would introduce additional parameters.

With these considerations in mind, we chose the 1D scheme as the basis of multiple models. We refer the reader to the works by (De Sanctis et al. 2005, 2010; Lasue et al. 2008) for details on the model, but we provide the main equations below for clarity. More sophisticated models exist to study the thermal evolution of cometary nuclei, considering several dimensions (Guilbert-Lepoutre et al. 2016), and refined descriptions of each thermophysical parameter (Davidsson 2021). However, because our purpose is to understand in detail the influence of the energy input, modulated by local topography and the global morphology of 67P’s nucleus, we need to keep our thermal evolution model relatively simple. Our model thus solves the heat diffusion equation: $${\rho }_{\text{bulk}}c\frac{\partial T}{\partial t}=\text{div}\left(\kappa \overrightarrow{\text{grad}}T\right)+S,$$(1)

where ρ_bulk [kg m⁻³] is the material’s bulk density, c [J kg⁻¹ K⁻¹] its heat capacity, κ [W m⁻¹ K⁻¹] its thermal conductivity, and S = Q_cr + Σ_α Q_α the energy sources and sinks. For our study, we take into account two such heat sources and sinks. The first one is the energy released upon crystallization of amorphous water ice, assuming it is exothermic: $$Q_{\text{cr}}=\lambda \left(T\right){\varrho }_{\text{am}}\text{Δ}{H}_{\text{ac}},$$(2)

where ϱ_am [kg m⁻³] is the mass of amorphous water ice per unit volume. The phase transition releases a latent heat ∆H_ac = 9×10⁴ J kg⁻¹ (Klinger 1981) at a rate of λ(T) = 1.05 × 10¹³ e^−5370/T s⁻¹, determined by Schmitt et al. (1989). The second is the energy loss (or gain) due to sublimation (or recondensation) of different ices present in the solid material. We assume a simple composition of dust and ice, with H₂O, CO, and CO₂ present as pure compounds in the initial icy matrix. For each ice, we have: $$Q_{\alpha }=-\psi \text{Δ}{H}_{\alpha }{q}_{\alpha },$$(3)

where ψ is the porosity; ∆H_α [J kg⁻¹] is the latent heat of sublimation of species α (H₂O, CO, or CO₂); and q_α is the related gas source term, which is obtained through mass conservation equations.

Assuming that sublimation of amorphous water ice is negligible, because the phase transition to crystalline water ice occurs first at lower temperatures, the set of mass balance equations may be written as: $$\frac{\partial {\varrho }_{\text{am}}}{\partial t}=-\lambda \left(T\right){\varrho }_{\text{am}},$$(4) $$\frac{\partial {\varrho }_{\text{cr}}}{\partial t}=\lambda \left(T\right){\varrho }_{\text{am}}-q_{{\text{H}}_2\text{O}}$$(5) $$\frac{\partial {\tilde{\varrho }}_{\alpha }}{\partial t}+\text{div}\overrightarrow{{\varphi }_{\alpha }}=q_{\alpha },$$(6)

where ϱ_am and ϱ_cr [kg m⁻³] are the mass per unit volume of amorphous and crystalline water ice, respectively, and ${\tilde{\varrho }}_{\alpha }$ [kg m⁻³] is the mass per unit volume of each gas species. We assume that the vapor and the solid phases are in local thermodynamic equilibrium, and the vapor phase behaves as an ideal gas (i.e., no interaction between species). Each gas flux $\overrightarrow{{\varphi }_{\alpha }}$ can thus be written separately, as: $$\overrightarrow{{\varphi }_{\alpha }}=-G_{\alpha }\overrightarrow{\text{grad}}{P}_{\alpha },$$(7)

with P_α being the partial pressure of each species, and G_α being a gas diffusion coefficient that generally depends on the structural parameters of the solid matrix (such as the porosity, the size of pores, or the tortuosity), and the temperature (see Prialnik et al. 2004 or Huebner et al. 2006 for details). For each volatile species, the gas source term can thus be written as: $$q_{\alpha }=\frac{1}{RT}\frac{\partial P_{\alpha }}{\partial t}-\text{div}\left(G_{\alpha }\overrightarrow{\text{grad}}{P}_{\alpha }\right),$$(8)

where R is the ideal gas constant.

2.1.2 Initial parameters

The composition and internal structure of cometary nuclei are generally poorly known. The Rosetta mission has, however, provided some crucial measurements for 67P, which are used as constraints in our model whenever possible. All the parameters included in the thermal evolution model have, nonetheless, not been measured, and we therefore make standard assumptions regarding the values of the unknown parameters (see Table 1, and Huebner et al. 2006 for details).

We note that some thermophysical parameters depend on one another. For example, the bulk density of a cometary nucleus can be written as: $${\rho }_{\text{bulk}}=\left(1-\psi \right){\left({\displaystyle \sum _i\frac{X_i}{{\rho }_i}}\right)}^{-1},$$(9)

where ψ is the porosity, X_i is the mass fraction of each individual component in the cometary material mixture, and ρ_i [kg m⁻³] is the corresponding solid density of each constituent. If we use the bulk density of 533 ± 6 kg m⁻³ measured by Rosetta for 67P’s nucleus (Pätzold et al. 2016), a degeneracy remains between the composition and the porosity to obtain this value. Moreover, we study processes affecting the ~100 m-surface layer, which might not have the same properties as the bulk of the nucleus deep inside. In our model, we thus chose the porosity and dust-to-ice mass ratio, which gives an associated bulk density on the same order of the observed bulk density. The composition and porosity also influence the value of thermal characteristics such as the thermal conductivity (i.e., $\kappa =f_{\psi }{f}_{\text{H}}\frac{{\displaystyle {\sum }_i{M}_i{\kappa }_i}}{{\displaystyle {\sum }_i{M}_i}}$, where f_ψ and f_H are respectively the correction factors to account for the porosity and the reduced contact between solid grains, also known as the Hertz factor, M_i is the mass per unit volume of each constituent i, and κ_i is their respective thermal conductivity) or the heat capacity (i.e., $c=\frac{{\displaystyle {\sum }_i{M}_i{c}_i}}{{\displaystyle {\sum }_i{M}_i}}$, where M_i is the mass per unit volume of each constituent i, and c_i is their respective heat capacity). Thus, we test several values of the most crucial characteristics, in order to assess their influence on the outcome of thermal evolution simulations. These parameters are: the initial porosity of the cometary material (Sect. 3.1); the dust-to-ice mass ratio (Sect. 3.2); the abundances of CO and CO₂ (Sect. 3.3); and the thickness of the dust mantle at the surface (Sect. 3.4).

2.2 Boundary conditions

2.2.1 Shape model

A shape model of 67P’s nucleus was reconstructed using the SPG technique on high-resolution images taken by the Rosetta/OSIRIS instrument (Preusker et al. 2015). The latest SHAP7 SPG shape model reaches a very high spatial resolution, with 44 million facets, reconstructed from 1500 OSIRIS’ NAC images (Preusker et al. 2017)¹. From this model, several lower-resolution models were derived: in this study, we use an SPG shape model composed of 124,938 facets². With this model, the typical average distance between two nodes of a facet is ~20 m. Local topography and roughness at smaller scale thus cannot be accounted for in our work.

2.2.2 Energy balance at the surface

The energy and mass conservation equations need to be constrained by boundary conditions. The energy equilibrium boundary condition at the surface is given for each facet by: $$\left(1-A_R\right)ℰ=\epsilon \sigma T^4+\kappa \frac{\partial T}{\partial r}+{\displaystyle \sum _{\alpha }{f}_{\alpha }\text{Δ}{H}_{\alpha }{Q}_{\alpha },}$$(10)

where A_R is the Bond albedo of the facet for which we compute the energy balance, ε is the emissivity, σ is the Stefan–Boltzmann constant, and T [K] is the surface equilibrium temperature. We allow for the presence of volatile species at the surface, so that sublimation is possible: f_α represents the fraction of the facet’s surface covered by these ices, and Q_α [kg m⁻² s⁻¹] is the corresponding sublimation rate. Finally, ℰ = E_⊙ + E_ir + E_VIS is the total energy flux received by a given facet, which takes into account the contributions (detailed below) from direct insolation E_⊙, and hence shadowing effects due to the complex global morphology of 67P’s nucleus, and from self-heating E_IR + E_VIS, namely the energy flux received by reflection and emission from neighboring facets in the visible and infrared, respectively.

Direct insolation is given by: $$E_{\odot }=\frac{F_{\odot }}{r_H^2}\cos \xi ,$$(11)

where F_⊙ [W m⁻²] is the solar flux at 1 au. The heliocentric distance r_H [au] and the local zenith angle ξ both vary with time. For each time step, we first retrieve the coordinates of the subsolar point using SPICE database kernels, which contain the information on both the rotation state of 67P’s nucleus and its orbital parameters. Then, the insolation geometry for each facet is computed with respect to the subsolar point’s coordinates. For facets located on the night side of the nucleus, we apply the following criterion: if cos ξ < 0 then E_⊙ = 0. To assess which facets are located in the shadow of global or local topographic features, we project the nodes of the shape model on a 2D plane normal to the zenith direction of the subsolar point. For each node, we compute its projected position along the normal direction, and test whether it is below an other facet: this node is then considered shadowed. When one of the three nodes of a facet is shadowed, we consider that the whole facet is shadowed.

Given the complex morphology of 67P’s nucleus, observed both on a global scale (two lobes) and on a local scale (e.g., El-Maarry et al. 2015), self-heating – the energy flux received by reflection or emission from neighboring facets – might be a significant additional source of energy. It is composed of two contributions, one from visible radiations reflected by mutually facing facets, E_VIS, and one from their thermal infrared emissions, E_IR. We note that the contribution from infrared radiations reflected by mutually facing facets is not taken into account, because it is always negligible compared to the other two self-heating contributions. The relative influence of self-heating thus depends on how facets from the shape model see each other: the self-viewing geometry is a function of the orientation of the mutual facing facets, both emitting and receiving. For the facet of interest, for which the energy balance is being computed, the visible contribution of self-heating can be written as: $$E_{\text{VIS}}={\displaystyle \sum _{\text{T}}{A}_{\text{T}}\frac{F_{\odot }}{r_H^2}\cos {\xi }_{\text{T}}\frac{S_{\text{T}}}{\pi }\frac{\cos {\zeta }_{\text{T}}\cos {\zeta }_{\text{R}}}{{\delta }_{\text{T}}^2}},$$(12)

where A_T is the Bond albedo of an emitting facet, ξ_T is its local zenith angle, S_T is its surface, ζ_T is the angle between the normal of the transmitter and the receiver facets, ζ_R is the angle between the normal of the receiving and the emitting facets, and δ_T is the distance between the two facets. The infrared contribution can be written as: $$E_{\text{IR}}={\displaystyle \sum _{\text{T}}\epsilon \sigma T_{\text{T}}^4\frac{S_{\text{T}}}{\pi }\frac{\cos {\zeta }_{\text{T}}\cos {\zeta }_{\text{R}}}{{\delta }_{\text{T}}^2}},$$(13)

In this equation, the surface temperature of each emitting facet is approximated using the energy balance: $$\left(1-A_{\text{T}}\right)\frac{F_{\odot }\cos {\xi }_{\text{T}}}{r_{\text{H}}^2}=\epsilon \sigma T_{\text{T}}^4,$$(14)

which accounts for direct insolation only, without any prerequisite knowledge of the importance of the self-heating contributions. When an emitting facet experiences night during a given time step, we set a minimum threshold of T_T = 20 K (various values have been tested and they do not result in any significant variation of the outcomes).

2.2.3 Time step and orbital evolution

Geometric calculations are performed with a cadence of one output every 8 min, over a full orbital revolution of the comet (i.e., ~6.44 yr). As a consequence, the thermal evolution model is run with a time step of 8 minutes, for ten full orbital revolutions. These ten revolutions represent the approximate time 67P has been evolving on its current orbit (Maquet 2015). However, before running the thermal evolution simulations for ten revolutions on the current orbit of 67P, we simulate an injection of the nucleus from the Kuiper Belt to the inner solar system. We use a standard multistage injection process, described by several successive orbits with semimajor axis and eccentricity values given in Table 2. These allow for the slow regression of ice sublimation fronts, and the amorphous-crystalline water ice boundary, below the surface, which mimic the thermal processing sustained by 67P prior to its current orbit (see, for example, Gkotsinas et al. 2022).

2.3 Selection of pits on the nucleus

The diversity of local morphological features at the surface of 67P has been recently reviewed by El-Maarry et al. (2019). Most circular depressions can be found in the northern hemisphere, where deep pits and steep cliffs are also observed. Pits in the southern hemisphere are scarcer, and typically wider and shallower than the ones found in the northern hemisphere. This dichotomy is explained by strong seasons affecting the nucleus, where the southern hemisphere sustains intense heating and erosion during the summer (Keller et al. 2015). For this study, we selected pits with different shapes and dimensions that can be representative of the different illumination conditions at the surface, on both hemispheres and on both lobes, with as much sampling in latitude as possible. We focus on large features, and therefore do not include cometary thermokarst depressions, for example (Bouquety et al. 2021a).

With these constraints in mind, we selected 30 pits: their positions on the surface, as well as morphological characteristics such as their approximate diameter and depth, are given in Table A.1. We note that not all features are circular or quasi-circular pits. Indeed, we also selected elongated pits and alcoves, as well as cliffs, in order to achieve our sampling goals and study possible evolutionary links between those features. We further note that some of the pits selected have shown activity, witnessed by Rosetta/OSIRIS (Vincent et al. 2015). For each pit, facets of the shape model were selected on the surrounding plateaus, the bottom, and the walls for a detailed study of each energy contribution (direct insolation, self-heating, and shadowing), and thermal evolution. One caveat of our method is that we do not account for any shape evolution. Indeed, it is impossible to know what these morphological structures looked like ten orbits ago, as we still do not know how, when, and through what process they were formed. Therefore, the erosion sustained at each time step is not used to modify the geometry of morphological structures. Instead, erosion after ten cometary orbits is assessed from the current shape of 67P’s nucleus, as observed by Rosetta.

All facets selected for our study can be localized on the 3D shape model (Fig. 2 top panel) and a 2D map of the highresolution shape model in an equidistant cylindrical projection (Fig. 2 bottom panel). We note that this map is based on a 12 million facet version of the SHAP7 shape model: we created a “rubber sheet” by putting each vertex point at an elevation proportional to its distance from the comet center above the plane of evenly spaced latitude and longitude. Shading was then realized through 3D rendering. The equidistant cylindrical projection cannot display the overhung areas, but we do not study any such feature here. We note that sophisticated map projections that do display the complete surface of the comet have been presented (e.g., Grieger 2019; Leon-Dasi et al. 2021).

Fig. 2 Display of the facets selected for the study of the 30 pits. Top: the location of the facets on the surface of 67P. The shape model presented is the SPG model composed of 124 938 facets (Preusker et al. 2017), which is used for the surface energy calculation. Bottom: the location of the facets on a 2D map of 67P, which is a projection of the high-resolution SPG shape model composed of 12 million facets.

Fig. 3 Pit selected for the study of the influence of initial parameters, and corresponding energy. A: the location of pit 5, and facets sampled on the plateaus, walls, and bottom. B: the energy received at the 15 facets over one complete orbit, averaged over a daily period window. The gray line marks the perihelion. C: the total quantity of energy integrated over one orbit (left) and the maximum reached during the perihelion (right).

3 Case of one pit: Assessing the influence of initial parameters

In this section, we study the case of one pit to understand the influence of each critical initial parameter on its evolution. The effect of these parameters on the outcomes of our thermal-evolution model will need to be kept in mind when discussing our results. To avoid any interference between parameters, each of them is studied independently of the others. This pit is located in the Seth region and highlighted in panel A of Fig. 3 (label 5 in Table A.1). It is on the big lobe’s northern hemisphere, away from the influence of shadowing by the small lobe, so as to avoid self-heating contributions due to the global shape of the nucleus. Hence, only self-heating due to the local topography of the depression can influence its evolution. Fifteen facets were selected at the bottom, on the walls, and on the plateau surrounding this pit. The energy received on the surface, which includes shadowing and self-heating contributions, is shown in panels B and C of Fig. 3.

Fig. 4 Progressive erosion sustained during ten full revolutions on 67P’s current orbit, for all facets of pit 5, and three values of the porosity: 60 (blue), 70 (red), and 80% (green). Vertical lines and numbers correspond to aphelion passages.

3.1 Influence of porosity

The Rosetta data allowed us to derive values for the internal porosity ranging from 75 to 85% (Hérique et al. 2016, for example). However, some areas of the surface appear to be consolidated material (El-Maarry et al. 2019), with a likely lower local porosity, although direct measurements have not been made. We thus tested three values for this parameter, in order to assess its influence on the outcomes of our thermal evolution simulations: 60, 70 and 80%. From these, we see that a higher porosity, ψ, induces a larger amount of erosion. Indeed, the total mass of eroded material is essentially the same in the different tests, driven by the total amount of energy received locally by each facet. The volume of this corresponding mass varies, however, increasing with an increasing porosity. As such, the extent of the erosion sustained after ten revolutions for ψ = 70% is (on average for all facets) ~30% higher than with ψ = 60%, and the erosion for ψ = 80% is ~50% higher than for ψ = 70% (see Figs. 4 and 5).

3.2 Influence of the dust-to-ice mass ratio

The bulk dust-to-ice mass ratio of a cometary nucleus is notoriously difficult to constrain, especially if it is inferred from the coma composition (Choukroun et al. 2020). For modeling the thermal evolution of comets, a value of one has historically been used, from the Giotto mission measurements at 1P/Halley (see Huebner et al. 2006, and references therein). It appears to be consistent with the Rosetta measurements for 67P (Choukroun et al. 2020). As seen in Figs. 6 and 7, we tested the effects of different dust-to-ice mass ratios, 0.5, one, and two, to assess its influence on the thermal evolution. As for porosity, the initial dust-to-ice mass ratio has a significant influence on the extent of the erosion sustained after ten full orbital revolutions. In fact, we see that erosion substantially increases with an increasing dust-to-ice mass ratio: it almost doubles when we double the ratio. This result is consistent with how the dust-to-ice mass ratio influences the value of thermophysical parameters, such as the thermal conductivity of the material. Indeed, as the dust-to-ice mass ratio increases, so does the thermal conductivity. Energy is thus transferred into deeper layers of the subsurface, leading to the sublimation of deep ice, and thus erosion at relatively higher depths.

Fig. 5 Erosion sustained after ten orbital revolutions for each facet of pit 5, and different values of porosity: 60 (A), 70 (B), and 80% (C).

Fig. 6 Progressive erosion sustained during ten full revolutions on 67P’s current orbit, for all facets of pit 5, and three values of the dust-to-ice mass ratio: 0.5 (blue), one (red), and two (green). Vertical lines correspond to aphelion passages.

3.3 Influence of the CO and CO₂ abundance

The abundance of CO and CO₂ in a cometary nucleus is also difficult to measure. Their bulk abundances are usually assessed from production rates measured in the coma, but the procedure is not straightforward. From modeling experiments (e.g., Prialnik 2006), we know that values integrated over a long period of time are more accurate than data obtained at a single moment on the orbit. Herny et al. (2021) showed that the nucleus of 67P can be considered uniform at the first order. Production rates of CO and CO₂ do vary significantly across the orbit, sometimes by several orders of magnitude (Fougere et al. 2016; Biver et al. 2019; Läuter et al. 2019; Combi et al. 2020). Besides, these volatiles are very sensitive to cumulative heating; hence, the uppermost surface layers are certainly depleted compared to the bulk values.

We tested the effect of the presence of CO and CO₂ by setting their initial mass fraction with respect to water to various values: 0% for both, 1% for both, and 5% CO with 15% CO₂. Adding such large amounts of CO and CO₂ to the initial inventory of volatiles, even if their sublimation boundary regresses below the surface during the multistage injection process, triggers some numerical instability for many facets. Fixing this numerical instability requires that we change the initial thermophysical parameters for these facets, which would defeat our purpose. Therefore, we do not compare the unstable facets further in this analysis (for instance, this is why in Fig. 8 we cannot display results for the 15 facets, as in the prior cases). We do, however, detail in the discussion the potential origin and implication of these facets’ numerical behavior. For facets that do not suffer from numerical instability, the sublimation of water ice remains the main driver for both activity and erosion. Most importantly, we note that all facets do not behave with the same pattern, in a departure from what was observed in previous tests, where all facets would follow the same trend: this effect can be observed in panel A of Fig. 8. Facets that receive the largest amount of energy integrated over one orbit are very active. However, they tend to build a dust layer at the surface after several perihelion passages, which prevents them from being active during subsequent perihelion passages (see panel B1 of Fig. 8). Facets receiving lower amounts of energy, on the other hand, are not active enough to build up such a dust layer. They do get a layer of dust at the surface, but of insufficient thickness to completely quench subsequent activity. Hence, they remain active for the ten full orbital revolutions (see panel B2 of Fig. 8). The dust production rate is, overall, similar for all facets remaining active throughout the ten orbital revolutions, whatever the amount of CO and CO₂ present in the icy phase. Facets for which the activity is quenched, however, stop emitting dust in the coma.

Overall, adding CO and CO₂ to the initial composition leads to a substantial complexity in the thermal simulations (numerical instabilities and unpredictable behavior), without significantly altering the outcomes when we simulate the thermal evolution for ten full orbital revolutions.

Fig. 7 Erosion sustained after ten orbital revolutions for each facet of pit 5, and different values of the dust-to-ice mass ratio: 0.5 (A), one (B), and two (C).

3.4 Thickness of the surface dust mantle

As highlighted by the results obtained in the previous tests, the thickness of the dust mantle at the surface may play a key role in the evolution of morphological features such as pits. The long-term survey of 67P’s nucleus by Rosetta has revealed that large smooth plains in the northern hemisphere are covered by a dust mantle, which originates from the southern hemisphere (Thomas et al. 2015a), following the ejection of dust particles around perihelion (Keller et al. 2015, 2017). The thickness of this mantle is unknown and almost certainly nonuniform (Hu et al. 2017; Davidsson et al. 2021). However, Davidsson et al. (2022) suggest that the dust mantle in the northern hemisphere may typically be thinner than 2 cm. Interestingly, Herny et al. (2021) noted that the assumption of the presence of a thin dust mantle was required in order to fit the Rosetta Orbiter Spectrometer for Ion and Neutral Analysis/Double Focusing Mass Spectrometer (ROSINA/DFMS) measurements, in particular to reproduce the patterns of volatile production rates. We have seen from the simulation outcomes discussed previously that a thin dust mantle naturally forms at the surface, without affecting the general results in terms of activity and erosion. Indeed, such a thin dust mantle is periodically produced and removed from the surface, after dust particles are dragged by escaping gas.

We further tested the influence of the thickness of a dust mantle, assuming that such a layer is initially present at the surface of the nucleus when the nucleus reaches its current orbit back in 1959 (Maquet 2015). Several values for the thickness were tested: 5 cm, 10 cm, 30 cm, 60 cm, and 1 m. We find that a 5cm-thick dust mantle is easily removed by cometary activity after the first perihelion passage. In previous tests, such very thin mantles did form and were removed, as described in the previous section, whenever the appropriate conditions were met. The presence of thin dust mantles on the surface of the facets is thus expected from the previous simulation results. As a consequence, adding an initial dust mantle of 5 cm gives the same simulation outcomes as if no dust mantle was considered. When the simulations start with a 10 cm-thick dust mantle, cometary activity is quenched for most facets of this pit (located in the northern hemisphere), as shown in Figs. 9 and 10. Only four facets remain active throughout the ten orbital revolutions. Their sustained activity might be attributed to two factors: they receive the most energy close to perihelion, which leads to a short period of strong activity that removes the dust mantle, or they receive the most energy integrated over one orbit, which allows them to eventually remove the dust mantle after several perihelion passages (Fig. 10 bottom panel). When the dust mantle thickness is larger than 10 cm (i.e., 30 cm, 60 cm, or 1 m), the activity remains quenched, and we do not observe any facets able to remove this layer during the ten orbital revolutions. We thus only show results for a thickness of 30 cm for comparison, in Fig. 10.

3.5 Set of uniform initial parameters

Whether cometary nuclei are homogeneous in composition or thermal and mechanical characteristics remains a matter of debate. Heterogeneities can be found at various spatial scales and might also be the result of evolution, due to initial nonuniform thermophysical properties, or simply a nonuniform insolation of the surface, due to the shape of the nucleus and seasonal variations (Guilbert-Lepoutre & Jewitt 2011). Therefore, it is entirely possible that different processes, related to nonuniform characteristics, are at the origin of the formation and evolution of sharp morphological features such as pits. However, in this study, we seek to quantify how energy received at the surface – through direct insolation, local and global self-heating, and shadowing effects – may give rise to morphological features of the spatial scale observed by Rosetta. Therefore, we seek a set of initial values for the thermophysical parameters, such that the same set can be used for all facets, on all our morphological features of interest, regardless of their location on the nucleus. These might not be representative of all the local conditions for all features studied, but are a good enough approximation to quantify local and global trends in the erosion rate. In the remainder of our study, we therefore consider a bulk porosity of 75%, and a dust-to-ice mass ratio of one. No CO or CO₂ was included in the ice mixture, to avoid numerical instabilities, but also because our tests show that they do not contribute to any significant change in the resulting erosion patterns, after ten orbital revolutions. Finally, no initial dust mantle at the surface was added. However, we stress that the formation of such a mantle is a natural consequence of cometary activity: the cyclic formation and destruction of such a layer of dust deposit is fully taken into account in our thermal evolution model. The previous sections detailed the effect of each of these critical parameters on the evolution outcomes for one specific pit. The following section will make use of the uniform set of parameters, to study 30 morphological features (circular and elongated pits, and alcoves) located across the surface of 67P.

Fig. 8 Influence of adding CO and CO2 to the ice composition. A: the erosion sustained after ten orbital revolutions for seven facets of pit 5, for various CO and CO2 abundances: no CO and CO2, 1% CO and 1% CO2, and 5% CO and 15% CO2 from left to right, respectively. B1 and B2: the activity patterns of two facets are given – H2O production rate on the left, CO2 and CO production rates in the middle column, and thickness of the dust layer on the right. The activity of facet 1 is quenched in the presence of CO and CO2. The behavior of facet 2 is given for comparison: it remains active in the presence of CO and CO2, preventing the formation of a dust layer (right plot).

4 Evolution of 30 morphological features across the surface of 67P

4.1 Energy received at the surface: General trends

To study the thermal processing of the 30 morphological features, a total of 380 facets were selected across the surface of 67P. For each facet, the energy input was computed, taking into account the effects of shadowing and self-heating as described in Sect. 2, and applied in Sect. 3. In Fig. 11 we show two quantities related to the energy input at the surface of 67P: the total energy integrated over one orbit, and the maximum energy flux received by each of the 380 facets. The maximum energy input is typically reached at perihelion for facets in the south, or just before perihelion for facets in the north. These two quantities were found to be essential to interpret the results of the thermal evolution model. Indeed, we see that the greatest amount of thermal processing – inducing substantial water ice sublimation and erosion – occurs during the perihelion passage, when the nucleus receives most of the direct solar energy on the southern hemisphere. As a result, the maximum energy is representative of this seasonal activity trend, and the maximum energy map does show the expected north-south dichotomy. However, we see from Fig.12 that a similar final amount of erosion can be achieved either by having a high amount of energy at or close to perihelion, or by having an increased amount of energy integrated over one orbit. However, we also see that for the same value of the integrated energy, facets that receive the bulk of the energy at perihelion tend to erode more (Fig. 12). Finally, these trends are not an absolute rule, and this justifies using a full subsurface thermophysical model, as we do below, rather than simply assessing the thermal behavior from surface energy balance maps.

Fig. 9 Progressive erosion sustained during ten full revolutions on 67P’s current orbit, for all facets of pit 5, and three values of initial dust mantle’s thickness: 0 cm (blue), 5 cm (red), and 10 cm (green). Vertical lines correspond to aphelion passages.

4.2 Effects of local topography and global shape

To the first order, the distribution of energy (integrated or maximum) is dominated by expected seasonal effects. We additionally see some variations for a given latitude and amongst facets related to one morphological feature. These can be generally attributed to shadowing and self-heating effects. Shadows are cast at the surface of 67P on a large scale (e.g., the neck area between the two lobes), as well as on a small, topographic scale (e.g., the bottom or part of the walls of deep circular pits). These can induce a significant decrease in the energy input, by as much as 70%, depending on the facets’ location and orientation. The effect of shadows can, however, be slightly offset by self-heating from neighboring facets (Fig. A.2). For most pits, self-heating contributes less than 20% of the total energy received at the surface. Thus, direct insolation dominates the energy input and self-heating is not the main activity driver. However, for several complex topographic configurations, where facets are not easily reached by direct insolation, self-heating can exceed the contribution from direct insolation. For these specific facets, the contribution of self-heating can reach more than 60% of the total energy received at the surface (Fig. A.2). They are typically located on the walls and at the bottom of deep circular pits. On a larger spatial scale, we also find such facets on alcoves close to the neck region, which are periodically in the shadow of the small lobe, and thus receive self-heating from it.

For the sake of completeness, we seek to quantify the relative contributions of the local topography versus the global morphology of the nucleus to the amount of self-heating. We thus compare the energy input for some facets of the shape model, and the energy input of the same facets when we numerically remove the small lobe from the shape model. This comparison is most informative for features 18 and 19 (also known as Seth_05 and Seth_04, respectively). These are two alcoves located close to the neck area, whose evolution is extremely affected by the presence of the small lobe. The integrated energy received over one orbit, with and without the small lobe in the shape model, is given in Fig. 13 (panels A and B). Facets on the alcoves receive up to 70% more energy when the small lobe is absent, due to the direct insolation reaching them. A detailed look at the various energy contributions informs us that the decrease in energy input from self-heating is not as significant as expected. For facets located at the bottom of those alcoves, direct insolation becomes the dominant source of energy, as expected, although the contribution of self-heating does not drop to zero. For one facet, there is even a slight increase in the self-heating contribution (~7%). This is due to the fact that surrounding facets receive much more direct insolation, and hence can transmit more energy. Overall, the contribution of the small lobe (vs. local topography) to the input of self-heating is not dominant (see Fig. 13). The small lobe contributes to up to ~ 22% of the total energy received by features 18 and 19. This is almost half of the total self-heating contribution for these alcoves, located in a region of the nucleus where the global contribution is maximum. However, in other regions, local topography is the major source of self-heating.

4.3 Thermal evolution simulations

The energy received at the surface of each facet, with the global distribution and trends described above, is the boundary condition for thermal evolution simulations performed over ten full orbital revolutions. This energy input was used to quantify the activity for each facet, for example, phase transition, gas production, dust mantling, and erosion. To keep our model relatively simple, we did two things: (1) we used a uniform set of initial parameters for each facet, as derived from Sect. 3; and (2) we did not account for any influence of shape evolution on the illumination conditions (i.e., erosion sustained at each time step was not used to modify the geometry of morphological structures). Instead, erosion after ten cometary orbits was assessed from the current shape of 67P’s nucleus, as observed by Rosetta. Global results, obtained for the 380 facets across 30 morphological features on the surface of 67P are represented in Fig. 14.

4.3.1 Latitudinal variations

We see that erosion at the surface is mostly correlated with the energy received at or close to perihelion. As a result, a stark contrast is observed between both northern and southern hemispheres. After ten orbital revolutions, erosion can reach up to about 77 m in the most active, southern regions in our study. In contrast, it does not exceed -30 m for most northern features. We see that facets directed toward the equator, while in the northern hemisphere, sustain enhanced erosion compared to other facets at the same latitude. For those, it can reach the same level of erosion as is seen in the southern hemisphere. As a consequence of the trends in the surface energy distribution described in the previous section, the pattern of latitudinal variations for erosion is clearly observed. Indeed, the amount of erosion after ten orbits decreases when facets are located closer to the north pole. In the northern hemisphere, facets sustaining the most erosion are those closest to the equator, or perpendicular to the equatorial plane, as they receive direct insolation around successive perihelion passages. Nonetheless, some of these first order latitudinal effects are mitigated in part, due to the complex topographic shape of 67P’s nucleus, which induces local self-heating. The important result for these general considerations is that the amount of erosion achieved after ten orbits (the assumed current period of time that the comet has spent in its current orbit) never reaches the observed dimensions of any of the observed morphological structures. For example, the smallest feature in our study (feature 12, also known as Ma’at_01) has a typical average size of ~130 m, and would sustain an increase of its diameter by only 10–15 m after ten orbits. The largest amount of erosion among our 380 facets remains below 80 m. Therefore, we confirm that pits cannot form by the progressive erosion of 67P’s surface.

Fig. 10 Influence of the presence of a dust mantle. A: the erosion sustained after ten revolutions with an initial dust mantle of 0, 5, 10 and 30 cm, left to right, respectively. B: the H2O production rate, the dust production, the progressive erosion, and the thickness of the dust layer, given for five facets, for an initial dust mantle of 10 cm.

Fig. 11 Energy flux received at the surface for all 380 facets, distributed over the 30 studied pits. A: the total quantity integrated over one complete orbit of 67P; B: the maximum reached during the perihelion passage.

4.3.2 Local variations

To the first order, the latitudinal pattern of erosion dominates. However, at the scale of each morphological structure, local trends appear similar across the surface. The first trend we observe is that erosion is generally more intense on the plateaus surrounding the pits when they are exposed to the Sun. In contrast, the bottoms of these pits do not sustain as much erosion, even after ten full orbital revolutions. This is especially true for circular pits with a high depth/diameter ratio (e.g., pits 1 or 2, also known as Seth_01 and Seth_02 with Seth_03 combined, respectively, and pit 12, also known as Ma’at_01). This general behavior tends to erase the local topography and leads to shallower features, such as those observed in the southern hemisphere. In general, the walls of the pits experience some differential processing, with erosion enhanced along a specific direction (e.g., features 1 or 5 in Fig. 15). This is directly related to the asymmetric distribution of the input energy, especially when some facets receive direct insolation while others mostly get a self-heating contribution. This suggests that, if we account for the shape evolution due to erosion, elongated pits are more thermally processed than small circular ones. As a consequence, our results are consistent with deep, circular, or quasi-circular pits, such as the pits labeled as 1 (Seth_01 on the big lobe) and 12 (Ma’at_01 on the small lobe) in our study, which are the least processed pits, or the best preserved under the current illumination conditions.

Fig. 12 Erosion sustained at each facet as a function of the energy they receive, integrated over one orbit. The dotted gray lines show the median of this energy for the large regions, i.e., the small and big lobes on the northern hemisphere and the southern hemisphere. The dashed gray lines show the median of erosion sustained by the facets in these regions. The color code provides an indication of the peak energy received at or close to perihelion.

5 Discussion

5.1 Local and global shape effects

Our results show that the local topography and the complex global shape of the nucleus can considerably impact the energy balance at the surface (Sect. 4.2). This is particularly true when considering the different sides of a given pit. As a result, some walls and bottoms of pits are not as easily reached by insolation as the corresponding exposed plateaus, making the onset of activity in the inner parts of these morphological features more difficult. It is thus necessary to take into account the effects of both shadowing and self-heating at the scale of these depressions. These processes are also important at the scale of 67P’s nucleus, because its specific bilobate shape leads to the neck region being highly shadowed during the northern day. While self-heating is found to be mostly negligible compared to direct insolation for most facets we studied, it can be an important energy source in some cases, especially at the bottom of pits and around the neck region, where direct insolation is limited. In such locations, the contribution of self-heating to the local energy balance can reach up to 60%. These results are consistent with earlier studies. For instance, Keller et al. (2015) showed that self-heating could reach 50% of the total energy received in some areas of the neck region. Macher et al. (2019) also showed that, even though the average contribution of self-heating in the regions they studied was evaluated to be 1% of the direct insolation, it can be enhanced in rough areas not reached by direct insolation. In these locations, it could reach as much as 50% of the direct insolation contribution. The important contribution of self-heating was also emphasized by Tosi et al. (2019), for deriving the temperature map at high spatial resolution (<15 m) from the Visible InfraRed Thermal Imaging Spectrometer (VIRTIS-M) data. The aforementioned studies were performed using various resolutions of 67P’s shape model, which suggests that shadowing and self-heating are important at all scales. Therefore, our results are not very sensitive to the choice of spatial resolution for the shape model: using the 125k-facets shape model, with an average distance between facets’ nodes of about 20 m (Marshall et al. 2018), allows morphological features to be sampled without increasing the computation time required if smaller facets are chosen. Overall, detailed knowledge of the energy balance at the surface on a local scale is thus a necessary condition to quantify the effect of thermally induced processes on the evolution of the cometary surface. However, as discussed below, we find that it is not sufficient to understand the evolution of the surface, since the energy input does not translate into phase transitions and erosion in a straightforward manner.

5.2 Nonuniform properties

We show that current illumination conditions cannot result in the formation of the deep circular pits with such characteristics as observed by Rosetta. In the southern regions, where sublimation-driven erosion is the most effective, erosion reaches ~80 m at best (Fig. 14). We now discuss how the choice of initial parameters used in our thermal evolution model may influence this outcome. For instance, an increased porosity could result in larger amounts of erosion, as much as 50% for facets that receive the most energy (Sect. 3.1). However, it is unlikely that the bulk material in the uppermost layers has a porosity greater than 75% (Ciarletti et al. 2015). An increased dust-to-ice mass ratio had a similar effect in our simulation outcomes (Sect. 3.2), although we identified that this was actually more due to the resulting increase in thermal conductivity than the composition itself. Therefore, local variations in composition or thermophysical properties could also induce different amounts of local erosion. Such local heterogeneities have indeed been identified at the surface of 67P, with a spatial scale of tens of meters, sometimes associated with the local exposure of volatile ices (e.g., Filacchione et al. 2016; Fornasier et al. 2016). On a global scale, differences between the small and the big lobes have been inferred from variations in their mechanical properties (El-Maarry et al. 2016), and physical characteristics. For instance, the small lobe has larger goose-bump features, fewer morphological changes, and less frequent and smaller frost areas than the big lobe (Fornasier et al. 2021). From these, the authors inferred that the small lobe might have a lower volatile content than the big lobe. Instead, we chose to apply a uniform set of initial parameters. Thus, our erosion rates could vary if we accounted for the actual heterogeneity of the nucleus. Based on the outcomes of simulations performed to select this set of initial parameters in Sect. 3, we can estimate that the final erosion would change by about 20% at most, due to local changes of porosity, composition, or thermal properties, as observed by the suite of instruments on board Rosetta. Nonetheless, our general trends, based on the relative erosion between plateaus and bottoms, and differential erosion, are not sensitive to these initial conditions. As a consequence, our quantitative study validates the qualitative trend suggested by Vincent et al. (2017) that sublimation-driven erosion leads to shallower and larger depressions, effectively erasing sharp geological features with time.

Fig. 13 Effects of the small lobe on the facing cliffs. A and B: the energy received at the surface of alcoves 18 and 19 integrated over one complete orbit with and without the small lobe in the shape model, respectively. C: the contribution of self-heating received from the small lobe only to the total energy received at the surface of structures 18 and 19.

Fig. 14 Erosion sustained after ten revolutions on 67P’s current orbit for all the 380 facets studied.

5.3 Dust mantle

The presence of a dust layer thicker than ~10 cm was able to quench the activity of most of the facets that we have studied (Sect. 3.4). If we consider that a thick (>10 cm) dust layer was initially present throughout the surface of 67P when it arrived on its current orbit in 1959, we would find erosion rates lower than those obtained in our simulations. It is interesting to note from Fig. 10 that, when the appropriate conditions are met, the activity of some facets is such that an initially thick dust mantle can be removed after several perihelion passages. In our simulations, thinner dust mantles are indeed periodically removed and formed as a direct consequence of ice sublimation and gas drag of dust particles. Evidence for such a cyclic formation and removal of dust with the seasons, and for fallback material, has been reported (e.g., Thomas et al. 2015a; Attree et al. 2019): a thickness of about 5 mm in northern regions has been reported (Herny et al. 2021). Through thermal evolution modeling, Davidsson et al. (2022) obtained a resulting dust mantle typically thinner than 2 cm. On a local scale, dust mantles may play a significant role. They would additionally be affected by the heterogeneous gravitational potential, impacting the local dust deposition at the surface. High-resolution observations by Rosetta/OSIRIS show that the bottom of deep, circular pits is relatively flat, and covered with a fine dust layer (e.g., features 1 and 2, also known as Seth_01, and Seth_02 and Seth_03, or feature 12, also known as Ma’at_01; Sierks et al. 2015). Some pits have boulders of various size on their floor, which Vincent et al. (2015) used as an indication of the erosion age of these structures. For instance, the authors suggested that the boulder-free floor of Ma’at_01 could represent the least eroded pit, while Ma’at_02 (feature 7 in our study) and Ma’at_03 would be increasingly eroded, with degraded walls and accumulated material within the pits. Our simulations cannot account for such effects. However, the degradation of walls after they are weakened and the accumulation of wall material at the bottom essentially lead to the same trend we found: pits become larger and shallower with time.

Fig. 15 Local examples of erosion achieved after ten orbits, highlighting differential erosion and flattening trends.

5.4 Active pits

Our study thus supports the hypothesis, initially made by Vincent et al. (2015), that the deep, circular pits are less processed (or better preserved) than the large or elongated ones. Interestingly, the more preserved features have been unambiguously revealed as the source of thin dust jets, arising from the edges of these depressions, which indicates that activity and erosion are currently occurring (Sierks et al. 2015). More generally, Vincent et al. (2015) identified two trends in the depth-to-diameter ratio (d/D) of pits at the surface of 67P: active pits have a high d/D (>0.3), while pits with no observed activity have a much smaller d/D.

From our results, we cannot exclude that large, relatively shallow pits could be active, as erosion is efficiently erasing the structures, especially in the southern regions (Sect. 4). Furthermore, these features typically receive high amounts of energy (integrated over an orbit, or at perihelion), such that adding moderately to highly volatile species triggered numerical instabilities in our simulations (Sect. 3.3). The sublimation of CO and CO₂ for these facets might actually trigger some outbursts, which no model can simulate, as the process is highly nonlinear. Indeed, these numerically unstable facets are typically found in areas of the pits that sustained the most erosion in previous tests (i.e., parts that received the most energy), either close to perihelion or integrated across the orbit. It is thus likely that these reflect bursts of activity driven by the sublimation of such species. However, the fact that no active outbursts were observed from these pits suggests that the sublimation fronts could actually be located deeper in the nucleus than in our model, after the insertion orbits. Moreover, the numerical instabilities were found within the first few orbital revolutions of 67P under current illumination conditions, and were thus not reflective of the time of the Rosetta observations.

For the best preserved structures, facets never experience such dramatic numerical behavior, yet Rosetta observations suggest some outbursts of activity. In our simulations, when volatile species are added, sublimation fronts slowly progress under the surface (Fig. A.4), yet continue to contribute to the activity. It is therefore very likely that CO and CO₂ remain close to the surface in these geological features, in accordance with their relatively unaltered nature. A further interesting aspect is that, when adding volatile species, facets that remain active (vs. those whose activity is quenched by progressive dust mantling) are those located at the edges of pits (Fig. 8). This corresponds to the observed activity of these pits (Vincent et al. 2015). Therefore, our results support the hypothesis that these morphological features are probably very well preserved, or are the least altered ones. Even without these additional volatiles, however, water ice is able to sublimate preferentially from the walls rather than the bottoms.

5.5 Implications for the evolution of pits

We have shown that cometary activity tends to erase surface features, so that deep, circular pits are likely the least processed morphological structures on the surface (Sect. 4.3). Clearly, these pits could not have been formed by sublimation-driven erosion. We have investigated very different illumination conditions across 67P’s surface. Under these conditions, the patterns of differential erosion, and the preference for eroding plateaus rather than bottoms of pits, are maintained. Therefore, we can extrapolate that different illumination conditions on a different orbit would have led to similar trends. Furthermore, even if the southern hemisphere is obviously more processed than the northern hemisphere, traces of larger depressions can be found, and there is no clear dependence of the distribution of depressions on latitude (Vincent et al. 2015). We can thus argue that pits were initially present on a global scale, and they likely evolved due to sublimation-driven erosion at various degrees on the surface of 67P.

Our results provide a quantitative confirmation for several studies since no quantification of the erosion through all the recent orbits has been performed before. Concerning the formation of pits, Ip et al. (2016) performed a morphological and dynamical study, by which they found that pits on JFCs were likely formed prior to acquiring their current orbital elements. Mousis et al. (2015) tested the formation of pits with three phase transitions (sublimation, amorphous water ice crystallization, and clathrate destabilization) and found that each of these processes would require a period of time much longer than the time spent by the comet in the inner solar system to form the observed pits. Guilbert-Lepoutre et al. (2016) also attested that it is very unlikely for 200-m pits to form under current illumination conditions. Such conditions are, however, prone to the formation of smaller-scale geological features, such as shallow depressions of several meters in depth, probably formed due to progressive seasonal erosion (Bouquety et al. 2021a,b). When it comes to the evolution of pits, Belton (2010) proposed an evolutionary sequence where pits are erased through cometary activity: initially found as acute depressions seen on 81P/Wild 2, they would progressively become shallower depressions as observed on 103P/Hartley 2, which is relatively older in terms of the sublimation process (Ip et al. 2016). Vincent et al. (2017) studied the global topography of comets observed by spacecrafts and reaffirmed this trend. This paper provides a quantification of the erosion rates sustained at the level of the pits during all the time that 67P spent as a JFC in the inner solar system, which vigorously reaffirms the previous studies, at least for 67P.

In the future, we will confront the trends established in our study by constraining the sublimation-driven erosion sustained by other cometary nuclei where pits have also been observed, in particular 103P/Hartley, (Syal et al. 2013), 81P/Wild 2, (Brownlee et al. 2004), and 9P/Tempel 1, (Thomas et al. 2013). When it comes to understanding the origin of pits, we argue that feature 1 (Seth_01), on the big lobe, and 12 (Ma’at_01), on the small lobe, are the least processed. Notwithstanding local heterogeneity giving rise to various pit sizes, these features are thus likely representative of pits as they were formed. This needs to be kept in mind when we seek to constrain the thermal or physical processes that carve these structures, and which remain to be identified: any process invoked needs to be able to excavate a significant volume of material in a quasi-circular shape.

6 Summary

We have investigated the erosion of morphological features at the surface of 67P/Churyumov-Gerasimenko, dominated by water ice-driven sublimation. We selected 380 facets of a medium resolution shape model of the nucleus, sampling 30 pits and alcoves across the surface. The energy balance at the surface was then computed with a high temporal resolution, and by including shadowing and self-heating contributions. We then applied a thermal evolution model to quantify the amount of erosion sustained after ten orbital revolutions under current illumination conditions.

Our study shows that a detailed knowledge of the energy balance at the surface on a local scale is a necessary condition to quantify the effect of thermally induced processes, but is not sufficient. Indeed, the energy input does not translate into phase transitions and erosion in a straightforward manner. Also, although seasons drive the global erosion trends, local topography can play a significant role in the final erosion state.

The erosional behavior on the surface revealed that morphological features such as pits and alcoves become larger and shallower with time: they are effectively erased through sustained cometary activity.

Finally, none of the surface structures we studied can be formed through progressive erosion. Pits Seth_01 and Ma’at_01 are among the least processed representatives of what pits would have looked like when they were formed, although the forming mechanism remains to be elucidated.

Acknowledgements

We thank the referee for their thorough review and comments which improved this paper. This study is part of a project that has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 802699). We gratefully acknowledge support from the PSMN (Pôle Scientifique de Modélisation Numérique) of the ENS de Lyon for the computing resources. We thank the European Space Research and Technology Centre (ESTEC) and the European Space Astronomy Centre (ESAC) faculty council for supporting this research. J.L. thanks the support by the CNES. A.B. thanks the Swedish National Space Agency (SNSA) grant 108/18 for its support.

Appendix A Additional information

Fig. A.1 Fraction of the energy input from self-heating relative to the total energy received. We highlight some examples where the self-heating contribution is significant.

Fig. A.2 Erosion (in meters) achieved after ten orbital revolutions, for all facets and morphological features studied. The blue box contains depressions located in the small lobe, the orange box contains the big lobe’s southern depressions, and the rest are located in the big lobe’s northern hemisphere.

Fig. A.3 Erosion sustained after ten orbital revolutions for all the structures we studied.

Fig. A.4 Subsurface retreat of the sublimation fronts of CO and CO2 for the two facets studied in Sect. 3.3. We display the two compositional cases: 1% CO and CO2, and 5% CO and 15% CO2. The gray vertical lines mark the perihelions.

References

All Tables

All Figures

	Fig. 1 Image from OSIRIS/NAC of a part of the Seth region on which we illustrate the type of depressions we study: pits and alcoves (half circular-pits, Vincent et al. 2015).
In the text

Fig. 2 Display of the facets selected for the study of the 30 pits. Top: the location of the facets on the surface of 67P. The shape model presented is the SPG model composed of 124 938 facets (Preusker et al. 2017), which is used for the surface energy calculation. Bottom: the location of the facets on a 2D map of 67P, which is a projection of the high-resolution SPG shape model composed of 12 million facets.

In the text

Fig. 3 Pit selected for the study of the influence of initial parameters, and corresponding energy. A: the location of pit 5, and facets sampled on the plateaus, walls, and bottom. B: the energy received at the 15 facets over one complete orbit, averaged over a daily period window. The gray line marks the perihelion. C: the total quantity of energy integrated over one orbit (left) and the maximum reached during the perihelion (right).

In the text

	Fig. 4 Progressive erosion sustained during ten full revolutions on 67P’s current orbit, for all facets of pit 5, and three values of the porosity: 60 (blue), 70 (red), and 80% (green). Vertical lines and numbers correspond to aphelion passages.
In the text

	Fig. 5 Erosion sustained after ten orbital revolutions for each facet of pit 5, and different values of porosity: 60 (A), 70 (B), and 80% (C).
In the text

	Fig. 6 Progressive erosion sustained during ten full revolutions on 67P’s current orbit, for all facets of pit 5, and three values of the dust-to-ice mass ratio: 0.5 (blue), one (red), and two (green). Vertical lines correspond to aphelion passages.
In the text

	Fig. 7 Erosion sustained after ten orbital revolutions for each facet of pit 5, and different values of the dust-to-ice mass ratio: 0.5 (A), one (B), and two (C).
In the text

Fig. 8 Influence of adding CO and CO2 to the ice composition. A: the erosion sustained after ten orbital revolutions for seven facets of pit 5, for various CO and CO2 abundances: no CO and CO2, 1% CO and 1% CO2, and 5% CO and 15% CO2 from left to right, respectively. B1 and B2: the activity patterns of two facets are given – H2O production rate on the left, CO2 and CO production rates in the middle column, and thickness of the dust layer on the right. The activity of facet 1 is quenched in the presence of CO and CO2. The behavior of facet 2 is given for comparison: it remains active in the presence of CO and CO2, preventing the formation of a dust layer (right plot).

In the text

	Fig. 9 Progressive erosion sustained during ten full revolutions on 67P’s current orbit, for all facets of pit 5, and three values of initial dust mantle’s thickness: 0 cm (blue), 5 cm (red), and 10 cm (green). Vertical lines correspond to aphelion passages.
In the text

	Fig. 10 Influence of the presence of a dust mantle. A: the erosion sustained after ten revolutions with an initial dust mantle of 0, 5, 10 and 30 cm, left to right, respectively. B: the H2O production rate, the dust production, the progressive erosion, and the thickness of the dust layer, given for five facets, for an initial dust mantle of 10 cm.
In the text

	Fig. 11 Energy flux received at the surface for all 380 facets, distributed over the 30 studied pits. A: the total quantity integrated over one complete orbit of 67P; B: the maximum reached during the perihelion passage.
In the text

Fig. 12 Erosion sustained at each facet as a function of the energy they receive, integrated over one orbit. The dotted gray lines show the median of this energy for the large regions, i.e., the small and big lobes on the northern hemisphere and the southern hemisphere. The dashed gray lines show the median of erosion sustained by the facets in these regions. The color code provides an indication of the peak energy received at or close to perihelion.

In the text

	Fig. 13 Effects of the small lobe on the facing cliffs. A and B: the energy received at the surface of alcoves 18 and 19 integrated over one complete orbit with and without the small lobe in the shape model, respectively. C: the contribution of self-heating received from the small lobe only to the total energy received at the surface of structures 18 and 19.
In the text

	Fig. 14 Erosion sustained after ten revolutions on 67P’s current orbit for all the 380 facets studied.
In the text

	Fig. 15 Local examples of erosion achieved after ten orbits, highlighting differential erosion and flattening trends.
In the text

	Fig. A.1 Fraction of the energy input from self-heating relative to the total energy received. We highlight some examples where the self-heating contribution is significant.
In the text

	Fig. A.2 Erosion (in meters) achieved after ten orbital revolutions, for all facets and morphological features studied. The blue box contains depressions located in the small lobe, the orange box contains the big lobe’s southern depressions, and the rest are located in the big lobe’s northern hemisphere.
In the text

	Fig. A.3 Erosion sustained after ten orbital revolutions for all the structures we studied.
In the text

	Fig. A.4 Subsurface retreat of the sublimation fronts of CO and CO2 for the two facets studied in Sect. 3.3. We display the two compositional cases: 1% CO and CO2, and 5% CO and 15% CO2. The gray vertical lines mark the perihelions.
In the text