Pits on Jupiter-family Comets and the Age of Cometary Surfaces

In 2001, Deep Space 1 explored comet 19P/Borrelly, revealing an elongated and extremely dark nucleus (Soderblom et al. 2002) with a highly variegated surface that can be divided into two terrain units: smooth and mottled terrain. A number of rounded depressions that could be identified as pits are visible in the mottled terrain down to a scale of ∼200 m. These pits have a similar size, which supports the fact that they might not be related to impact craters but rather the result of activity-related mechanisms (Vincent et al. 2015a). While the smooth terrain is located in the sunward direction, the mottled terrain appears largely inactive, as it is not associated with sources of gas and dust. The best spatial resolution obtained from Deep Space 1 is of the order of 50 m pixel^–1, compared to the average 10–15 m pixel^–1 obtained by the subsequent flyby missions described below. As a consequence, it is not possible to exploit the terrain model derived from the flyby images and sample pits in any useful way for this study.

The Stardust images have revealed the presence of pits on the surface of 81P/Wild 2 (hereafter 81P). They vary in both size and shape, some structures reaching several dozens of meters and up to ∼2 km in diameter. Brownlee et al. (2004) identified two types of pits: circular and irregular-shaped. Circular pits further exhibit two primary morphologies: pit-halo and flat-floored. Their origin has been presumed to be linked with impacts, possibly combined with sublimation and ablation processes. Indeed, hypervelocity impact experiments have successfully replicated pit-halo and flat-floored craters by impacting resin-coated sand with different degrees of porosity (Brownlee et al. 2004). In that framework, the lack of small impact structures (with sizes <0.5 km) would be attributed to surface erosion or a limited number of impactors within the corresponding size range. Additionally, 81P exhibits noncircular depressions that have been assumed to be formed by a combination of sublimation, mass wasting, and ablation processes (Brownlee et al. 2004).

Observations by Deep Impact and Stardust/NExT revealed that 9P/Tempel 1 (hereafter 9P) has a very pitted surface (Belton et al. 2013), with 380 pits ranging in diameter from tens to hundreds of meters (up to ∼900 m) and a depth of up to 25 m. Two of these depressions are considered as plausible impact craters (Thomas et al. 2007). Belton et al. (2013) inferred that JFCs would enter the inner solar system lacking "primitive" craters (i.e., formed through an intense, early collisional bombardment), and that most of these pits would likely result from outbursts of cometary activity. Indeed, they suggested that outbursts could account for the formation of 96% of them, and the process could contribute to a significant portion of total nucleus mass loss, in addition to sublimation. Finally, Belton et al. (2013) proposed that a few acute depressions may have resulted from sinkhole collapse because the expected formation timescale for these surface structures substantially exceeds the corresponding sublimation timescale.

Similarly to 9P, the surface of comet 103P/Hartley 2 (hereafter 103P) displays depressions indicative of a formation process different than impacts. Bruck Syal et al. (2013) proposed that most surface structures, including circular depressions or pits, could be the products of evolving jets arising from vents, active during several orbits. In this framework, surface material located on the periphery of a vent could fall into pits or cracks during periods of low activity, eventually leading to shallower structures. This would bring warmer material into contact with the colder, icy material located at the bottom of the vent. This process could also apply to material tumbling from scarps and ridges, both at the surface of 103P and 9P (Farnham et al. 2013). The relation between jets and pits was supported by Thomas et al. (2013b), who additionally investigated the hypothesis of collapsing subsurface cavities.

Taken all together, these observations suggest that pits may be ubiquitous on cometary surfaces, and that a link with cometary activity may exist. Moreover, the pits observed on 9P, 81P, 103P, and 67P display some morphological similarities in shape and dimensions (Vincent et al. 2015a; Ip et al. 2016). However, some of the pits observed on these comet nuclei exhibit a lower depth-to-diameter (d/D) ratio compared to 67P. Active pits on 67P have an average d/D of ∼0.73, while inactive pits have a shallower d/D of ∼0.26 (Vincent et al. 2015a). The observed average d/D is ∼0.2 for 81P (Kirk et al. 2005; Vincent et al. 2015a) and ∼0.1 for 9P (Thomas et al. 2013a). Also, the aspect of comet 103P's surface resembles that of 9P, exhibiting a smooth appearance with no evident deep pits (Ip et al. 2016).

