# Topographically Influenced Evolution of Large-scale Changes in Comet 67P/Churyumov-Gerasimenko's Imhotep Region

Abhinav S. Jindal, Samuel P. D. Birch, Alexander G. Hayes, Orkan M.

Umurhan, Raphael Marschall, Jason M. Soderblom, Jean-Baptiste Vincent, Dennis Bodewits

► To cite this version:

Abhinav S. Jindal, Samuel P. D. Birch, Alexander G. Hayes, Orkan M. Umurhan, Raphael Marschall, et al.

HAL Id: insu-03869674

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

Submitted on 24 Nov 2022

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

# Topographically Influenced Evolution of Large-scale Changes in Comet 67P/Churyumov–Gerasimenko’s Imhotep Region

Abhinav S. Jindal<fnref n="1" /> ID, Samuel P. D. Birch<fnref n="2" /> ID, Alexander G. Hayes<fnref n="1" /> ID, Orkan M. Umurhan<fnref n="1" />,³,<fnref n="4" /> ID, Raphael Marschall<fnref n="5" /> ID, Jason M. Soderblom<fnref n="2" /> ID, Jean-Baptiste Vincent<fnref n="6" /> ID, and Dennis Bodewits<fnref n="7" /> ID

Received 2022 April 22; revised 2022 June 1; accepted 2022 July 2; published 2022 August 16

## Abstract

Large portions of comet 67P/Churyumov–Gerasimenko’s northern hemisphere are blanketed by fallback material consisting of centimeter-sized particles termed the smooth terrains. Observations from the Rosetta mission show that the most drastic transient changes during 67P’s 2015 perihelion passage occurred within a subset of these deposits. However, we still do not understand the processes driving these changes, limiting our overall understanding of how comets evolve over both seasonal and multiorbit timescales. Herein we provide a complete documentation of scarp-driven activity on 67P’s largest smooth terrain deposit, a highly active portion of the Imhotep region that is the southernmost of all smooth terrain basins on 67P. We also present a thermal model that accurately predicts when and where scarpes originate during the course of the observed activity. Assuming a uniform surface composition, our model shows that activity is heavily controlled by local topography rather than the presence of ice-enhanced hot spots on the surface. Scarps within the smooth terrain deposits in central Imhotep also exhibit a peculiar behavior, where three scarpes originate from the same location but at different times and migrate in opposite directions. This behavior indicates that the landscape retains a memory of previous cycles of erosion and deposition, reflected by the depth of the volatile-rich layer. Future work will need to couple our thermal model with a landscape evolution model in order to explain the complete dynamic evolution of these terrains.

Unified Astronomy Thesaurus concepts: Comet nuclei (2160); Small Solar System bodies (1469); Comets (280)

## 1. Introduction

The Rosetta mission to comet 67P/Churyumov–Gerasimenko (hereafter 67P) has opened the doors to a new era of small-body geology by allowing us to directly interrogate the processes that drive cometary surface evolution. With the high-resolution imaging and topographic data collected by Rosetta, we finally have the necessary inputs to apply terrestrial geomorphic principles so as to understand the coupled climate–landscape evolution of these small bodies.

Rosetta data revealed that the northern hemisphere of 67P is dominated by vast, topographically smooth deposits of centimeter-to decimeter-sized particles, termed the smooth terrains (Birch et al. 2017). The obliquity and eccentricity of 67P’s orbit combine such that the net solar insolation during southern summer, which approximately corresponds with perihelion, is orders of magnitude greater than the net insolation during northern summer (Hu et al. 2017). As a result, a greater number of particles are liberated from the nucleus during southern summer (Thomas et al. 2015). Liberated particles follow ballistic trajectories, where the smallest grains escape to the coma, but a large fraction of centimeter-sized or larger particles fall back on to the surface, eventually accumulating into the cold, gravitational lows in the northern (winter) hemisphere, forming the smooth terrains (Lai et al. 2017).

Between the ejection process and the long flight time across the nucleus for the fallback material (Davidsson et al. 2021), many of the most volatile ices (like CO and CO₂) are likely lost. The exposure of these particles, however, is such that water ice likely remains embedded within them (De Sanctis et al. 2015; Fornasier et al. 2016; Davidsson et al. 2021). This is supported by Rosetta’s mass spectrometer measurements, which show an increase in water production in the northern hemisphere and far higher CO and CO₂ outgassing rates in the south (Gasc et al. 2017). Together, this implies that the smooth terrains comprise water ice–rich particles deposited during each perihelion passage that blanket the surface with a layer of still-unknown thickness.

The survival of ice within the individual grains that build up the smooth terrains means that on 67P’s subsequent orbit, erosion of this ice-rich sediment layer can occur on approach to the Sun. This retained ice is the driver for the many rapid changes observed within smooth terrains (Groussin et al. 2015; El-Maarry et al. 2017; Hu et al. 2017; Pajola et al. 2017a; Birch et al. 2019). While changes in consolidated terrains appear to occur on multiorbit timescales (as indicated by the small number of changes observed in consolidated terrains during the Rosetta mission and previous works such as Steckloff & Samarasinha 2018), the smooth terrains are evolving on much shorter timescales. As such, it is essential to understand the evolution of these terrains in order to decipher the overall surface evolution of the comet. Moreover, the smooth terrains are not unique to 67P and are also prevalent on other observed Jupiter family comets, like 19P/Borrelly (Soderblom et al. 2002; Britt et al. 2004), 9P/Tempel 1 (Veverka et al. 2013;

<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. 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">¹Cornell Center for Astrophysics and Planetary Sciences, Cornell University, Ithaca, NY 14853, USA; asj59@cornell.edu</footnote>

<footnote n="2">²Department of Earth, Atmospheric and Planetary Science, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
³SETI Institute, 189 Bernardo Way, Mountain View, CA 94043, USA</footnote>

<footnote n="4">⁴NASA Ames Research Center, Mail Stop 245-3, Moffett Field, CA 94035, USA</footnote>

<footnote n="5">⁵CNRS, Observatoire de la Côte d'Azur, Laboratoire J.-L. Lagrange, Nice, France</footnote>

<footnote n="6">⁶Deutsches Zentrum fur Luft-und Raumfahrt (DLR), Institut für Planetenforschung, Berlin, Germany</footnote>

<footnote n="7">⁷Physics Department, Leach Science Center, Auburn University, Auburn, AL 36849, USA</footnote>
Thomas et al. 2013a), and 103P/Hartley 2 (Thomas et al. 2013b), suggesting a similar origin. Understanding the processes capable of modifying smooth terrains may also provide critical context for future comet sample return missions that target them (e.g., the Comet Astrobiology Exploration SAmple Return; Squyres et al. 2018).

Unlike previous data sets, data from Rosetta's Optical, Spectroscopic, and Infrared Remote Imaging System (OSIRIS; Keller et al. 2007) allow access to the spatial scales where the processes responsible for eroding the smooth terrains can be directly studied. The temporal coverage of the Rosetta mission (2+ yr) also permits us to study decimeter-scale changes and modifications of the smooth terrains. To date, however, only a few studies have used these data to document and/or model surface changes on 67P (Groussin et al. 2015; Pajola et al. 2017a; El-Maarry et al. 2017; Hu et al. 2017; Birch et al. 2019), with each study documenting activity only over a portion of the Rosetta mission and on terrains largely in 67P's northern hemisphere and not representative of 67P more broadly. As such, though we know that most changes likely occur within 67P's smooth terrains, we do not understand how the observed changes relate to each other, how representative any set of changes are, or how they contribute to the overall erosional and depositional history of 67P.

To resolve this, we choose the Imhotep region to document the full complexity of the evolution of smooth terrains and couple it with a thermal model to explain the observations, building upon previous studies (Groussin et al. 2015; El-Maarry et al. 2017) that documented a subset of activity. The Imhotep region is a vast sedimentary basin on 67P's large lobe (Figure 1), home to the largest deposit of smooth terrains on the comet (Auger et al. 2015). Unlike all other smooth terrains on 67P, Imhotep's equatorial location ensures that it never experiences polar winter and will therefore be active throughout perihelion. This results in greater activity relative to other areas on the comet, including multiple still-unexplained outbursts (Knollenberg et al. 2016; Agarwal et al. 2017), though such outbursts were not observed within the central smooth terrain deposits that our work focuses on. Despite this high activity, Imhotep still serves as a sink for sediment transported from across the comet (Davidsson et al. 2021), suggesting that the influx of material still outpaces activity from the region over multiyear timescales. This complex behavior of erosional and depositional processes acting simultaneously, therefore, makes Imhotep an ideal location to study the evolution of smooth terrains on comet surfaces. Moreover, we observe this region over the entirety of the Rosetta mission, documenting exactly when large-scale changes occur and when the region remains mostly inactive.

## 2. Migrating Scarps in the Imhotep Basin

We utilize image data from the OSIRIS narrow-angle camera (NAC) on board the Rosetta spacecraft to document and measure all large-scale changes within the Imhotep region over the course of the Rosetta mission. We found changes in the landscape from 2015 May to December, with no observed activity from either 2014 August through 2015 April or 2016 January until the end of the Rosetta mission (2016 August). While hundreds of images for Imhotep were acquired throughout the Rosetta mission, we rejected the majority of them due to poor image resolution, viewing geometries hindering the observation of the scarps, or no changes being observed to take place. This left us with a total of

24 images for our analysis (see Table 1 in Appendix A), 22 of which were acquired between 2015 May and 2016 January. These images were projected in a common reference frame using the shapeViewer software (Vincent 2018) and subsequently imported as layers into ArcGIS to allow more precise georeferencing and measurements of migration rates and directions. The migration rates of individual scarps, of course, are an average across the time interval between the two images of interest (see Table 2 in Appendix B). We note that activity may occur over far shorter time periods (e.g., minutes to hours) with long intervals in between of limited/no activity. The rates we document are therefore the minimum rates of migration, with the relative rates between scarps of most interest. For the average migration direction, because the scarps are observed to expand radially, we trace rays from the origin point of the scarp to points along the scarp boundary. Averaging these angles provides an average direction of a given scarp.

### 2.1. Observations of Imhotep during 67P's Perihelion Passage

The majority of Imhotep's smooth terrain deposits (hereafter the "central" deposit, termed Imhotep "a" in Thomas et al. 2018) are concentrated between 7°N and 41°S, within a gravitational low and surrounded by steep (up to ~700 m of relief) topography (Figure 1). There are additional isolated deposits of smooth terrains within the surrounding alcoves and other smaller deposits scattered across the broader region. The central deposit retains the majority of the region's smooth terrains, covering an area of 0.8 km² and gently sloping toward the east (Auger et al. 2015). There is little obvious texture to the smooth terrains in this region (Figure 1), suggesting that it is made up of material far below Rosetta's resolution and likely made of centimeter-sized granular particles, like all other smooth terrains observed up close (Mottola et al. 2015; Pajola et al. 2017b).

