© A. Bouquety et al. 2022

1 Introduction

The Rosetta mission of the European Space Agency (ESA) provided detailed data of the surface of the nucleus of comet 67P/Churyumov-Gerasimenko (hereafter 67P). In particular, the Optical Spectroscopic and Infrared Remote Imaging System (OSIRIS; Keller et al. 2007) returned high-resolution images of the nucleus and its surface. The analysis of these images revealed an unexpected morphological diversity (Groussin et al. 2015b; Thomas et al. 2015b; El-Maarry et al. 2019), and the nucleus has been divided into 19 main geomorphological regions (El-Maarry et al. 2015; Thomas et al. 2015b) that are further divided into 71 subregions (Thomas et al. 2018). This large number shows the heterogeneity of 67P nucleus surface and reveals different geological processes and morphologies, depending on the terrain properties and composition. The terrains on 67P are divided into two categories: the consolidated terrains, and the nonconsolidated terrains (Thomas et al. 2015b; El-Maarry et al. 2019).

The consolidated terrains correspond to the exposed surface that appears indurated and rough. All terrains with a gravitational slopes exceeding 50° are consolidated (Groussin et al. 2015a), which includes all the cliffs. The erosion of cliffs induces landslides (Lucchetti et al. 2019), which are mass-wasting events called cliff collapses that create boulders (Pajola et al. 2015) at their feet. When they are located on terrains with a gravitational slope <50°, the cohesion of consolidated materials allows the formation of various morphologies, such as pits (Thomas et al. 2015b), fractures at different scales (El-Maarry et al. 2015; Matonti et al. 2019), strata (Massironi et al. 2015), and meter-size thermal contraction crack polygons (Auger et al. 2018). The presence of water ice in the subsurface and the increasing sublimation rate close to perihelion have been proposed to explain the observed cliff collapses, pits, or contraction crack polygons (Groussin et al. 2015b; Pajola et al. 2016; Vincent et al. 2016; El-Maarry et al. 2017; Hu et al. 2017; Birch et al. 2019).

The nonconsolidated terrains are characterized by a thin layer of materials that is a few meters thick and covers the underlying topography (Thomas et al. 2018). These units were called fine-particle deposit units (hereafter FPDs) by Lee et al. (2016) and Giacomini et al. (2016), and they are defined as units that are covered by a material with an apparently smooth and homogeneous texture at the meter scale. The origin of FPDs is the accumulation of materials that is caused by the seasonal redistribution of materials throughout the nucleus surface induced by activity (Keller et al. 2015; Thomas et al. 2015a). Several depressions have been observed on FPDs (Fig. 1). These depressions are topographic hollows with an arched shape that has a size of meters to tens of meters, and a depth of several meters (Groussin et al. 2015b; Pajola et al. 2016; Vincent et al. 2016; El-Maarry et al. 2017; Hu et al. 2017; Birch et al. 2019; Bouquety et al. 2021). Some of these depressions grew by several meters during perihelion approach; for example, depressions in the Hapi region grew and migrated by 11 m over 12 days (Birch et al. 2019). However, the origin of these depressions remains unclear, and two hypotheses are currently proposed. The first hypothesis is that they are located next to a cliff collapse and therefore correspond to a mass-wasting event (Pajola et al. 2016; Vincent et al. 2016; El-Maarry et al. 2017), as demonstrated by new boulders at the feet of some cliffs (Vincent et al. 2015; Bouquety et al. 2021). The second hypothesis is that they result from the erosion of dust-covered terrains, with dust dragged away by outgassing (Keller et al. 2015; Groussin et al. 2015b; Shi et al. 2016; Hu et al. 2017; Birch et al. 2019); this erosion process could be a transient event induced by the perihelion passage. For both hypotheses, the process at the origin of the depression formation remains to be fully understood.

In a previous study (Bouquety et al. 2021), we performed a comparative morphometrical analysis of two of these depressions that are located on a cliff edge in the Ash region. We found that these two depressions grew by several meters during perihelion passage, with two preferential growths: a first growth close to the cliff, associated with the apparition of new boulders at its foot, which we interpreted as triggered by a combination of water-ice sublimation and a landslide, and a second growth on the opposite side of the cliff, which we interpreted as controlled by water-ice sublimation alone. Again, we did not fully constrain the process at the origin of the growth of the two depressions, but only mentioned the association with a landslide and water-ice sublimation.

In this paper, we aim to understand the formation and evolution of depressions on 67P better. We apply our comparative morphometrical analysis to a larger sample of depressions (131 instead of 2) that are distributed throughout the nucleus of 67P. To better understand the mechanism responsible for their origin and evolution, and in particular, the role of water-ice sublimation and/or a landslide, we compared these depressions on 67P to possible analogs on Earth and Mars.

Fig. 1 NAC images of depressions in the Ma’at (a) and in the Anuket region (b).

2 Toward an interplanetary comparison

The geological diversity of telluric planets (the exogeodiversity; Bétard & Peulvast 2019), satellites, and small bodies shows that depression-like morphologies are common in the entire solar system. A depression is defined as a topographic hole whose bottom altitude is lower than that of the neighboring area, and whose slopes converge toward it (Hargitai & Kereszturi 2015). Even though these morphological characteristics are common to all depressions, the process that formed and shaped the different types of depression is strongly linked to the geological context. We list several examples of possible geological contexts below to illustrate their diversity (Hargitai & Kereszturi 2015):

1. Depressions in a volcanic context (a caldeira, a lava lake, etc.) are caused by the erosion induced by volcanic activity (Fig. 2a);

2. Impact craters are depressions resulting from the collision of an object with a surface (Fig. 2b);

3. Mass-wasting depressions are created by landslides and/or tectonic activity (rockfall, avalanches, etc.; Fig. 2c);

4. Glacial depressions such as glacial cirques or kettle holes are shaped by the direct or indirect erosion of a glacier (Fig. 2d);

5. Sinkholes are linked to the erosion of liquid water in a karstic terrain (Fig. 2e);

6. Finally, permafrost depressions (thermokarst depressions, thermokarst lakes, etc.) are formed by the thawing of permafrost, that is, a surface that remains at or below 0 °C for at least two subsequent years (Harris 1988; French 2007; Fig. 2f).

Fig. 2 Examples of planetary depressions in different contexts. (a) Okmok caldera (Alaska). (b) Martian impact crater in Lunae planum (Colour and Stereo Surface Imaging System, CaSSIS). (c) Landslide in Quebec. (d) Kettle holes in Wrangell-Saint-Élie national parc (Alaska). (e) Sinkhole in Russia. (f) Thermokarst lake in the Siberian permafrost. (g) Collapsed pingo (Canada). (h) Retrogressive thaw slump in the Siberian permafrost.

2.1 Cometary environment