In this work, we want to understand how cometary activity may modify these surface structures and whether signatures of their formation process can be inferred from their expected evolution through sustained activity. We thus apply the same method as used to study the evolution of pits at the surface of 67P (Benseguane et al. 2022) to quantify the amount of erosion sustained by pits at the surface of 9P, 81P, and 103P under their current illumination conditions. We note that 19P/Borrelly cannot be included in this study, as there is no shape model of its nucleus with a sufficient spatial resolution to be used. We are interested in two quantities: the erosion sustained during each orbital revolution (i.e., the erosion per orbit in the following) and the erosion sustained as a result of multiple revolutions. We summarize in Section 3 the different steps of our method. Results for each comet are presented in Section 4 and discussed in Section 5.

For each comet we study, a high-resolution shape model is key to apply our surface energy model. We thus use the highest spatial resolution available for each nucleus, so to capture the effects of both their global shape, and local topography, pits in particular:

• 81P—derived from the Stardust Navcam images by Farnham et al. (2005), ⁶
• 9P—derived from the images obtained by the Deep Impact and Stardust missions by Farnham & Thomas (2013a), ⁷
• 103P—derived from images obtained by the EPOXI mission by Farnham & Thomas (2013b). ⁸

Shape models for comets 81P, 9P, and 103P are not as spatially resolved as those available for comet 67P. This lower quality (for our purpose) is to be expected, as these were derived from limited observations during flyby missions, whereas 67P was escorted during the 2 yr of an orbiting mission. Again, we note that we had to exclude 19P from our study, as no available shape model has a good enough resolution for our study. More specifically, the shape models of 103P and 9P do not always reach the spatial resolution to unambiguously show features as deep as the pits observed on corresponding surface images. For example, the southern hemisphere of 81P was not observed during the Stardust flyby; missing data have been completed by a smooth ellipsoid of revolution based on the average surface ellipsoid of the comet.

It is important to mention that the smoother appearance of depressions on the surfaces of the studied comets is not solely due to the comparatively lower resolution of their shape models relative to 67P but is also because some of the pits observed on these comets exhibit a lower depth-to-diameter (d/D) ratio compared to 67P, as described in Section 2.5. Despite these limitations, we are able to identify a number of pits and alcoves on the shape models. We select facets similar to those of 67P, located on the plateaus, walls, and bottoms of each structure, to the best of our ability and the local resolution of these shape models.

We select a minimum of 10 pits of each nucleus, located across all latitudes of each nucleus. Indeed, we show that seasonal effects dominate the global erosion trends on 67P (Benseguane et al. 2022), so covering the entire nucleus is important in this study too. We thus apply selection criteria similar to 67P. First, we aim to sample latitudes as much as possible, to assess the influence of seasonal mechanisms. Second, we focus on large pits (rather than small and shallow depressions, as seen on the smooth terrains of 67P), characterized by steep walls and flat bottoms, that have sizes ranging from tens to hundreds of meters. Indeed, pits ranging from ∼150 m to ∼1 km exhibit a size–frequency distribution that is similar to 67P and those observed on 9P and 81P (Vincent et al. 2015a; Ip et al. 2016). Effectively, we exclude smaller thermokarstic features from this study, as we did for 67P, because their formation and evolution appear to be different than the large and deep circular pits (see Benseguane et al. 2022 and references therein). On the corresponding shape models, we select multiple facets on different sides of each pit (plateaus, bottoms, and walls). When we select a pit on the images that is not easily identified on the corresponding shape models, we pick facets with appropriate latitude and longitude. In these limited cases, our modeling outputs would serve the purpose of constraining the seasonal trends and explore possible evolution scenarios for the corresponding surface features. For that same purpose, though we mostly select facets located in the northern hemisphere of 81P, we also select a group of facets in its southern hemisphere, even if there is no direct evidence for pits. This will allow us to compare erosion rates across the entire nucleus to form a more complete picture. For that purpose, we not only select several pits across the entire surface of 103P but also add single facets randomly distributed across the surface in order to test the effects of this nucleus's unique elongated shape and complex rotational properties. For each facet of the shape models, we then compute the thermal environment, including self-heating and shadowing, either by neighboring facets or due to the complex global morphology of the nucleus as described below.

The surface energy and thermal evolution models used in this study are described in Benseguane et al. (2022). We provide a brief summary below. The total energy ${ \mathcal E }$ received by each facet selected for this study is the sum of different contributions: direct insolation E_⊙ (accounting for shadowing effects) and self-heating, i.e., the energy received by reflection and emission from neighboring facets in the visible E_VIS and infrared E_IR. Each contribution is given below. Direct insolation is given by