Rosetta's first observations of Imhotep revealed two prominent scarps within the central deposit (scarps "h" and "n" in Figure 2(a)) and one additional scarp in one of the alcoves on the western end of Imhotep (scarp "m" in Figure 2(a)). These scarps existed in the first Rosetta images of 67P obtained in late 2014 and did not change until sometime after 2015 May 24 (Figure 2(a)). Groussin et al. (2015) documented the first changes observed to occur in this region, the growth of a few large scarps between 2015 May 24 and July 11. Our work builds upon theirs, documenting numerous previously undocumented changes that continue throughout perihelion, until 2015 December 6. Throughout this period where we observe changes, the Sun always rises over the eastern end of Imhotep, with the first light falling upon the large cliffs (~400 m high) on the western margin of the basin (Figure 1).

The first phase of changes occurred between 2015 June 3 and June 27. We identify six unique scarps, "a," "b," "d," "e," "c," and "f" (Figures 2(b) and (c)), in the central deposit, each of which began at topographic discontinuities like boulders ("e," "f," and "c"), alcove edges ("b"), and preexisting scarp boundaries ("a" and "d"). Scarps "a," "b," and "c" are the first to appear on 2015 June 5, followed shortly after by scarp "d" on June 13. Scarp "d" appears from the same location as scarp "a" but migrates in the opposite direction (Figures 2(c) and 3(a)). Scarps "e" and "f" are then observed on June 27, though we are uncertain of their source or initiation date due to these
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Overview of the Imhotep region. (a) The Imhotep region (highlighted in yellow) lies on the large lobe of 67P and has the equator (dashed green line) passing through its northern portion. (b) Imhotep consists of a central deposit of smooth terrains bounded by high-standing topography on all sides. (c) Shown is the topography surrounding Imhotep as viewed from the north of Imhotep. Source: ESA—European Space Agency (https://creativecommons.org/licenses/by-sa/3.0/igo), (d) Shown is a high-resolution view of two large boulders within the smooth terrain deposits. The granular nature of the deposit is evident by the speckle pattern in the smooth terrains at such resolutions (0.2 m pixel⁻¹).</figcaption></figure>

regions being within shadows in the only intervening image acquired on June 13. Given their significant size by the time we first observe them, however, we favor a formation closer to June 5. The eastward-moving scarps ("a," "b," and "e") move at moderate time-averaged rates of ~23 cm hr⁻¹, while the southward-moving scarp "c" and westward-moving scarp "d" move significantly slower, at ~11 and ~14 cm hr⁻¹, respectively (Figure 4(b)). Details as to how the rates were calculated are provided in Appendix B.

By June 27, the initial six scarps expanded across most of the northern half of the central Imhotep basin, consistent with the findings of Groussin et al. (2015). However, Groussin et al. (2015) stopped their analysis after this date, missing a wealth of additional changes, specifically, two additional scarps on this date, with an additional circular pit emerging on July 2. Scarp "g" begins near two large boulders and migrates at ~25 cm hr⁻¹, marking the most southern migration yet observed. Scarp "h," which serves as the initial seed point for both scarps "a" and "d" (Section 3), begins migrating westward, moving at the slowest average rate observed, ~7 cm hr⁻¹ (Figure 4(b)). Finally, one small pit ("i" in Figure 2(d)) expands radially at an

average rate of ~22 cm hr⁻¹ within this region, eventually intersecting with both scarps "f" and "b." By July 11, the scarps and pit all appear to have eroded an entire overlaying layer of smooth terrains, revealing deeper and visually identical smooth terrains, but remain restricted northward of 21°S latitude.

As 67P approached its perihelion on August 13, activity systematically progressed southward. From July 11 to 31, scarp "h" continues expanding westward at a slow average speed of ~7 cm hr⁻¹. Similarly, scarp "g" expands eastward with an average speed of ~26 cm hr⁻¹. Two new scarps also appear on July 11. Scarp "j" originates near the easternmost edge of scarp "h," migrating toward the west at an average rate of ~18 cm hr⁻¹. Scarp "k" emerges from what appears to be a radially expanding pit, near the same region as scarp "c," migrating southward at a much slower average rate of ~4 cm hr⁻¹. The material that this scarp migrates through also appears to be freshly deposited, as that material just recently buried the region that scarp "c" swept through 2 weeks prior on June 27. Scarp "l" also appears from the edge of an alcove on July 15, at nearly the same exact location as scarp "b," and just like scarp "b," it migrates eastward until July 29 at an average speed of
~21 cm hr⁻¹. On July 26, there is yet more activity further south, as scarps “m” and “n” mobilize, exposing underlying bright material. Scarp “m” migrates toward the southwest at an average rate of ~25 cm hr⁻¹, while scarp “n” migrates eastward at an average rate of ~37 cm hr⁻¹. At the same time, scarp “o” forms the base of another alcove (Figure 2(f)) and radially expands eastward at ~45 cm hr⁻¹ until it intersects with scarp “g” on August 6.

We also observe continuous deposition of material in the northern portion of the central deposit during this same period of scarp migration. Scarps “a,” “b,” “c,” “d,” “e,” and “f” all stop migrating around 2015 July 11, after which they appear more muted, with less prominent shadows along their boundaries, suggesting some degree of burial. As image resolution and illumination conditions did not change significantly, this suggests that the more diffuse appearance of features in the northern half of the basin may result from mass sediment deposition. Deposition of material also continues for multiple weeks, with scarp “a” being the first one to get completely buried by July 11. This is followed shortly after by the complete burial of scarps “b,” “c,” and “d” by July 26, pit “i” by July 31, and scarp “e” by August 6 (1 week before perihelion). The region containing scarp “f” is in shadow from July 26 to November 29, making it impossible to determine its fate. Finally, while scarp “h” is not completely buried under sediment, it ceases to migrate by July 31 and appears significantly more muted by September 4. All of these scarps are located north of 14°S, implying a deposition event taking place in the northern portion of Imhotep, while scarp-driven erosion simultaneously persists further south.

After perihelion, we are limited in the quality and number of images available to observe changes on the surface. However, the images available throughout this time still show numerous large-scale changes. The most prominent observation is of a mass deposition event across the entire basin in the weeks following perihelion (August 23–September 4). Specifically, some features that were actively migrating in both the southern and northern extents of the basin (e.g., scarps “l” and “o”) have largely disappeared. Scarp “o” gets completely buried by August 23, followed by scarp “l” on September 4. Scarps “g” and “j” also appear more muted in their appearance during this time, but they continue to migrate, overlapping with one another on August 5, until scarp “j” finally overtakes scarp “g” by August 23. Both of these scarps continue their migration until September 4 (scarp “g”) and October 26 (scarp “j”). Scarps “m” and “n” also continue their migration postperihelion through the most southern portion of the basin, finally coming to a standstill by September 4 and August 23, respectively. A new scarp, “p,” also appears in the same alcove as scarp “e” on August 23; however, unlike “e,” “p” moves in the opposite, westward direction at ~6 cm hr⁻¹ (Figure 2(g), 3). Finally, scarp “k” continues its expansion throughout perihelion, progressively moving southward, strikingly similar to scarp “c.”

By October 26, scarps “m” and “n” are no longer visible, suggesting that deposition of new sediment now outpaces erosion across nearly the entire basin. Scarp “p” still undergoes marginal migration, ceasing any noticeable activity sometime between October 26 and November 29. Finally, scarp “k” stops moving on 2015 December 6, marking the end of all obvious, large-scale activity within the Imhotep basin.

# 3. Thermal Modeling and Scarp Initiation

## 3.1. Model Setup

The changes we observe are driven by thermal processes. To quantify this, we leverage a thermal model developed by Umurhan et al. (2022) but applied to 67P. We calculate the net energy input (direct + reradiated) for each facet of the SHAP4S shape model (Preusker et al. 2015), downsampled to 50,000 facets. Though there is a more updated shape model for 67P (SHAP7; Preusker et al. 2017), the additional images inputted into the newer model largely cover latitudes beyond our region of interest and therefore do not affect our calculations. The 50,000 facet version of the SHAP4S model is further cropped to only include the larger lobe of the comet, such that we only account for facets of the shape model that may reradiate energy onto Imhotep. This leaves us with an approximately 16,000 facet shape model, making our calculations far less computationally expensive.

The amount of energy directly incident on each facet, *i*, depends upon the heliocentric distance of the comet and the subsolar latitude. The subsolar latitude for the duration of the Rosetta mission is calculated from SPICE, and the heliocentric distance is obtained from JPL Horizons. Given the uncertain nature of cometary orbits through time, for simplicity, the subsolar latitude values for multiple orbits are obtained by duplicating the SPICE values and using a cubic interpolation to fill in the gaps. This is reasonable because 67P’s spin axis was rather stable throughout the Rosetta mission (Kramer & Läuter 2019).

We assume a two-layer structure for the surface, where a semi-infinite volatile-rich layer is blanketed by a more ice-depleted top layer. This assumption stems from the top layer losing some of its volatiles over the course of 67P’s rotation, with the thickness of this layer being the diurnal skin depth. Additional layers could be added in future runs, at the cost of computational time. Finally, the surface is divided into “cliffs” and “plains” based on their gravitational slopes (Marschall et al. 2017), where the cliffs have a thermal inertia factor 1.5x greater (Groussin et al. 2019) than the plains (Figure 5). Both the smooth and consolidated terrains are modeled assuming the same albedo and composition. Other than the thermal inertia, the surface is assumed to be uniform across the entire region, contrary to previous works (e.g., Agarwal et al. 2017).

For our purposes, we adopt the time-asymptotic fast Fourier transform thermal modeling framework found in Umurhan et al. (2022). We use a two-step approach, first to assess the daily averaged insolation received by every surface element, followed by an estimate of the total sublimation on any given day. Calculating the total energy input to each facet on the shape model involves calculating both the directly incident and reradiated energy. To calculate the reradiated energy on every facet, we first find all facets that are visible to a given facet (Figure 5(a)) and calculate the angles between the surface normals of these facets to derive an apparent luminosity. The details for this step are found in Umurhan et al. (2022). Next, we obtain the thermal solution by splitting the thermodynamic response of every facet within the region of interest into a day phase and night phase and applying appropriate boundary conditions at the interface of the ice-rich and ice-depleted layers and at the surface. The boundary conditions at the interface of the two layers include ensuring continuity of the temperature and heat flux across the interface for both the day and night phases. We assume that sublimation occurs at the interface
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Georeferenced images of the Imhotep region. (a) Preexisting scarps “h,” “m,” and “n” can be observed on 2015 May 24 before activity begins in the Imhotep region. (b) By 2015 June 5, new scarps “a,” “b,” and “c” emerge in the northern portions of the basin. (c) New scarps “d,” “e,” “f,” and “g” are observed on 2015 June 27. (d) New scarp “i” is observed on 2015 July 2. (e) New scarps “j” and “k” emerge on 2015 July 11. At the same time, the boundaries of the scarps in the northern portions of the basin appear more muted, indicating sediment deposition. (f) By August 1, scarp activity reaches the southern end of the basin, while deposition continues toward the north. New scarps “l” and “o” emerge during this time, while preexisting scarps “m” and “n” also begin their migration. (g) New scarp “p” is first observed on 2015 August 23. Deposition of material spreads across the basin during this time, with several scarps completely disappearing under newly deposited sediment. (h) After 2015 October 27, scarps “g” and “j” no longer migrate, whereas scarps “p” and “k” continue migrating until 2015 November 29 and December 6, respectively. (i) Scarps “g,” “j,” “k,” and “p” are still visible on 2016 January 23 and remain visible for the entirety of the Rosetta mission. These scarps likely act as seed points for the next cycle of erosion as the comet approaches perihelion again.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Average direction of migration for scarps within the Imhotep basin. (a) Scarps initiating between 2015 May 24 and June 13. (b) Scarps initiating between 2015 June 13 and July 1. (c) Scarps initiating between 2015 July 1 and August 13.</figcaption></figure>

of the two layers only during the daytime; hence, balancing the heat flux during the daytime includes an irreversible consumption of energy to drive sublimation. Sublimation during nighttime does not need to be considered, since the extremely low thermal inertia of cometary surfaces means they do not store much heat. We provide full details on the implementation of this step in Appendix C.

Broadly, the depth-dependent thermal energetics are determined over the course of one full orbit in which, for every surface element i, we solve the standard one-dimensional heat equation,

$$ \rho C_p \frac{\partial T_i}{\partial t} = \frac{\partial}{\partial z} K_i \frac{\partial T_i}{\partial z}, \qquad (1) $$

for the temperature $T_i(z, t)$ subject to the boundary condition that the thermal flux goes to zero at sufficient depth,

$$ K_i \frac{\partial T_i}{\partial z} \bigg|_{z \to -\infty} = 0. \qquad (2) $$

We assess the balance of the daily averaged incoming solar radiation received at the surface, $(1-A)f_{\odot,i}(t)$ plus the amount of daily averaged radiation locally received from the surrounding terrain visible to surface element i. This latter radiation is denoted here as $F_{ij}(T_j, \dots)$, where j is the subset of all facets that are visible to facet i, as alluded to above (Figure 5).

The total received radiation is balanced by losses to space via blackbody radiation plus thermal conduction into the subsurface and energy consumed in sublimation, $F_s$. Thus, the boundary condition at the surface ($z=0$) is given by

$$ (1 - A)f_{\odot,i}(t) + F_{ij}(T_j, \dots) = \varepsilon \sigma T_i^4 \bigg|_{z=0} + K_i \frac{\partial T_i}{\partial z} \bigg|_{z=0} + F_s. \qquad (3) $$

In developing the orbital timescale solution, $F_s$ is set to zero. The $F_s$, from which a mass-loss rate is calculated, is then treated as a first-order correction and calculated on the basis of a diurnal submodel. The model then assesses the total sublimation based on the effective thermal energy available to the very top layers of the surface (i.e., the total incoming radiation minus the losses conducted into the subsurface)—an approach that is sufficiently satisfactory, since the thermal skin depth is much smaller than the corresponding seasonal skin depth (see also Table 3 in Appendix C).

For simplicity, we only assume H₂O sublimation for our model, which is a reasonable assumption, since most of the material in the Imhotep region is fallback material depleted in CO₂ and other highly volatile ices (De Sanctis et al. 2015; Fornasier et al. 2016). We also assume a uniform distribution of material on the surface; i.e., all regions have an equal amount of ices available for sublimation, and any difference in the sublimation rates is a direct consequence of excess reradiated energy from the surrounding topography.

## 3.2. Correlation of Model Predictions with Observations

To compare the predictions from our model with the observed activity, we divide our study of the activity in Imhotep into five phases throughout 2015: (1) May 24–June 5, (2) June 5–June 27, (3) June 27–July 11, (4) July 11–August 6, and (5) August 6–August 23. While activity continues beyond August 23, no new scarps are observed; hence, we only compare predictions from
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Duration of migration and average migration rates for scarps within the Imhotep basin. Colors are the same as in Figure 3. (a) The solid bars represent the duration for which scarps migrate, with the dashed line indicating the uncertainty in this duration (due to the cadence of the images). The two stars connected by a dotted line indicate the period within which a given scarp is buried under sediment. The red vertical dashed line indicates perihelion passage for 67P. Scarps are arranged from top to bottom in decreasing order of latitude. (b) The range of average migration rates for each scarp is shown, with outliers denoted by single points. Scarps are arranged from left to right in order of when they start migrating.</figcaption></figure>

our model up to this point. For each phase, we calculate an integrated mass-loss rate (i.e., the sum of the daily mass-loss rates) starting 4 days prior to the start date of the phase up to the end date of the phase. The additional 4 days are included to account for the thermal history of the landscape, which will have an influence on the initiation of scarp migration. We then threshold the integrated mass-loss rate in each phase at the 80th percentile to find the facets that are driving the majority of activity in each phase (Figure 6) and are therefore most likely to exceed our unknown “critical volume” (CV; Section 4.3) necessary to modify the surface. We note, however, that this calculation of the mass-loss rate is based on a static landscape, using the shape model’s topography as an input. Hence, regions of high mass-loss rates in our model predict the seed points for scarp activity and not necessarily their subsequent evolution, as the topography clearly changes as the landscape evolves. As topography heavily influences where activity occurs, the destruction/creation of new topography by activity itself is a feedback our model cannot yet account for.

Our model shows that the most equatorial parts of Imhotep are the first to warm up, matching our observations (Figure 6). As the subsolar point shifts southward, the locations of the peak temperatures and mass-loss rates follow. Observed activity during phase 1 starts with the emergence of scarp “a” on June 3, followed by scarps “b” and “c” on June 5. Each of these scarps originates from a region for which our model predicts high mass-loss rates (Figure 6) and keeps expanding into phase 2. Phase 2 involves the commencement of migration of scarp “h” and the emergence of four new scarps: “d,” “e,” “g,” and “f.” While scarps “h,” “d,” and “g” emerge from regions of high mass-loss rate in phase 2, scarps “e” and “f” emerge from regions where the mass-loss rate lies below the 80th percentile for this phase (Figure 6). However, scarp “e” is formed sometime between May 24 and June 27, while scarp “f” is formed between June 5 and June 27. Accordingly, it is not possible to put a precise date on when these scarps formed, as

they are in shadow in the Rosetta images obtained in this time period. Therefore, both scarp “e” and scarp “f” form in regions of high mass-loss rate during phase 1, which may be more representative of their actual initiation date. Phase 3 includes the emergence of two new scarps: “a” and “k.” While scarp “j” does not directly emerge from a location of high mass-loss rate in phase 3, it is location is correlated with a region with high mass loss over the previous phase (Figure 6). Scarp “k” seems to be an anomaly, since it emerges in a region with relatively low mass loss in both phase 2 and phase 3. However, it is important to note that scarp “c” previously migrated through this region, eroding off a layer of unknown thickness. Hence, our two-layer thermal model assumption may not completely hold in this case. Phase 4 involves the commencement of migration of preexisting scarps “m” and “n.” Both of these scarps lie in regions with high mass loss in the period when they initiate activity (Figure 6). Scarps “o” and “l” also emerge from regions of high mass loss during phase 4. Finally, scarp “p” also emerges from a region of high mass loss during phase 5.

While it is not possible to predict an exact point of origin for the scarps, we can get a reliable estimate by assuming the scarps expand radially. So the scarps may not originate exactly from regions of high mass loss predicted by our model, but they do appear to at least emerge in close proximity to such regions. Based on these correlations, we conclude that our thermal model accurately predicts sublimation-driven scarp migration within the smooth terrains of Imhotep. The fact that we have not assumed any ice-enhanced hot spots on the surface while setting up our model indicates that scarp activity is largely driven by influences from local topography. Any differences in composition would then be due to activity postemplacement. Regions for which our model predicts a high mass flux but no activity is observed may indicate areas where the two-layer assumption for our model fails; volatiles might lie deeper below the surface (i.e., there may be a larger dust cover)
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Thermal model setup. (a) The surface is divided into cliffs and plains, and a shape analysis is performed to calculate the total thermal input for every facet. A receiver facet (indicated in red) within a small portion of Imhotep on the SHAP4S model is highlighted in this example; this calculation is repeated for all facets on the shape. (b) We assume a two-layered structure under each facet where a volatile-rich layer (dark blue) is blanketed by a more ice-depleted (though not ice-free) top layer (light blue).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Correlations between regions with high mass-loss rates and observed activity. (a)–(e) Georeferenced images (the date the image was acquired can be found on the top left of each panel) of the Imhotep region with labeled scarps as seen in Figure 2. (f)–(j) Same images as in panels (a)–(e). The regions shaded in yellow are regions of the highest integrated mass loss (80th percentile or higher) within the selected period of time (f) 2015 May 24–June 5, (g) 2015 June 1–June 27, (h) 2015 June 23–July 11, (i) 2015 July 7–August 6, and (j) 2015 August 2–August 23). Stars indicate the points of origin for scarps originating within the given time period, with the colors of the stars the same as the colors of the scarp labels in panels (a)–(e).</figcaption></figure>

