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

Aurélie Guilbert-Lepoutre, Selma Benseguane, Laurine Martinien, Jérémie

Lasue, Sébastien Besse, Björn Grieger, Arnaud Beth

► To cite this version:

Aurélie Guilbert-Lepoutre, Selma Benseguane, Laurine Martinien, Jérémie Lasue, Sébastien Besse, et al.. Pits on Jupiter-family Comets and the Age of Cometary Surfaces. The Planetary Science Journal, 2023, 4, 10.3847/PSJ/ad083a. insu-04473131

HAL Id: insu-04473131

https://insu.hal.science/insu-04473131

Submitted on 23 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Pits on Jupiter-family Comets and the Age of Cometary Surfaces

Aurélie Guilbert-Lepoutre, Selma Benseguane, Laurine Martinien, Jérémie

Lasue, Sébastien Besse, Björn Grieger, Arnaud Beth

► To cite this version:

Aurélie Guilbert-Lepoutre, Selma Benseguane, Laurine Martinien, Jérémie Lasue, Sébastien Besse, et al.. Pits on Jupiter-family Comets and the Age of Cometary Surfaces. The Planetary Science Journal, 2023, 4, 10.3847/PSJ/ad083a. insu-04473131

HAL Id: insu-04473131

https://insu.hal.science/insu-04473131v1

Submitted on 23 Feb 2024

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
OPEN ACCESS

<figure><img src="image_2.png" /></figure>

CrossMark

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

Aurélie Guilbert-Lepoutre<fnref n="1" />, Selma Benseguane<fnref n="1" />, Laurine Martinien<fnref n="1" />, Jérémie Lasue<fnref n="2" />, Sébastien Besse<fnref n="3" />, Björn Grieger<fnref n="4" />, and Arnaud Beth<fnref n="5" />

Received 2023 August 29; revised 2023 October 13; accepted 2023 October 29; published 2023 November 20

## Abstract

Large and deep depressions, also known as pits, are observed at the surface of all Jupiter-family comets (JFCs) imaged by spacecraft missions. They offer the opportunity to glimpse the subsurface characteristics of comet nuclei and study the complex interplay between surface structures and cometary activity. This work investigates the evolution of pits at the surface of 81P/Wild 2, 9P/Tempel 1, and 103P/Hartley 2, in continuation of the work by Benseguane et al. on 67P/Churyumov-Gerasimenko. Pits are selected across the surface of each nucleus, and high-resolution shape models are used to compute the energy they receive. A thermal evolution model is applied to constrain how cometary activity sustained under current illumination conditions could modify them. Similar to what was found for 67P, we show that erosion resulting from water-driven activity is primarily controlled by seasonal patterns that are unique to each comet as a consequence of their shape and rotational properties. However, progressive erosion sustained after multiple perihelion passages is not able to carve any of the observed pits. Instead, cometary activity tends to erase sharp morphological features; they become wider and shallower over time. Our results reinforce the evolutionary sequence evidenced from independent measurables to transform “young” cometary surfaces, with sharp surface topography prone to outbursts, into “old” cometary surfaces. Finally, we suggest that the mechanism at the origin of the pits on JFCs should be able to carve these structures in a region of the solar system where water ice does not sublimate; the Centaur phase thus appears critical to understand JFC surface properties.

*Unified Astronomy Thesaurus concepts: Short period comets (1452); Comet nuclei (2160); Comet surfaces (2161); Theoretical models (2107)*

## 1. Introduction

The surfaces of comet nuclei display a diversity of morphological features (terraces, fractures, or boulders, for example; Massironi et al. 2015; El-Maarry et al. 2019) that have a complex interplay with cometary activity and the cycle of material across the nucleus (e.g., see Pajola et al. 2022, for a review). Surface depressions offer the opportunity to connect the subsurface properties of comet nuclei with the thermal physical processes actively shaping them (e.g., Vincent et al. 2015b; Benseguane et al. 2022; Davidsson et al. 2022). Shallow depressions (a few meters deep) generally observed on smooth terrains might be seasonal in nature, shaped by sublimation-driven activity orbit after orbit (Groussin et al. 2015; Vincent et al. 2016; El-Maarry et al. 2017; Birch et al. 2019; Bouquety et al. 2021; Davidsson et al. 2022). However, larger depressions, typically tens to several hundreds of meters deep, are also observed (Vincent et al. 2015a; El-Maarry et al. 2019) and cannot be linked with such seasonal activity (Benseguane et al. 2022). These structures, also known as pits, have been observed on all of the Jupiter-family comets (JFCs) imaged by spacecraft (Pajola et al. 2022). At the surface of 67P/Churyumov-Gerasimenko (hereafter 67P), they are mostly present on the northern hemisphere (Leon-Dasi et al.

2021). We recently studied their evolution through erosion as a result of thermally driven water-ice sublimation under current illumination conditions (Benseguane et al. 2022). We showed that through sustained cometary activity, erosion tends to erase such sharp morphological features. This result has implications for the origin of the pits, as the modeled pit evolution suggests that none of the structures observed at the surface of 67P could be formed through progressive erosion on a typical JFC orbit. This needs to be consolidated or refuted through the study of the evolution of similar structures seen on other comets.

## 2. Characteristics of the Pits Observed on Spacecraft Targets