To define the types of depression that might be observable on 67P, we must take the specificity of the cometary environment into account. First, 67P is a Jupiter-family comet (JFC) from the Kuiper Belt (Levison & Duncan 1997). It therefore formed and remained in a very cold environment (<60 K; Davidsson et al. 2016) for most of its lifetime, and was only recently (typically 100 kyr; Levison & Duncan 1997) injected into the inner Solar System, where it experienced higher insulation conditions. However, due to its low thermal inertia (Groussin et al. 2019), the interior did not experience significant temperature changes since its formation and remained very cold (Davidsson et al. 2016). Only the upper first meters are affected by thermal heating, with surface temperatures increasing up to 400 K at perihelion for 67P, but with nightside temperatures still below the sublimation temperature of water ice (<180 K). This very cold, icy, environment allowed the Visible, Infrared and Thermal Imaging Spectrometer (VIRTIS) instrument on board the Rosetta spacecraft to detect water ice (Capaccioni et al. 2015; Pommerol et al. 2015; Barucci et al. 2016) and even CO₂ ice (Filacchione et al. 2016) on its surface.

Next, 67P is an active body that releases dust and gas during its perihelion passage (Sierks et al. 2015; Vincent et al. 2016; Lai et al. 2019). This activity results from the sublimation of ices in the subsurface, mainly H₂O, CO₂, and CO (Bockelée-Morvan et al. 2016; Läuter et al. 2018; Biver et al. 2019; Haack et al. 2021). The occurrence of outgassing events increases when the comet approaches the Sun because the solar insulation increases. At each perihelion passage, the comet loses material (Pätzold et al. 2018), leading to the apparition of new morphological features or changes in existing ones (Groussin et al. 2015b; El-Maarry et al. 2017; Bouquety et al. 2021).

Because the ice lies close to the surface and because of the morphological changes induced by the sublimation of ice on 67P, we decided to focus on permafrost depressions. More precisely, we will compare the depressions on 67P with permafrost depressions on Earth and Mars, which are well documented in the literature.

2.2 Permafrost depressions on Earth and Mars

On Earth and Mars, permafrost is often associated with ice-related landforms. Several of these morphologies can be identified as depressions. In order to select the most appropriated ones, which will become our baseline for the comparison with 67P, we list here the terrestrial and Martian permafrost depressions identified in the literature, and compare their characteristics with the two depressions studied in Bouquety et al. (2021). Four main types of permafrost depressions exist on Earth and Mars (French 2007; Hargitai & Kereszturi 2015).

1. The collapsed pingo. A pingo is a perennially frozen ice core mound (nonglacial) that is primarily formed by the injection of water. It often shows radial cracking. The melting of the frozen ice core induces the slumping of the pingo sides and then its collapse, which appears as a circular or arcuate ridge of material surrounding the depression (Harris 1988; French 2007; Fig. 2g). We did not select the collapsed pingo for the comparison with 67P because the depressions on 67P are rimless and more elongated than circular, and moreover, we did not detect any pingo (mound) on 67P.

2. The retrogressive thaw slump (RTS; Fig. 2h). The RTS occurs only in thawing permafrost and creates a scarp through the thawing of the ice-rich permafrost. The fluid creates a debris lobe characterized by falls, slumps, slides, and flows. We did not select the RTS for the comparison with 67P because 67P lacks fluids and liquids, and moreover, the majority of depressions on 67P are not associated with deposits (Sect. 4).

3. The thermokarst depression. The terrestrial arctic lowlands landscape exhibits a multitude of depression-like forms, called thermokarst. The rising temperature induced by climate changes modifies the thermal equilibrium of the ground and leads to the melting of the ice in the permafrost (French 2007; Soare et al. 2008). Consequently, the sediments within the permafrost are no longer linked by the ice, which weakens the permafrost and ultimately induces its subsidence and the creation of a thermokarst depression. The permafrost water-ice volume must be 30–50% for the subsidence to occur (Costard & Kargel 1995; French 2007; Soare et al. 2008; Séjourné et al. 2011). On Earth, the melting of water ice causes liquid water to fill the depression, therefore it is called a thermokarst lake. When water evaporates from the lake, a dry depression remains that is called an alas (Fig. 3a). On Earth, alases can only be formed by sublimation in very peculiar areas (Fairbanks tunnel, Alaska; Douglas & Mellon 2019). On Mars, thermokasrt depressions have also been observed, and two hypotheses have been proposed for their formation, either the melting of ground-ice followed by its evaporation (Costard & Kargel 1995; Soare et al. 2008), or the sublimation of ground-ice (Morgenstern et al. 2007; Lefort et al. 2009). Since the current Martian environment does not allow stable liquid water at its surface, sublimation is the most accepted hypothesis today (Morgenstern et al. 2007; Ulrich et al. 2010; Séjourné et al. 2011, 2019; Costard et al. 2021). As on Earth, the sublimation of the permafrost water ice leads to a subsidence. On Mars, thermokarst depressions are called scalloped depressions, and they are analogs to thermokarst depressions on Earth (Costard & Kargel 1995; Morgenstern et al. 2007; Soare et al. 2008; Ulrich et al. 2010; Séjourné et al. 2011, 2019; Costard et al. 2021).

Finally, because of their similarities in terms of context, morphologies, and environment, we decided to compare the depressions on 67P with the alases (i.e., dry thermokarst depressions formed by melting and evaporation) on Earth and with the scalloped depressions (i.e., thermokarst depressions formed by sublimation) on Mars. The lack of liquid water in alases (in contrast to thermokarst lakes) makes it easier to obtain topographic information such as the depth. More precisely, we focused on alases in the Alaskan arctic coastal plain (ACP) on Earth and on scalloped depressions in the Utopia Planitia region on Mars. We selected these areas because data with similar spatial resolution are available that can be compared with data of 67P (images and topography). Moreover, the selected areas on Earth and Mars are representative of alases and scalloped depressions found on these planets.

2.3 Alases on Earth: The arctic coastal plain

The ACP exhibits tens of thousands of thermokarst lakes and alases formed in a thick (>10m) permafrost (Fig. 3a; French 2007; Arp et al. 2011; Hinkel et al. 2012; Ulrich et al. 2017). The alases have sizes in the range 0.1–30km, with an elongated to circular shape. They are characterized by inward terraces that formed by periodic evaporation or drainage, and by the absence of rims. They have a relatively flat floor and shallow depth (1 to 15 m), with steep walls (20° on average), and can present a slope asymmetry (French 2007; Arp et al. 2011; Hinkel et al. 2012; Ulrich et al. 2017).

2.4 Scalloped depressions on Mars: Utopia planitia

Utopia planitia is one of the three main topographic basins in the northern hemisphere of Mars (Tanaka et al. 2014). This area is covered by an ice-rich deposit that is tens of meter thick and forms a permafrost (Costard & Kargel 1995; Morgenstern et al. 2007; Soare et al. 2008; Séjourné et al. 2012, 2019). This area presents several periglacial morphologies, such as polygons or scalloped depressions (Fig. 3b, Séjourné et al. 2011, 2019). The scalloped depressions result from ground-ice degradation (Costard & Kargel 1995; Soare et al. 2008; Morgenstern et al. 2007; Lefort et al. 2009; Séjourné et al. 2011).

The scalloped depressions in Utopia planitia have a circular to elliptical shape and extend over meters to kilometers. They are up to several tens of meters deep (Costard & Kargel 1995; Lefort et al. 2009; Morgenstern et al. 2007; Séjourné et al. 2011; Soare et al. 2008; Ulrich et al. 2010). These depressions are characterized by an N-S asymmetric profile with curvilinear steps. The pole-facing slopes are steeper (5°–20°) than the equator-facing slopes (2°–5°) (Morgenstern et al. 2007; Levy et al. 2010; Séjourné et al. 2011). This asymmetry is the result of depressions shaped by the obliquity-driven insulation (Morgenstern et al. 2007; Séjourné et al. 2011).