than we initially assume for our calculations, which may result
from activity in the region (both erosion and deposition)
disturbing the distribution of the material.

To further validate our model, we compare the total water
production predicted from Imhotep by our model to that
predicted by a direct simulation Monte Carlo model called
UltraSPARTS (Wu & Lian 2003; Wu et al. 2004; Wu & Tseng 2005; Su 2013; Marschall et al. 2016, 2020). Ultra-SPARTS<fnref n="8" /> simulates the gas flow from the surface into the inner gas coma and provides a direct connection between local gas densities in the coma and the gas flux at the surface. The gas production rate of the comet can be constrained by comparing the simulated gas densities with the in situ measurements from the Rosetta Orbiter Spectrometer for Ion and Neutral Analysis (ROSINA). Using this method, Marschall et al. (2020) determined the global gas production rate of the comet for the duration of the Rosetta mission. They assumed that there is no regional heterogeneity of the activity distribution. Although 67P does have well-known regional heterogeneities (e.g., Marschall et al. 2016; Fougere et al. 2016; Marschall et al. 2019; Combi et al. 2020), they do not play a significant role in determining the global gas production rate (Marschall et al. 2020). Here we use the global gas production rate calculated by Marschall et al. (2020) to calculate the total mass loss from Imhotep. During the period between 2015 May 20 and 2016 January 15, the contribution from the Imhotep region alone amounted to a volatile mass loss of (3.3 ± 0.6) × 10⁸ kg. The total water production from Imhotep during this period predicted by our model is 3.6 × 10⁸ kg. The two values agree within the errors, thus providing an important independent check on our model through coma measurements. Moreover, the parameters used in our model predict a thermal diffusivity of 3.5 × 10⁻<fnref n="8" /> m² s⁻¹, which is consistent with the thermal diffusivity of other comets (Festou et al. 2004; Steckloff et al. 2021; Lisse et al. 2022).

# 4. Topography and Evolution of Imhotep's Triple-scarp System

## 4.1. Unique Observations in the Central Deposit of Imhotep