### 2.1. 19P/Borrelly

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⁻<fnref n="1" />, compared to the average

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ LGL TPE, UMR 5276 CNRS, Université Lyon 1, ENS, Villeurbanne, France</footnote>

<footnote n="2">² IRAP, UMR 5277 CNRS, Université de Toulouse 3, CNES, Toulouse, France</footnote>

<footnote n="3">³ European Space Agency (ESA), European Space Astronomy Centre (ESAC), Villanueva de la Cañada, Spain</footnote>

<footnote n="4">⁴ Aurora Technology B.V. for the European Space Agency, ESAC, Madrid, Spain</footnote>

<footnote n="5">⁵ Department of Physics, Imperial College London, Prince Consort Road, London SW7 2AZ, UK</footnote>
OPEN ACCESS

<figure><img src="image_2.png" /></figure>

CrossMark

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

Aurélie Guilbert-Lepoutre<fnref n="1" />, Selma Benseguane<fnref n="1" />, Laurine Martinien<fnref n="1" />, Jérémie Lasue<fnref n="2" />, Sébastien Besse<fnref n="3" />, Björn Grieger<fnref n="4" />, and Arnaud Beth<fnref n="5" />

Received 2023 August 29; revised 2023 October 13; accepted 2023 October 29; published 2023 November 20

## Abstract

Large and deep depressions, also known as pits, are observed at the surface of all Jupiter-family comets (JFCs) imaged by spacecraft missions. They offer the opportunity to glimpse the subsurface characteristics of comet nuclei and study the complex interplay between surface structures and cometary activity. This work investigates the evolution of pits at the surface of 81P/Wild 2, 9P/Tempel 1, and 103P/Hartley 2, in continuation of the work by Benseguane et al. on 67P/Churyumov-Gerasimenko. Pits are selected across the surface of each nucleus, and high-resolution shape models are used to compute the energy they receive. A thermal evolution model is applied to constrain how cometary activity sustained under current illumination conditions could modify them. Similar to what was found for 67P, we show that erosion resulting from water-driven activity is primarily controlled by seasonal patterns that are unique to each comet as a consequence of their shape and rotational properties. However, progressive erosion sustained after multiple perihelion passages is not able to carve any of the observed pits. Instead, cometary activity tends to erase sharp morphological features; they become wider and shallower over time. Our results reinforce the evolutionary sequence evidenced from independent measurables to transform “young” cometary surfaces, with sharp surface topography prone to outbursts, into “old” cometary surfaces. Finally, we suggest that the mechanism at the origin of the pits on JFCs should be able to carve these structures in a region of the solar system where water ice does not sublimate; the Centaur phase thus appears critical to understand JFC surface properties.

*Unified Astronomy Thesaurus concepts: Short period comets (1452); Comet nuclei (2160); Comet surfaces (2161); Theoretical models (2107)*

## 1. Introduction

The surfaces of comet nuclei display a diversity of morphological features (terraces, fractures, or boulders, for example; Massironi et al. 2015; El-Maarry et al. 2019) that have a complex interplay with cometary activity and the cycle of material across the nucleus (e.g., see Pajola et al. 2022, for a review). Surface depressions offer the opportunity to connect the subsurface properties of comet nuclei with the thermal-physical processes actively shaping them (e.g., Vincent et al. 2015b; Benseguane et al. 2022; Davidsson et al. 2022). Shallow depressions (a few meters deep) generally observed on smooth terrains might be seasonal in nature, shaped by sublimation-driven activity orbit after orbit (Groussin et al. 2015; Vincent et al. 2016; El-Maarry et al. 2017; Birch et al. 2019; Bouquety et al. 2021; Davidsson et al. 2022). However, larger depressions, typically tens to several hundreds of meters deep, are also observed (Vincent et al. 2015a; El-Maarry et al. 2019) and cannot be linked with such seasonal activity (Benseguane et al. 2022). These structures, also known as pits, have been observed on all of the Jupiter-family comets (JFCs) imaged by spacecraft (Pajola et al. 2022). At the surface of 67P/Churyumov-Gerasimenko (hereafter 67P), they are mostly present on the northern hemisphere (Leon-Dasi et al.

2021). We recently studied their evolution through erosion as a result of thermally driven water-ice sublimation under current illumination conditions (Benseguane et al. 2022). We showed that through sustained cometary activity, erosion tends to erase such sharp morphological features. This result has implications for the origin of the pits, as the modeled pit evolution suggests that none of the structures observed at the surface of 67P could be formed through progressive erosion on a typical JFC orbit. This needs to be consolidated or refuted through the study of the evolution of similar structures seen on other comets.

## 2. Characteristics of the Pits Observed on Spacecraft Targets

### 2.1. 19P/Borrelly

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⁻<fnref n="1" />, compared to the average

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ LGL TPE, UMR 5276 CNRS, Université Lyon 1, ENS, Villeurbanne, France</footnote>

<footnote n="2">² IRAP, UMR 5277 CNRS, Université de Toulouse 3, CNES, Toulouse, France</footnote>

<footnote n="3">³ European Space Agency (ESA), European Space Astronomy Centre (ESAC), Villanueva de la Cañada, Spain</footnote>

<footnote n="4">⁴ Aurora Technology B.V. for the European Space Agency, ESAC, Madrid, Spain</footnote>

<footnote n="5">⁵ Department of Physics, Imperial College London, Prince Consort Road, London SW7 2AZ, UK</footnote>
10–15 m pixel⁻¹ 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.