$$\begin{eqnarray}&&{E}_{\odot }=\displaystyle \frac{{F}_{\odot }}{{r}_{H}^{2}}\,\cos \xi ,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{E}_{\mathrm{VIS}}=\displaystyle \sum _{T}\,\,{{ \mathcal A }}_{T}\,\displaystyle \frac{{F}_{\odot }}{{r}_{H}^{2}}\cos {\xi }_{T}\,\displaystyle \frac{{S}_{T}}{\pi }\,\displaystyle \frac{\cos {\zeta }_{T}\,\cos {\zeta }_{R}}{{\delta }_{T}^{2}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{E}_{\mathrm{IR}}=\displaystyle \sum _{T}\,\varepsilon \sigma {T}_{T}^{4}\,\displaystyle \frac{{S}_{T}}{\pi }\,\,\displaystyle \frac{\cos {\zeta }_{T}\cos {\zeta }_{R}}{{\delta }_{T}^{2}},\end{eqnarray} \tag{ 3 }$$

where F_⊙ [W m⁻²] is the solar flux at 1 au, ${{ \mathcal A }}_{T}$ is the Bond albedo of an emitting facet (${{ \mathcal A }}_{T}$ = 0.06 for all calculations), ξ_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 emitting facets, and δ_T is the distance between the two facets. The value of the emissivity ε is 0.95 for all calculations. T_T is computed by considering 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. The heliocentric distance r_H [au] and the local zenith angle ξ both vary with time, as described in the following subsection.

The total energy ${ \mathcal E }$ is used in the surface boundary condition of a 1D thermal evolution model. This condition is given by

$$\begin{eqnarray}&&(1-{{ \mathcal A }}_{R})\,{ \mathcal E }=\varepsilon \sigma {T}^{4}+\kappa \displaystyle \frac{\partial T}{\partial r}+{f}_{{{\rm{H}}}_{2}{\rm{O}}}\,{\rm{\Delta }}{H}_{{{\rm{H}}}_{2}{\rm{O}}}\,{Q}_{{{\rm{H}}}_{2}{\rm{O}}},\end{eqnarray} \tag{ 4 }$$

where ${{ \mathcal A }}_{R}$ is the Bond albedo (with a value of 0.06 in all calculations) of the facet for which we compute the energy balance, ${ \mathcal E }={E}_{\odot }+{E}_{\mathrm{VIS}}+{E}_{\mathrm{IR}}$ is the total energy received at its surface, T [K] is the surface equilibrium temperature, ${f}_{{{\rm{H}}}_{2}{\rm{O}}}$ is the fraction of the facet covered by water ice, ${\rm{\Delta }}{H}_{{{\rm{H}}}_{2}{\rm{O}}}$ is the latent heat of water-ice sublimation, and ${Q}_{{{\rm{H}}}_{2}{\rm{O}}}$ [kg m⁻² s⁻¹] is the corresponding sublimation rate. We aim at constraining how the patterns of energy received at the surface (diurnal but most significantly seasonal) influence the activity of each nucleus and the erosion of its surface features. Thus, we assume that the thermal and physical characteristics are the same for each comet so that the contributions of varying parameters may be removed. In order to make our results comparable to our study of pits on 67P, we use the same set of initial parameters, a hypothesis with consequences discussed in Section 5. For instance, we assume that the material is a simple mixture of two components, water ice and dust, with a mass fraction ratio of 1 and a porosity of 75%. The resulting thermal conductivity is reduced by a Hertz factor of 0.005 to account for the limited contact between grains in this porous structure (e.g., see Guilbert-Lepoutre et al. 2023, for further description of this parameter).

Benseguane et al. (2022) showed the influence of each of these parameters on the erosion rates sustained on 67P, so we do not repeat these here, since the effects are the same. The chosen values for the initial parameters are listed and justified in Benseguane et al. (2022). Since the behavior of cometary material in our model depends on the energy received at the surface of each facet and its heating rate, we can extend the conclusions of Benseguane et al. (2022) on the influence of these parameters for this work. We note that including CO and CO₂ in the volatile mixture was not altering the evolution trends in any significant manner, since the most significant source of erosion was the sublimation of water ice. Therefore, we do not add these species in our mixture. The thermal evolution model includes the usual features such as heat and gas diffusion, phase transitions for water ice (crystallization and sublimation), drag of dust particles by the vapor phase, and formation of a dust mantle at the surface (Lasue et al. 2008).

The surface energy model is computed with a time step of 8 minutes for each comet. This allows one to achieve a good description of the diurnal patterns of heating and the resulting activity for any combination of spin state and shape. We note that the extreme members in that respect are 67P, studied by Benseguane et al. (2022), and 103P, which has a complex rotation state (Belton et al. 2013; Knight et al. 2015). For each time step, we first retrieve the coordinates of the subsolar point using SPICE kernels available on the WebGeocalc platform, ^{9 , 10} which contain the information on each nucleus's rotation state, pole orientation, and orbital parameters. Then, the insolation geometry for each facet is computed with respect to these subsolar point coordinates. As a result of 103P's complex rotation state, the SPICE kernels are probably only valid for the duration of the EPOXI flyby; extrapolations before and after the flyby duration may not be accurate. The impact of this complex spin state will be discussed with the results. We run our thermal evolution calculations to study the impact of current illumination conditions. However, each comet has "acquired" its current orbit following a distinct orbital evolution and has been holding it for a different period of time since its latest orbital change. We thus use a different number x of orbital revolutions for each comet to reflect their recent past history, based on backward dynamical integrations performed by Ip et al. (2016): x = 6 orbits for 81P, x = 13 orbits for 9P, and x = 20 orbits for 103P. In the following, we present both the erosion sustained during each orbital revolution, which allows us to compare comets, and the total erosion sustained under current illumination conditions after x revolutions, which allows us to assess whether erosion driven by the sublimation of water ice could carve the pits as they are observed on each individual comet.

The spin state of 81P is such that the subsolar point crosses a large range of latitudes near perihelion (from −60° to 60°; see Figure 1). As a result, the total amount of energy per orbit received by each facet selected in our study is relatively uniform across the surface. The slight asymmetry between the pre- and post-perihelion latitudes and corresponding heliocentric distances results in the southern hemisphere receiving, on average, almost twice as much energy as the northern hemisphere. As for 67P (Benseguane et al. 2022), the latitudinal effects (i.e., seasonal) dominate the energy distribution at the surface of 81P, but the local shape can also play a key role. Indeed, we note that only a few facets located on the walls of some pits receive lower amounts of energy in total, about half the maximum amount received by others, as a result of shadowing effects from neighboring facets. The contribution of self-heating to the total energy is relatively small, but it can account for up to 30% of the total energy in the shadowed regions, where direct insolation is weak. It should be noted, however, that because 81P's pits are quite large (∼2 km for the largest one; Brownlee et al. 2004), most facets are ultimately exposed to direct insolation.

Zoom In Zoom Out Reset image size

The erosion resulting from the sublimation of water ice is, consequently, relatively uniform for all selected facets in the north (Figure 1), amounting to 4–5 m per orbital revolution and ∼15–30 m at most after six orbital revolutions (Table 1). We note here that the southern hemisphere may erode more than the northern hemisphere. Indeed, the facets we considered in the south sustain a maximum erosion of the order of 40 m after six revolutions. Nonetheless, the lack of actual images and the corresponding poor quality of the shape model for this region of the nucleus do not allow us to draw any further conclusions. These results imply that erosion driven by water sublimation is likely not the primary process responsible for the formation of the large pits of diameters up to 2 km studied here. We note that in general, the pits on 81P are large enough for facets located at their bottoms to behave similarly to facets located on the surrounding plateaus. Both "parallel" planes thus erode in the same way, so that the depth of these features should be expected to remain relatively constant with time.

Comet 9P is, in the context of this study, the "simplest" comet. Due to its low obliquity, the largest amounts of energy are received at perihelion by facets located in the 0°–20° latitude range. The northern hemisphere receives slightly more energy than the southern hemisphere, and both patterns of peak and total energy correlate well with the latitude of the subsolar point (see Figure 2). In this regard, seasonal patterns resulting from a combination of shape and rotational properties also dominate the energy patterns, in a way even more obvious than for 67P or 81P. We note local differences for facets selected on each pit, which suggest patterns of differential erosion similar to those obtained for 67P (Benseguane et al. 2022). Erosion caused by water-driven outgassing is larger than for 81P, and because we consider a total of 13 orbits, it can eventually become substantial in equatorial regions. However, it never exceeds 100 m in total (Table 1), whereas the observed dimensions of the pits can reach several hundred meters across for the largest ones (Thomas et al. 2013b).

Zoom In Zoom Out Reset image size

Comet 103P is in non-principal-axis rotation with apparently changing component rotation periods (Belton et al. 2013; Knight et al. 2015). Its rotational properties are, in fact, so complex that the SPICE kernels have a limited range of validity around the EPOXI flyby of the nucleus, and propagating coordinates for the subsolar points to the whole orbit is not necessarily possible. This is nonetheless the best we can do at this point to assess the influence of activity on the evolution of surface features. We must keep this effect in mind to interpret our results. The nucleus spins in an excited long-axis mode, with its rotational angular momentum per unit mass and rotational energy per unit mass slowly decreasing, while the degree of excitation in the spin increases through perihelion passage (Belton et al. 2013). To further complicate the picture, the nucleus has a very elongated shape.

These characteristics are reflected in the complex distribution of energy received by the nucleus's surface (Figure 3), which exhibits not only a latitudinal trend (as observed for the other comet nuclei) but also strong variations across longitude, especially around the equator region. Overall, equatorial regions and nearby northern latitudes receive a substantial amount of energy around perihelion, while the southern and extreme western equatorial regions receive less energy during this period. The contribution of self-heating to the total energy on 103P is minimal, accounting for less than 10% of the total energy. This contribution is extremely low compared to 67P or 81P. This is primarily due to the pits on 103P having a low d/D (depth-to-diameter) ratio compared to 67P or 81P (see Table 1). Additionally, the low spatial resolution of the shape model may play a role in limiting the ability to effectively reproduce shadowing and self-heating effects on a scale smaller than 10 m. The two effects might compensate for each other; however, we can calculate that with an additional 10% of surface energy due to self-heating, the final erosion would be enhanced by ∼15% for the most eroded facets.

Zoom In Zoom Out Reset image size

Erosion on 103P is strongly correlated with the peak of energy received at or close to perihelion (see Figure 3). However, because of 103P's complex rotation, almost all facets are ultimately exposed at perihelion. This correlation can be attributed to the fact that erosion on 103P occurs predominantly during brief periods of intense heating, outside of which the energy is insufficient to cause facets to erode. This is in contrast to 81P or 9P, where energy is more consistently distributed throughout the entire near-perihelion passage, instead of occurring in brief peaks. At perihelion, 103P is also the comet closest to the Sun, at a distance of 1.05 au, compared to 81P (1.59 au), 9P (1.54 au), or even 67P (1.24 au). As a result, 103P is the comet that sustains the most erosion after 20 revolutions under current illumination conditions (Table 1).

For all comets considered in this work, cometary activity driven by the sublimation of water ice and the resulting erosion are primarily controlled by direct insolation and thus display strong seasonal patterns. Indeed, the thermal processing of each facet considered in this study depends not only on the peak and total energy it receives but also on how energy is delivered as a function of time to produce subtle effects that appear unique to each comet nucleus, in particular its shape and rotational properties, as shown in Figure 4. For comets 9P, 103P, and 81P, no global shape effects are observed due to the lack of substantial shape irregularities on a global scale, unlike the acute shape of 67P. The unique shape of comet 67P, for instance, makes the pits located near the neck of the nucleus more susceptible to shadowing effects by the smaller lobe (Benseguane et al. 2022). However, we found that local shape effects (i.e., linked to local topography) can be significant at the scale of a given pit. Self-heating can contribute important fractions to the total energy in the deep pits and steep cliffs of 81P, for instance, accounting for 30% of the total energy input. In contrast, it is minimal for 9P and 103P (<10%), where surface features are wider and shallower.

Zoom In Zoom Out Reset image size

Facets located at the bottom and on the wall of circular depressions can be affected by shadowing, compared to exposed plateau facets. Consequently, if these depressions are deep enough, they tend to become shallower over time due to water-driven erosion. This is due to the combination of two effects. First, plateaus tend to erode more than bottoms (being more exposed to direct insolation), so that pits tend to become shallower with time. In addition, walls sustain some differential erosion; over several perihelion passages, pits also become wider over time. When this trend is not observed, it is because the corresponding pits are already large or shallow, i.e., not deep enough compared to their diameter, as seen in most of the pits on 9P and 103P, or the pits are large enough that the bottom facets are directly exposed to the Sun, as in the case of the large pits of 81P. While 81P is the least eroded nucleus in our study, our results imply that most of its largest pits are already wide enough to prevent any further change in depth. We note here that additional processing of cliffs can occur following their collapses, as observed for 67P (Vincent et al. 2016; Pajola et al. 2017; El-Maarry et al. 2019). By this mechanism (not accounted for in our model), the filling of pits with debris material would tend to further reduce their depth. Overall, we suggest that pits can reach a depth-to-diameter (d/D) ratio that seems to prevent any further change due to erosion alone. Future studies could explore critical d/D thresholds in different illumination conditions, determining when pits are able to evolve or when their morphology becomes relatively fixed as a function of material properties. Ip et al. (2016) suggested that the d/D ratios of the large pits are mostly within the range of 0.1–0.3. In comparison, active pits studied by Vincent et al. (2015a) have a large d/D ratio (>0.3) and a small diameter (<300 m). This statistical result coupled with backward dynamical integrations suggests that large circular depressions could have grown from small and deep ones via erosive mass wasting of the surrounding areas (Cheng et al. 2013; Vincent et al. 2015a, 2016).

Overall, we find that sharp depressions are likely erased with time as a result of sustained cometary activity.

Most significantly, erosion sustained after the multiple perihelion passages is not able to carve large depressions with the observed size and shape on any of the comets we studied. Of course, some limitations arise from our methodology, most notably from the assumed uniform thermal and physical characteristics, for all pits we have studied. Local heterogeneities (in composition, albedo, or thermal properties, for example) could actually enhance the local erosion computed in our simulations. Significant deviations from our results can only be achieved with extreme values for the initial parameters we have considered; for example, a combination of 70% surface water ice with a porosity larger than 90% can double the amount of erosion. As a result, within the range of plausible parameters (see Benseguane et al. 2022, for a review), erosion could be increased by up to 20%–30% at most; this does not affect our general trends or the general conclusion that cometary activity tends to erase sharp surface features.

Comet 103P is an extreme example of the effects described above. In our simulations, the small northern lobe of 103P's nucleus is very active and experiences the most erosion as a result of its preferential exposure to the Sun at perihelion. The northern lobe was indeed observed to be active during the EPOXI flyby. More precisely, jets were clustered in the rough topography of the small northern lobe and mid- to northern part of the big lobe (A'Hearn et al. 2011). Nevertheless, the observed high activity may have been enhanced by an abundance of volatile species in specific regions of the small lobe not accounted for in our model. A'Hearn et al. (2011) indeed determined that different species were being ejected from the different parts of the nucleus, with H₂O vapor coming primarily from the waist and CO₂, H₂O ice, and organics coming primarily from the top of the small lobe. Taking into account a higher abundance of volatile species in the small lobe within our thermal evolution calculations results in increased erosion compared to the erosion yielded by a homogeneous nucleus assumption. This would further emphasize the contrast with the southern big lobe. EPOXI additionally revealed distinct terrains on the nucleus, with a smooth "waist" connecting two rougher lobes (A'Hearn et al. 2011; Knight & Schleicher 2013; Thomas et al. 2013b). No significant difference was noted in the concentration or appearance of the pits between the two lobes, though. In light of our results, we could argue that this might be caused by erosion being sustained in a similar manner across the two lobes. We mentioned that due to 103P's complex rotation state, the kernels used to derive the latitude of the subsolar point might not be accurate outside of the EPOXI flyby duration. Of interest to our results, the key effect of this spin state is primarily to expose the whole nucleus at perihelion, where peak energy fluxes are received at all latitudes to efficiently trigger water-ice sublimation and erosion (Figure 3). This occurs within the validity range of the SPICE kernels, so ultimately, our results should not be severely affected.