Scarps “h,” “a,” and “d” in the central deposit of Imhotep are unlike any other observed within Imhotep or elsewhere on 67P. Specifically, while some scarps initiate at the same location as other scarps (e.g., scarps “b” and “l,” “e” and “p,” and “p” and “k”), they are always separated by multiple weeks in time, often with a depositional event in the intervening time. For example, scarps “b” and “l” both start from the same location near a cliff on the western margin of the central deposit, migrate eastward away from the cliff, and are separated in time by multiple weeks with a clear depositional event in the interim period.

Though scarps “h,” “a,” and “d” originate from the same location within the central deposit, they all begin their migration at slightly different times with no clear deposition occurring in between start dates. They also all migrate in different directions and at different rates. This behavior is unexpected given that all three scarps start from the same exact point and should thus have similar heat inputs. We suspect that additional physical processes may need to be considered.

Specifically, scarp “a” marks the beginning of activity in Imhotep, radially expanding from the preexisting scarp “h” on June 5. Scarp “a” migrates toward the east at a rate consistent with other scarps in the region at the time (Figure 4(b)). Scarp “d” begins from nearly the same location as scarp “a,” delayed by about 1 week (June 13), expanding in the opposite direction and at a somewhat reduced rate (Figures 3(a) and 4(b)). This is

then followed by the migration of scarp “h” itself, which begins moving nearly 3 weeks after scarp “a” (June 27) and moves the slowest of all three scarps, also in a westward direction (Figures 3(b) and 4(b)).

To further analyze the three-scarp system and its surrounding terrain, we generated photoclinometric digital terrain models (DTMs) of the area using the techniques outlined in Tang et al. (2019). The process involves extracting a low-resolution topographic map from the shape model for our region of interest and generating a synthetic image (shading) corresponding to this topographic map (under specified viewing geometries) by assuming uniform photometric properties over the surface. This is followed by iteratively adjusting the topographic map to maximize the agreement between the synthetic and original images, in this case, one that was acquired on June 27. The outcome is a high-resolution DTM that matches the pixel scale of the original image (Figure 7(a)).

Using the derived DTM, we calculate the height of each scarp in the region on 2015 June 27. Scarp “h” shows the greatest depth of ~3–4 m, followed by scarp “a” with depths of ~1–1.5 m and then scarp “d” at ~0.8 m (Figures 7(a) and (b)). These results suggest that scarp migration rates are not principally controlled by scarp depth. While the shallower of the two westward-migrating scarps moves faster, both of the westward-moving scarps still migrate slower than the eastward scarp.

## 4.2. Mass-loss Rate

We define an instantaneous volume loss rate per unit width (V) along the scarp as the migration rate (v) multiplied by the scarp height (H) (as defined in Figure 7(b)). This rate relates the instantaneous scarp retreat rate for each scarp, a relevant quantity for future modeling on the dynamics of these features. We find that the instantaneous volume loss for scarp “a” is on the order of 0.24–0.36 m³ m⁻¹ hr and scarp “h” is 0.18–0.24 m³ m⁻¹ hr, with scarp “d” removing the smallest amount per unit time, 0.08 m³ m⁻¹ hr. While scarps “h” and “a” may have similar volume loss rates, scarp “d” is significantly less efficient, despite being similar in height to scarp “a.”

We can also integrate these rates along the entire scarp boundary (L), as measured on June 27, to get an instantaneous mass-loss rate for each scarp. In this order-of-magnitude exercise, we assume that the scarp depths remain largely constant along each scarp boundary and adopt a substrate material density to match that of the bulk nucleus (ρ ≃ 535 kg m⁻³; Pätzold et al. 2018), as we lack knowledge of the deposit density. We then calculate an instantaneous mass-loss rate (dust and ice) of

$$ML = \rho * v * H * L. \quad (4)$$

We calculate mass-loss rates of 11.2–16.9 kg s⁻¹ for scarp “a,” 6–8 kg s⁻¹ for scarp “h,” and 1.7 kg s⁻¹ for scarp “d.” This shows that the discrepancies in total mass lost from the surface between the scarps is even more dramatic than the volume loss per unit width, with scarp “a” being responsible for more mass loss than the other two scarps combined.

We interpret the discrepancy between V and ML calculated for the scarps to result from geometric effects related to how the features evolve. Specifically, the scarps all appear to expand radially, often from a single point (Section 2). As such, the scarp boundary length scales as L ~ v, with the mass-loss rate therefore scaling as ML ~ v². This implies that faster-moving

<footnote n="8">⁸ Commercialized by http://plasmati.com.tw/.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Topography of the triple-scarp system in Imhotep. (a) DTM for the region encompassing the triple scarp on 2015 June 27. (b) Vertically exaggerated line profile (taken along the dashed line in panel (a)) showing the topography of the triple-scarp system. The slope along the scarp face is ~10°. The stars along the profile indicate corresponding points in panel (a). The average velocity over the duration of migration (v) and the height of each scarp (H) are also indicated. Note that the ice-depleted layer is far thinner than the scarp depth (see Figure 10).</figcaption></figure>

scarps are exposing more of the active scarp wall (Birch et al. 2019) to erosion at any given time, thereby releasing more total mass per unit time. This positive feedback continues until a given scarp has swept through an entire region, consuming all available sublimable ice.

### 4.3. Hypothesis and Interpretations

We suspect that these differences result from four possible factors: (1) differences in local illumination conditions, (2) spatial heterogeneities in the leftover ice volume from 67P's previous apparition and migration of scarps, (3) stratigraphic effects, and/or (4) local small-scale slope instabilities.

First, using our thermal model, we calculate the difference in the incoming solar radiation within two boxes of the same size, one covering the region through which scarp “a” migrates and the other the region through which scarps “h” and “d” migrate (Figure 8). Our calculations reveal that the area through which scarp “a” migrates (cyan box) receives only marginally more solar insolation (about $2 \text{ W m}^{-2}$ more), averaged over the course of the day, than the region through which scarps “h” and “d” migrate (pink box). This small difference in solar illumination is likely insufficient to cause the large time differences in the start of migration of the three scarps.

Second, the difference in the onset of activity could also point toward a difference in the initial thickness of the ice-depleted layer on either side of scarp “h” due to ongoing evolution and the landscape’s “memory.” Specifically, we interpret the initial appearance of scarp “h” (and, by extension, scarps “m” and “n”) as a remnant, final erosional event from the comet’s previous orbit, analogous to the final scarps we

observed to migrate during 67P’s current passage (e.g., “k,” “j,” and “p”). As such, the western portion of scarp “h” would have escaped any postemplacement erosion and been exposed on the surface for a greater time, up to many months, compared to the eastern portion, which would have lost overlying material as the scarp moved through the region. If this is correct, the western side would have spent a longer period of time on the surface than the eastern side, allowing it to accumulate a thicker ice-depleted upper layer. While we assume that ice exists at all depths (Figure 10; Section 4), we hypothesize here that the sublimation of some CV of ice is necessary to begin the positive feedback cycle of scarp erosion (this is corroborated on other comets by past works such as Belton & Melosh 2009; Blum et al. 2014). Specifically, assuming that ice sublimation is from within the particles themselves, once a sufficiently large volume of particles is activated (denoted by the white line reaching the yellow line in Figure 9), the entire volume itself can be mobilized; the exact details of this process are beyond the scope of this work. This instantaneous removal of material, exposure of previously hidden ice-rich materials, and generation of new topography then begins a positive feedback cycle, whereby the scarps cascade through neighboring regions until an entire layer is removed or new sediment is deposited onto the region. Once a layer is removed (or a new layer is added), the process resets (i.e., the white line moves upward in Figure 9). This may explain why multiple scarps occur in the same region (e.g., scarps “b” and “l”) but could also further delay the onset of activity. Therefore, this “memory” that the landscape has due to past/ongoing erosion and deposition results in a dynamic
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Variability of direct solar input within the triple-scarp region. (a) Cropped image of the triple-scarp region taken on 2014 September 5. The region remains unaltered until 2015 June 3, when scarp “a” starts migrating. (b) Modeled direct solar input within the triple-scarp region on 2015 June 1. The cyan box covers the region through which scarp “a” migrates; the pink box covers the area through which scarps “d” and “h” migrate.</figcaption></figure>

evolution of the surface, where thermal energy is continuously battling to activate a CV. Under this framework, a shallower depth to activate a CV on the eastern side (i.e., the yellow line is not as deep in Figure 9) could allow scarp “a” to initiate before scarps “d” or “h,” beginning the cascade process earlier.

Third, the slow migration rate of scarp “h” may also point to stratigraphic differences in the volume of ice within the local subsurface. Specifically, scarps “d” and “a” may only be removing an entire seasonally processed layer of material (Figure 10), the thickness ($L_{seasonal}$) of which is simply the thermal skin depth, given by Macher et al. (2019),

$$L_{seasonal} = \sqrt{\frac{kP_{comet}}{\pi\rho c_{p}}}, \qquad (5)$$

which is related to the density ($\rho$), thermal conductivity ($k$), and heat capacity ($c_p$) of the material and time that the material has been exposed on the surface ($P_{comet}$). For values estimated at 67P, this amounts to $L_{seasonal} \simeq 1.5$ m if we assume that such sediment was deposited at perihelion and reactivated the following orbit. This depth is similar to the scale of scarps “a” and “d.” With less overall ice (denoted by shades of blue in Figure 10) required to sublimate to remove a given volume of material (both ice and dust; Birch et al. 2019) after the scarp migration process has been initiated, scarps “d” and “a” will therefore migrate faster. Scarp “h,” meanwhile, is sublimating through deeper, more ice-rich layers (darkest blue layers in Figure 10), consuming more energy per unit distance as the scarps migrate. Since it must sublimate through greater and greater ice volumes, scarp “h” would therefore migrate at the slowest rate of all three scarps (Birch et al. 2019).

Finally, scarp “d” could have initiated as a result of a slope instability (loosely analogous to pore pressure-driven shallow landsliding on Earth; e.g., Montgomery & Dietrich 1994) along the boundary of the initial scarp “h.” This would instantaneously remove the ice-depleted upper layer on the western side of scarp “h,” exposing the underlying CV layer of ice,

thereby initiating the migration process. This could also result in more refractory material accumulating at the base of scarp “h,” which would further delay its initiation and slow its migration (Figure 7(b)). Scarp “d” would therefore be a side effect of the large initial height of scarp “h,” where we would otherwise expect two scarps radially expanding away from each other.

Through this analysis, we conclude that all of the above processes are occurring, though not equally for all scarps. Scarp “d” could have been aided by a slope instability and then had its migration occur in a similar manner to scarp “a,” though slower, due to both a slightly more processed upper layer and slightly cooler temperatures. Similarly, scarp “h” can be slowed by higher ice fractions, with its onset similarly delayed by less favorable illumination conditions and greater depths to the CV ice layer.

## 5. Implications for Cometary Surface Evolution