## 2.2. 81P/Wild 2

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).

## 2.3. 9P/Tempel 1

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.

## 2.4. 103P/Hartley 2

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.

## 2.5. Ensemble Properties

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.

# 3. Methods

## 3.1. Shape Models for Comets 81P, 9P, and 103P

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),<fnref n="6" />

9P—derived from the images obtained by the Deep Impact and Stardust missions by Farnham & Thomas (2013a),<fnref n="7" />

103P—derived from images obtained by the EPOXI mission by Farnham & Thomas (2013b).<fnref n="8" />

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

<footnote n="6">⁶ https://pdssbn.astro.umd.edu/holdings/sdu-c-navcam-5-wild2-shape-model-v2.1/dataset.shtml</footnote>

<footnote n="7">⁷ https://pdssbn.astro.umd.edu/holdings/dif-c-hriv_its_mri-5-tempel1-shape-v2.0/dataset.shtml</footnote>

<footnote n="8">⁸ https://pdssbn.astro.umd.edu/holdings/dif-c-hriv_mri-5-hartley2-shape-v1.0/dataset.shtml</footnote>
10–15 m pixel⁻¹ 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.

## 2.2. 81P/Wild 2

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).

## 2.3. 9P/Tempel 1

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.

## 2.4. 103P/Hartley 2

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.

## 2.5. Ensemble Properties

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.

# 3. Methods

## 3.1. Shape Models for Comets 81P, 9P, and 103P

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),<fnref n="6" />

9P—derived from the images obtained by the Deep Impact and Stardust missions by Farnham & Thomas (2013a),<fnref n="7" />

103P—derived from images obtained by the EPOXI mission by Farnham & Thomas (2013b).<fnref n="8" />

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

<footnote n="6">⁶ https://pdssbn.astro.umd.edu/holdings/sdu-c-navcam-5-wild2-shape-model-v2.1/dataset.shtml</footnote>

<footnote n="7">⁷ https://pdssbn.astro.umd.edu/holdings/dif-c-hriv_its_mri-5-tempel1-shape-v2.0/dataset.shtml</footnote>

<footnote n="8">⁸ https://pdssbn.astro.umd.edu/holdings/dif-c-hriv_mri-5-hartley2-shape-v1.0/dataset.shtml</footnote>
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.

## 3.2. Selection of Pits on Each Comet

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.