Overall, we find that progressive erosion, driven by the sublimation of water ice under current illumination conditions, is not able to form pits at the surface of JFCs because (a) the total erosion, even after calculations taking into account several orbital revolutions of orbits, remains lower than the observed dimensions of the pits, and (b) it tends to erase sharp features, which become shallower and wider with time.

If sharp features are indeed erased by erosion, driven in our simulations by the sublimation of water ice, then as a corollary, we can infer that the deepest, most circular pits are likely the most primitive or best-preserved pits. From the results of our thermal evolution model, including those obtained for 67P by Benseguane et al. (2022), we can "rank" the primitiveness of these surface structures observed on these four comet nuclei. Comet 81P would have the least processed pits (or best-preserved), followed by the northern hemisphere of 67P, then its southern hemisphere. Comet 9P overlaps in the thermal processing space with 67P, with its southern hemisphere being relatively unprocessed compared to its equatorial region. Finally, 103P is by far the most thermally processed JFC we have studied, due to a combination of aggravating factors: smaller perihelion distance, larger number of orbital revolutions close to the Sun, and complex rotation that leads to the relatively uniform processing of the entire surface.

Interestingly, Vincent et al. (2017) performed a statistical analysis of the distribution of large-scale topographic features on 67P and found that cliff height correlates with surface erosion rates and follows a power law with an average cumulative power index of –1.69. They suggested that topography could be used to trace a comet's erosional history. In this framework, large and sharp cliffs would characterize primordial surfaces, while eroded surfaces would display smaller blocks (e.g., boulders, pebbles, and dust). The power-law index of the corresponding topography cumulative height distribution could indicate how primitive a comet nucleus is. They performed the same statistical analysis of surface features observed on 81P, 9P, and 103P (see Table 2 of Vincent et al. 2017). They found that 67P and 81P would have experienced similar degrees of erosion, while comets such as 9P and 103P would be more eroded, in agreement with the suspected past dynamical histories for each comet (see Figure 11 of Vincent et al. 2017). They concluded that a comet recently entering the inner solar system would have a p-index of topographic height of around −1.5. Older comets show larger power indices, up to about −2.3.