Our complete documentation of changes within the Imhotep region of comet 67P highlights the fact that the surface is rapidly cycling from being net erosional to net depositional. Specifically, scarp activity starts in Imhotep on 2015 June 5, 2 months before perihelion, with the initial month of activity localized between 3° N and 14°S. This is followed by deposition of sediment north of 14°S in the weeks leading up to perihelion, as evidenced by the partial or complete burial of scarps in this region. Simultaneously, activity continues to spread toward the more southern parts of the central deposit, with several new scarps mobilizing as far as 35°S. After perihelion, scarp-driven erosion continues across all latitudes spanned by the central deposit, with deposition of material observed in tandem, often within a few hundred meters of eroding scarps. Therefore, this complicated behavior at Imhotep suggests finer-scale dynamics, both spatially and temporally, than described in previous studies (Thomas et al. 2015; El-Maarry et al. 2017; Keller et al. 2017; Birch et al. 2019), and that to understand the overall evolution of any given smooth terrain deposit, we need to
<figure><img src="image_12.png" /></figure>

**Figure 9.** The propagation of a thermal heat wave into the smooth terrains. The dotted line indicates the thermal heat wave, while the yellow line indicates the depth the thermal wave needs to propagate to in order to initiate the positive feedback cycle of scarp migration. Once the thermal wave reaches this depth, the yellow line resets, and the cycle needs to repeat in order to reinitiate scarp migration.

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

**Figure 10.** Hypothesized stratigraphy of ices within the smooth terrains. The particles closest to the surface are the most ice-depleted, and the ice content increases with depth. The figure shows three discrete subsurface ice concentrations, while the actual distribution will smoothly vary. The dashed orange, red, and green lines represent the progression of scarps "a," "d," and "h," respectively.

observe activity at a finer temporal scale than just before/after perihelion.

Our observations also suggest that the fallback particles themselves are icy, consistent with previous work (Davidsson et al. 2021), and that activity is likely not primarily controlled

by surface and/or subsurface heterogeneities resulting from 67P's early evolution (Agarwal et al. 2017). First, we observe that portions of the Imhotep basin that receive sediment fallback also serve as points of origin for scarps. Second, we observe that yet more smooth terrains underlie regions where
erosion recently took place. With no underlying, putatively ice-
rich bedrock, this suggests that ice must be retained within the
particles themselves as they land in the region. Any subsequent
activity is due to the sublimation of those embedded ices, the
location and timing of which primarily depend on the
surrounding topography.

Specifically, this topographic control is revealed by our thermal
model, which reliably predicts when and where activity begins in
the region and matches independent observations from ROSINA.
Even though we assume a uniform composition across the
surface, our model still predicts “hot spots” of increased activity,
the location of which are primarily controlled by influences from
local topography. We also observe that all scarps originate from
topographic discontinuities, further suggesting a topographic
influence on when and where activity begins.

However, while the influences from topography can explain
a majority of the observed dynamics, the unique behavior of the
triple-scarp system in Imhotep suggests further complexity. To
explain the triple scarp’s evolution, we hypothesize that the
landscapes may have a “memory.” We propose that, in addition
to lateral and vertical heterogeneities in the local subsurface
volatile content that formed during the comet’s formation,
variations may also reflect previous cycles of erosion/
deposition. As deposition and erosion occur simultaneously,
whether a given location is capable of removing material
depends on the topography, whether it has received a new layer
of material, and how long that material resides within a given
hot spot. Therefore, to initiate the positive feedback cycle of
scarp erosion, the sublimation of a CV of ice is required, which
in turn depends upon how the volatile content changes with
depth and whether the surface lies within one of the
topographically controlled hot spots at any given time.

Future work should couple our thermal modeling with
landscape evolution so as to explain the dynamic evolution of
these granular deposits. With smooth terrains observed on all other

comets for which we have resolved surface images (Sunshine et al.
2016), we can use the understanding and tools we have developed
here at the Imhotep region of 67P to understand the broader
evolution of cometary smooth terrains. Finally, as the smooth
terrains are likely targets of upcoming sample return missions, our
work greatly benefits such concepts by constraining the evolution
and context of any returned materials.

In summary, our work demonstrates that local-scale activity
and sediment transport are equally important to the global-scale
south-to-north (Keller et al. 2017) transport of sediment.
Further, these local effects are not necessarily related to
variations in near-surface volatile contents (Oklay et al. 2016;
Fornasier et al. 2016; Agarwal et al. 2017). Instead, reradiation
from topography alone can create local enhancements in
activity. Meanwhile, lateral and vertical spatial variations in the
local ice content are heavily influenced by the comet’s
evolution, which we term the “memory.” Activity from 67P’s
previous orbit has a significant effect on the evolution of the
surface on the subsequent passage. Together, these local-scale
changes must be accounted for so as to understand the physical
and chemical evolution of cometary surfaces.

This research was supported by the Rosetta Data Analysis
Program (No. 80NSSC19K1307) and the Heising-Simons
Foundation (51 Pegasi b Fellowship to S.P.D.B). This research
has made use of the scientific software shapeViewer (www.
comet-toolbox.com). All OSIRIS image data presented in this
paper were downloaded and are freely available on ESA's
Planetary Science Archive (https://archives.esac.esa.int/psa)
and NASA's Planetary Science Data System (https://pds.nasa.
gov). Topographic and spectrophotometric data that support the
plots within this paper and other findings of this study are
available in the appendices or from the corresponding author
upon request.
## Appendix A
List of Images Used for Analysis

Table 1 contains a list of 26 Rosetta images directly used in this work. Twenty-four of these images (acquired between 2014 September–2016 July) were used to make direct measurements of the scarp migration rates and directions.

<table><caption>**Table 1**
OSIRIS NAC Images Used for Analysis</caption><thead><tr><th>Image ID</th><th>Time of Image Acquisition (MM/DD/YYYY hh:mm:ss)</th><th>Target-center Distance (km)</th><th>Pixel Scale (m pixel<sup>-1</sup>)</th><th>Figure</th></tr></thead><tbody><tr><td>N20140803T052116161ID30F82</td><td>08/03/2014 05:22:24</td><td>333.4</td><td>6.2</td><td>Figure 1(a)</td></tr><tr><td>N20140903T034422640ID30F22</td><td>09/03/2014 03:45:32</td><td>67.2</td><td>1.3</td><td>Figure 8</td></tr><tr><td>N20150524T013440826ID30F22</td><td>05/24/2015 01:35:57</td><td>165.7</td><td>3.1</td><td>Figures 2(a), 6(a) and (f)</td></tr><tr><td>N20150605T070306479ID20F23</td><td>06/05/2015 07:04:23</td><td>206.1</td><td>3.9</td><td>Figure 2(b)</td></tr><tr><td>N20150613T134156317ID20F22</td><td>06/13/2015 13:43:13</td><td>200.3</td><td>3.7</td><td>NA</td></tr><tr><td>N20150627T132707811ID20F22</td><td>06/27/2015 13:28:25</td><td>194.2</td><td>3.6</td><td>Figures 2(c), 6(b) and (g)</td></tr><tr><td>N20150701T035546768ID20F22</td><td>07/01/2015 03:57:03</td><td>159.9</td><td>3</td><td>NA</td></tr><tr><td>N20150702T053556341ID20F22</td><td>07/02/2015 05:37:13</td><td>161.1</td><td>3</td><td>Figure 2(d)</td></tr><tr><td>N20150711T124159342ID20F23</td><td>07/11/2015 12:43:16</td><td>161.4</td><td>3</td><td>Figures 2(e), 6(c) and (h)</td></tr><tr><td>N20150715T151802760ID30F22</td><td>07/15/2015 15:19:20</td><td>164.8</td><td>3.1</td><td>NA</td></tr><tr><td>N20150726T234311743ID20F22</td><td>07/26/2015 23:44:29</td><td>169</td><td>3.2</td><td>NA</td></tr><tr><td>N20150729T130610758ID20F22</td><td>07/29/2015 13:07:28</td><td>187</td><td>3.5</td><td>NA</td></tr><tr><td>N20150731T081942769ID20F22</td><td>07/31/2015 08:21:00</td><td>194.8</td><td>3.6</td><td>NA</td></tr><tr><td>N20150801T212507780ID30F22</td><td>08/01/2015 21:26:25</td><td>209.4</td><td>3.9</td><td>Figures 2(f), 6(d) and (i)</td></tr><tr><td>N20150805T120150769ID30F22</td><td>08/05/2015 12:03:08</td><td>255.3</td><td>4.8</td><td>NA</td></tr><tr><td>N20150806T122006031ID20F23</td><td>08/06/2015 12:21:24</td><td>267.8</td><td>5</td><td>NA</td></tr><tr><td>N20150823T170124554ID20F22</td><td>08/23/2015 17:02:43</td><td>339.8</td><td>6.4</td><td>Figure 2(g)</td></tr><tr><td>N20150904T053045720ID20F22</td><td>09/04/2015 05:32:04</td><td>397</td><td>7.4</td><td>NA</td></tr><tr><td>N20150915T175211532ID20F22</td><td>09/15/2015 17:53:30</td><td>331</td><td>6.2</td><td>NA</td></tr><tr><td>N20151026T125738791ID20F22</td><td>10/26/2015 12:58:59</td><td>312.7</td><td>5.9</td><td>Figure 2(h)</td></tr><tr><td>N20151104T213925829ID20F22</td><td>11/04/2015 21:40:46</td><td>243</td><td>4.5</td><td>NA</td></tr><tr><td>N20151129T023654724ID30F22</td><td>11/29/2015 02:38:16</td><td>122.9</td><td>2.3</td><td>NA</td></tr><tr><td>N20151206T211312655ID30F22</td><td>12/06/2015 21:14:34</td><td>99.6</td><td>1.9</td><td>NA</td></tr><tr><td>N20160123T234521190ID30F22</td><td>01/23/2016 23:46:44</td><td>74.7</td><td>1.4</td><td>Figures 2(i), 6(e) and (j)</td></tr><tr><td>N20160609T193002728ID30F22</td><td>06/09/2016 19:31:29</td><td>30.6</td><td>0.6</td><td>NA</td></tr><tr><td>N20160730T120845524ID30F22</td><td>07/30/2016 12:10:13</td><td>9</td><td>0.2</td><td>Figure 1(d)</td></tr></tbody></table>

Note. An OSIRIS WAC image, W20140905T062922570ID30F13, was used for Figure 1(b). Some images used for analysis do not appear in the figures.
## Appendix B
### Calculating the Average Migration Rates for Scarps

As Rosetta did not image the surface continuously, we can only observe the progression of scarps intermittently whenever images were acquired. This means that we do not know precisely the instantaneous migration rate. Rather, we calculate an average migration rate using the total distance traveled by a scarp between two images where its position demonstrated a measurable change. We note that this only gives us a minimum rate of migration for the scarps, since they may be inactive for long intervals between any two given images. However, such

rates are of value, in particular when compared scarp-to-scarp. It is these relative rates that are of most interest in our work.

Table 2 shows the calculation of the average migration rates for each individual scarp. The mean distance covered is calculated by taking the average of three independent measurements of the distance covered by a scarp relative to its previous observed position. While the migration rate column also calculates an average migration rate for the scarp between two images, a net average migration rate is also calculated for each scarp from the total distance it covered over its period of activity. These rates are noted in the main text.