## 3.3. Surface Energy and Thermal Evolution Model

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_\odot$ (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

$$E_{\odot} = \frac{F_{\odot}}{r_H^2} \cos \xi, \qquad (1)$$

$$E_{VIS} = \sum_{T} A_{T} \frac{F_{\odot}}{r_{H}^{2}} \cos \xi_{T} \frac{S_{T}}{\pi} \frac{\cos \zeta_{T} \cos \zeta_{R}}{\delta_{T}^{2}}, \qquad (2)$$

$$E_{IR} = \sum_{T} \varepsilon \sigma T_{T}^{4} \frac{S_{T}}{\pi} \frac{\cos \zeta_{T} \cos \zeta_{R}}{\delta_{T}^{2}}, \qquad (3)$$

where $F_\odot$ [W m⁻²] is the solar flux at 1 au, $A_T$ is the Bond albedo of an emitting facet ($A_T = 0.06$ for all calculations), $\xi_T$ is its local zenith angle, $S_T$ is its surface, $\zeta_T$ is the angle between the normal of the transmitter and the receiver facets, $\zeta_R$ is the angle between the normal of the receiving and emitting facets, and $\delta_T$ is the distance between the two facets. The value of the emissivity $\varepsilon$ 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 $\xi$ 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

$$(1 - A_R) \mathcal{E} = \varepsilon \sigma T^4 + \kappa \frac{\partial T}{\partial r} + f_{H_2O} \Delta H_{H_2O} Q_{H_2O}, \qquad (4)$$

where $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_{VIS} + E_{IR}$ is the total energy received at its surface, $T$ [K] is the surface equilibrium temperature, $f_{H_2O}$ is the fraction of the facet covered by water ice, $\Delta H_{H_2O}$ is the latent heat of water-ice sublimation, and $Q_{H_2O}$ [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
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.

## 3.2. Selection of Pits on Each Comet

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.

## 3.3. Surface Energy and Thermal Evolution Model

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_\odot$ (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

$$E_{\odot} = \frac{F_{\odot}}{r_H^2} \cos \xi, \qquad (1)$$

$$E_{VIS} = \sum_{T} A_{T} \frac{F_{\odot}}{r_{H}^{2}} \cos \xi_{T} \frac{S_{T}}{\pi} \frac{\cos \zeta_{T} \cos \zeta_{R}}{\delta_{T}^{2}}, \qquad (2)$$

$$E_{IR} = \sum_{T} \varepsilon \sigma T_{T}^{4} \frac{S_{T}}{\pi} \frac{\cos \zeta_{T} \cos \zeta_{R}}{\delta_{T}^{2}}, \qquad (3)$$

where $F_\odot$ [W m⁻²] is the solar flux at 1 au, $A_T$ is the Bond albedo of an emitting facet ($A_T = 0.06$ for all calculations), $\xi_T$ is its local zenith angle, $S_T$ is its surface, $\zeta_T$ is the angle between the normal of the transmitter and the receiver facets, $\zeta_R$ is the angle between the normal of the receiving and emitting facets, and $\delta_T$ is the distance between the two facets. The value of the emissivity $\varepsilon$ 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 $\xi$ 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

$$(1 - A_R) \mathcal{E} = \varepsilon \sigma T^4 + \kappa \frac{\partial T}{\partial r} + f_{H_2O} \Delta H_{H_2O} Q_{H_2O}, \qquad (4)$$

where $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_{VIS} + E_{IR}$ is the total energy received at its surface, $T$ [K] is the surface equilibrium temperature, $f_{H_2O}$ is the fraction of the facet covered by water ice, $\Delta H_{H_2O}$ is the latent heat of water-ice sublimation, and $Q_{H_2O}$ [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).

### 3.4. Orbital Considerations for Each Comet

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,<fnref n="9" />,<fnref n="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.

## 4. Thermal Processing of Pits

### 4.1. 81P/Wild 2

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

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Left: erosion of each selected facet on the surface of 81P as a function of latitude and energy peak. Right: latitude of the subsolar point as a function of heliocentric distance.</figcaption></figure>

<table><caption>**Table 1**
Summary of Results, Including Data for 67P from Benseguane et al. (2022): Number of Orbital Revolutions, Perihelion Distance *q*, Maximum Erosion, and Comparison to the Pits' Diameter and Depth</caption><thead><tr><th>Comet</th><th>No. of Orbits</th><th><i>q</i><br>(au)</th><th>Max.<br>Erosion<br>(m)</th><th>Pit<br>Diameter<br>(m)</th><th>Pit<br>Depth<br>(m)</th></tr></thead><tbody><tr><td>67P</td><td>10</td><td>1.24</td><td>77.23</td><td>100s</td><td>10s–100s</td></tr><tr><td>81P</td><td>6</td><td>1.59</td><td>28.25<sup>a</sup></td><td>100s–1000s</td><td>10s–100s</td></tr><tr><td>9P</td><td>13</td><td>1.54</td><td>83.81</td><td>10s–100s</td><td>10s</td></tr><tr><td>103P</td><td>20</td><td>1.06</td><td>265.10</td><td>10s–100s</td><td>10s</td></tr></tbody></table>

**Note.**

ª For 81P, the maximum is given for facets corresponding to observed regions of the surface, excluding the southern hemisphere.

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.

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

<footnote n="9">⁹ https://wgc.jpl.nasa.gov:8443/webgeocalc/#NewCalculation</footnote>

<footnote n="10">¹⁰ http://spice.esac.esa.int/webgeocalc/#NewCalculation</footnote>
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).

### 3.4. Orbital Considerations for Each Comet

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,<fnref n="9" />,<fnref n="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.

## 4. Thermal Processing of Pits

### 4.1. 81P/Wild 2

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

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Left: erosion of each selected facet on the surface of 81P as a function of latitude and energy peak. Right: latitude of the subsolar point as a function of heliocentric distance.</figcaption></figure>

<table><caption>**Table 1**
Summary of Results, Including Data for 67P from Benseguane et al. (2022): Number of Orbital Revolutions, Perihelion Distance *q*, Maximum Erosion, and Comparison to the Pits' Diameter and Depth</caption><thead><tr><th>Comet</th><th>No. of Orbits</th><th><i>q</i><br>(au)</th><th>Max.<br>Erosion<br>(m)</th><th>Pit<br>Diameter<br>(m)</th><th>Pit<br>Depth<br>(m)</th></tr></thead><tbody><tr><td>67P</td><td>10</td><td>1.24</td><td>77.23</td><td>100s</td><td>10s–100s</td></tr><tr><td>81P</td><td>6</td><td>1.59</td><td>28.25<sup>a</sup></td><td>100s–1000s</td><td>10s–100s</td></tr><tr><td>9P</td><td>13</td><td>1.54</td><td>83.81</td><td>10s–100s</td><td>10s</td></tr><tr><td>103P</td><td>20</td><td>1.06</td><td>265.10</td><td>10s–100s</td><td>10s</td></tr></tbody></table>

**Note.**

ª For 81P, the maximum is given for facets corresponding to observed regions of the surface, excluding the southern hemisphere.

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.

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

<footnote n="9">⁹ https://wgc.jpl.nasa.gov:8443/webgeocalc/#NewCalculation</footnote>

<footnote n="10">¹⁰ http://spice.esac.esa.int/webgeocalc/#NewCalculation</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Left: erosion of each selected facet on the surface of 9P as a function of latitude and energy peak. Right: latitude of the subsolar point as a function of heliocentric distance.</figcaption></figure>

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.

## 4.2. 9P/Tempel 1

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).

## 4.3. 103P/Hartley 2

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

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Left: erosion of each selected facet on the surface of 103P as a function of latitude and energy peak. Right: latitude of the subsolar point as a function of heliocentric distance.</figcaption></figure>

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.

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
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Left: erosion of each selected facet on the surface of 9P as a function of latitude and energy peak. Right: latitude of the subsolar point as a function of heliocentric distance.</figcaption></figure>

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.

## 4.2. 9P/Tempel 1

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).

## 4.3. 103P/Hartley 2

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

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Left: erosion of each selected facet on the surface of 103P as a function of latitude and energy peak. Right: latitude of the subsolar point as a function of heliocentric distance.</figcaption></figure>

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.

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
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Erosion per orbit as a function of total energy integrated over one orbit. The color code and the increasing size of the symbols give the (increasing) peak energy, usually received close to perihelion.</figcaption></figure>

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).