Fig. 3 Examples of measured thermokarst depressions. (a) Alases in Alaska on Earth (digital orthophoto quadrangle, DOQ, 5mpixel−1). (b) Scalloped terrain in Utopia planitia on Mars (HiRISE image, 25 cm pixel−1). (c and e) Associated digital elevation nodels. (d and f) Slope maps. (g and h) Topographical profiles along the length direction (in blue in panels d and f). The pink lines show the width measurement.

3 Data and method

3.1 Comparative morphometrical analysis

In order to constrain the origin of the depressions on 67P, we performed a morphometrical comparison between these depressions, the alases on Earth, and the scalloped depressions on Mars. The main objective of this method is to identify geometrical and morphological similarities and differences between these morphologies. We already successfully applied this technique to distinguish between glacial and fluvial erosion on Mars (Bouquety et al. 2019, 2020) and to quantify the changes observed in two depressions on 67P (Bouquety et al. 2021).

The first step of this analysis is to establish the list of geometric parameters that can all be measured on the three studied bodies (67P, Earth, and Mars). We followed the standard approach of Ulrich et al. (2010); Séjourné et al. (2011); Morgenstern et al. (2011), and Niu et al. (2014), and selected ten geometric parameters for each depression. We list them in Table 1. This set of parameters provides important constraints to study the depressions morphology. Several parameters are related to the size of the depression (length, width, depth, perimeter, area). The elongation is related to the shape of the depression and to the circularity, which is a good proxy to evaluate the deviation of an object shape to a perfect circle (Morgenstern et al. 2011). Finally, the slope is a key parameter for the local context. We used the ArcGIS software to perform all the measurement.

To validate our morphometrical results, we used a statistical Mann–Whitney test for each parameter and/or depression. This rank-based test allows validating whether two studied populations belong (or do not belong) to the same statistical sample. The test is characterized by its p-value: a p-value <0.01 indicates that the two populations are statically different. This test is particularly relevant for comparing depressions on Earth, Mars, and 67P.

3.1.1 Earth and Mars