Based on these results, Kokotanekova et al. (2018) proposed that the phase function–albedo correlation they had previously found (Kokotanekova et al. 2017) might be explained by the erosion of pits and rough surface topography. Based on their hypothesis, rough surfaces with steep phase functions would gradually evolve toward smoother terrains with decreased phase function coefficients. Our results stem from a distinct method, providing a physical model to both empirical studies (Vincent et al. 2017; Kokotanekova et al. 2018); they point to the same evolutionary sequence. Assuming that this overall interpretation is correct, a decreasing phase function coefficient would provide a useful observable to characterize the level of erosion of a cometary surface.

With a very different prism, Kelley et al. (2021) examined several outbursts observed on comet 46P/Wirtanen and found that the mass estimates were similar to or an order of magnitude larger than the mini-outbursts observed at comets 9P and 67P. They hypothesized that mini-outbursts on comets could be associated with steep terrain features, like cliffs and scarps, based on observations linking such mini-outbursts of comet 67P to these terrain features, and even their collapse (Vincent et al. 2019). Based on this assumption, they analyzed the outburst frequencies of comets 67P, 9P, 46P, and 103P. They suggested that the observed differences may be related to distinct surface terrains. Comets 67P and 9P displayed significantly higher outburst frequencies compared to 46P and 103P. Interestingly, comet 46P would appear as an evolutionary intermediate between 103P (very processed) and 9P (moderately processed) in terms of surface topography and erosion implied from the work presented here. Indeed, this comet has performed a number of orbital revolutions since its discovery, similar to 9P, with a perihelion distance decreasing from ∼1.6 au (i.e., similar to 9P) to ∼1.05 au (i.e., similar to 103P; Krolikowska & Sitarski 1996). This is in agreement with the processing sequence proposed by Vincent et al. (2017) and Kokotanekova et al. (2018).