# 5. Discussion

## 5.1. Erosion and the Evolution of Pits

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.

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.

## 5.2. The Case of 103P/Hartley 2

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
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Erosion per orbit as a function of total energy integrated over one orbit. The color code and the increasing size of the symbols give the (increasing) peak energy, usually received close to perihelion.</figcaption></figure>

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).

# 5. Discussion

## 5.1. Erosion and the Evolution of Pits

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.

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.

## 5.2. The Case of 103P/Hartley 2

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.

### 5.3. Implications for the Aging of Cometary Surfaces

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
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.

### 5.3. Implications for the Aging of Cometary Surfaces

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
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Illustration of the evolutionary sequence between young (81P) and old (103P) cometary surfaces. Total erosion is calculated for each facet after a given number of orbital revolutions unique to each comet nucleus (10 for 67P, 13 for 9P, 6 for 81P, and 20 for 103P) to reflect their evolution under current illumination conditions. The color code (and increasing symbol size) gives the (increasing) peak energy received by each facet, typically around perihelion. Sequences suggested from the cumulative distribution of surface roughness and topography by Vincent et al. (2017) and the rate of mini-outbursts by Kelley et al. (2021; excluding 81P) are overlaid; the vertical spacing of the comets on these corresponding scales is qualitative, to serve an illustrative purpose.</figcaption></figure>

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.

## 5.4. On the Origin of Pits

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
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Illustration of the evolutionary sequence between young (81P) and old (103P) cometary surfaces. Total erosion is calculated for each facet after a given number of orbital revolutions unique to each comet nucleus (10 for 67P, 13 for 9P, 6 for 81P, and 20 for 103P) to reflect their evolution under current illumination conditions. The color code (and increasing symbol size) gives the (increasing) peak energy received by each facet, typically around perihelion. Sequences suggested from the cumulative distribution of surface roughness and topography by Vincent et al. (2017) and the rate of mini-outbursts by Kelley et al. (2021; excluding 81P) are overlaid; the vertical spacing of the comets on these corresponding scales is qualitative, to serve an illustrative purpose.</figcaption></figure>

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.

### 5.4. On the Origin of Pits

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).

## 5.5. Perspectives

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.

## 6. Summary

We investigate the evolution of pits on the surface of JFCs visited by space missions, i.e., 81P/Wild 2, 9P/Tempel 1, and

103P/Hartley 2, by applying the same method as for 67P/Churyumov-Gerasimenko (Benseguane et al. 2022). On each comet shape model, we select facets to sample at least 10 pits across the surface, distributed at all latitudes. The energy balance at the surface is then computed by including shadowing and self-heating contributions and used as a boundary condition of a 1D thermal evolution model to quantify the amount of erosion sustained after a number of orbital revolutions. This number is selected for each comet to correspond to the number of perihelion passages on the current orbit: 6 orbits for 81P, 13 orbits for 9P, and 20 orbits for 103P. We find the following.

1. Similar to what was found for 67P (Benseguane et al. 2022), erosion resulting from water-driven activity is primarily controlled by direct insolation. Strong seasonal patterns thus arise. However, our results suggest that erosion depends not only on the peak and total energy the surface 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.

2. Progressive erosion sustained after multiple perihelion passages is not able to carve large depressions of the observed size and shape on any of the comets we studied.

3. Cometary activity tends to erase sharp morphological features; they become wider and shallower over time.

4. Because the same patterns hold for four comet nuclei, our results can reinforce the evolutionary sequence evidenced from independent measurables, such as surface topography and roughness (Vincent et al. 2017), the phase function coefficient (Kokotanekova et al. 2018), or the rate of observed mini-outbursts (Kelley et al. 2021), that transform “young” cometary surfaces with sharp surface topography prone to outbursts into “old” cometary surfaces.

5. We suggest that the mechanism at the origin of the pits on JFCs should be able to carve these features in a region of the solar system where water ice does not sublimate; the Centaur phase of JFCs thus appears critical to understanding their surface properties.

## Acknowledgments

We warmly thank Marc Costa Sitja, Alfredo Escalante Lopez, Dave Schleicher, and Boris Semenov for their assistance with the SPICE kernels, especially relating to 103P/Hartley 2. We also thank the anonymous referee, whose comments helped improve this manuscript. This work benefited from discussions with the ISSI team #504 “The life cycle of comets” led by Rosita Kokotanekova. This project 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. A.B. was supported by grants ST/S000364/1 and ST/W001071/1.

## ORCID iDs

Aurélie Guilbert-Lepoutre https://orcid.org/0000-0003-2354-0766

Selma Benseguane https://orcid.org/0009-0006-6704-9388
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).

## 5.5. Perspectives

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.

## 6. Summary

We investigate the evolution of pits on the surface of JFCs visited by space missions, i.e., 81P/Wild 2, 9P/Tempel 1, and

103P/Hartley 2, by applying the same method as for 67P/Churyumov-Gerasimenko (Benseguane et al. 2022). On each comet shape model, we select facets to sample at least 10 pits across the surface, distributed at all latitudes. The energy balance at the surface is then computed by including shadowing and self-heating contributions and used as a boundary condition of a 1D thermal evolution model to quantify the amount of erosion sustained after a number of orbital revolutions. This number is selected for each comet to correspond to the number of perihelion passages on the current orbit: 6 orbits for 81P, 13 orbits for 9P, and 20 orbits for 103P. We find the following.