<table><caption>**Table 2**
Measurements of Scarp Migration Rates</caption><thead><tr><th>Scarp</th><th>Time of Image Acquisition (MM/DD/YYYY hh:mm)</th><th>Mean Distance Covered (m)</th><th>Time Elapsed (days)</th><th>Migration Rate (cm hr<sup>-1</sup>)</th><th>Average Migration Rate (cm hr<sup>-1</sup>)</th></tr></thead><tbody><tr><td rowspan="6">Scarp "a"</td><td>05/24/2015 01:35</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>06/05/2015 07:04</td><td>52</td><td>12.23</td><td>18</td><td>...</td></tr><tr><td>06/13/2015 13:43</td><td>35</td><td>8.28</td><td>17</td><td>23</td></tr><tr><td>06/27/2015 13:28</td><td>81</td><td>13.99</td><td>24</td><td>...</td></tr><tr><td>07/01/2015 03:57</td><td>27</td><td>3.6</td><td>31</td><td>...</td></tr><tr><td>07/02/2015 05:37</td><td>20</td><td>1.07</td><td>80</td><td>...</td></tr><tr><td rowspan="4">Scarp "b"</td><td>06/13/2015 13:43</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>06/27/2015 13:28</td><td>93</td><td>13.99</td><td>28</td><td>24</td></tr><tr><td>07/01/2015 03:57</td><td>27</td><td>3.6</td><td>31</td><td>...</td></tr><tr><td>07/11/2015 12:43</td><td>40</td><td>10.36</td><td>16</td><td>...</td></tr><tr><td rowspan="3">Scarp "c"</td><td>05/24/2015 01:35</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>06/05/2015 07:04</td><td>19</td><td>12.23</td><td>6</td><td>10</td></tr><tr><td>06/27/2015 13:28</td><td>68</td><td>22.27</td><td>13</td><td>...</td></tr><tr><td rowspan="6">Scarp "d"</td><td>06/05/2015 07:04</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>06/13/2015 13:43</td><td>16</td><td>8.28</td><td>8</td><td>...</td></tr><tr><td>06/27/2015 13:28</td><td>35</td><td>13.99</td><td>10</td><td>14</td></tr><tr><td>07/01/2015 03:57</td><td>14</td><td>3.6</td><td>16</td><td>...</td></tr><tr><td>07/02/2015 05:37</td><td>24</td><td>1.07</td><td>94</td><td>...</td></tr><tr><td>07/11/2015 12:43</td><td>33</td><td>9.3</td><td>15</td><td>...</td></tr><tr><td rowspan="3">Scarp "e"</td><td>06/27/2015 13:28</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>07/02/2015 05:37</td><td>43</td><td>4.67</td><td>39</td><td>22</td></tr><tr><td>07/11/2015 12:43</td><td>32</td><td>9.3</td><td>14</td><td>...</td></tr><tr><td rowspan="6">Scarp "g"</td><td>07/02/2015 05:37</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>07/11/2015 12:43</td><td>51</td><td>9.3</td><td>23</td><td>...</td></tr><tr><td>07/15/2015 15:19</td><td>28</td><td>4.11</td><td>28</td><td>25</td></tr><tr><td>07/26/2015 23:44</td><td>56</td><td>11.35</td><td>20</td><td>...</td></tr><tr><td>07/29/2015 13:07</td><td>29</td><td>2.56</td><td>48</td><td>...</td></tr><tr><td>07/31/2015 08:21</td><td>14</td><td>1.8</td><td>33</td><td>...</td></tr><tr><td rowspan="4">Scarp "h"</td><td>06/13/2015 13:43</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>06/27/2015 13:28</td><td>20</td><td>13.99</td><td>6</td><td>7</td></tr><tr><td>07/11/2015 12:43</td><td>29</td><td>13.97</td><td>9</td><td>...</td></tr><tr><td>07/26/2015 23:44</td><td>27</td><td>15.46</td><td>7</td><td>...</td></tr><tr><td rowspan="3">Scarp "i"</td><td>07/02/2015 05:37</td><td>...</td><td>...</td><td>...</td><td>22</td></tr><tr><td>07/11/2015 12:43</td><td>43</td><td>9.3</td><td>19</td><td>...</td></tr><tr><td>07/15/2015 15:19</td><td>28</td><td>4.11</td><td>28</td><td>...</td></tr><tr><td rowspan="6">Scarp "j"</td><td>07/11/2015 12:43</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>07/15/2015 15:19</td><td>25</td><td>4.11</td><td>25</td><td>...</td></tr><tr><td>07/26/2015 23:44</td><td>92</td><td>11.35</td><td>34</td><td>...</td></tr><tr><td>07/31/2015 08:21</td><td>45</td><td>4.36</td><td>43</td><td>18</td></tr><tr><td>08/05/2015 12:03</td><td>75</td><td>5.15</td><td>60</td><td>...</td></tr><tr><td>08/23/2015 04:39</td><td>108</td><td>17.69</td><td>26</td><td>...</td></tr><tr><td></td><td>09/04/2015 04:42</td><td>36</td><td>12</td><td>12</td><td>...</td></tr></tbody></table>
<table><caption>**Table 2**
(Continued)</caption><thead><tr><th>Scarp</th><th>Time of Image Acquisition (MM/DD/YYYY hh:mm)</th><th>Mean Distance Covered (m)</th><th>Time Elapsed (days)</th><th>Migration Rate (cm hr<sup>-1</sup>)</th><th>Average Migration Rate (cm hr<sup>-1</sup>)</th></tr></thead><tbody><tr><td></td><td>10/26/2015 11:44</td><td>76</td><td>52.29</td><td>6</td><td>...</td></tr><tr><td rowspan="6">Scarp "k"</td><td>07/11/2015 12:43</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>07/26/2015 23:44</td><td>38</td><td>15.46</td><td>10</td><td>...</td></tr><tr><td>07/31/2015 08:21</td><td>22</td><td>4.36</td><td>21</td><td>4</td></tr><tr><td>09/15/2015 17:53</td><td>17</td><td>46.4</td><td>1</td><td>...</td></tr><tr><td>10/26/2015 11:44</td><td>29</td><td>40.74</td><td>3</td><td>...</td></tr><tr><td>11/29/2015 02:38</td><td>40</td><td>33.62</td><td>5</td><td>...</td></tr><tr><td rowspan="3">Scarp "l"</td><td>07/11/2015 12:43</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>07/15/2015 15:19</td><td>19</td><td>4.11</td><td>19</td><td>21</td></tr><tr><td>07/26/2015 23:44</td><td>60</td><td>11.35</td><td>22</td><td>...</td></tr><tr><td rowspan="4">Scarp "m"</td><td>07/26/2015 23:44</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>07/29/2015 13:07</td><td>13</td><td>2.56</td><td>22</td><td>25</td></tr><tr><td>08/01/2015 21:26</td><td>25</td><td>3.35</td><td>31</td><td>...</td></tr><tr><td>08/05/2015 12:03</td><td>19</td><td>3.61</td><td>22</td><td>...</td></tr><tr><td rowspan="4">Scarp "n"</td><td>07/26/2015 23:44</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>07/29/2015 13:07</td><td>21</td><td>2.56</td><td>34</td><td>37</td></tr><tr><td>08/01/2015 21:26</td><td>31</td><td>3.35</td><td>39</td><td>...</td></tr><tr><td>08/05/2015 12:03</td><td>32</td><td>3.61</td><td>37</td><td>...</td></tr><tr><td rowspan="4">Scarp "o"</td><td>07/26/2015 23:44</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>07/29/2015 13:07</td><td>21</td><td>2.56</td><td>34</td><td>...</td></tr><tr><td>07/31/2015 08:21</td><td>21</td><td>1.8</td><td>48</td><td>45</td></tr><tr><td>08/01/2015 21:26</td><td>16</td><td>1.54</td><td>43</td><td>...</td></tr><tr><td rowspan="3">Scarp "p"</td><td>08/23/2015 04:39</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>09/15/2015 17:53</td><td>71</td><td>23.55</td><td>12</td><td>6</td></tr><tr><td>11/04/2015 21:40</td><td>36</td><td>50.16</td><td>3</td><td>...</td></tr></tbody></table>

## Appendix C
### Details of Simple Thermodynamic Model for Diurnal Activity

We follow the framework described in Umurhan et al. (2022). We split the thermodynamic response of a localized layer into a day phase and a night phase, where all corresponding variables are respectively denoted hereafter with “d” and “n” subscripts. As discussed in the main text, we assume that a given subsurface layer is divided up into a relatively low volatile ice-containing top layer, with corresponding thermal inertial $\mathcal{I}$, sitting atop a semi-infinite high volatile ice-containing bottom layer with $\mathcal{I}_b$ located a depth $h$ beneath the surface $z=0$. Given the definition of the thermal inertia, $\mathcal{I} = \sqrt{\rho C_p K}$, we assume that the only property distinguishing these two layers is their respective thermal conductivities $K$. As such, we assume that the thermal conductivity of the bottom layer, $K_b$, is a fraction $\xi$ of the top layer’s conductivity, $K$, i.e., $K_b = \xi K$. It then follows that $\mathcal{I}_b = \sqrt{\xi} \cdot \mathcal{I}$.

The subsurface temperature profile during night is modeled by

$$T_n = T_{\text{int}} + \begin{cases} (T_{\text{sn}} - B_n - T_{\text{int}})e^{kz} + B_n e^{-kz}; & -h < z < 0; \\ C_n e^{k_b(z+h)}; & z < -h; \end{cases} \quad (\text{C1})$$

where $T_{\text{int}}$, $T_{\text{sn}}$, $B_n$, and $C_n$ are unknown constants to be determined via satisfaction of certain matching and boundary conditions. The variable $T_{\text{sn}}$ denotes the nighttime surface temperature. The thermal skin depths in the top and bottom layers are $\ell$, $\ell_b$, respectively, and are related to the corresponding wavenumbers by

$$\frac{2\pi}{\ell} \equiv k = \sqrt{\frac{\rho C_p \omega}{K}}; \quad \frac{2\pi}{\ell_b} \equiv k_b = \sqrt{\frac{\rho C_p \omega}{K_b}}, \quad (\text{C2})$$

where $\omega$ is the diurnal frequency $2\pi/P_{\text{day}}$, in which $P_{\text{day}} \approx 12.4$ hr. It follows that $k_b = k/\sqrt{\xi}$. By construction, the solution for $T_n$ expressed in Equation (C1) limits to $T_{\text{int}}$ for $z \to -\infty$, which we interpret as representing the deep interior temperature, as yet unknown. The dayside temperature solution follows the same form:

$$T_d = T_{\text{int}} + \begin{cases} (T_{\text{sd}} - B_d - T_{\text{int}})e^{kz} + B_d e^{-kz}; & -h < z < 0; \\ C_d e^{k_b(z+h)}; & z < -h. \end{cases} \quad (\text{C3})$$