The results we present here agree with the evolutionary sequence proposed from independent observables (Vincent et al. 2017; Kokotanekova et al. 2018; Kelley et al. 2021). With a distinct method, we can provide a physical framework for this evolutionary sequence that transforms "young" cometary surfaces, with sharp surface topography prone to spark mini-outbursts, into "old" cometary surfaces that are eroded and do not experience as many mini-outbursts, as summarized in Figure 5. By "young," what we mean here is that a comet nucleus's surface has undergone relatively little modification resulting from water-driven activity, although the nucleus is located in a region where water is efficiently sublimating (typically with q < 2.5 au). This is a consequence of each nucleus's unique past dynamical history and rotational properties. In this study, we have only skimmed over the most recent influence of each comet's past orbital evolution by accounting for a number of perihelion passages, which amounts to ∼40–130 yr of thermal processing at most (for 81P and 103P, respectively). However, the dynamical evolution since the comet nuclei left the outer solar system reservoirs to reach the orbit on which they are currently observed is much longer and more complex and entails some thermal processing that is not accounted for in our work (e.g., Gkotsinas et al. 2022). On the other hand, "old" refers to a surface that has significantly changed as a result of thermal processing, leading to water-driven activity and substantial surface erosion.

Zoom In Zoom Out Reset image size