1. Similar to what was found for 67P (Benseguane et al. 2022), erosion resulting from water-driven activity is primarily controlled by direct insolation. Strong seasonal patterns thus arise. However, our results suggest that erosion depends not only on the peak and total energy the surface 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.

2. Progressive erosion sustained after multiple perihelion passages is not able to carve large depressions of the observed size and shape on any of the comets we studied.

3. Cometary activity tends to erase sharp morphological features; they become wider and shallower over time.

4. Because the same patterns hold for four comet nuclei, our results can reinforce the evolutionary sequence evidenced from independent measurables, such as surface topography and roughness (Vincent et al. 2017), the phase function coefficient (Kokotanekova et al. 2018), or the rate of observed mini-outbursts (Kelley et al. 2021), that transform “young” cometary surfaces with sharp surface topography prone to outbursts into “old” cometary surfaces.

5. We suggest that the mechanism at the origin of the pits on JFCs should be able to carve these features in a region of the solar system where water ice does not sublimate; the Centaur phase of JFCs thus appears critical to understanding their surface properties.

## Acknowledgments

We warmly thank Marc Costa Sitja, Alfredo Escalante Lopez, Dave Schleicher, and Boris Semenov for their assistance with the SPICE kernels, especially relating to 103P/Hartley 2. We also thank the anonymous referee, whose comments helped improve this manuscript. This work benefited from discussions with the ISSI team #504 “The life cycle of comets” led by Rosita Kokotanekova. This project 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. A.B. was supported by grants ST/S000364/1 and ST/W001071/1.

## ORCID iDs

Aurélie Guilbert-Lepoutre https://orcid.org/0000-0003-2354-0766

Selma Benseguane https://orcid.org/0009-0006-6704-9388
Jérémie Lasue https://orcid.org/0000-0001-9082-4457

Arnaud Beth https://orcid.org/0000-0002-5644-2022

## References

A'Hearn, M. F., Belton, M. J., Delamere, W. A., et al. 2011, *Sci*, **332**, 1396

Belton, M. J. S., Thomas, P., Li, J.-Y., et al. 2013, *Icar*, **222**, 595

Benseguane, S., Guilbert-Lepoutre, A., Lasue, J., et al. 2022, *A&A*, **668**, A132

Betzler, A. S. 2023, *MNRAS*, **523**, 3678

Birch, S. P. D., Hayes, A. G., Umurhan, O. M., et al. 2019, *GeoRL*, **46**, 12794

Bouquety, A., Jorda, L., Groussin, O., et al. 2021, *A&A*, **649**, A82

Brownlee, D. E., Horz, F., Newburn, R. L., et al. 2004, *Sci*, **304**, 1764

Bruck Syal, M., Schultz, P. H., Sunshine, J. M., et al. 2013, *Icar*, **222**, 610

Cheng, A. F., Lisse, C. M., & A'Hearn, M. 2013, *Icar*, **222**, 808

Clements, T. D., & Fernandez, Y. 2021, *AJ*, **161**, 73

Davidsson, B. J. 2021, *MNRAS*, **505**, 5654

Davidsson, B. J., Schloerb, F. P., Fornasier, S., et al. 2022, *MNRAS*, **516**, 6009

El-Maarry, M. R., Groussin, O., Keller, H. U., et al. 2019, *SSRv*, **215**, 36

El-Maarry, M. R., Groussin, O., Thomas, N., et al. 2017, *Sci*, **355**, 1392

Farnham, T., Duxbury, T., & L. J.-Y. 2005, SHAPE MODELS OF COMET WILD 2, SDU-C-NAVCAM-5-WILD2-SHAPE-MODEL-V2.1, NASA Planetary Data System, https://pds.nasa.gov/ds-view/pds/viewDataset.jsp?dsid=SDU-C-NAVCAM-5-WILD2-SHAPE-MODEL-V2.1

Farnham, T., & Thomas, P. 2013a, PLATE SHAPE MODEL OF COMET 9P/TEMPEL 1 V2.0, DIF-C-HRIV/ITS/MRI-5-TEMPEL1-SHAPE-MODEL-V2.0, NASA Planetary Data System, https://pds.nasa.gov/ds-view/pds/viewDataset.jsp?dsid=DIF-C-HRIV%2FITS%2FMRI-5-TEMPEL1-SHAPE-V2.0

Farnham, T. L., Bodewits, D., Li, J.-Y., et al. 2013, *Icar*, **222**, 540

Farnham, T. L., & Thomas, P. C. 2013b, Plate Shape Model of Comet 103P/Hartley 2 V1.0, DIF-C-HRIV/MRI-5-HARTLEY2-SHAPE-V1.0, NASA Planetary Data System, https://pds.nasa.gov/ds-view/pds/viewDataset.jsp?dsid=DIF-C-HRIV%2FMRI-5-HARTLEY2-SHAPE-V1.0

Gkotsinas, A., Guilbert-Lepoutre, A., Raymond, S. N., & Nesvorny, D. 2022, *ApJ*, **928**, 43

Groussin, O., Sierks, H., Barbieri, C., et al. 2015, *A&A*, **583**, A36

Guilbert-Lepoutre, A. 2012, *AJ*, **144**, 97