In a similar fashion, the unknown $T_{\text{sd}}$ denotes the daytime surface temperature, and $T_d$ similarly limits to $T_{\text{int}}$ for $z \to -\infty$. The solutions across the layer interfaces at $z = -h$ must satisfy the continuity of temperature, as well as a differential flux
<table><caption>**Table 3**
Various Quantities and Parameters and Their Adopted Values Where Appropriate</caption><thead><tr><th>Symbol</th><th>Definition</th><th>Formula<sup>a</sup></th><th>Typical Value<sup>a</sup></th></tr></thead><tbody><tr><td>ω</td><td>Diurnal frequency of 67P</td><td></td><td>1.407 × 10<sup>-4</sup> s<sup>-1</sup></td></tr><tr><td>ρ</td><td>Mass density of surface materials</td><td></td><td>535 kg m<sup>-3</sup></td></tr><tr><td>C<sub>p</sub></td><td>Specific heat at constant pressure</td><td></td><td>500 J K<sup>-1</sup> kg<sup>-1</sup></td></tr><tr><td>K</td><td>Conductivity of refractory materials</td><td></td><td></td></tr><tr><td>K<sub>b</sub></td><td>Conductivity of matrix ice and scarps</td><td></td><td></td></tr><tr><td>ξ</td><td>Ratio of conductivities</td><td>K<sub>b</sub>/K</td><td>2.25</td></tr><tr><td>I</td><td>Thermal inertia of plains</td><td>√ρC<sub>p</sub>K</td><td>50 tiu</td></tr><tr><td>I<sub>b</sub></td><td>Thermal inertia of scarps and ice</td><td>√ξ · I</td><td>80 tiu</td></tr><tr><td>h</td><td>Thickness of refractory cover</td><td></td><td>~0.01 m</td></tr><tr><td>k</td><td>Thermal decay wavenumber</td><td>k = √ρC<sub>p</sub>ω/K</td><td></td></tr><tr><td>l</td><td>Diurnal skin depth</td><td>l = 2π/k</td><td>l ~ 0.02 m</td></tr><tr><td>T<sub>sn</sub></td><td>Nighttime surface temperature</td><td></td><td></td></tr><tr><td>T<sub>sd</sub></td><td>Daytime surface temperature</td><td></td><td></td></tr><tr><td>T<sub>noon</sub></td><td>Noontime temperature of sublimating layer</td><td></td><td></td></tr><tr><td>T<sub>0</sub></td><td>Reference temperature</td><td></td><td>220, 240 K<sup>b</sup></td></tr><tr><td>P<sub>0</sub></td><td>Reference H<sub>2</sub>O vapor pressure at T = T<sub>0</sub></td><td></td><td>3, 30 Pa<sup>b</sup></td></tr><tr><td>f<sub>comb</sub></td><td>Cosmic microwave background flux</td><td></td><td>3 × 10<sup>-6</sup> W m<sup>-2</sup></td></tr><tr><td>f<sub>☉</sub></td><td>Impingent daytime solar flux</td><td></td><td>100–1000 W m<sup>-2</sup></td></tr><tr><td>f<sub>bg</sub></td><td>Background radiative flux</td><td></td><td>≥f<sub>comb</sub></td></tr><tr><td>A</td><td>Albedo</td><td></td><td>~0.1</td></tr><tr><td>ε</td><td>Emissivity</td><td></td><td>~0.9</td></tr><tr><td>L</td><td>Enthalpy of sublimation for H<sub>2</sub>O</td><td></td><td>2.83 × 10<sup>6</sup> J kg<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>k</td><td>Boltzmann constant</td><td></td><td>1.38 × 10<sup>-23</sup> J K<sup>-1</sup></td></tr><tr><td>m<sub>H2O</sub></td><td>Mass of H<sub>2</sub>O atom</td><td></td><td>2.98 × 10<sup>-26</sup> kg</td></tr><tr><td>T<sub>a</sub></td><td>Activation temperature</td><td>Lm<sub>H2O</sub>/k</td><td>6208 K</td></tr><tr><td>σ</td><td>Stefan–Boltzmann constant</td><td></td><td>5.67 × 10<sup>-8</sup> W K<sup>-4</sup> m<sup>-2</sup></td></tr><tr><td>T<sub>b</sub></td><td>Daytime brightness temperature</td><td>[ (1-A)f<sub>☉</sub> ]<sup>1/4</sup> / (I√ωT<sub>b</sub>)</td><td></td></tr><tr><td>Γ</td><td>Spencer number</td><td></td><td></td></tr></tbody></table>

**Notes.**

a Where applicable.

b Based on fit to Fray & Schmitt (2009) valid in range of 180 K < T < 260 K.

condition that takes into account the amount of energy consumed in sublimating H<sub>2</sub>O ice from the interfacial surface. We quantify the latter effect using the formulation found in Hobley et al. (2018), where by exploiting the vapor pressure's Arrhenius temperature dependence, it was shown that the diurnally averaged sublimation—the majority of which occurs in the moments around noon with peak noontime temperature T<sub>noon</sub>—can be analytically approximated in terms of an averaged sublimated gas mass flux (Σ̇ in units of kg m<sup>-2</sup> s<sup>-1</sup>) given by

$$ \Sigma = \frac{P_{\text{vap}}(T_{\text{noon}})}{2\pi\sqrt{L}}, \qquad (\text{C4}) $$

where L is the enthalpy of sublimation, here for H<sub>2</sub>O ice, in which L ≈ 2.83 × 10<sup>6</sup> J kg<sup>-1</sup> K<sup>-1</sup>, and P<sub>vap</sub>(T) is H<sub>2</sub>O's T-dependent vapor pressure. We analytically represent P<sub>vap</sub> in terms of the formula

$$ P_{\text{vap}} \approx P_0 \exp \left[ \frac{T_a}{T_0} - \frac{T_a}{T_{\text{noon}}} \right]; \quad T_a \equiv \frac{Lm_{\text{H2O}}}{k}, \qquad (\text{C5}) $$

where k is the Boltzmann constant, and T<sub>a</sub> is the activation temperature. To the degree of accuracy that we require, the formula in Equation (C5) satisfactorily quantifies H<sub>2</sub>O's vapor pressure, particularly in the temperature range we expect to be

significant for during Imhotep's observed sublimation activity, the range 200 K < T < 250 K. As such, we select P<sub>0</sub> ≈ 30(3) Pa, which is H<sub>2</sub>O's vapor pressure at T<sub>0</sub> = 240(220) K (Fray & Schmitt 2009). Finally, the diurnally averaged rate of energy consumed by H<sub>2</sub>O sublimation per unit area, denoted by F<sub>s</sub>, is given by

$$ F_s = L\Sigma = \frac{\sqrt{L}}{2\pi} P_{\text{vap}}(T_{\text{noon}}). \qquad (\text{C6}) $$

The solution procedure involves enforcing matching conditions across the interface at z = -h. We start with the evening-side profile by requiring that the temperatures are equal when approaching the interface from above and below, i.e.,

$$ T_n|_{z \to -h^+} = T_n|_{z \to -h^-}, \qquad (\text{C7}) $$

together with a heat flux condition,

$$ K \partial_z T_n|_{z \to -h^+} - K_b \partial_z T_n|_{z \to -h^-} = 0. \qquad (\text{C8}) $$

Similarly, we require the continuity of the temperature field at z = -h for the dayside solution,

$$ T_d|_{z \to -h^+} = T_d|_{z \to -h^-}. \qquad (\text{C9}) $$

However, the dayside flux condition differs from the nightside requirement by taking into account the irreversible consumption
of energy to drive sublimation from the interfacial layer,

$$K \partial_z T_d |_{z \to -h^+} - K_b \partial_z T_d |_{z \to -h^-} = F_s = \frac{\sqrt{L}}{2\pi} P_{\text{vap}}(T_{\text{noon}}). \quad (\text{C10})$$

The final condition requires that the total energy consumed by sublimation over the course of the day is reflected in the net energy conducted into and out of the surface ($z=0$). In other words,

$$K \partial_z T_d |_{z \to 0} + K \partial_z T_n |_{z \to 0} = F_s. \quad (\text{C11})$$

The simultaneous satisfaction of these five connection conditions (Equations (C7)–(C11)) determines expressions for five of the unknown coefficients, i.e., $C_n$, $B_n$, $B_d$, $C_d$, and $T_{\text{int}}$, in terms of the three remaining unknowns, $T_{\text{sn}}$, $T_{\text{sd}}$, and $T_{\text{noon}}$:

$$B_n = \frac{(T_{\text{int}} - T_{\text{sn}})(\sqrt{\xi} - 1)}{1 - \sqrt{\xi} + (1 + \sqrt{\xi})e^{2hk}},$$

$$B_d = \frac{(T_{\text{int}} - T_{\text{sd}})(\sqrt{\xi} - 1) - e^{hk}\phi_d}{1 - \sqrt{\xi} + (1 + \sqrt{\xi})e^{2hk}},$$

$$C_n = - \frac{T_{\text{int}} - T_{\text{sn}}}{\cosh hk + \sqrt{\xi} \sinh hk},$$

$$C_d = - \frac{T_{\text{int}} - T_{\text{sd}} + \phi_d \sinh hk}{\cosh hk + \sqrt{\xi} \sinh hk},$$

with

$$T_{\text{int}} = \frac{T_{\text{sd}} + T_{\text{sn}}}{2} - \frac{1}{2} \left( \frac{\cosh hk + \sqrt{\xi} \sinh hk - 1}{\sqrt{\xi} \cosh hk + \sinhk} \right) \phi_d, \quad (\text{C12})$$

where

$$\phi_d(T_{\text{noon}}) \equiv \frac{F_s}{kK} = \frac{\sqrt{L}P_0}{2\pi\mathcal{I}\sqrt{\omega}} \exp\left[\frac{T_a}{T_0} - \frac{T_a}{T_{\text{noon}}}\right]. \quad (\text{C13})$$

The remaining unknowns are assessed by imposing the energy balance boundary condition at the surface $z=0$. This amounts to saying that the incoming received radiation is balanced by blackbody losses to space and thermal conduction into/out of the interior. We assume that the material has an albedo $A$ and surface emissivity of $\varepsilon$. Thus, on the daytime surface, we have

$$(1 - A) f_{\odot} = \varepsilon \sigma T_{\text{sd}}^4 + K \partial_z T_d |_{z \to 0}, \quad (\text{C14})$$

where $f_{\odot}$ represents the incoming solar radiation at the comet's instantaneous orbital location, and $(1-A)f_{\odot}$ denotes the total amount actually absorbed. The thermal surface conductive flux on the dayside may be expressed utilizing the soluti on for $T_d$ thus far developed. After some algebra (facilitated by

Mathematica), the above expression may be reexpressed as

$$(1 - A)f_{\odot} = \varepsilon \sigma T_{\text{sd}}^4 + \mathcal{I} \sqrt{\omega} \left[ \frac{T_{\text{sd}} - T_{\text{sn}}}{2} \Xi(hk) + \frac{\phi_d(T_{\text{noon}})}{2} (1 + \Upsilon(hk)) \right], \quad (\text{C15})$$

where

$$\Xi(hk) \equiv \frac{\sqrt{\xi} \cosh hk + \sinhk}{\cosh hk + \sqrt{\xi} \sinhk},$$

$$\Upsilon(hk) \equiv \frac{1}{\cosh hk + \sqrt{\xi} \sinhk}. \quad (\text{C16})$$

We note that $\Upsilon > 0$ and $\Xi > 0$, and their limiting forms with respect to $0 < hk$ are