Various scenarios have been proposed in the literature to explain the origin of pits (see Benseguane et al. 2022, for a review of 67P, and the Introduction of this study). Ip et al. (2016) found that such features with steep walls and flat bottoms, with sizes between 150 m and 1 km on 67P, have the same size–frequency distribution as those on 81P and 9P. Therefore, it seems reasonable to assume that the pits observed at the surface of JFCs share a similar origin, so that their morphological characteristics appear similar. These steep walls and flat bottoms make them different from the bowl-shaped impact craters found on the Moon or asteroids (Brownlee et al. 2004; Ip et al. 2016). Although we cannot exclude that some pits may remain associated with impact events, these could instead be considered as a signature of some process related to cometary activity rather than the result of collisions.

In light of the discussion above, we would like to highlight the conclusions of Belton et al. (2013), who suggested that most pits on the surface of 9P would likely be the most common surface features related to outbursts of activity. Additionally, Pozuelos et al. (2014) suggested that cometary outbursts could be the origin of the pits observed on 81P. From our results, key aspects of how water-driven cometary activity fuels the evolution of pits need to be recalled. First, we see that sharp features tend to be erased, as they become wider and shallower with time. Second, latitudinal effects are so strong that patterns of differential erosion tend to elongate initially circular features. If we assume that pits formed as cylindrical structures, erosion with time would lead to elongated features departing from this initial morphology. In order to carve deep, almost circular structures, it is therefore crucial that water ice does not sublimate whenever pits are formed. Therefore, it appears that pits might have been formed before JFC nuclei crossed the water snow line, a suggestion made by Ip et al. (2016), who studied the past dynamical history of these comets.