Guilbert-Lepoutre, A., Davidsson, B. J. R., Scheeres, D. J., & Ciarletti, V. 2023, Comets III (Tucson, AZ: Univ. Arizona Press), in press

Ip, W.-H., Lai, I.-L., Lee, J.-C., et al. 2016, *A&A*, **591**, A132

Kareta, T., Sharkey, B., Noonan, J., et al. 2019, *AJ*, **158**, 255

Kelley, M. S., Farnham, T. L., Li, J.-Y., et al. 2021, *PSJ*, **2**, 131

Kirk, R. L., Duxbury, T. C., Horz, F., et al. 2005, *LPSC*, **36**, 2244

Knight, M. M., Mueller, B. E. A., Samarasinha, N. H., & Schleicher, D. G. 2015, *AJ*, **150**, 22

Knight, M. M., & Schleicher, D. G. 2013, *Icar*, **222**, 691

Kokotanekova, R., Snodgrass, C., Lacerda, P., et al. 2017, *MNRAS*, **471**, 2974

Kokotanekova, R., Snodgrass, C., Lacerda, P., et al. 2018, *MNRAS*, **479**, 4665

Krolikowska, M., & Sitarski, G. 1996, *A&A*, **310**, 992

Lasue, J., De Sanctis, M. C., Coradini, A., et al. 2008, *P&SS*, **56**, 1977

Leon-Dasi, M., Besse, S., Grieger, B., & Küppers, M. 2021, *A&A*, **652**, A52

Lin, Z.-Y. 2023, *PASJ*, **75**, 462

Massironi, M., Simioni, E., Marzari, F., et al. 2015, *Natur*, **526**, 402

Miles, R. 2016, *Icar*, **272**, 356

Pajola, M., Höfner, S., Vincent, J. B., et al. 2017, *NatAs*, **1**, 0092

Pajola, M., Vincent, J.-B., El-Maarry, M. R., & Lucchetti, A. 2022, Comets III (Tucson, AZ: Univ. Arizona Press), in press

Pozuelos, F. J., Moreno, F., Aceituno, F., et al. 2014, *A&A*, **571**, A64

Rousselot, P., Korsun, P. P., Kulyk, I., Guilbert-Lepoutre, A., & Petit, J.-M. 2016, *MNRAS*, **462**, S432

Rousselot, P., Kryszczyńska, A., Bartczak, P., et al. 2021, *MNRAS*, **507**, 3444

Snodgrass, C., & Jones, G. H. 2019, *NatCo*, **10**, 5418

Soderblom, L. A., Becker, T. L., Bennett, G., et al. 2002, *Sci*, **296**, 1087

Thomas, P., A'Hearn, M., Belton, M. J. S., et al. 2013a, *Icar*, **222**, 453

Thomas, P. C., A'Hearn, M. F., Veverka, J., et al. 2013b, *Icar*, **222**, 550

Thomas, P. C., Veverka, J., Belton, M. J., et al. 2007, *Icar*, **187**, 4

Vincent, J.-B., Bodewits, D., Besse, S., et al. 2015a, *Natur*, **523**, 63

Vincent, J.-B., Farnham, T., Kührt, E., et al. 2019, *SSRv*, **215**, 30

Vincent, J.-B., Hviid, S. F., Mottola, S., et al. 2017, *MNRAS*, **469**, S329

Vincent, J.-B., Oklay, N., Marchi, S., Höfner, S., & Sierks, H. 2015b, *P&SS*, **107**, 53

Vincent, J.-B., Oklay, N., Pajola, M., et al. 2016, *A&A*, **587**, A14

Wierzchos, K., & Womack, M. 2020, *AJ*, **159**, 136
Jérémie Lasue https://orcid.org/0000-0001-9082-4457

Arnaud Beth https://orcid.org/0000-0002-5644-2022

## References

A'Hearn, M. F., Belton, M. J., Delamere, W. A., et al. 2011, *Sci*, **332**, 1396

Belton, M. J. S., Thomas, P., Li, J.-Y., et al. 2013, *Icar*, **222**, 595

Benseguane, S., Guilbert-Lepoutre, A., Lasue, J., et al. 2022, *A&A*, **668**, A132

Betzler, A. S. 2023, *MNRAS*, **523**, 3678

Birch, S. P. D., Hayes, A. G., Umurhan, O. M., et al. 2019, *GeoRL*, **46**, 12794

Bouquety, A., Jorda, L., Groussin, O., et al. 2021, *A&A*, **649**, A82

Brownlee, D. E., Horz, F., Newburn, R. L., et al. 2004, *Sci*, **304**, 1764

Bruck Syal, M., Schultz, P. H., Sunshine, J. M., et al. 2013, *Icar*, **222**, 610

Cheng, A. F., Lisse, C. M., & A'Hearn, M. 2013, *Icar*, **222**, 808

Clements, T. D., & Fernandez, Y. 2021, *AJ*, **161**, 73

Davidsson, B. J. 2021, *MNRAS*, **505**, 5654

Davidsson, B. J., Schloerb, F. P., Fornasier, S., et al. 2022, *MNRAS*, **516**, 6009

El-Maarry, M. R., Groussin, O., Keller, H. U., et al. 2019, *SSRv*, **215**, 36

El-Maarry, M. R., Groussin, O., Thomas, N., et al. 2017, *Sci*, **355**, 1392