In the ACP region on Earth, we performed our measurements on alases using the USGC-digital orthophoto quadrangle (5 mpixel⁻¹, Fig. 3a) and the associated digital elevation model (DEM) from the 3D elevation program (5 m pixel⁻¹, Fig. 3c, Earth Resources Observation And Science (EROS) Center (2017).

In the Utopia planitia region of Mars, we performed our measurements on scalloped depressions using the High Resolution Imaging Science Experiment (25 cmpixel⁻¹, Fig. 3b; McEwen et al. 2007) on the Mars Reconnaissance Orbiter (MRO) spacecraft, and its associated DEM (Fig. 3e).

For Earth and Mars, we used the ArcGIS slope tool to construct the slope map from the DEM (Figs. 3d,f). These data are georeferenced, therefore we can directly extract the value of each geometric parameter of Table 1 from the shapefiles we drew in ArcGIS.

3.1.2 67P

To identify the depressions on 67P, we used the Narrow Angle Camera (NAC) images from the OSIRIS instrument (Keller et al. 2007). We established three criteria, based on Bouquety et al. (2021), to validate our identification:

1. The edges of the depression must be continuous and easy to follow over tens of meters;

2. The texture inside the edges must be different from that of the surrounding terrain. Usually, the texture inside a depression appears less rough than outside (Fig. 4a);

3. The depression must be visible in at least three NAC images with different illumination conditions to remove possible artifacts.

It was not always possible to follow the depression edges to form a closed shape. In these cases, we followed the edges as far as we could, and then closed the shapefile with a straight line.

From the thousands of NAC images acquired during the Rosetta mission, we selected 180 images that are distributed among all regions. We selected these images to cover a wide range of spatial resolution and provide the global (typically 1 mpixel⁻¹) and local (typically 0.1 mpixel⁻¹) context (see Fig. 1 in Bouquety et al. 2021).

The NAC data set offers both stereo coverage and varying illumination conditions that were used to derive a high-resolution DEM of the nucleus via the stereophotoclinometry method, or SPC (Gaskell et al. 2008; Jorda et al. 2016a). We preferred the SPC model over the photogrammetric model (Preusker et al. 2017) for this study because the gentle slopes of the depressions are well-suited for the SPC reconstruction method.

For each NAC image, we extracted the corresponding DEM (1 m pixel^-1), the spacecraft–comet distance (Fig. 4b), and the emission angle (Fig. 4c). Then, we calculated the gravitational heights (GH, Fig. 4d) and the gravitational slopes (GS, Fig. 4e) derived from the DEM.

In contrast to Earth and Mars, the images on 67P are not georeferenced. Thus, we cannot measure the depression geometry without first calculating the pixel scale of each image from the spacecraft position. Moreover, in a given image, the pixel scale varies from one pixel to the other due to local topographic variations. To accurately compute the pixel scale, we extracted the spacecraft-to-comet (S/C) distance and the emission angle of each pixel inside each depression using the ArcGIS extract-by-mask tool. Then, we calculated the mean distance d and emission angle e of each depression by averaging the value for each pixel. The final pixel scale was calculated as p ~ d × IFOV/cos e, where IFOV = 18.82 μrad for the NAC.

The uncertainty on each geometric parameter measurement is twofold. First, there is the uncertainty on the pixel scale, which is the root mean square (rms) value of d and 1/cos e for each depression (Jorda et al. 2016a; Bouquety et al. 2021). Next, there is an error coming from the operator using ArcGIS. We estimated this error to ±2 pixels for the length, width, and perimeter. The total error is the root sum squared of the pixel scale error and the operator error. Overall, the typical error of individual measurement (i.e., for each depression) is 4% for the length and width, 4% for the perimeter, 11% for the area, 2% for the elongation, and 10% for the circularity. From the accuracy of the DEM, we estimate the error on the height to be 0.1 m and the error on a gravitational slope to be 5 deg (Groussin et al. 2015a; Jorda et al. 2016a).

In contrast to Bouquety et al. (2021), we did not select images before and after perihelion for each depression, but only made our measurements at a single time. Following each depression over the duration of the Rosetta mission to track potential morphological changes goes beyond the scope of this paper. This is a limit of our method, but as shown by Bouquety et al. (2021), most changes are expected to be small, typically a few meters, which is comparable to our measurement error bar.

3.2 Thermal analysis

In order to study the thermal environment of each depression on 67P, we used a thermal model to compute the maximum temperature (K), the total erosion (m), and the total energy received from the Sun (Jm⁻²) in each depression over one revolution.

The thermal model is a surface model that takes the solar insulation, the surface thermal emission, and the sublimation of water ice into account. Cast shadows are accounted for in the model. The surface temperature was computed for two cases, first assuming no water ice to compute the maximum temperature, and then assuming pure water ice to compute the maximum erosion. The total energy received from the Sun was independent of the surface composition and identical for the two cases. The surface energy balance was described according to Groussin et al. (2004) for case one and case two by Eqs. (4) and (5), respectively, ignoring the η parameter (i.e., η = 1). While the temperature and energy received from the Sun are robust values that can be used quantitatively, the erosion is complex, depends on many parameters, and is highly variable across the surface. The calculated value of erosion should therefore only be used qualitatively to compare different areas, but not quantitatively.

The thermal model was applied to the shape model derived from the SPC method (Jorda et al. 2016a). We used the version with 12.6 million facets, which provides a resolution exceeding 1000 facets (on average) for each depression. To avoid computational time issue resulting from this large number of facets, the thermal model did not take heat conductivity into account, nor the self-heating resulting from facets seeing each other.

The thermal model take into account 19 different heliocentric distances over one revolution of 67P, distributed unevenly to ensure a constant decrease in solar energy. With this distribution of heliocentric distances, the time step increased from 30 days at perihelion to 300 days at aphelion. Smaller time steps would probably increase the accuracy of our calculations, but we are limited by CPU time for this 12.6 million facet shape model. At each heliocentric distance, the thermal model was run over one rotation of 67P, with a time step of 1/18th of the rotation period of 12.4 h. At each time step and for each facet, we computed the surface temperature and the solar energy (assuming no water ice) and the erosion (assuming pure water ice). Solar energy and erosion were integrated over one complete revolution. Figure 5 illustrates the results for one depression. We plot the maximum temperature (K), total erosion (m), and total energy received from the Sun (J m⁻²) in the area over one orbit.

To facilitate the thermal analysis and smooth the local thermal variations, we divided each depression into three units, corresponding to the bottom of the depression (the flat area inside the depression edges), the slope break, and the outer border of the depression (Figs. 5d, e). For each unit, we computed the average value of the maximum temperature, total erosion, and total energy received from the Sun.

Fig. 4 Dataset derived from the SPC shape model that we used to measure the geometric parameters of depressions on 67P. (a) NAC image of depressions located in the Ma’at region. (b) Distance between the camera and the surface and (c) emission angle, both draped on an NAC image. (d) Gravitational heights (GH) draped on an NAC image that we used to extract the topographical profile in panel f. (e) Gravitational slopes (GS) draped on an NAC image. (f) Topographic profile along the blue line in panel d, extracted from GH.

Fig. 5 Results from the thermal analysis of a depression in the Ma’at region. (a) Maximum temperature (K), (b) total erosion (m), (c) total energy received from the Sun (Jm−2), (d) depression subdivision into three units, and (e) topographical profile showing the subdivisions.

4 Results of the morphometrical analysis

4.1 Database of geometric parameters

We manually measured the geometric parameters of a total of 381 depressions: 200 scallop depressions on Mars, 50 alases on Earth, and 131 depressions on 67P. We extracted 10 geometric parameters for each depression, that is, a total of 3810 parameters. Table 2 summarizes the results, and Fig. 6 illustrates the statistical distribution of the measured parameters.

4.2 Results for 67P

In the 180 NAC images, we identified 131 depressions that met our criteria (Sect. 3.1.2). The first result is that the depressions are exclusively located in FPDs, without a preferential orientation. We calculated the density of a depression in each region by dividing the number of depressions by the region area (Fig. 7). The Seth, Imothep, Hapi, Aten, and Aswan regions have the highest depression density with more than 5 depressions km^-2 (Fig. 7). This result is consistent with previous studies that indicated a higher fraction of FPDs in these regions (Thomas et al. 2015b; Lee et al. 2016; Giacomini et al. 2016; El-Maarry et al. 2019). However, there is no general trend because the Babi and Bes regions, for example, have few depressions, althougha significant fraction of their surface is covered by FPDs.

The depressions present a wide range of morphologies, depending on their local topographical context (Figs. 8a,b), their local distribution (Figs. 8c, d), and their association (or disconnection) with neighboring structures (Figs. 8e, f):

1. Topographical context: 58% of the depressions are located in a flat terrain (mean slope 3.8°; Fig. 8a), where the slope is null or constant, that is, without a slope break. The other 42% are all located at a cliff edge (mean slope 11°; Fig. 8b), with a sharp break in slope.

2. Local distribution: 71% of the depressions are single, isolated (Fig. 8c), that is, without a neighboring depression within a distance of 100 m. Some of these isolated depressions are located in a small FPD unit and are surrounded by a terrain that is dominated by exposed consolidated materials. Conversely, 29% of the depressions are organized in groups, less than 100 m from each other, and form a depression field (Fig. 8d). In these fields, most of the depressions are smaller (<50 m in length and width) than the isolated depressions, and they seem to merge together to form a unique, larger depression (Fig. 8b).

3. Associated structures: Even if most (>86%) of the depressions have no associated structures, 6% of them are associated with thermal contraction crack polygons (Fig. 8e), 8% with boulders at the cliff feet (Fig. 8f), and none with both (polygons and boulders). The depressions associated with polygons are all isolated and always located in areas surrounded by exposed consolidated materials. The depressions associated with boulders are all located at a cliff edge; the boulders are located several meters down, at the cliff feet, and probably come from previous landslides triggered by cliff collapse.

The geometry of the depressions on 67P reveals a large diversity (Table 2 and Fig. 6). Depressions have a length of 55 ± 41 m, a width of 81 ± 63 m, and an area of 4731 ± 3002m², with a few extreme values. Remarkably, all the measured depressions follow the same perimeter-over-area trend, regardless of their topographic context, local distribution, or associated structures (Fig. 9). Moreover, small depressions can be deeper and/or steeper than large ones, that is, the size of a depression and its other geometric parameters are not correlated. This apparently random spatial distribution across the surface allows us to consider all the depressions as a single statistical sample, whose formation and evolution is likely driven by the same process. This process is not related to the topographic context, the local distribution, or the association with neighboring structures, but depends on additional factors such as possibly the FPD physic-ochemical properties and/or the thermal environment, as we discuss below.

Fig. 6 Box and whisker plot for all parameters on the three bodies. The black lines at each end correspond to the maximum and minimum value. The box represents the interquartile range; Q1 and Q3 are the bottom and top of this box. The median value is represented by the black lines inside the box. Abbreviations are shown in Table 2. Note the logarithmic scale in panel a and the linear scale in panel b.

Fig. 7 Density map of depressions in the different regions of the nucleus of 67P.

5 Interplanetary morphometrical comparison

5.1 Comparison between Earth and Mars: analogous depressions

Following Morgenstern et al. (2011), we plot in Fig. 10, the perimeter versus area for the alases and scalloped depressions. The scalloped depressions on Mars follow the same trend as alases on Earth. The scalloped depressions have a circularity index of 0.77 ± 0.15, similar to the circular index of 0.83 ± 0.10 for alases. As a comparison, we also add in Fig. 10 the perimeter versus area for 119 randomly measured Martian craters, which are known to be among the most circular features in the Solar system (C ~ 1). The trend followed by craters is very different from the trend of alases or scalloped depressions. Finally, we performed a Mann–Whitney statistical test. This confirmed that the alases on Earth and the scalloped depressions on Mars probably belong to the same statistical sample (Table 3), demonstrating the already known geological analogy between these two types of thermokarst depressions (Costard & Kargel 1995; Lefort et al. 2009; Morgenstern et al. 2007; Séjourné et al. 2011; Soare et al. 2008; Ulrich et al. 2010).

5.2 Comparison between 67P, Earth, and Mars

Now, we compare the database of geometric parameters built for the depressions on 67P, whose origin has to be determined, with those built for alases on Earth and scalloped depressions on Mars.

5.2.1 Size, area, perimeter, and circularity

Table 2 shows the comparison between the alases on Earth, the scalloped depressions on Mars, and the depressions on 67P. The depressions on 67P are ten times smaller than the alases and four times smaller than the scalloped depressions on average. Their circularity index is 28% smaller than that of the alases and 22% smaller than that of the scalloped depressions, but the statistical distributions of the circularity index overlap at 1σ. In a perimeter versus area plot as in Fig. 10, the depressions on 67P follow the same trend as alases and scalloped depressions. This trend differs from that followed by the Martian craters. These results indicate a morphological analogy between the depressions on 67P, the alases on Earth, and the scalloped depressions on Mars. This analogy is reinforced by a Mann–Whitney test. The test indicates that all these depressions could indeed come from the same population (Table 3), excepted for the Martian craters, as expected.

5.2.2 Slope asymmetry

As indicated in Sect. 2, the scalloped depressions on Mars are characterized by a slope asymmetry (Morgenstern et al. 2007; Séjourné et al. 2011). This characteristic is indeed visible in the scalloped terrains we selected on Mars (Fig. 11b) and also on the alases on Earth (Fig. 11a). Our morphometrical analysis reveals that the depressions on 67P also exhibit a slope asymmetry (Fig. 11c), with no preferential orientation. In fact, 98% of the measured depressions on 67P have a steeper slope upstream, from 6° to 19°, with a mean value of 14°, and a more gentle down-slope, from 0.2° to 4°, with a mean value of 2.4°. 2% of the depressions do not show this slope asymmetry. They are small depressions (<200m²), for which the spatial resolution of the gravitational slopes is insufficient to detect minor slope changes.

The mean slope values measured for the 67P depressions are comparable to those characterizing thermokarst depressions on Earth and Mars, which are in the range 5°−20° for the steep slope and 2°−5° for the gentle slope. The slope asymmetry observed in the 67P depressions is therefore similar to the asymmetry observed on alases and scalloped depressions. This additionally reinforces the morphological analogy between these three types of depression.

Fig. 8 NAC images showing different examples of depression context. The white arrows show the depressions. (a) Depression in a flat terrain in the Anuket region. (b) Depression at a cliff edge in the Aswan region. (c) Isolated depression in the Ash region. (d) Depression field in the Imhotep region. (e) Depressions associated with thermal contraction crack polygons in the Anhur region. (f) Depressions associated with boulders in the Imhotep region.

Fig. 9 Area vs. perimeter for the measured depressions on 67P. The differently colored squares correspond to the different depression contexts (see text for details).

Fig. 10 Area vs. perimeter for the thermokarst lake measured in the ACP (circles), the scalloped depressions (triangles), Martian craters (dia-monds),and depressions on 67P(square). The blue colorscale represents the circularityof the measured features.

Fig. 11 Topographic profiles showing the asymmetric slopes of the studied depressions. (a) Circular depression in the ACP on Earth. (b) Scalloped depression in Utopia planitia on Mars. (c) Depression in the Seth region of 67P. The dashed line indicates the location of the topographic profile, and the green arrows are markers.

Fig. 12 Depth value of the measured depressions and for other terrestrial thermkarstic depressions in Siberia (Ulrich et al. 2010), the Lena Delta (Morgenstern et al. 2011), and in the Qinghai-Tibet engineering corridor (Niu et al. 2014).

5.2.3 Depth

We compared the depression depth for 67P, Earth, and Mars (Fig. 12). The depressions on 67P are 30% shallower than the scalloped depressions on Mars on average, but the two statistical distributions widely overlap at 1σ. This is confirmed by the Maan-Whitney test. The test indicates that the statistical difference between the two populations is not significant, with a p-value of 0.18 (i.e., > 0.01).

Nevertheless, the maximum depression depth on 67P (>15 m) is three times larger than the maximum depth of alases on Earth (<5 m). To determine whether this difference is real or a bias resulting from the selected alases, we compared our results with other depth measurements of thermokarst depressions on Earth. Ulrich et al. (2010) studied the alases in the Siberian ice-rich permafrost and obtained a maximum value of 30 m for their depth. Morgenstern et al. (2011) measured 2327 thermokarst depressions in the Yedoma landscape of the Lena Delta and obtained a depth of between 4 m and 35 m. Finally, Niu etal. (2014) studied 2163 thermokarst depressions in the Qinghai-Tibet engineering corridor. Their depth lies between 0.71 m and 2.95 m. A wide range of depths for alases on Earth clearly exists. The depth depends on the location of the depression., Some depths are comparable to the depth measured for the 67P depressions.

Overall, this comparison reveals a similar depth of a few meters for the depressions on 67P and the scalloped depressions on Mars. This is comparable to the depth of some (but not all) alases on Earth.

5.2.4 Thermokarst depressions on 67P

Our morphometrical analysis reveals that the depressions on 67P share strong similarities with alases on Earth and scalloped depressions on Mars:

1. Their statistical distributions of circularity index overlap at 1σ;

2. In a perimeter versus area plot, the depressions on 67P follow the same trend as alases and scalloped depressions, and a Mann-Whitney test further confirms this point;

3. They all show a slope asymmetry along their length profiles, with a steep up-slope and a gentle down-slope. Moreover, the slope values are comparable;

4. The depressions are all a few meters deep.

Moreover, because they all share a common cold and icy environment, we consider alases on Earth and scalloped depressions on Mars as analogs to 67P depressions. This analogy allows us to constrain the erosive agent better that shaped the depressions on 67P, and to highlight the fact that the thermokarst is the main process. We propose to call this type of depression on cometary nuclei thermokarst depression. It was formed by thermokarst (i.e., permafrost thawing) triggered by sublimation. This analogy and its implications are discussed in detail in Sect. 7.

6 Thermal environment of the depressions

Our thermal model (Sect. 3.2) allowed us to compute the maximum temperature, the total erosion, and the total energy received from the Sun during one revolution for each depression. More precisely, we computed the average value of each of the above quantities for the three parts of each depression, that is, the bottom of the depression, the slope break, and the outer border (Figs. 5d, e). Figure 13 illustrates the statistical distribution of the measured thermal parameters for the different regions of 67P.

The regions in the southern hemisphere and from the equator have a higher average temperature (>340K) and stronger average erosion (>5 m) than those in the northern hemisphere (<300K, <3m). The only exception to this general trend is the Aten region, whose average temperature and erosion are similar to regions in the northern hemisphere, although Aten is located close to the equator. This exception simply reflects the fact that Aten is spread over a wide range of latitudes, from 0° to 45°N (El-Maarry et al. 2019), and that the depressions in this region are all located in its northern part and therefore correspond to a thermal environment in the northern hemisphere.

This north–south dichotomy is not observed for the total energy received from the Sun (Fig. 13c). This is caused by the seasonal variations on 67P. The northern hemisphere is mainly illuminated far from the Sun (at aphelion), that is, it receives less energy, but for a very long time (years), while the southern hemisphere is mainly illuminated close to the Sun (at perihelion), that is, it receives more energy, but for a far shorter time (months). Many regions also vary strongly because they are spread in latitude. The only exceptions are Anuket, with a high energy and small variations, and Anhur, with a low and almost constant energy. It reflects their peculiar position, which is well centered on the equator, where the energy is maximum for Anuket, and in the neck, where projected shadows often hide the Sun for Anhur.

In order to study the thermal behavior of the depressions in more detail, we display the percentage of the highest values for each parameter according to the area of the three depressions (Fig. 14). The slope break is the key location. Temperature, erosion, and energy received from the Sun reach their highest values for most (63–78%) of the depressions. In the area of the slope break, most of the changes induced by sublimation occur. Three quarters of these depressions are located in a flat terrain. The bottom of a depression has the highest values especially for the temperature (37%) and the erosion (22%). This is caused by an edge effect: the limit between the area of the slope break and the bottom is poorly defined (insufficient spatial resolution); the area of the slope break contaminates the bottom and artificially increases its maximum temperature and erosion.

Overall, this analysis demonstrates that the characteristic asymmetric slope that is observed in all cometary thermokarst depressions has a thermal origin. As on Mars, the slope break is the preferential erosional area in which sublimation occurs (Séjourné et al. 2011). Most depressions follow the same thermal behavior: their slope break occurs in the area of highest temperature and erosion. This allows us to conclude that they are potentially still active today.

7 Discussion

7.1 Comparison with other hypotheses

The subsidence hypothesis is a new possible scenario for the formation of depressions. In this subsection, we discuss this new scenario with respect to previously proposed scenarii, namely cliff collapse (Vincent et al. 2015; Pajola et al. 2016; El-Maarry et al. 2017; Bouquety et al. 2021) and dust ejection dragged by outgassing (Keller et al. 2015; Groussin et al. 2015b; Shi et al. 2016; Hu et al. 2017; Birch et al. 2019).

In our previous study (Bouquety et al. 2021), we found that the depressions were shaped by sublimation and/or a landslide. In this large-scale study, we conclude that only 8% of the 131 measured depressions on 67P are associated with boulders, that is, are located at a cliff edge with boulders at its feet. 92% of the depressions must be explained by a different process. Altogether, the small fraction of depressions associated with boulders, the analogy with Earth and Mars depressions, and the fact that most (58%) depressions are found in flat terrains with no possible mass-wasting (mean slope of 3.8°) allow us to conclude that subsidence triggered by sublimation is the main process that created and/or modified these depressions. We propose that landslides are secondary events that occur as a response to the destabiliza-tion of the soil generated by the thermokarst phenomenon (i.e., subsidence).

The general scenario for dust erosion via sublimation has been proposed by several authors to explain the transient activity of depressions formed in FPDs (Keller et al. 2015; Groussin et al. 2015b; Shi et al. 2016; Hu et al. 2017; Birch et al. 2019). In particular, Birch et al. (2019) developed a model that explains the scarp migration and activity observed in the Hapi smooth terrains. Based on dust jets coming from the location of depressions (Shi et al. 2016) and on the fact that water sublimation is the main process that drives the depression evolution, together with the assumption that the FPD is a mixture of dust and water ice, the model of Birch et al. (2019) demonstrates that during perihelion approach, the sublimation of water ice leads to a dust jet that anhydrates the FPD, ejecting the non-ice material from the depressions, excepted for the largest boulders (>1 m). This process leads to the propagation of a sublimation front. This explains the depressions growth (i.e., the scarp migration).

Our thermokarst scenario and the model for migrating scarps developed by Birch et al. (2019) agree that the evolution of depressions results from the propagation of a sublimation front, as also proposed in our previous paper on this topic (Bouquety et al. 2021). However, they do not agree on the process that initiated the formation of the depression. For Birch et al. (2019), it is the same sublimation process, with the radiation coming from neighboring boulders or cliffs, which “act as seed points for an instability that propagates across the entire plain”. For our thermokarst scenario, the subsidence resulting from the permafrost thawing triggered by sublimation is the key process. The advantage of the subsidence is that it works on all terrains covered by a FPD, independently of the neighboring morphologies and topography. The subsidence creates the necessary slopes for the sublimation front to propagate and the depression to grow, following the evolution model of Birch et al. (2019), for example.

7.2 Periglacial system in cometary permafrost

We identified 131 cometary thermokarst depressions on 67P. They cover ≈1% of the nucleus surface and are exclusively located in the FPD layer. These morphologies are analogous to the thermokarst depressions on Earth and Mars and reveal permafrost. On Earth and Mars, the thermokarst depressions are often associated with other morphologies related to permafrost, such as thermal contraction polygons, and they form a periglacial landscape.

On 67P, Auger et al. (2018) identified more than 3600 thermal contraction polygons thatare similarto those found on Earth and Mars. The polygons cover1.5% of the cometary surface and are formed by the diurnal temperature variations in a consolidated layer of water ice. In our study, we found that 6% of the cometary thermokarst depressions are associated with polygons (Figs. 15c and d). This association can only be detected in specific areas in which FPDs from the thermokarst depression is in contact with exposed consolidated material, but it might indeed hold for a larger fraction of depressions in which the polygons are hidden below FPDs. On Earth and Mars, depressions and polygons can be found in the same material, and polygons can be observed around and inside depressions. This association of morphologies forms a periglacial landscape on 67P.

By analogy with Earth and Mars, the FPD layer of a few meters on 67P is therefore cometary permafrost, composed of an ice-rich mixture of dust and ice, as confirmed by spectroscopic observations (Groussin et al. 2015b; Shi et al. 2016; El-Maarry et al. 2017; Birch et al. 2019). The observed thermokarst depressions are caused by the thawing and subsidence of this cometary permafrost.

Fig. 13 Box and whisker plot of the thermal parameters for the different regions of 67P and for the different part of each depression: BD for the bottom of the depression, SB for the slope break, and D.o.B for the outer border of the depression. (a) Temperature (K), (b) total erosion (m), and (c) total energy received from the Sun (109 J/m2). The background color is an indication of the latitude: [−90°; −30°] in blue, [−30°; +30°] in green, and [+30°; +90°] in yellow.

Fig. 14 Histogram illustrating in which part of the depressions (bottom, slope break, or outer part), the temperature, erosion, and energy are highest. The slope break is the key location. Temperature, erosion, and energy received from the Sun reach their highest value for most depressions.

Fig. 15 Example of a periglacial system on (a) the ACP on Earth, (b) Utopia planitia on Mars, (c) the Anhur region on 67P, and (d) the Ma’at region on 67P. In each case, the thermokarst depressions are associated with thermal contraction polygons.

7.3 Composition and properties of the permafrost or FPD layer

The subsidence of permafrost provides an important constraint on its composition because it only occurs when the permafrost is supersaturated, that is, when the volume of ice exceeds the volume of pores (Costard & Kargel 1995; French 2007; Morgenstern et al. 2007; Soare et al. 2008; Levy et al. 2010; Ulrich et al. 2010; Séjourné et al. 2011). On Earth, supersaturated permafrost is estimated to have a volume fraction of water ice between 50% and 80% (French 2007). On Mars, the permafrost in Utopia planitia contains a volume fraction of water ice of 50%–80% (Séjourné et al. 2019), measured by the SHARAD (shallow radar) instrument of the Mars Reconnaissance Orbiter mission (MRO). A similar volume fraction of water ice of 50–80% in the FPD layer of 67P is therefore likely a good approximation.

This value can be set in perspective using the work of Haack et al. (2021), who simulated the cometary surface erosion in the laboratory. Their experiments revealed that sublimation morphologies such as the depressions studied here are strongly dependent on the ratio of water ice to dust mass and on the insulation direction. They estimated that a volume fraction of water ice between 25% and 50% is more susceptible to producing sublimation-driven morphologies. Experiments by Poch et al. (2016a,b), and Kaufmann & Hagermann (2018) also show that for a high volume fraction of water ice ≥66%, a network of organic filaments can form. This limits the subsidence. It is difficult to reconcile this laboratory experiment at the centimeter scale with our observed features at the meter scale, but overall, the above discussion indicates a volume fraction of water ice of about 50% in the FPD layer of 67P. When we assume a density of refractory material of 2 (e.g., organics) to 3.5 times (e.g., silicates) higher than that of water ice, this corresponds to a dust-to-water mass ratio in the range 2–3.5 for this external FPD layer. Our value is at the lower end of the very wide range found for the refractory-to-ice mass ratio on 67P, which increases from lower than 1–15, depending on the method used to derive it (Choukroun et al. 2020).

Remarkably, none of the observed thermokarst depressions on 67P exposes the consolidated material underlying the FPD. The depression floor is always cover by FPDs, which seems to be thick enough for the depression to develop into it. The depth of the depression therefore provides a lower limit of the FPD thickness, as is also the case for the scalloped depressions on Mars (Séjourné et al. 2012). Figure 16 displays a map of the FPD minimum thickness measured in the different regions of 67P. The average value is 4.8 m, with some variations across the surface, from 2.5 m to 8 m. From the analysis of an impact crater in an FPD layer in Ash, Nébouy et al. (2022) also found a minimum thickness of 5 m that agrees with these measurements. This is thicker than the value predicted by models, which is from 5 mm to 2 m (Hu et al. 2017; Oklay et al. 2017), but it is similar to the estimates of El-Maarry et al. (2017), which lie between 1 and 5 m. The largest depression depths are found in the southern hemisphere (Bes and Atum), which is consistent with the fact that this part of the nucleus is most strongly eroded according to thermal models (Keller et al. 2015). However, even if the southern hemisphere receives more energy and is more eroded, the outgassing drags dust particles that are either ejected from the nucleus or are redeposited with a ballistic trajectory in the northern plains (Keller et al. 2015). This redistribution of dust particles from south to north can fill the northern depressions, decrease their depth, and even likely erase the smallest depressions. Thus, the measured depth depressions in the northern hemisphere may be underestimated.

8 Formation and evolution of thermokarst depressions induced by sublimation

Despite the morphometrical similarities between the thermokarst depressions observed on 67P, Earth and Mars, a notable difference has to be underlined. While for 67P and Mars, the subsidence is likely to be triggered by sublimation, on Earth, it is induced by the thawing of the ice permafrost that creates liquid water. This leads to a thermokarst lake. Nevertheless, our comparative morphometrical analysis demonstrated that they can be considered geologically analogous, even if the process at the origin of the subsidence is different (sublimation versus melting).

We now propose the following scenario for the formation and evolution of cometary thermokarst depressions (Fig. 17):

(a) Initial state. The FPD layer is composed of a mixture of dust and water ice. The amount of water ice is not homogeneous in this layer. This initial state is characterized by a relatively smooth surface, without observable topographic hollows;

(b) Subsidence induced by sublimation. Approaching perihelion, the solar energy received by the surface increases, leading to an increase in water ice sublimation, which in turn decreases the concentration of water ice in the FPD layer. This increases the porosity and destabilizes the surface, to the point where it triggers a subsidence. Because cliff walls are thought to be a source of cometary activity (Vincent et al. 2016), water ice can concentrate at the cliff tops (e.g., due to sublimation or recondensation cycles in the subwall layer; Sanctis et al. 2015). This in turn might favor the apparition of subsidence above the cliff walls. This effect is more pronounced when the water ice content is higher, and it leads to the apparition of topographic hollows at the surface. These hollows could also correspond to underlying topographic heterogeneities. During our processing of the NAC images, we observed several of these topographic hollows, which we were unable to classify as depressions because they did not meet all our criteria (Sect. 3.1.2). However, because most of these hollows are located close to existing depressions, they might correspond to the first stage of a thermokarst depression;

(c) Asymmetrical growth of the depression. The orientation of the edges of the topographic hollows is different from that of its bottom. When the comet approaches the Sun, the side of the hollow facing the Sun sublimates more strongly and becomes the preferential direction of erosion, with a retreat of the sublimation front. The hollow becomes a growing depression. The retreat of the sublimation front is an erosion process by which exposed water ice sublimates and drags dust with it, resulting in a retreat of the scarp. The subsidence might also play a key role here, in addition to dust dragged by gas, and it might accelerate the overall depression growth process;

(d) Depression evolution. The depression continues to grow, mainly on the side that faces the dominant insulation. Close to a cliff, a collapse can occur and create new deposits at the cliff feet. In the long term, adjacent depressions can merge to form a larger depression, and they can also become inactive if the FPD layer is desiccated or if a change in insulation conditions is caused by a modification of the pole orientation or of the comet orbit.

Fig. 16 Global map of the minimum thickness of the FPD layer based on the depressions depth.

9 Conclusions

We studied the morphology and thermal environment of 131 depressions on the nucleus of comet 67P in detail. We compared these depressions with 50 thermokarst depressions (alases) on Earth, and 200 scalloped terrains on Mars. From this comparison, we came to the conclusion that these three types of depressions are geological analogs. In particular, they have similar morphological properties in terms of area versus perimeter evolution, depth, and asymmetric slope. They all show a slope break on one side, which is the direction of preferential growth (i.e., erosion). They are all shaped by a thermokarstic process (i.e., permafrost thawing) triggered by either sublimation (for 67P and Mars) or melting (Earth). We propose to call these depressions on 67P cometary thermokarst depressions.

We compared our results with previous hypotheses about the formation of these depressions. The low number of depressions associated with boulders (8%) seems to exclude landslide as the main driver of their formation. Landslides are more likely secondary events that occur as a response to the desta-bilization of the soil generated by the thermokarst phenomenon (i.e., the subsidence). Our study agrees with the dust-jet erosion model developed in the literature (e.g., Birch et al. 2019), which explains the growth of depressions by propagation of a sublimation front. We strongly favor the creation of depressions by a subsidence process, however, rather than by ejection of dust particles from outgassing.

On 67P, the thermokarst depressions are exclusively located in FPDs terrains. By analogy with the Earth and Mars, this FPD layer of a few meters thickness can be considered cometary permafrost and should contain a volume fraction of water ice of about 50% for subsidence to occur. The depressions on 67P are distributed throughout the nucleus, with no preferential distribution for a particular hemisphere or region. They are either single and isolated (71%) or are grouped in clusters (29%). 58% of the depressions are located in flat terrains (mean slope 3.8°), and the others are located at cliff edges. Their formation and evolution are driven by the Sun; the side facing the dominant insulation is the preferential direction of erosion. In the long term, adjacent cometary thermokarst depressions can merge to form a larger depression, and they can also become inactive if there is a change in insulation conditions caused by a modification in pole orientation or in the comet orbit. Figure 17 illustrates the scenario for the formation and evolution of cometary thermokarst depressions.

Our results contribute to a better understanding of the periglacial system on comet 67P, with morphologies associated with permafrost and driven by sublimation. The analogy with Earth and Mars is of upmost interest because it demonstrates that similar morphologies can exist in very different gravitational environments. In a near future, we could aim to predict the evolution of these structures, and others, to have a better view of the erosion of cometary nuclei.

Fig. 17 Proposed scenario for the formation and evolution of cometary thermokarst depressions (see Sect. 7.3 for details).

Acknowledgements

OSIRIS was built by a consortium of the Max-Planck-Institut fur Sonnensystemforschung, Katlenburg-Lindau, Germany; CISAS University of Padova, Italy; the Laboratoire d’Astrophysique de Marseille, France; the Instituto de Astrofisica de Andalucia, CSIC, Granada, Spain; the Research and Scientific Support Department of the ESA, Noordwijk, Netherlands; the Instituto Nacional de Tecnica Aeroespacial, Madrid, Spain; the Universidad Politechnica de Madrid, Spain; the Department of Physics and Astronomy of Uppsala University, Sweden; and the Institut fur Datentechnik und Kommunikationsnetze der Technischen Universitat Braunschweig, Germany. The support of the national funding agencies of Germany (DLR), France (CNES), Italy (ASI), Spain (MEC), Sweden (SNSB), and the ESA Technical Directorate is gratefully acknowledged. We thank the Rosetta Science Operations Centre and the Rosetta Mission Operations Centre for the successful rendezvous with comet 67P/Churyumov-Gerasimenko. A. Bouquety, O. Groussin, and L. Jorda, acknowledge the support from CNES.

Appendix A List of NAC images

References

All Tables

All Figures

	Fig. 1 NAC images of depressions in the Ma’at (a) and in the Anuket region (b).
In the text

Fig. 2 Examples of planetary depressions in different contexts. (a) Okmok caldera (Alaska). (b) Martian impact crater in Lunae planum (Colour and Stereo Surface Imaging System, CaSSIS). (c) Landslide in Quebec. (d) Kettle holes in Wrangell-Saint-Élie national parc (Alaska). (e) Sinkhole in Russia. (f) Thermokarst lake in the Siberian permafrost. (g) Collapsed pingo (Canada). (h) Retrogressive thaw slump in the Siberian permafrost.

In the text

Fig. 3 Examples of measured thermokarst depressions. (a) Alases in Alaska on Earth (digital orthophoto quadrangle, DOQ, 5mpixel−1). (b) Scalloped terrain in Utopia planitia on Mars (HiRISE image, 25 cm pixel−1). (c and e) Associated digital elevation nodels. (d and f) Slope maps. (g and h) Topographical profiles along the length direction (in blue in panels d and f). The pink lines show the width measurement.

In the text

Fig. 4 Dataset derived from the SPC shape model that we used to measure the geometric parameters of depressions on 67P. (a) NAC image of depressions located in the Ma’at region. (b) Distance between the camera and the surface and (c) emission angle, both draped on an NAC image. (d) Gravitational heights (GH) draped on an NAC image that we used to extract the topographical profile in panel f. (e) Gravitational slopes (GS) draped on an NAC image. (f) Topographic profile along the blue line in panel d, extracted from GH.

In the text

	Fig. 5 Results from the thermal analysis of a depression in the Ma’at region. (a) Maximum temperature (K), (b) total erosion (m), (c) total energy received from the Sun (Jm−2), (d) depression subdivision into three units, and (e) topographical profile showing the subdivisions.
In the text

Fig. 6 Box and whisker plot for all parameters on the three bodies. The black lines at each end correspond to the maximum and minimum value. The box represents the interquartile range; Q1 and Q3 are the bottom and top of this box. The median value is represented by the black lines inside the box. Abbreviations are shown in Table 2. Note the logarithmic scale in panel a and the linear scale in panel b.

In the text

	Fig. 7 Density map of depressions in the different regions of the nucleus of 67P.
In the text

Fig. 8 NAC images showing different examples of depression context. The white arrows show the depressions. (a) Depression in a flat terrain in the Anuket region. (b) Depression at a cliff edge in the Aswan region. (c) Isolated depression in the Ash region. (d) Depression field in the Imhotep region. (e) Depressions associated with thermal contraction crack polygons in the Anhur region. (f) Depressions associated with boulders in the Imhotep region.

In the text

	Fig. 9 Area vs. perimeter for the measured depressions on 67P. The differently colored squares correspond to the different depression contexts (see text for details).
In the text

	Fig. 10 Area vs. perimeter for the thermokarst lake measured in the ACP (circles), the scalloped depressions (triangles), Martian craters (dia-monds),and depressions on 67P(square). The blue colorscale represents the circularityof the measured features.
In the text

	Fig. 11 Topographic profiles showing the asymmetric slopes of the studied depressions. (a) Circular depression in the ACP on Earth. (b) Scalloped depression in Utopia planitia on Mars. (c) Depression in the Seth region of 67P. The dashed line indicates the location of the topographic profile, and the green arrows are markers.
In the text

	Fig. 12 Depth value of the measured depressions and for other terrestrial thermkarstic depressions in Siberia (Ulrich et al. 2010), the Lena Delta (Morgenstern et al. 2011), and in the Qinghai-Tibet engineering corridor (Niu et al. 2014).
In the text

Fig. 13 Box and whisker plot of the thermal parameters for the different regions of 67P and for the different part of each depression: BD for the bottom of the depression, SB for the slope break, and D.o.B for the outer border of the depression. (a) Temperature (K), (b) total erosion (m), and (c) total energy received from the Sun (109 J/m2). The background color is an indication of the latitude: [−90°; −30°] in blue, [−30°; +30°] in green, and [+30°; +90°] in yellow.

In the text

	Fig. 14 Histogram illustrating in which part of the depressions (bottom, slope break, or outer part), the temperature, erosion, and energy are highest. The slope break is the key location. Temperature, erosion, and energy received from the Sun reach their highest value for most depressions.
In the text

	Fig. 15 Example of a periglacial system on (a) the ACP on Earth, (b) Utopia planitia on Mars, (c) the Anhur region on 67P, and (d) the Ma’at region on 67P. In each case, the thermokarst depressions are associated with thermal contraction polygons.
In the text

	Fig. 16 Global map of the minimum thickness of the FPD layer based on the depressions depth.
In the text

	Fig. 17 Proposed scenario for the formation and evolution of cometary thermokarst depressions (see Sect. 7.3 for details).
In the text