The morphological characteristics of the least processed pits we could identify imply a formation scenario where a rather explosive mechanism, able to carve a large amount of cometary material in a short period of time, occurs in a region of the solar system where water ice is not sublimating, or the freshly formed features would become progressively elongated and eroded. The Centaur phase experienced by each JFC (e.g., Gkotsinas et al. 2022) may be key here to understanding the origin of such surface features. In the giant planet region, several phase transitions occur that can lead to cometary activity among Centaurs (e.g., crystallization of amorphous water ice or CO₂ sublimation or segregation; Guilbert-Lepoutre 2012; Davidsson 2021, respectively). Sudden thermally induced events such as clathrate destabilization and the crystallization of amorphous ice could lead to outbursts of activity (Miles 2016; Wierzchos & Womack 2020), potentially leading to the formation of pits. Furthermore, several Centaurs are prone to recurrent, sporadic outbursts of activity, like 29P/Schwassmann–Wachmann 1 (Wierzchos & Womack 2020; Clements & Fernandez 2021; Betzler 2023; Lin 2023) or 174P/Echeclus (Kareta et al. 2019; Rousselot et al. 2016, 2021).

Our hypothesis for the formation of pits could be tested in the future by the Comet Interceptor mission (Snodgrass & Jones 2019). This mission is designed to encounter a Dynamically New Coment (DNC) which typically experiences only limited processing in the giant planet region in comparison to JFCs. Observing no pits could imply that their origin is linked to a process that exclusively affects JFCs. Since the key difference between these two populations rests mainly on their orbital evolution, the formation mechanism of pits should be sought there, resulting from the JFCs' dynamical evolution. Alternatively, a limited number of pits could be observed, which we might attribute to thermal processing prior to the flyby, possibly on the inbound part of the orbit, or the early processing prior to the ejection of the nucleus in the Oort Cloud. Comparisons with the characteristics of the pits observed on JFCs (depth, diameter, and location with respect to the subsolar point, for example) would help to pinpoint the origin of these surface features. Finally, observing as many pits and as much sharp topography at the surface of such a pristine comet nucleus as on JFCs would suggest that these are signatures of mechanisms at play during the earliest stages of comet formation, rather than the signature of processes at play during the Centaur phase of JFCs. Indeed, the implication would be that such rough surface topography would be common to all comet nuclei before water sublimation sets in, regardless of their subsequent orbital evolution.

Evidence for the evolutionary sequence provided by this work and the prior studies by Vincent et al. (2017), Kokotanekova et al. (2018), and Kelley et al. (2021) highlights the importance of more space- and ground-based observations of comet nuclei. In particular, the best way to verify the validity of this sequence is to (a) increase the number of comets and (b) use multiple independent techniques to cross-check the resulting sequences. We advocate that programs targeting JFCs at various stages of their evolution will be of primary importance in advancing our understanding of these objects. In particular, space missions toward active or outbursting Centaurs (such as 29P/Schwassmann–Wachmann 1) would certainly prove instructive. Missions toward and targeted telescope observations of less evolved Centaurs, especially those currently orbiting beyond Saturn and less processed from a statistical point of view (Gkotsinas et al. 2022), would also correct the blind spot we currently have in our understanding of the evolution from the outer solar system to the JFC population.