Farnham, T., Duxbury, T., & L. J.-Y. 2005, SHAPE MODELS OF COMET WILD 2, SDU-C-NAVCAM-5-WILD2-SHAPE-MODEL-V2.1, NASA Planetary Data System, https://pds.nasa.gov/ds-view/pds/viewDataset.jsp?dsid=SDU-C-NAVCAM-5-WILD2-SHAPE-MODEL-V2.1

Farnham, T., & Thomas, P. 2013a, PLATE SHAPE MODEL OF COMET 9P/TEMPEL 1 V2.0, DIF-C-HRIV/ITS/MRI-5-TEMPEL1-SHAPE-MODEL-V2.0, NASA Planetary Data System, https://pds.nasa.gov/ds-view/pds/viewDataset.jsp?dsid=DIF-C-HRIV%2FITS%2FMRI-5-TEMPEL1-SHAPE-V2.0

Farnham, T. L., Bodewits, D., Li, J.-Y., et al. 2013, *Icar*, **222**, 540

Farnham, T. L., & Thomas, P. C. 2013b, Plate Shape Model of Comet 103P/Hartley 2 V1.0, DIF-C-HRIV/MRI-5-HARTLEY2-SHAPE-V1.0, NASA Planetary Data System, https://pds.nasa.gov/ds-view/pds/viewDataset.jsp?dsid=DIF-C-HRIV%2FMRI-5-HARTLEY2-SHAPE-V1.0

Gkotsinas, A., Guilbert-Lepoutre, A., Raymond, S. N., & Nesvorny, D. 2022, *ApJ*, **928**, 43

Groussin, O., Sierks, H., Barbieri, C., et al. 2015, *A&A*, **583**, A36

Guilbert-Lepoutre, A. 2012, *AJ*, **144**, 97

Guilbert-Lepoutre, A., Davidsson, B. J. R., Scheeres, D. J., & Ciarletti, V. 2023, Comets III (Tucson, AZ: Univ. Arizona Press), in press

Ip, W.-H., Lai, I.-L., Lee, J.-C., et al. 2016, *A&A*, **591**, A132

Kareta, T., Sharkey, B., Noonan, J., et al. 2019, *AJ*, **158**, 255

Kelley, M. S., Farnham, T. L., Li, J.-Y., et al. 2021, *PSJ*, **2**, 131

Kirk, R. L., Duxbury, T. C., Horz, F., et al. 2005, *LPSC*, **36**, 2244

Knight, M. M., Mueller, B. E. A., Samarasinha, N. H., & Schleicher, D. G. 2015, *AJ*, **150**, 22

Knight, M. M., & Schleicher, D. G. 2013, *Icar*, **222**, 691

Kokotanekova, R., Snodgrass, C., Lacerda, P., et al. 2017, *MNRAS*, **471**, 2974

Kokotanekova, R., Snodgrass, C., Lacerda, P., et al. 2018, *MNRAS*, **479**, 4665

Krolikowska, M., & Sitarski, G. 1996, *A&A*, **310**, 992

Lasue, J., De Sanctis, M. C., Coradini, A., et al. 2008, *P&SS*, **56**, 1977

Leon-Dasi, M., Besse, S., Grieger, B., & Küppers, M. 2021, *A&A*, **652**, A52

Lin, Z.-Y. 2023, *PASJ*, **75**, 462

Massironi, M., Simioni, E., Marzari, F., et al. 2015, *Natur*, **526**, 402

Miles, R. 2016, *Icar*, **272**, 356

Pajola, M., Höfner, S., Vincent, J. B., et al. 2017, *NatAs*, **1**, 0092

Pajola, M., Vincent, J.-B., El-Maarry, M. R., & Lucchetti, A. 2022, Comets III (Tucson, AZ: Univ. Arizona Press), in press

Pozuelos, F. J., Moreno, F., Aceituno, F., et al. 2014, *A&A*, **571**, A64

Rousselot, P., Korsun, P. P., Kulyk, I., Guilbert-Lepoutre, A., & Petit, J.-M. 2016, *MNRAS*, **462**, S432

Rousselot, P., Kryszczyńska, A., Bartczak, P., et al. 2021, *MNRAS*, **507**, 3444

Snodgrass, C., & Jones, G. H. 2019, *NatCo*, **10**, 5418

Soderblom, L. A., Becker, T. L., Bennett, G., et al. 2002, *Sci*, **296**, 1087

Thomas, P., A'Hearn, M., Belton, M. J. S., et al. 2013a, *Icar*, **222**, 453

Thomas, P. C., A'Hearn, M. F., Veverka, J., et al. 2013b, *Icar*, **222**, 550

Thomas, P. C., Veverka, J., Belton, M. J., et al. 2007, *Icar*, **187**, 4

Vincent, J.-B., Bodewits, D., Besse, S., et al. 2015a, *Natur*, **523**, 63

Vincent, J.-B., Farnham, T., Kührt, E., et al. 2019, *SSRv*, **215**, 30

Vincent, J.-B., Hviid, S. F., Mottola, S., et al. 2017, *MNRAS*, **469**, S329

Vincent, J.-B., Oklay, N., Marchi, S., Höfner, S., & Sierks, H. 2015b, *P&SS*, **107**, 53

Vincent, J.-B., Oklay, N., Pajola, M., et al. 2016, *A&A*, **587**, A14

Wierzchos, K., & Womack, M. 2020, *AJ*, **159**, 136