$$\lim_{hk \to 0} \Upsilon(hk) = 1 - \sqrt{\xi}hk + \mathcal{O}(hk)^2,$$

$$\lim_{hk \to \infty} \Upsilon(hk) \approx \frac{2e^{-hk}}{1 + \sqrt{\xi}} \quad (\text{C17})$$

and

$$\lim_{hk \to 0} \Xi(hk) = \sqrt{\xi} + \mathcal{O}(hk),$$

$$\lim_{hk \to \infty} \Xi(hk) \approx 1. \quad (\text{C18})$$

Similarly, the energy balance statement for the nightside is

$$(1 - A) f_{\text{bg}} = \varepsilon \sigma T_{\text{sn}}^4 + K \partial_z T_n |_{z \to 0}, \quad (\text{C19})$$

where $f_{\text{bg}}$ is the impudent background flux received on the nightside surface. This will depend upon how much of the local surface is exposed to other parts of the landscape. At minimum, $f_{\text{bg}} = f_{\text{cmb}}$, i.e., the flux coming from the cosmic microwave background. Similar rewriting of the nighttime surface flux condition in terms of the solutions developed reveals

$$(1 - A) f_{\text{bg}} = \varepsilon \sigma T_{\text{sn}}^4 + \mathcal{I} \sqrt{\omega} \times \left[ - \frac{T_{\text{sd}} - T_{\text{sn}}}{2} \Xi + \frac{\phi_d(T)}{2} (1 - \Upsilon) \right]. \quad (\text{C20})$$

Finally, we must establish the value of $T_{\text{noon}}$, which is determined by evaluating $T_d$ at $z = -h$, which, after some more algebraic rearranging, turns into the relationship

$$T_{\text{noon}} + \Theta(hk) \phi_d(T_{\text{noon}}) = \frac{1}{2} (T_{\text{sd}} + T_{\text{sn}} + (T_{\text{sd}} - T_{\text{sn}})\Upsilon), \quad (\text{C21})$$

with

$$\Theta(hk) \equiv \frac{(3 + \xi)\cosh 2hk - 4\cosh hk + 4\sqrt{\xi}(\sinh 2hk - \sinhk) + 1 - \xi}{4\sqrt{\xi}\cosh 2hk + 2(1 + \xi)\sinhk}. \quad (\text{C22})$$
The various limiting forms for Θ are

$$ 
\begin{aligned}
\lim_{hk \to 0} \Theta &= hk + \mathcal{O}(hk)^2, \\
\lim_{hk \to \infty} \Theta &= \frac{3 + \sqrt{\xi}}{2 + 2\sqrt{\xi}} + \mathcal{O}(e^{-hk}).
\end{aligned}
\quad (C23) $$

We also observe that the prefactor

$$ \frac{\sqrt{L}P_0}{2\pi I \sqrt{\omega}} \quad (C24) $$

has the units of a temperature. We can use the relationship in Equation (C21) to rewrite the two surface boundary conditions free of the Arrhenius form $\phi_d$. As such, Equation (C15) now appears as

$$ (1 - A)f_{\odot} = \varepsilon \sigma T_{sd}^4 + \frac{1}{2} I \sqrt{\omega} \left[ (T_{sd} - T_{sn})\Xi + T_{sd} \frac{(1 + \Upsilon)^2}{2\Theta} + T_{sn} \frac{1 - \Upsilon^2}{2\Theta} - \frac{1 + \Upsilon}{\Theta} T_{noon} \right], \quad (C25) $$

and similarly, Equation (C20) as rewritten reads

$$ (1 - A)f_{bg} = \varepsilon \sigma T_{sn}^4 + \frac{1}{2} I \sqrt{\omega} \left[ -(T_{sd} - T_{sn})\Xi + T_{sd} \frac{1 - \Upsilon^2}{2\Theta} + T_{sn} \frac{(1 - \Upsilon)^2}{2\Theta} - \frac{1 - \Upsilon}{\Theta} T_{noon} \right]. \quad (C26) $$

The complete solution of this system that determines $T_{sd}$, $T_{sn}$, $T_{noon}$ amounts to the simultaneous solution of the three nonlinear algebraic relationships found in Equations (C21), (C25), and (C26), where $\phi_d(T_{noon})$ is given in Equation (C13), together with the functions $\Xi$, $\Upsilon$, and $\Theta$ expressed in Equations (C16) and (C22).

We observe an important limiting form. When the layer depth approaches the surface, i.e., $h \to 0$, it follows that $K\partial_z T_d|_{z \to 0} \approx K\partial_z T_d|_{z \to -h^+}$. It also stands to reason that $T_{noon} \to T_{sd}$ in this limit as well. Then, based on the jump condition in Equation (C10), it follows that

$$ \lim_{h \to 0^+} K \partial_z T_d \bigg|_{z \to 0} \to K_b \partial_z T_d \bigg|_{z \to 0^-} + \frac{\sqrt{L}}{2\pi} P_{vap}(T_{sd}). \quad (C27) $$

Thus, the energy balance condition (Equation (C15)) is limited to the familiar form

$$ (1 - A)f_{\odot} = \varepsilon \sigma T_{sd}^4 + K_b \partial_z T_d|_{z \to 0} + \frac{\sqrt{L}}{2\pi} P_{vap}(T_{sd}), \quad (C28) $$

where we have now dropped the “±” superscripting. In other words, the last term on the right-hand side of the above relationship accounts for the direct loss of energy to space due to material sublimation. In this case, there is no need to solve Equation (C21) for $T_{noon}$.

ORCID iDs

Abhinav S. Jindal https://orcid.org/0000-0002-1158-3446

Samuel P. D. Birch https://orcid.org/0000-0002-4578-1694

Alexander G. Hayes https://orcid.org/0000-0001-6397-2630

Orkan M. Umurhan https://orcid.org/0000-0001-5372-4254

Raphael Marschall https://orcid.org/0000-0002-0362-0403

Jason M. Soderblom https://orcid.org/0000-0003-3715-6407

Jean-Baptiste Vincent https://orcid.org/0000-0001-6575-3079

Dennis Bodewits https://orcid.org/0000-0002-2668-7248

References

Agarwal, J., Della Corte, V., Feldman, P. D., et al. 2017, MNRAS, 469, s606

Auger, A.-T., Groussin, O., Jorda, L., et al. 2015, A&A, 583, A35

Belton, M. J., & Melosh, J. 2009, Icar, 200, 280

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

Birch, S. P. D., Tang, Y., Hayes, A. G., et al. 2017, MNRAS, 469, S50

Blum, J., Gundlach, B., Mühle, S., & Trigo-Rodriguez, J. 2014, Icar, 235, 156

Britt, D., Boice, D., Buratti, B., et al. 2004, Icar, 167, 45

Combi, M., Shou, Y., Fougere, N., et al. 2020, Icar, 335, 113421

Davidsson, B. J., Birch, S., Blake, G. A., et al. 2021, Icar, 354, 114004

De Sanctis, M. C., Capaccioni, F., Ciarniello, M., et al. 2015, Natur, 525, 500

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

Festou, M. C., Keller, H. U., & Weaver, H. A. (ed.) 2004, Comets II (Tucson, AZ: Univ. Arizona Press)

Fornasier, S., Mottola, S., Keller, H. U., et al. 2016, Sci, 354, 1566

Fougere, N., Altwegg, K., Berthelier, J. J., et al. 2016, A&A, 588, A134

Fray, N., & Schmitt, B. 2009, P&SS, 57, 2053

Gasc, S., Altwegg, K., Balsiger, H., et al. 2017, MNRAS, 469, S108

Groussin, O., Attree, N., Brouet, Y., et al. 2019, SSRv, 215, 29

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

Hobley, D. E. J., Moore, J. M., Howard, A. D., & Umurhan, O. M. 2018, NatGe, 11, 901

Hu, X., Shi, X., Sierks, H., et al. 2017, A&A, 604, A114

Keller, H. U., Barbieri, C., Lamy, P., et al. 2007, SSRv, 128, 433

Keller, H. U., Mottola, S., Hviid, S. F., et al. 2017, MNRAS, 469, S357

Knollenberg, J., Lin, Z. Y., Hviid, S. F., et al. 2016, A&A, 596, A89

Kramer, T., & Läuter, M. 2019, A&A, 630, A4

Lai, I.-L., Ip, W.-H., Su, C.-C., et al. 2017, MNRAS, 462, S533

Lisse, C. M., Gladstone, G. R., Young, L. A., et al. 2022, PSJ, 3, 112

Macher, W., Kömle, N., Skorov, Y., et al. 2019, A&A, 630, A12

Marschall, R., Markkanen, J., Gerig, S.-B., et al. 2020, FrP, 8, 227

Marschall, R., Mottola, S., Su, C. C., et al. 2017, A&A, 605, A112

Marschall, R., Rezac, L., Kappel, D., et al. 2019, Icar, 328, 104

Marschall, R., Su, C. C., Liao, Y., et al. 2016, A&A, 589, A90

Montgomery, D. R., & Dietrich, W. E. 1994, WRR, 30, 1153

Mottola, S., Arnold, G., Grothues, H.-G., et al. 2015, Sci, 349, aab0232

Oklay, N., Vincent, J. B., Fornasier, S., et al. 2016, A&A, 586, A80

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

Pajola, M., Lucchetti, A., Fulle, M., et al. 2017b, MNRAS, 469, S636

Preusker, F., Scholten, F., Matz, K.-D., et al. 2015, A&A, 583, A33

Preusker, F., Scholten, F., Matz, K.-D., et al. 2017, A&A, 607, L1

Pätzold, M., Andert, T. P., Hahn, M., et al. 2018, MNRAS, 483, 2337

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

Squyres, S. W., Nakamura-Messenger, K., Mitchell, D. F., et al. 2018, LPSC, 49, 1332

Steckloff, J. K., Lisse, C. M., Safrit, T. K., et al. 2021, Icar, 356, 113998

Steckloff, J. K., & Samarasinha, N. H. 2018, Icar, 312, 172

Su, C. C. 2013, PhD thesis, National Chiao Tung Univ., Taiwan

Sunshine, J. M., Thomas, N., El-Maarry, M. R., & Farnham, T. L. 2016, JGRE, 121, 2194

Tang, Y., Birch, S. P. D., Hayes, A. G., et al. 2019, A&A, 630, A10

Thomas, N., El Maarry, M., Theologou, P., et al. 2018, P&SS, 164, 19

Thomas, N., Davidsson, B., El-Maarry, M. R., et al. 2015, A&A, 583, A17

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

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

Umurhan, O. M., Grundy, W. M., Bird, M. K., et al. 2022, PSJ, 3, 110

Veverka, J., Klaasen, K., A'Hearn, M., et al. 2013, Icar, 222, 424

Vincent, J. B. 2018, LPSC, 49, 1281

Wu, J.-S., & Lian, Y. 2003, CF, 32, 1133

Wu, J.-S., & Tseng, K.-C. 2005, IJNME, 63, 37

Wu, J.-S., Tseng, K.-C., & Wu, F.-Y. 2004, CoPhC, 162, 166