© The Authors 2025

1 Introduction

Comets are among the most primitive remnants in our Solar System and they are considered key objects for studying its early stages.. However, even the comets that are currently the subject of observations are not genuinely pristine and they have been, instead, continuously altered from their origins due to collisional processes (e.g., Stern 1995), temporal or periodic outbursts (e.g., Hartmann 1993), or the loss of volatiles (e.g., Prialnik et al. 2004). Since these evolutionary processes leave traces on the cometary surfaces, it is important to explore the geological features of comets.

As of December 2024, six periodic comets have been explored by space projects (1P/Halley by Giotto; Keller et al. 1987; 9P/Tempel 1 by Deep Impact; A’Hearn et al. 2005 and Stardust-NexT; Thomas et al. 2013a; 19P/Borrelly by Deep Space 1; Soderblom et al. 2002; 67P/Churyumov-Gerasimenko by Rosetta; Thomas et al. 2015b; 81P/Wild 2 by Stardust; Brownlee et al. 2004; and 103P/Hartley 2 by EPOXI; Thomas et al. 2013b; hereafter, these comets are referred to as 1P, 9P, 19P, 67P, 81P, and 103P). From images taken by onboard cameras, the comet nuclei showed remarkable diversity in shape. This observational evidence may imply that their origins and modification processes were complex and multiphase, and that they had occurred in various directions (Cheng et al. 2013). The highresolution images also revealed diverse geological substructures on their surfaces, such as quasicircular depressions (Brownlee et al. 2004; Belton et al. 2013; Vincent et al. 2015a), terraces (Basilevsky & Keller 2006), mesas (Bruck Syal et al. 2013), and pinnacles (Basilevsky & Keller 2006; Basilevsky et al. 2017). These various substructures are thought to have distinct origins linked to the physical properties, heterogeneity, or seasonal evolution of their parent bodies (for a comprehensive review, see El-Maarry et al. 2019).

Quasicircular depressions are prevalent substructures found in all resolved cometary surfaces. Distinct from ordinary impact craters with a bowl-like morphology (Melosh 1980), cometary depressions commonly show flat bottoms and steep walls without rims, implying that they were heavily modified from the primitive state (Basilevsky & Keller 2006; Cheng et al. 2013). As a result, there is little consensus regarding their initial formation, although several hypotheses have been proposed. For example, Brownlee et al. (2004) first identified depressions on the 81P surface and suggested that they might have resulted from direct impact cratering. However, subsequent investigations have consistently indicated that the commonly observed morphology of the cometary depressions, featuring flat floors and steep walls, cannot be attributed solely to the impact process (Basilevsky & Keller 2006; Holsapple & Housen 2007; Vincent et al. 2015b). Belton et al. (2013) proposed an alternative idea, namely, “mini-outburst” events that were observed during the Stardust-NexT encounter of 9P could account for the formation of approximately 96% of the observed pits. After Rosetta discovered similar features on 67P, Vincent et al. (2015a) stated that sinkhole collapse could have caused the formation of some of the large and deep circular depressions. This necessitated pre-existing subsurface cavities of similar depths (d ∼ 200 m) and a numerical approach was also used to investigate this scenario. Considering the thermal history of 67P, Mousis et al. (2015) conducted a numerical simulation of the sublimation of crystalline ices. These authors concluded that such large cavities could only form if a dust mantle was not present on the surface, leading to the assertion that current illumination conditions alone cannot fully account for the large depressions observed on the comets. Mousis et al. (2015) argued alternatively that other thermal processes, such as clathrate destabilization or crystallization of amorphous ices, can explain such cavities of the observed depressions. In this paper, we use the geological term “depressions” to describe the excavated quasicircular features seen on comets, which have also been referred to as “craters” (e.g., Vincent et al. 2015b) or “pits” (e.g., Benseguane et al. 2022) in a number of works in the literature.

Although the question of the formation process remains open, recent investigations by Rosetta demonstrated a shared evolutionary trajectory of depressions on 67P, influenced by cometary activity. The major result of these studies has been that sublimation-driven erosion continually obliterates depressions by preferentially expanding their lateral dimensions through cliff retreatment, rather than by deepening them (Marschall et al. 2016; Attree et al. 2018; Benseguane et al. 2022; Guilbert-Lepoutre et al. 2023). Episodic cliff collapses from sudden outbursts or jets are also attributed to this expansion (Vincent et al. 2016a,b; Attree et al. 2018).

From a broader perspective, this tendency implies a common surface evolutionary sequence that predominantly smooth comets throughout sublimation activity. For example, Vincent et al. (2017) compared the general trend of surface topography for four Jupiter-family comets (JFCs; 9P, 67P, 81P, and 103P) and found that large cliff is predominantly found on the younger surface (67P and 81P). Similar evolutionary sequences have been proposed by ground-based observation, such as the correlation of surface ages with the phase function and albedo (Kokotanekova et al. 2017, 2018) or the frequency of minioutburst (Kelley et al. 2021).

Recent investigations have revealed the clear evolutionary path of 67P depressions under the influence of sublimation activity: the depressions are enlarged and eventually smoothed. However, little is known about the evolution of surface depressions under the influence of sublimation activity for comets other than 67P. If the hypothesis of “depression enlargement” holds for other JFCs beyond 67P, it is possible to broaden our investigation and revisit other JFCs that have already been explored. Among them, 81P particularly draws our attention because nearly the entire explored side of the nucleus was covered with quasicircular depressions (Brownlee et al. 2004). This abnormally high depression density is not fully explained by a single formation process.

Specifically, the numerical approaches undertaken by Guilbert-Lepoutre et al. (2023) determined the erosion rate of depressions on three JFCs (9P, 81P, and 103P), suggesting that water sublimation-driven erosion plays a role in depression enlargement and contributing to approximately a few percent of the overall excavated features. Nevertheless, there is still a lack of observational constraints to support this numerical approach. Guilbert-Lepoutre et al. (2023) estimated the erosion rates of JFCs based on the surface properties constrained by the Rosetta measurement for 67P (Benseguane et al. 2022). However, the surface active fraction of 81P (~10%; Sekanina 2003) is expected to be an order of magnitude greater than that of 67P (~1%; Marschall et al. 2016), which may imply ten times faster surface erosion. Indeed, the typical depression sizes on 81P were nearly an order of magnitude larger than those on other JFCs, exceeding the empirical saturation limit of impact crater densities (Basilevsky & Keller 2006). This paper attempts to determine whether this diversity in depression is related to the past activity of 81P. Notably, the 81P depressions had similar traits to those of other JFCs in terms of morphology and power slope of the size distribution, implying common formation mechanisms among these depressions (Thomas et al. 2013a; Ip et al. 2016).

81P is a typical JFC (a = 3.45 au, e = 0.537, i = 3.24°, and P = 6.4 yr) first discovered in the 1978 perihelion approach (Sekanina & Yeomans 1985). It is thought to have been injected into the current JFC orbit by a very close encounter with Jupiter (~0.0061 au) in 1974. This event dramatically altered the perihelion distance from q = 5 au to q = 1.69 au, dragging the body inside the water activity region. Since then, 81P has remained in its current orbit for five periods of perihelion (1978, 1984, 1990, 1997, and 2003) until the Stardust encounter in 2004. Although there were slight changes in the perihelion distance between 1984 and 1990 (q = 1.7 au → q = 1.6 au) due to a moderate Jupiter encounter in 1988, the orbit of 81P during this period was relatively stable without significant variations in orbit or activity level (Sekanina 2003).

The primary objective of our study is to constrain the potential “excavated volume” of the 81P depressions during its 30-year JFC phase (1974–2004) under the observational consideration of past activities. We studied the water and dust environment when 81P was inside the water snowline (r_h ≲ 3 au) based on previous ground-based observations. Given the inherent limitations stemming from incomplete in situ data on 81P activities, we addressed this gap through several models of water and dust with appropriate assumptions. We then compared the results with the observed 81P depression sizes and volumes. This comparison enables us to constrain the upper-limit contribution of 81P activity to the overall degree of depression enlargement.

In Sect. 2, we describe the methods and results to estimate the volumes and sizes of 81P depressions. Then, we present the water (Sect. 3) and dust (Sect. 4) environment of 81P to estimate the total excavated volume during its 30-year orbit and we compare the volume with that of the depressions. We discuss the validity, reliability, and novelty of the results in Sect. 5. We summarize this paper in Sect. 6.

2 Volume of depressions

This section describes the methods and results for the depression volume calculation. First, we describe the method used to estimate depression volumes using Stardust images (Sects. 2.1–2.2). Then, we report our results of the depression volume estimates (Sect. 2.3).

2.1 Shape model

Images taken by the Stardust onboard camera Navcam (Newburn et al. 2003) revealed approximately half of the 81P surface by capturing solar phase angle coverage ranging from −72° to + 103°. While these images revealed a predominant presence of quasicircular depressions, the shady region on the opposite side remains unexplored (Brownlee et al. 2004). In this study, we use the terms “explored” and “unexplored” by Stardust to distinguish the two regions. The overall three-dimensional (3D) configuration of 81P closely conforms to a triaxial ellipsoid (Duxbury et al. 2004), with semiaxes measuring 2.607 km × 2.002 km × 1.350 km (Kirk et al. 2005). The highest resolution achievable from the Navcam images was 14 m pixel⁻¹ (Brownlee et al. 2004). Farnham et al. (2005) constructed a shape model based on these observed images. The model comprises a total of 17 513 facets, and the unexplored side is artificially extrapolated with the fit ellipsoid (Fig. 1a).

Fig. 1 Schematic process to calculate depression volumes in this study. (a) Shape model of 81P used in this study. The unexplored side is extrapolated to the fit ellipsoid. (b) Elevation map of 81P. (c) Three defined regions of the explored side: left foot & right foot (orange), Rahe (red), and Hemenway (green). (d) Depressions extracted by the edge detection technique (same as the outermost contours in Fig. 2).

2.2 Derivation of the depression volume

In this subsection, we outline the methodology for estimating the total volume of the 81P depressions. The term “volume” may carry ambiguity because these depressions are not entirely closed systems but are connected to neighboring topographies We operationally define the volume as the surface “degree of the cavity,” delineated by the deviation from the smooth zero-level potential. This definition aligns with the conventional framework employed to express the volume of impact craters (Melosh 1980) Importantly, these depressions are presumed to have been formed through excavation processes, enabling us to gauge the extent of surface excavation on 81P through this metric.

We determined the relative elevation of the facets by measuring the shortest distance between the point of each facet and the fit ellipsoid (Fig. 1b). We then used an edge detection technique to detect depressions and delineate borders. The Canny edge detector, a widely employed 2D edge detection operator, was utilized for this purpose (Canny 1986). In principle, this detector identifies sharp boundaries in 2D images marked by rapid elevation changes. To apply the detector to the 81P shape model, we utilized 2D projected elevation maps of the original shape model To minimize the distortion resulting from the 2D projection of the 3D ellipsoid, we divided the explored side into three distinc regions: left foot & right foot, Rahe, and Hemenway (Fig. 1c).

The detector not only delineated the closed rims of the depressions but also highlighted various characteristic linear and curved features, potentially arising from the surface roughness of the comet’s nucleus (Brownlee et al. 2004). To match the detector output with only the depression borders, we visually determined that the negative elevation contours were well aligned with the detector output. This process yielded “reference” zero-level elevations for these three regions, which were −50 m for the left foot & right foot, −80 m for Rahe, and −40 m for Hemenway. Considering the facets beneath these contour levels as the depressions (Fig. 1d), we computed the volume of the j-th depression as follows: $$V_j={\displaystyle \sum _i{d}_i}{A}_i\cos {\theta }_i,$$(1)

where d_i is the relative depth of the i-th facet from the reference contour level, A_i is the surface area of the facet, and θ_i is the angle between the normal direction of the facet and the direction from the facet to the shortest point of the reference ellipsoid. The associated error of depths was considered by choosing the different points (three nodes and centers) on the facets, potentially containing a few meters of uncertainty for each facet (Farnham et al. 2005).

We also estimated the dimensions (depth and diameter) of the depressions. We defined the depth of the j-th depression (d_j) as the distance between the deepest facet in the cavity and the reference contour level. Defining the depression diameter (D) was more complicated due to the irregular depression morphology. We derived two types of diameters, D_min and D_max , which were defined as the diameters of the maximum circle enclosed by the depression borders and the minimum circle enclosing the depression borders, respectively. These definitions complement the previous definition of connecting two selected points on the opposite side of the borders (Kirk et al. 2005).

2.3 Results

Using the edge detection technique for the 81P shape model, we detected a total of 15 depressions on the explored side. This detection is listed in Table 1 and shown in Fig. 2. Averages in depression volume were calculated by using the center of the facets, while the associated errors were calculated by using the three nodes of the facets. Because some depressions partially overlapped, we derived the total volumes for the overlapping depressions. In these cases, we distinguished each depression based on visual inspection.

We measured the depth (d) and diameters (D_min and D_max) of each depression. The largest depressions in terms of diameter are located in the right foot (#1) and Rahe (#9). These two depressions have D > 1000 m. Due to the irregularity of the surface, some depressions exhibit a factor of two difference between D_min and D_max (see depressions #2 and #5 in Table 1). The average depth-to-diameter ratios are d/D_min ~ 0.2 and d/D_max ~ 0.12, which are in line with previous estimates (d/D ~ 0.21; Kirk et al. 2005).

We derived the total volume of these 15 depressions as ${0.291}_{-0.048}^{+0.025} {\text{km}}^3$. The three largest depressions (right foot, left foot, and Rahe; Brownlee et al. 2004) occupy most of this total volume. The left foot and Rahe are each likely to consist of two depressions. Five depressions (Depressions #1, #3, #4, #9, and #10) encompass −90% of the total depression volume. The right foot and Rahe regions, in particular, encompass almost ~80% of the total depression volume and contain the two largest depressions (depressions #1 and #9). Adopting the 81P nucleus bulk density (ρ_bulk = 600 kg m⁻³ ; Davidsson & Gutiérrez 2006), the total volume of these depressions corresponds to a mass of $M_{\text{depression }}={1.74}_{-0.29}^{+0.15}\times {10}^{11} \text{kg}$.

Fig. 2 Depressions were detected in three regions: Left foot & right foot, Rahe, and Hemenway. The black solid lines show the results of the edge detection technique. The red contours indicate the relative elevation labeled with depth (m). The outermost contours of each region represent the borders of the detected depressions. The x- and y-axes of each region are projected distances (km). In each plot, the horizontal and vertical axes indicate arbitrary Cartesian coordinates on the triaxial ellipsoid surface.

3 Water environments

In this section, we describe the water environments of 81P. We present our model to calculate the water sublimation in Sect. 3.1 and our results in Sect. 3.2.

3.1 Water sublimation model

In this subsection, we outline our 81P water sublimation model, which combines numerical simulations and observational constraints. The model enables us to acquire the full-phase water production curve throughout the orbit. The comet water production trend is inherently first-order dependent on the illuminated cross-section, and this dependency can result in asymmetric production before and after perihelion passage, particularly when the rotational pole is oblique (Marshall et al. 2019).

For this reason, we estimated the water sublimation environment, taking into account its known orbit, shape model (Sect. 2.1), and spin axis orientation (I, Φ) = (55°, 150°), where I is the obliquity and Φ the longitude of the pole orientation (Sekanina et al. 2004). Since the orbit remained stable for 30 years after the encounter with Jupiter in 1974, we did not consider the small variation in the orbit (Sekanina 2003). Additionally, we did not consider the nongravitational effect of the 81P spin state due to its small change (Gutiérrez & Davidsson 2007).

We simulated the water sublimation environment following the frameworks established in previous works (e.g., Steckloff et al. 2015; Marshall et al. 2019). We postulated a cometary surface where sublimation activity occurs instantaneously (i.e., zero thermal inertia) and is fully homogeneous across the nucleus surface (i.e., pure water ice surface). Under these assumptions, the energy balance among solar radiation, thermal emission, and water ice sublimation on the nucleus surface is given by $$\frac{S_{\odot }\left(1-A_{\text{H}}\right)}{r_{\text{h}}^2}\cos {\theta }_i=\epsilon \sigma T_i^4+Z\left(T_i\right)L_{\text{ice}},$$(2)

where the subscript i indicates the i-th facet of the shape model, S_⊙ = 1367 W m⁻² is the solar constant equivalent to the solar flux at 1 au, A_H = 0.01 is the bolometric Bond albedo (Lamy et al. 2007), r_h is the heliocentric distance in units of au (practically dimensionless in the equation), θ_i is the solar incidence angle of the facet, ε = 0.9 is the emissivity, σ is the Stefan- Boltzmann constant, T_i is the sublimation temperature of the facet, and L_ice = 2.6 × 10⁶ J kg⁻¹ is the latent heat of water ice.

In Eq. (2), Z(T) (kg s⁻¹ m⁻²) is the temperature-dependent water sublimation rate. It is given by Hertz–Knudsen equation (Steiner & Koemle 1991): $$Z(T)=\frac{p_{\text{evp }}(T)}{\sqrt{2\pi m_{{\text{H}}_2\text{O}}{k}_{\text{B}}T}},$$(3)

where p_evp(T) is the evaporation gas pressure obtained by solving the Clausius–Clapeyron equations (Keller & Thomas 1989), $m_{{\text{H}}_2\text{O}}$ is the water molecular mass, and k_B is the Stefan– Boltzmann constant. Combining Eqs. (2) and (3) enables us to calculate the equilibrium sublimation temperature and the corresponding sublimation rate on the facets at a given solar incident flux. Then, we estimated the total water production rate ($Q_{{\text{H}}_2\text{O}}$ ; molecule s⁻¹ ) on the entire nucleus surface at a given heliocentric distance by the following equation: $$m_{{\text{H}}_2\text{O}}{Q}_{{\text{H}}_2\text{O}}\left(r_{\text{h}}\right)={\displaystyle \sum _i{Z}_i}\left(r_{\text{h}}\right)A_i.$$(4)

We computed the full-phase water production rates (Eq. (4)) for 200 evenly divided true anomalies. The water sublimation rate for every orbital step was averaged daily by rotating the subsolar point over one comet day at 25 uniform intervals. For each computation, we considered only the solar-facing facets (θ_i < 90° ). We did not consider the change in the heliocentric distance during one comet rotation (∆r_h ~ 10⁻⁴ au per comet rotation). We ignored the effects of shadowing and self-heating by other facets.

We confined our simulations to r_h < 3 au not only because the water sublimation drops significantly around r_h = 3 au but also because the simplified approach described above loses validity beyond this distance. Beyond this line (r_h > 3 au), the role of heat conduction becomes nonnegligible (Steckloff et al. 2015), and the primary sublimation species shift from water to other supervolatiles, such as CO or CO₂ (Davidsson et al. 2022). Neither of these factors were incorporated into our formulations. While recognizing the potential limitations of these simplifications, we justified our approach as suitable for our specific objective, considering that water sublimation-driven erosion is the predominant mechanism for surface mass loss (Kossacki & Czechowski 2018).

Nevertheless, the areal water production rates are highly overestimated under the assumptions above (zero thermal inertia and a pure water ice surface) because a large fraction of the surface could be covered with an inert dust mantle. In pursuit of a more realistic representation, we modified the second assumption by considering that the surface comprises well-mixed pure water ice and refractories. This involved introducing a scaling parameter f into the right-hand side of Eq. (4), that is, $$m_{{\text{H}}_2\text{O}}{Q}_{{\text{H}}_2\text{O}}\left(r_{\text{h}}\right)=f{\displaystyle \sum _i{Z}_i}\left(r_{\text{h}}\right)A_i(0\le f\le 1),$$(5)

where f is often referred to as an “effective active fraction” and is interpreted as the fraction of pure water ice on the surface; it is used to scale the idealized free water sublimation rate to the real rate (Keller & Thomas 1989). The value of f does not convey a physical meaning for characterizing real cometary surfaces, as the expected sublimation process is more intricate (e.g., involving sublimation from subsurface ice).

We constrained the possible range of f by considering the observed water production rate data and their standard errors. We used the water production rates during the two perihelion passages in 1997 and 2010, when extensive observations were conducted due to favorable observation conditions (Fink et al. 1999). These data sets are summarized in Table A.1. Since single-night water production measurements often have large errors, we binned the data points at 50-day intervals and used the average values within each bin. Referring to the variance-defined weights, we calculated the weighted mean water production rate in each interval as $$Q_{50}=\frac{{\displaystyle {\sum }_i^N{w}_i}{Q}_i}{{\displaystyle {\sum }_i^N{w}_i}},$$(6)

where N is the number of water production data points within each bin, Q_i is the individual measurement of the water production rate within the interval, $w_i=1/{\sigma }_i^2$ is the weight defined by the variance, and σ_i is the standard error of the water production rate from the individual measurements. Here, the subscript H₂O is omitted for simplicity. Similarly, the weighted standard error in each interval is calculated as $${\sigma }_{50}=\sqrt{\frac{{\displaystyle {\sum }_i^N{w}_i}{\left(Q_i-Q_{50}\right)}^2}{\frac{N}{N-1}{\displaystyle {\sum }_i^N{w}_i}}}.$$(7)

The water production rates calculated using Eqs. (6) and (7) are included in Table A.1. We set the lower and upper limits of the effective active fraction ( f_min and f_max ) according to the values enclosing all the binned data and their standard errors (see Fig. 3).

Specifically, the effective active fraction, f , is a parameter that varies with the heliocentric distance and with the nucleus coordinates if the surface is heterogeneous. In general, f tends to decrease at a larger heliocentric distance (Marschall et al. 2020) and shows large surface heterogeneity (Attree et al. 2019). We did not rigorously include this effect in our study.

Fig. 3 Seasonal water production rate of 81P scaled with observational data. The red, green, blue, magenta, and cyan points show the observational data from previous studies, and the black-filled circles are the average data within the 50-day bin (Table A.1). The gray and blue-filled regions represent the possible ranges of the water production rates from the entire surface (gray) and from only the explored side (blue), respectively, which is scaled by the average data. The dotted line shows the dust escape rates adopted from Pozuelos et al. (2014a). The black arrow indicates the day on which Stardust encountered 81P.

3.2 Results

We simulated seasonal water production rates under simplified conditions: a homogeneous surface with zero thermal inertia. Here, we considered the orbital motion, the shape model, and the measured spin orientation. We averaged the individual observational data with 50-day intervals and scaled the simulated water production curve into the binned data. Here, the range of the scaling factor ( f ) was determined by encompassing the whole binned data and their errors. This criterion produces the same order of water activity trends with observational data, while mitigating the large scatters of individual observations.

Figure 3 (gray-filled region) presents the resultant seasonal water production rate. It reveals asymmetric water production with respect to the perihelion, peaking approximately 40 days before the perihelion passage. This computational result is in line with the observed peak of the 81P water production rate (40–60 days before the perihelion; de Val-Borro et al. 2010). The overall trend of the water production roughly follows the dust escape rate (dotted line), which was reported by Pozuelos et al. (2014a). The peak of dust escape rate also occurs 40 days before the perihelion.

We find that f is in the range of 0.08 to 0.16 (the lower and upper ends of filled regions in Fig. 3 correspond to f = 0.08 and 0.16, respectively). We derive the peak water production rates for the minimum and maximum cases as $Q_{{\text{H}}_2\text{O}}\sim 5.6\times {10}^{27}$ molecule s⁻¹ and $Q_{{\text{H}}_2\text{O}}\sim 11.5\times {10}^{27}$ molecule s⁻¹, respectively, and the total water loss during one comet orbit as 4.6− 9.5 x 10⁹ kg. When we limit the calculation to the explored side (blue-filled region), we find that the surface lost 1.3−2.7 x 10⁹ kg of water during one comet orbit (29% of the total water loss).

The orbit-integrated water production in the southern hemisphere is two to three times greater than in the northern hemisphere (Fig. 4). The tilted pole orientation to the orbital plane can explain this latitudinal asymmetry. The illuminated area and subsequent water production on the inbound trajectory are twice as large as those on the outbound trajectory due to the tilted pole (de Val-Borro et al. 2010). During the period around the perihelion passage, the subsolar latitude gradually moves from the southern to the northern side.

As a result of asymmetric water production, the explored side (mostly the northern hemisphere) contributes to a total activity of ~30% during one cometary orbit around the Sun, while the same area consists of ~40% of the entire surface (22 of 56 km²). The Stardust spacecraft flew by 81P approximately 116 days after its perihelion passage (black arrow in Fig. 3), making it possible to observe the northern hemisphere. Our calculation result for the water production rate suggests that this comet’s activity had already been reduced by a factor of two with respect to its peak and that it was observed in the less active northern hemisphere (opposite the active southern hemisphere).

The explored side in Fig. 4 shows a heterogeneous water sublimation among the facets, potentially depending on the slopes. This may accelerate the lateral expansion of the depressions rather than deepening the bottoms, as observed in 67P (Attree et al. 2018; Benseguane et al. 2022; Guilbert-Lepoutre et al. 2023). We further discuss the heterogeneous erosion of the explored side in Sect. 5.8.

Fig. 4 Orbit-integrated water production rates of 81P. (left) Surface distribution of water production. The left and right figures preferentially show the explored northern and the unexplored southern side, respectively. (right) Latitudinal-dependent water production. The filled circles and the error bars show average sublimation rates and standard deviations of the facets inside the latitudinal bin. The dotted line shows the fraction of the explored side with respect to the entire surface inside the latitudinal bin.

4 Dust environments

In this section, we describe the dust environments of 81P. We present our dust ejection model in Sect. 4.1 and the results in Sect. 4.2. Finally, we summarize the results of depressions, water sublimation, and dust ejection of 81P and derive the depression excavation rate in Sect. 4.3.

4.1 Dust ejection model

In this subsection, we outline our model for estimating the amount of dust ejected from the nucleus surface. Dust particles can be ejected from the nucleus surface primarily due to the gas lifting force generated by the expansion of sublimating volatiles (Grun et al. 1989). In a general examination of the interplay among the lifting force, nucleus gravity, and dust particle properties, we identify three distinct cases of the dust ejection process:

• (i)

  Dust particles that are ejected from the surface and eventually escape from the comet Hill sphere (Ṁ_escape),

• (ii)

  Dust particles that are ejected from the surface but subsequently fall back to the surface (Ṁ_fallback),

• (iii)

  Dust particles that remain on the surface without lifting.

Here, we define Ṁ_escape (kg s⁻¹) as the dust escape rate and Ṁ_fallback (kg s⁻¹) as the dust fallback rate. Although it is not critical to consider the dust particles in case (iii) in this section, they should be taken into account to fully describe the nucleus surface environment.

Given the water production environment that determines the lifting force, the differences among these three cases depend solely on the dust particle mass (or size) if we assume homogeneous and spherical dust particles (m_d = 4πa³ρ_d/3, where a and ρ_d denote the dust radius and the bulk mass density, respectively). Here, we set ρ_d = 600 kg m⁻³ for dust particles from 81P (Niimi et al. 2012).

Assuming that the size distribution of the dust production follows a single-power law, the number of dust particles ejected in a unit of time within the size range of a_min < a < a_max is expressed as $$n_{\text{ej}}(a)=k_{\text{p}}{\left(\frac{a}{a_0}\right)}^{-b_{\text{d}}},\left(a_{\min }<a<a_{\max }\right),$$(8)

where k_p(s⁻¹) is a normalization parameter, b_d is the power index of the size distribution, and a₀ = 1 cm is the reference dust radius. We fixed the minimum radius of the dust particles to a_min = 1 µm (Pozuelos et al. 2014a). amax represents the “maximum liftable radius,” which is determined by the balance between the gas drag force and the nucleus gravity at the surface. Assuming adiabatic spherical expanding gas, an analytical solution for a_max can be expressed as (Zakharov et al. 2018) $$a_{\max }\left(r_{\text{h}}\right)=\frac{3m_{{\text{H}}_2\text{O}}\text{Z}\left(r_{\text{h}}\right)R_{\text{N}}^2{C}_{\text{D}}\sqrt{k_{\text{B}}{T}^{*}\gamma /m_{{\text{H}}_2\text{O}}}}{8{\rho }_{\text{d}}{M}_{\text{N}}G},$$(9)

where γ = 4/3 is the specific heat ratio of water, G is the gravitational constant, R_N = 2 km is the effective nucleus radius, and M_N = 2.3 x 10¹³ kg is the nucleus mass (ρ_bulk = 600 kg m⁻³; Davidsson & Gutiérrez 2006). T^* is the gas temperature on the sonic surface, which is related to the sublimation temperature (T in Eq. (2)) by T^* = 0.8T (Lukianov & Khanlarov 2000). C_D is the gas drag coefficient; for a spherical dust particle, it can be expressed as follows (Bird 1994): $$C_{\text{D}}(s)=\frac{2s^2+1}{s^3\sqrt{\pi }}\exp \left(-s^2\right)+\frac{4s^4+4s^2-1}{2s^4}\mathrm{erf}(\text{s})+\frac{2\sqrt{\pi }}{3 \text{s}},$$(10)

where s is the speed ratio between the expanding gas and the gas thermal velocity. The adiabatic solution yields $s=\sqrt{\gamma /2}$ at the surface. Notably, a_max depends on the local water sublimation rate (Z) and therefore depends on the heliocentric distance along the comet orbit.

The sum of Ṁ_escape and Ṁ_fallback is calculated using n_ej as $${\dot{M}}_{\text{escape }}+{\dot{M}}_{\text{fallback }}={\displaystyle {\int }_{a_{\min }}^{a_{\max }}{m}_{\text{d}}}(a)\left|\frac{d{n}_{\text{ej}}(a)}{da}\right|da.$$(11)

Eq. (11) is obtained based on mass flux conservation, indicating that dust either escapes or falls back once ejected. We did not account for the disintegration of particles, which could alter the dust size distribution over time (Davidsson 2024). We separated Eq. (11) as follows: $$\begin{array}{l}{\dot{M}}_{\text{escape }}={\displaystyle {\int }_{a_{\min }}^{a_{\text{escmax }}}{m}_{\text{d}}}(a)\left|\frac{d{n}_{\text{ej}}(a)}{da}\right|da\\ =\frac{4\pi }{3}\frac{b_{\text{d}}}{3-b_{\text{d}}}{\rho }_{\text{d}}{k}_{\text{p}}\left[{\left(\frac{a_{\text{escmax}}}{a_0}\right)}^{-b_{\text{d}}}{a}_{\text{escmax}}^3-{\left(\frac{a_{\min }}{a_0}\right)}^{-b_{\text{d}}}{a}_{{\min }^3}\right],\end{array}$$(12)

and $$\begin{array}{l}{\dot{M}}_{\text{fallback }}={\displaystyle {\int }_{a_{\text{escmax }}}^{a_{\max }}{m}_{\text{d}}}(a)\left|\frac{d{n}_{\text{ej}}(a)}{da}\right|da\\ =\frac{4\pi }{3}\frac{b_{\text{d}}}{3-b_{\text{d}}}{\rho }_{\text{d}}{k}_{\text{p}}\left[{\left(\frac{a_{\max }}{a_0}\right)}^{-b_{\text{d}}}{a}_{{\max }^3}-{\left(\frac{a_{\text{escmax}}}{a_0}\right)}^{-b_{\text{d}}}{a}_{\text{escmax}}^3\right],\end{array}$$(13)

where a_escmax represents the “maximum escape radius,” which distinguishes between escaping and falling back dust particles. Accordingly, a_escmax is a key parameter for determining Ṁ_escape and Ṁ_fallback . Here, Ṁ_escape is a quantity that can be measured from ground-based observations of dust comas and dust tails. In fact, Pozuelos et al. (2014a) conducted intensive observations of the dust coma and tail of 81P and derived the time-dependent Ṁ_escape . In contrast, Ṁ_fallback cannot be determined directly from ground-based observations but can be predicted using a possible range of a_escmax , as described below.

To determine a_escmax , we utilized the 81P dust ejection model constrained by Pozuelos et al. (2014a). During the 2010 apparition, the authors conducted a comprehensive study on the evolution of the optical coma brightness and the dust tail morphology of 81P across a broad heliocentric distance range both before (r_h < 2.7 au) and after (r_h < 3.2 au) the perihelion passage. Pozuelos et al. (2014a) used ~300 ground-based observation data points obtained through 2010 apparition, which were mostly observed by amateur astronomical association Cometas-Obs (to find the complete dataset, refer to their Fig. 4). The authors fitted the data with Monte Carlo dust tail modeling (Moreno 2009) to analyze the dust tail morphology in terms of the initial dust ejection parameters (see Pozuelos et al. 2014b, for additional information on the usage of the model). This approach enabled them to constrain the seasonal variations in the dust escape rate (Ṁ_escape), the dust velocity along the normal direction (v_d) of a₀ = 1 cm particles at R = 20R_N , where R is the cometocentric distance, and the power index of the dust size distribution (b_d). Fig. 5 presents a reproduction of their dust modeling, adopted from Fig. 3 in Pozuelos et al. (2014a). We excluded two minor transient outbursts that occurred in 2010 (Pozuelos et al. 2014a) to characterize the general trend of the activity. The effects of these cometary outbursts on the depressions are discussed separately in Sect. 5.10.1.

Following the works of Pozuelos et al. (2014a), we considered dust velocity (v_d) of a₀ = 1 cm particles at R = 20R_N. At this distance, the gas drag force on the dust particles is deemed negligible, and v_d reaches the terminal velocity (Zakharov et al. 2018; Pätzold et al. 2019). We further assumed that v_d can be determined by the dust size. Accordingly, the size-dependent v_d is expressed as $$v_{\text{d}}\left(a,R=20R_{\text{N}}\right)=v_{\text{d}}\left(a=a_0,R=20R_{\text{N}}\right){\left(\frac{a}{a_0}\right)}^{-{\alpha }_{\text{v}}},$$(14)

where α_v is the power index of the size-dependent velocity distribution at a cometocentric distance of R = 20R_N (i.e., when the dust reaches the terminal velocity). Wallis (1982) derived α_v ~ 0.5 via an analytical approach. In-situ measurements obtained using GIADA onboard Rosetta yielded α_v = 0.96 ± 0.54 at 5R_N ≲ R ≲ 15R_N from 67P (Della Corte et al. 2015). Given the large uncertainty associated with these in situ measurements, we considered α_v = 1.0 ± 0.5.

Calculating Eq. (14) enables us to determine a_escmax by comparing v_d with the escape velocity $v_{\text{esc}}(R)=\sqrt{2G{M}_{\text{N}}/R}\text{\hspace{0.17em}\hspace{0.17em}at\hspace{0.17em}\hspace{0.17em}}R=20R_{\text{N}}$. In Fig. 6, we illustrate examples of how the a_escmax values are determined in two cases at small and large heliocentric distances. In the left column, the escape velocity at R = 20R_N (horizontal dotted line) is fixed to v_esc = 0.23 m s⁻¹, and the a_max (vertical dotted line) is determined by the water production rate (Eq. (9)). Then, a_escmax is determined based on the intersection points between v_d and v_esc or v_d and a_max. As a result, the particles escape if the dust velocity is smaller than the escape velocity (v_d < v_esc) and fall back if the dust velocity is larger (v_d > v_esc). As seen from the right column of Fig. 6, α_v is the key parameter for determining a_escmax . Considering the full range of α_v (i.e., 0.5 ≤ α_v ≤ 1.5), there is a factor of three difference in a_escmax at the small heliocentric distance. This tendency decreases when the comet is far from the Sun because weak solar insolation at a large heliocentric distance narrows the gap between v_d(a₀ = 1 cm) and v_esc.

We did not incorporate additional forces, such as the solar radiative force or solar gravity, which are deemed negligible in this region in the short time period under consideration here (Marschall et al. 2020). Therefore, we can simplify the analysis of dust motion in this regime by focusing solely on the nucleus gravity.

Fig. 5 Heliocentric dependencies of (left) the dust ejection rate, (center) the dust velocity of a0 = 1 cm-sized particles at a cometocentric distance of R = 20RN, and (right) the power index of the dust size distribution. The plots are adopted from Fig. 3 in Pozuelos et al. (2014a), while we excluded the data relevant to the two small outburst events observed in the 2010 apparition.

Fig. 6 Dust velocity (vd) at R = 20RN and maximum escape radius (aescmax) in our dust ejection model. The upper and lower panels represent examples at small and large heliocentric distances, respectively. The left column shows the size-dependent dust velocity at a cometocentric distance of R = 20RN . The horizontal dotted lines indicate the escape velocity (vesc) at the cometocentric distance. The vertical dotted lines indicate the maximum liftable radius (amax) at the heliocentric distance. The intersection points with horizontal dotted and colored lines correspond to the maximum escape radius at the given αv values. The right column shows the determined maximum escape radius as a function of αv. The horizontal dotted lines represent the maximum liftable radius, showing the upper limit of the maximum escape radius at given heliocentric distances.

4.2 Results

As expressed in Eq. (13), the amount of fallback debris inherently depends on the two edges of the fallback particle sizes, the maximum liftable size (a_max) and the maximum escape size (a_escmax), under the assumption of a single power law distribution. We believe that a_max can be determined solely by the balance of the local water sublimation environment and the nucleus surface gravity. We then calculated a_escmax by assuming the power index of the particle velocity distribution (α_v) at R = 20R_N and using the observational constraints for the dust velocity along the normal direction (v_d) of a₀ = 1 cm particles at R = 20R_N obtained by Pozuelos et al. (2014a).

Figure 7 shows the a_escmax values with respect to the distance from the Sun (with the perihelion distance subtracted, r_h-r_h,peri). As predicted from Fig. 6, a_escmax increases as α_v decreases. As an upper limit of a_escmax , we plot the maximum liftable radius (a_max) in Fig. 7. a_escmax and a_max are derived by independent methods: a_escmax from the dust tail model in Pozuelos et al. (2014a) and the dust velocity distribution and a_max from the water sublimation model. In principle, the condition a_escmax = a_max implies that all the particles ejected from the composite surface will escape the comet, so there is no fallback material (see the interval of the integral in Eq. (13)).

Figure 8 shows the dust fallback rate (Ṁ_fallback) with respect to the heliocentric distance and different α_v values. For comparison, we plot the dust escape rate (Ṁ_escape) from Pozuelos et al. (2014a). Due to the negative correlation between a_escmax and αv , Ṁ_fallback increases as α_v increases. Here, Ṁ_fallback approaches zero for small α_v (≲ 1). Otherwise, for a large α_v (≳ 1), Ṁ_fallback with respect to the heliocentric distance is proportional to Ṁ_escape. The normalization coefficient (k_p) in Eqs. (12) and (13) analytically explains this proportionality. A larger Ṁ_escape value leads to a larger kp , which simultaneously drives up Ṁ_escape and Ṁ_fallback . More intuitively, an increase in dust production (Ṁ_escape + Ṁ_fallback) increases both Ṁ_escape and Ṁ_fallback.

We integrated Ṁ_fallback through one orbital revolution around the Sun. According to the result, $M_{\text{fallback}}={4.06}_{-4.05}^{+14.84}\times {10}^9$ kg, where we obtain the average at α_v = 1.0, upper limit at α_v = 1.5, and lower limit at α_v = 0.5. Importantly, Ṁ_fallback approaches zero for a small α_v. We compare M_fallback with the dust escape during one comet orbit (M_escape) and obtain $M_{\text{fallback}}/M_{\text{escape}}={0.36}_{-0.36}^{+1.32}.$

In principle, v_esc and a_max also depend on nucleus mass (M_N), which has large uncertainty due to the poor constraint of the 81P bulk density (ρ_bulk = 200–800 kg m⁻³; Davidsson & Gutiérrez 2006; Szutowicz et al. 2008; Sosa & Fernández 2009). Here we fixed ρ_bulk = 600 kg m⁻³ and corresponding nucleus mass M_N = 2.3 × 10¹³ kg. The effect of nucleus bulk density on the fallback mass is discussed in Sect. 5.7.

Fig. 7 Maximum escape radius (aescmax) at different heliocentric distances and αv values. The dotted line represents the maximum liftable radius and indicates the upper limit of the maximum escape radius; thus, the colored lines cannot exceed this limit.

Fig. 8 Dust fallback rate (Ṁfallback) with respect to the distance from the Sun and the αv values. The black solid line represents the dust mass escape rate adopted from Pozuelos et al. (2014a).

4.3 Excavation rate of 81P depressions

In this subsection, we summarize our results and compare the depression mass (M_depression ; Sect. 2) with the water sublimation ($Q_{{\text{H}}_2\text{O}}$; Sect. 3) and dust ejection (M_escape + M_fallback; Sect. 4) to estimate the “excavation rate” of the depressions. We express the mass excavation rate (kg s⁻¹) of the nucleus surface as follows: $${\dot{M}}_{\text{excavate }}=m_{{\text{H}}_2\text{O}}{Q}_{{\text{H}}_2\text{O}}+{\dot{M}}_{\text{escape }}+{\dot{M}}_{\text{fallback }}.$$(15)

If we only consider the explored side of 81P, we can revise the above equation to $${\dot{M}}_{\text{excavate }}=m_{{\text{H}}_2\text{O}}{\left[Q_{{\text{H}}_2\text{O}}\right]}_{\text{explored }}+\left({\dot{M}}_{\text{escape }}+{\dot{M}}_{\text{fallback }}\right)\frac{{\left[Q_{{\text{H}}_2\text{O}}\right]}_{\text{explored }}}{Q_{{\text{H}}_2\text{O}}},$$(16)

where ${\left[Q_{{\text{H}}_2\text{O}}\right]}_{\text{explored}}$ is the water production rate on the explored side (the blue-filled region in Fig. 3). The second term of Eq. (16) suggests that water gas lifts up dust in proportion to its sublimation activity.

Figure 9 shows the results for Ṁ_excavate from the explored side (only three cases, namely, α_v = 0.5,1.0, 1.5, are presented in the plot). In the case of M_N = 2.3×10¹³ kg (ρ_bulk = 600 kg m⁻³; Davidsson & Gutiérrez 2006), we constrain the total excavated mass from the explored side during one integrated orbit by M_excavate = 3.7−9.1 × 10⁹ kg, where the lower limit occurs when α_v = 0.5 and f = 0.08 and the upper limit occurs when α_v = 1.5 and f = 0.16 (i.e., dependent on the dust fallback and water production). The estimated maximum rate of excavated mass is Ṁ_excavate ~ 250 kg s⁻¹, occurring ~20 days after the perihelion. Approximately 70% of the surface material is excavated after the perihelion due to asymmetric water production.

Comparing M_depression with M_excavation from the explored side enables us to determine the excavation rate of the 81P depressions. The result yields M_excavation/M_depression = 0.02–0.06, implying 2–6% of the total depressions mass (or volume) on the explored side are excavated every comet orbit around the Sun. From the Jupiter encounter in 1974 to the Stardust encounter in 2004, 81P revolved five times, maintaining its constant JFC orbit and the activity level (Sekanina 2003). Accordingly, we conclude that the total excavated mass of the 81P depressions due to the JFC-phase activity is up to 30% of the total depression mass.

We summarize our findings in Table 2. Throughout the estimation, we fixed the nucleus bulk density with ρ_bulk = 600 kg m⁻³. The effect of bulk density on our results is discussed in Sect. 5.7. In addition, the estimate is an ideal upper limit of the excavation rate of the depressions, assuming all the water sublimation occurs inside the depressions. We further discuss the heterogeneity of the local excavation in Sect. 5.8, based on a consideration of the 81P topography and the solar illumination conditions.

Fig. 9 Mass excavation rate (Ṁexcavate) from the explored side of 81P. Three cases of αv are presented. The filled range of each plot is attributed to the uncertainty in the water production rate.

5 Discussion

In this section, we discuss the validity, reliability, and novelty of these results by comparing them with previous research. We first discuss the novelty of this study in Sect. 5.1. We discuss the consistency of 81P activity level during 30 years in Sect. 5.2, the reliability and interpretability of our results on the number of depressions in Sect. 5.3, the water environment in Sect. 5.4, and fallback debris in Sect. 5.5. Then, we describe our derivation of the dust-to-ice mass ratio range of 81P under a few assumptions in Sect. 5.6. We discuss the effect of nucleus bulk density and the surface heterogeneity on our results in Sects. 5.7 and 5.8, respectively. Finally, we discuss the implications of our results for the evolution and origin of 81P depressions in Sects. 5.9–5.10.

Throughout Sects. 5.3–5.6, we discuss our results on the depression count, water environment, fallback debris, and dust-to-ice mass ratio of 81P by comparing them with those for 67P. Both comets share similar orbital traits and histories, having comparable semimajor axes (a = 3.45 au for 81P and a = 3.46 au for 67P), perihelion distances (q = 1.69 au for 81P and q = 1.21 au for 67P), obliquity of the rotation pole axes (I = 55° for 81P and I = 52° for 81P), and years of injection into the current JFC orbit (in 1974 for 81P and in 1959 for 67P). This comparative analysis enhances the understanding of the characteristics of 81P by leveraging insights gained from the well-studied comet 67P.

5.1 Novelty of this study

We conducted unique research by integrating observational evidence from 81P (Pozuelos et al. 2014a) into the analysis of depression enlargement, addressing the potential limitations of previous research. In fact, Guilbert-Lepoutre et al. (2023) examined the depression enlargement of various JFCs. However, their approach relies on thermal parameters determined in 67P research; thus, it may overlook differences in the surface conditions of 81P and 67P. In contrast, our study incorporates observational constraints from Pozuelos et al. (2014a) when analyzing the enlargement of 81P depressions, representing the first attempt at treating comets other than 67P. Although our estimates of the dust environment of 81P (including escaping and fallback dust debris) may not be as accurate as in situ measurements of 67P (Pätzold et al. 2019; Marschall et al. 2020), we provide reasonable results considering the broad range of parameters (f, α_v), as well as strong observational evidence from Pozuelos et al. (2014a), who observed 81P across a wide range of orbital locations.

5.2 Long-term activity of 81P

Throughout the procedures described in Sects. 3–4, we used ground-based observational data on water and dust production rates to estimate the seasonal activity of 81P. For the water production rates (Sect. 3), we relied on data from the 1997 and 2010 apparitions, when the visibility of the comet was favorable (Table A.1). Additionally, as described in Section 4.1, our dust ejection model was developed using observational data from Pozuelos et al. (2014a), who collected data across a wide range of heliocentric distances during the 2010 apparition.

We assumed that 81P’s activity remained consistent across its orbits after it was injected into its current orbit in 1974. From 1974 to 2010, data on water and dust production are available only for the 1997 and 2010 apparitions. However, visual magnitude data for earlier apparitions were analyzed by Sekanina (2003). According to his study, the brightness of the comet was consistent during the first two apparitions (1978 and 1984) but dimmed by about one magnitude in the third apparition (1990). This decrease coincided with an increase in the perihelion distance, from q = 1.49 au to q = 1.58 au, after a moderate encounter with Jupiter in 1987. However, Sekanina (2003) noted that this change in brightness may not be directly related to the change in the perihelion distance, as the brightness returned to previous levels during the fourth apparition (1997), even though the perihelion distance remained unchanged. The comet was not observed during the fifth apparition (2003), but its activity was likely similar to previous levels, as the activity observed in the fourth (1997) and sixth (2010) apparitions was consistent (de Val-Borro et al. 2010).

In summary, the overall activity level of 81P remained stable for the six apparitions from 1974 to 2010, despite a temporary brightness decrease in the 1990 apparition. Thus, we confidently extended our water and dust activity models based on the 1997 and 2010 data to earlier 81P orbits.

5.3 Reliability of our depression count

We detected a total of 15 depressions. This number is less than the 23 depressions detected by Kirk et al. (2005) in the same region. It is not clear why there is a discrepancy. We utilized the edge detection algorithm, while Kirk et al. (2005) detected depressions by visual inspection of a false color elevation map, which made it difficult for us to trace their results. However, we ensured the inclusion of all the large depressions (D ≳ 500 m) in our counting. Kirk et al. (2005) estimated about 13 depressions with D > 500 m, while we detected 13 depressions with D ≳ 400 m (Table 1). On the other hand, small depressions with D ≲ 400 m may have led to the discrepancy in counts between our study (2 depressions) and Kirk et al. (2005) (about 10 depressions). Nevertheless, the volume of each depression with D < 400 m corresponds to <1% of our estimates (< 0.001 km³/0.290 km³). Therefore, we expect that uncounted small depressions contribute only to 10% of the total volume of depression. As seen in Sect. 2.3, the five largest depressions encompass the majority (90%) of the total volume of depression. Therefore, the discrepancy in the depression counting between this study and Kirk et al. (2005) would not influence our main result, as it has minor contributions to the total volume of depression.

Although the total number of detected depressions differs between Kirk et al. (2005) and this study, this study shows good agreement with Kirk et al. (2005) regarding the depression diameters. They derived the dimensions of the depressions, with depths ranging from 50 m to 500 m and diameters ranging from 250 m to 2500 m, and obtained a depth-to-diameter ratio of d/D ~ 0.2. We derive d/D ~ 0.12-0.2 from the depression diameter range of 100-1000 m. We also find that the d/D ratio in our research is consistent with those of 9P and 67P, most of which had ratios of d/D ~ 0.1 (9P, Thomas et al. 2013a) and d/D = 0.1–0.2 (67P, Ip et al. 2016), respectively. Despite the small number of depression samples in our research, interestingly, the d/D ratio is also similar to those of impact craters observed on other airless bodies, such as lunar craters (d/D = 0.8–0.17; Wu et al. 2022) and Ryugu (d/D = 0.09 ± 0.02; Noguchi et al. 2021). This similarity in the d/D ratios may imply that we cannot rule out the impact origins of depressions on JFCs, although they have a distinct flat bottom and steep cliff morphology. The morphological trait suggests that the JFC depressions may have undergone significant modification from their original shapes (Basilevsky & Keller 2006; Cheng et al. 2013). Further discussion on the impact origin of these depressions is provided in Sect. 5.10.2.

5.4 Asymmetric water production

Our water sublimation model successfully reproduces the observed asymmetry in water production for the perihelion. For example, the combined observational data (depicted as black dots in Fig. 4) yield $Q_{{\text{H}}_2\text{O}}=(9.3\pm 2.0)\times {10}^{27}$ molecule s⁻¹ 25 days before perihelion and $Q_{{\text{H}}_2\text{O}}=(5.6\pm 0.5)\times {10}^{27}$ molecule s⁻¹ 25 days after perihelion. Our model closely matches these values to the accuracy of the observation, predicting $Q_{{\text{H}}_2\text{O}}=8.4\times {10}^{27}$ molecule s⁻¹ and $Q_{{\text{H}}_2\text{O}}=6.9\times {10}^{27}$ molecule s⁻¹ on the same days, assuming an average effective active fraction (f = 0.12). While our model exhibits slightly smaller variations, it broadly aligns with the observational data, considering the associated error bars. The observed asymmetry has often been attributed to the seasonal effect of two confined jets positioned near the north pole (latitude >75°) and the southern hemisphere (Farnham & Schleicher 2005). However, the exclusion of this effect in our model suggests that the asymmetry arises primarily from the shape of the nucleus, coupled with its high obliquity (I = 55°). This result is consistent with earlier discussions regarding the asymmetric water production trend of JFCs (including 81P) (Marshall et al. 2019).

The integrated water production over one cometary orbit indicates that the water production in the southern hemisphere is two to three times greater than that in the northern hemisphere. A critical aspect of this finding is that the explored side represents a relatively less eroded region. Similar disparities were reported in 67P with a high obliquity (I = 52°), where the southern side showed a factor of four greater erosion rate than did the northern side due to intensive solar insolation in the southern hemisphere near the perihelion (Keller et al. 2015; Thomas et al. 2015a).

Our model estimates the effective active fraction within the range of 0.08 < f < 0.16. This value is consistent with that of a previous study (f = 0.081; Jewitt 2021). By comparison, 67P exhibits f ~ 0.01 under similar assumptions (Marschall et al. 2016) and f ~ 0.06 under more realistic surface conditions (Keller et al. 2015). Therefore, the surface of 81P appears more active than those of typical JFCs (Sekanina 2003; Pozuelos et al. 2014a).

5.5 Amount of fallback debris

Our investigation provides a constraint on the total fallback debris mass of M_fallback < 2.0 × 10¹⁰ kg during one comet orbit (Table 3). Assuming a homogeneous distribution of the debris over the entire surface (~56 km²), we infer that the debris would produce a dust layer of ≲50 cm thickness per orbit around the Sun. Here, we posit that the bulk density of this dust layer equals the nucleus bulk density (ρ_bulk = 600 kg m⁻³; Davidsson & Gutiérrez 2006). If we further assume that the debris accumulates solely on regions with low gravitational slope (<30°), the layer thickness is estimated to double (≲1 m). Here, we postulate that the low gravitational slope region constitutes 50% of the entire surface based on the shape model of the explored side. The slope criterion considers the angle of repose of the debris (Marschall et al. 2017).

This is the first attempt to constrain the amount of fallback debris other than for 67P (Pätzold et al. 2019; Marschall et al. 2020; Davidsson et al. 2021). Despite the simplicity of our modeling compared to previous work, our estimate for the fallback debris aligns well with that of 67P. We constrain the ratio between fallback and escape to $M_{\text{fallback}}/\left(M_{\text{escape}}+m_{{\text{H}}_2\text{O}}{Q}_{{\text{H}}_2\text{O}}\right)\sim 1.2$ or less, which is similar to the 1.8–4.3 range found for 67P (Pätzold et al. 2019). The estimated thickness of the fallback debris layer also corresponds to the findings for 67P: Davidsson et al. (2021) estimated decimeter-scale dust deposits (0.1–1 m) per comet orbit depending on the location in the nucleus, and Marschall et al. (2020) suggested a debris layer of ≲40 cm per comet orbit over the smooth regions on 67P.

However, the amount of fallback debris is strongly affected by α_v and we cannot decide which value between 0.5 < α_v < 1.5 is suitable for 81P. To retrieve the reliability of our dust model, we calculated the fallback debris of 67P adopting our model. The result is presented in the last column of Table 3. Overall, 67P shows approximately five times larger fallback debris compared to 81P, at the same α_v. Comparing our results with previous research, a sophisticated dust ejection model of Marschall et al. (2020) estimated the amount of fallback debris of 67P during one comet orbit as ${6.8}_{-6.8}^{+11}\times {10}^8$ kg. Our result yields a similar amount of debris when α_v falls between 0.6 and 0.7. If this value can be adaptable to 81P, the estimated fallback debris is less than < 0.1 cm of deposits to the surface. However, due to the lack of observational constraints, we cannot decide which value is correct between this large uncertainty.

Our results suggest that 81P activity would add ≲5 m to the debris layer during the 30-year JFC phase. However, the actual distribution of the fallback debris may differ significantly from our simplified assumption. The spatial distribution of the debris may be heterogeneous, as observed on 67P, where much of the debris is transported from the southern hemisphere to the northern hemisphere during the polar night on the northern side near the perihelion (Keller et al. 2015; Thomas et al. 2015a; Lai et al. 2016). Local debris removal due to self-cleaning effects or intensive outbursts is also nonnegligible (Davidsson et al. 2021). Nevertheless, as pointed out by Davidsson et al. (2021), similar mass transport (from the active south to the quiescent north) is unlikely to occur on 81P. Admittedly, our estimation is less sophisticated than those of well-studied 67P models. Nonetheless, we suggest that the actual fallback layer on 81P likely reached a few meters.

Fig. 10 Dust-to-ice mass ratio with respect to the different heliocentric distances and αv values. The dust mass is the sum of the escaping and fallback masses.

5.6 Dust-to-ice mass ratio

This research provides insights into a crucial parameter, namely, the mass ratio between dust and water ice, which is a fundamental characteristic of icy bodies. Determining this ratio is challenging. The conventional approach involves simultaneously measuring the production rates of escaping dust (A fρ; A’Hearn et al. 1984; Fink & Rubin 2012) and water $\left(Q_{{\text{H}}_2\text{O}}\right)$. However, this approach likely introduces bias, as the ratio varies with the heliocentric distance, with comae becoming more dusty at large distances (A’Hearn et al. 1995). Similarly, the average dust-to-ice mass ratio along the integrated orbit is likely biased.

As an alternative approach, including the dust fallback rate in the calculation would mitigate the potential bias. Our estimation yields a dust-to-ice mass ratio $\left(\left(M_{\text{fallback}}+M_{\text{escape}}\right)/m_{{\text{H}}_2\text{O}}{Q}_{{\text{H}}_2\text{O}}\right)$ between 1 and 7. A similar approach suggests a similar value (<2.5) for 67P (Marschall et al. 2020). However, as shown in Fig. 10, this ratio depends on the orbital location (especially the heliocentric distance), yielding values smaller than unity at large heliocentric distances and values substantially larger than unity near the peak production level (i.e., near the perihelion).

We attribute the heliocentric distance dependency of the dust-to-ice mass ratio in Fig. 10 to two primary factors. First, dust ejection from the surface is contingent upon the gas sublimation environment, which determines the total dust amount (M_fallback + M_escape). The maximum liftable size (a_max) is a function of the local gas sublimation rate ($Q_{{\text{H}}_2\text{O}}$; Eq. (9)). At a large heliocentric distance, a_max would be small, resulting in most of the remaining dust not being lifted. Second, the fallback "dust" (Ṁ_fallback) consists of not only refractory but also volatile ice. Evidence of this was observed during the Rosetta mission, where smooth terrains consisting of fallback debris frequently showed ice sublimation activity (Hu et al. 2017). Furthermore, Davidsson et al. (2021) conducted a thermophysical simulation of ejected dust and concluded that decimeter-sized particles (the typical upper limit for fallback debris) retain >95% of water ice during the first 12 hours post-ejection, with even cm-sized particles preserving >40% of initial water ice. This could potentially introduce a bias in Fig. 10, overestimating the refractory dust mass, as sublimation from these ejected particles would still contribute to the total water production rate.

Moreover, our model does not rigorously account for the long-term evolution of the dust-to-ice mass ratio. Fallback dust deposits are known to suppress activity over time (Attree et al. 2018, 2023), leading to a gradual reduction in activity levels and an increase in the dust-to-ice mass ratio near the surface. This retroactive effect, which is absent from our model, could introduce further bias. The magnitude of this bias depends on the amount of fallback debris generated per comet orbit (see Sect. 5.5). If the dust layer deposits are less than <1 cm per orbit (α_v ≲ 1), the dust-to-ice mass ratio will evolve slowly, and the current values will remain close to the previous estimates. However, if the dust layer deposits exceed >10 cm per orbit (α_v ≳ 1), the ratio will increase rapidly, introducing substantial bias into current estimates.

In short, it would be difficult to derive the true surface dust-to-ice mass ratio without comprehensive in situ measurements, as were obtained by the Rosetta mission. Nevertheless, we suggest that the dust-to-ice mass ratio of 81P near the perihelion (i.e., the peak value in Fig. 10) could reflect a close approximation to the true ratio. This approach mitigates the biases caused by the factors mentioned above. First, a_max increases near the perihelion due to intensive gas sublimation, incorporating most of the surface dust components into the mass calculation. Our estimation indicates a maximum liftable dust size of a_max ∼ 50 cm near the perihelion. Second, the volatile loss of the lifted dust would be maximized because of the intensive solar radiation and the longer flight time within the coma due to the high ejection velocity. As a result, we suggest that the surface dust-to-ice mass ratio of 81P falls between 2 and 14 depending on the water production efficiency (f) and the velocity distribution of the ejected dust (α_v). This result implies a dusty surface environment of 81P, which is consistent with that of 67P (3–7.5; Choukroun et al. 2020).

5.7 Effect of bulk density

Throughout this study, we have assumed a bulk density for 81P of ρ_bulk = 600 kg m⁻³. However, estimates of the bulk density of 81P widely vary between different studies. Davidsson & Gutiérrez (2006) proposed a bulk density of ρ_bulk < 600–800 kg m⁻³, while Szutowicz et al. (2008) suggested ρ_bulk = 400 ± 200 kg m⁻³ , and Sosa & Fernández (2009) estimated ${\rho }_{\text{bulk}}={300}_{-300}^{+500}$ kg m⁻³. As a result, this introduces a large uncertainty not only in the estimation of depression mass (M_depression) but also in the total nucleus mass. This, in turn, affects key parameters in our dust ejection model, particularly in calculations of the maximum escape radius (a_escmax), the maximum liftable particle size (a_max), and the amount of fallback debris (M_fallback).

The effect of the bulk density on our findings is summarized in Table 4. The estimated bulk density of 81P from previous studies lies between ρ_bulk = 200–800 kg m⁻³, therefore, we considered four different values, ρ_bulk = 200, 400, 600, and 800 kg m⁻³ , and calculated the corresponding M_depression and M_fallback . It is clear that M_depression increases linearly with the bulk density. In contrast, M_fallback decreases as the bulk density increases due to the reduction in a_max , which limits the contribution of large dust particles to the total dust mass. The central value of M_fallback is determined when α_v = 1.0, with the lower and upper limits defined by α_v = 0.5 and 1.5, respectively. As noted in Sect. 4.2, no fallback material is present in our dust model when α_v = 0.5. In addition, in all cases, M_fallback is highly sensitive to the choice of α_v.

Variations in bulk density affect both M_depression and excavated mass (M_excavate) from the explored side. Consequently, a wide range of “depression excavation rates” (M_excavate/M_depression) is possible, ranging from 1% to 28% of the depressions excavated per comet orbit. If we exclude the cases of very low bulk density (ρ_bulk = 200 kg m⁻³), “typical” M_excavate/M_depression falls between 1% and 11% per comet orbit. Although we cannot completely rule out the case of low bulk density, we conclude that the total volume excavated from the 81P depressions during the JFC phase is less than ≲50% of the depression volume (considering that 81P completed five revolutions around the Sun for 30 years).

This estimate could be an ideal (i.e., oversimplified) upper limit of the depression excavation rates, assuming that all water sublimation occurs inside the depressions. We further discuss the excavated volume of the depression in Sect. 5.8, considering the heterogeneous activity of the nucleus.

Fig. 11 Water sublimation model from different activity models. (left & center) Orbit-integrated water sublimation of 81P’s explored side. The left shows a fully homogeneous model, and the center shows a depression-active model, respectively. (right) Time-dependent water production rates of two models.

5.8 Effect of activity location on the nucleus

In Sect. 4.3, we compare the depression volume with the surface excavated volume to estimate the depression excavation rate during the 30-year JFC phase. Our results constrain an upper limit for the depression excavation rate, with M_excavate/M_depression < 6% per orbit and <30% over the 30-year JFC phase. This notably assumes that all excavations occurred “inside” the depressions. Admittedly, our study does not rigorously define the exact location of the excavation on the nucleus surface, suggesting that not all sublimation-driven activity necessarily contributes to depression enlargement. Therefore, our estimate represents an upper limit for depression enlargement.

We assumed homogeneous activity across the surface, modeling water sublimation with a single scaling factor f (Eq. (5)). However, this assumption overlooks the potential for a heterogeneous activity across the nucleus, where different surface regions may have different f values, as found on the 67P’s surface (i.e., the f value features a regional dependency; Marschall et al. 2020; Attree et al. 2023). Identifying these variations is challenging for 81P. In Fig. 11, we present two extreme activity models: one assuming homogeneously active regions (our default assumption) and the other assuming only facets within the depressions are active (the "depression-active" model). In the latter case, the faces outside the depressions are assumed to have f = 0, and the f values inside the depressions are determined to fit the observational data. However, no significant differences are observed between the models. This highlights a limitation of our research because the available ground-based observational data are insufficient to precisely pinpoint the locations of activity on the nucleus.

Two primary factors contribute to this limitation. First, the depressions on 81P are evenly distributed on the explored side, resulting in minimal variation in overall activity trends. Second, the explored side of 81P contributes only a small fraction of total water production, while the unexplored southern hemisphere plays a more substantial role, especially near perihelion.

We further explored whether the degree of erosion could be different inside and outside the depressions. Two main factors likely drive heterogeneous erosion: solar energy input and gravitational slope. Solar energy input is the primary driver of heterogeneous activity (Kossacki & Jasiak 2019; Benseguane et al. 2022) and is influenced by the oblique rotation pole (Marshall et al. 2019) and local topography (Kossacki & Jasiak 2019; Benseguane et al. 2022). For 81P, these effects are pronounced due to its large obliquity (I = 55°; Sekanina et al. 2004) and rough topography (Brownlee et al. 2004; Vincent et al. 2017). Gravitational slope also plays a role in determining the thickness of local refractory dust deposits, which can hinder solar heat transfer and reduce local sublimation activity. As shown by Hu et al. (2017), even a 10 mm thick dust layer can reduce activity to 20% of that on clean surfaces (without dust deposits). Our calculations suggest that fallback debris could accumulate to thicknesses exceeding 1 m in low-sloped regions (see Sect. 5.5), potentially suppressing subsurface sublimation.

However, according to our calculation, there are little differences inside and outside the depressions, in terms of the total solar energy input and the fraction of high gravitational slope (>30°). The average solar energy input is similar for the facets inside (3.4 GJ m⁻² orbit⁻¹) and outside (3.3 GJ m⁻² orbit⁻¹) depressions. Moreover, the area fraction with a gravitational slope larger than >30° is the same for inside and outside (55%). In conclusion, the regions inside and outside the depressions do not exhibit significant differences in surface excavation, indicating nearly homogeneous erosion across the surface. Based on their surface area fraction, depressions are estimated to contribute approximately 40% of the total excavation. Accordingly, the estimated depression excavation rate is M_total/M_dep < 3% per orbit and less than 15% over the 30-year JFC phase.

5.9 Implications and potential caveats of depression enlargement

Our estimate of the excavated volume (in comparison with the depression volume) implies that the depressions on 81P barely expanded in diameter over its 30 years as a JFC. We calculate the excavation volume over 30 years, compare it with the depression diameters and volumes, and find that the diameters have increased by only approximately 1% of the original diameters (2–20 m). Our estimate is roughly consistent with the maximum erosion of the depressions of 81P determined in a previous study (28 m; Guilbert-Lepoutre et al. 2023). This consistency supports the results of this previous numerical study, although their approach (a numerical approach based on thermal evolution) was different from ours (mostly based on the observational evidence of dust and gas ejections).

However, we recognize that our constraint does not fully account for the potential enlargement of the depressions. For example, cliff collapse can enlarge the depression diameter without a significant volume loss by accumulating collapsed chunks on depression bottoms. Evidence of cliff collapse and a subsequent increase in the number of boulders below the cliff was observed for 67P (Groussin et al. 2015; Pajola et al. 2015; Davidsson 2024). Some observational evidence on 81P, such as the talus found in the left foot, suggests that cliff collapses might have occurred in past orbits (Brownlee et al. 2004). However, while such accumulation processes can increase the depression diameter, they do not result in volume excavation. Therefore, our analysis of depression “enlargement” pertains to depression volume but may not fully encompass changes in depression size due to these unknown enlargement mechanisms.

5.10 Implications for the origins of 81P depressions

Because sublimation-driven activity contributes less to depression enlargement according to our research, there should be other factors involved in this phenomenon. Here, we discuss two potential mechanisms, namely, cometary outbursts (Sect. 5.10.1) and impact cratering (Sect. 5.10.2), as alternative origins of depressions on 81P and other JFCs.

5.10.1 Outbursts

Outbursts are observed across various comets, displaying a wide range of magnitudes and frequencies. However, their contribution to the formation of depressions on 81P remains uncertain. Pozuelos et al. (2014a) observed two outbursts during the 2010 apparition, estimating the total mass ejected from these events as M_ej ~ 9.2 × 10⁸ kg M_ej ~ 3.0 × 10⁸ kg. Assuming these outbursts excavated surface mass (with p_bulk = 600 kg m⁻³; Davidsson & Gutiérrez 2006) and that all the excavated mass escaped from the nucleus, the corresponding excavated volume is approximately V ~ 0.001 km³ for each outburst. The total excavated volume from these two outbursts is thus less than 1% of the total depression volume on the explored side (0.288 km³). Obviously, the outbursts observed by Pozuelos et al. (2014a) in the 2010 apparition are not associated with the depressions observed during the Stardust encounter in 2004. Therefore, it is unlikely that the outbursts strongly contributed to the formation of depressions within the 30-year JFC orbit. Although we definitively exclude the possibility that 81P underwent large-scale outbursts after it became a JFC to create depressions on the surface, the possibility of outbursts at large heliocentric distances cannot be ruled out (Gronkowski & Wesołowski 2016).

Fig. 12 Cumulative depression size distributions of three JFCs (9P, 67P, and 81P) and one Kuiper Belt object Arrokoth. The data used were obtained from the following studies: 9P: Thomas et al. (2013a); 67P: Ip et al. (2016); 81P: Kirk et al. (2005); and Arrokoth: Spencer et al. (2020).

5.10.2 Impact cratering

Impact cratering is one of the most significant surface modification processes for Solar System bodies, including comets. Given the short dynamic lifetime of JFCs (< 10⁴ years; Levison & Duncan 1997), it is unlikely that impact cratering formed dense depressions during the JFC orbit (Vincent et al. 2015b). Similarly, the contribution of impact cratering during the Centaurs (i.e., pre-JFC) phase is negligible compared to the Kuiper Belt, given their short dynamic lifetime and low impact probability (Durda & Stern 2000).

Ip et al. (2016) investigated depressions on three JFCs (9P, 67P, and 81P) and find that they share a similar power-law slope in the depression size distribution, implying that the depressions commonly formed before they became JFCs. These slopes are also consistent with the impact craters on Arrokoth, the Kuiper Belt object (Fig. 12). This may imply that the depressions on the three JFCs formed during the long-term residence in their reservoir. However, direct comparisons of size distributions between JFCs and Arrokoth should be treated with caution, since JFCs likely originated from different subpopulations with Arrokoth (Nesvorný et al. 2017). Furthermore, the surface densities of the JFC depressions are one to two orders of magnitude higher than those of the Arrokoth depressions. To account for this excess, 10 to 100 times more frequent impacts must have occurred compared to the expected impact rates on Kuiper Belt objects (Dell’Oro et al. 2013; Abedin et al. 2021). It is not certain whether such a diversity of crater densities occurs in the Kuiper Belt region, while incorporating a similar size distribution for the different subpopulations. Further investigation is required to shed light on this question.

6 Conclusions and summary

We studied the water and dust environment of 81P (water escape, dust escape, and dust fallback) and estimated the excavation rates of depressions, comparing them with the depression volume. We conjecture that the depressions were preserved during the 30-year JFC phase. We summarize our main results below.

1. The total excavated mass (volume) on the explored side of 81P is 3.7–9.1 × 10⁹ kg (0.006-0.015 km³) during one comet orbit around the Sun. This constrains the maximum excavation rate of the depression (0.288 km³) to within 5% per one cometary orbit and less than 30% during the 30-year JFC phase;

2. We constrained the total amount of fallback debris to M_fallback < 2.0 × 10¹⁰ kg during one cometary orbit around the Sun. The expected thickness of the debris is ≲50 cm in one comet orbit, assuming a homogeneous distribution of the debris. This result is consistent with the literature findings for 67P;

3. We constrained the dust-to-ice mass ratio of the 81P surface to between 2 and 14, depending on the water production efficiency ( f ) and the velocity distribution of the ejected dust (α_v). This result is consistent with the findings for 67P.

Acknowledgements

This research was supported by a National Research Foundation of Korea (NRF) grant funded by the Korean government (MEST) (No. 2023R1A2C1006180). We also appreciate the anonymous reviewer for their careful reading and insightful comments.

Appendix A Additional materials

References

All Tables

All Figures

Fig. 1 Schematic process to calculate depression volumes in this study. (a) Shape model of 81P used in this study. The unexplored side is extrapolated to the fit ellipsoid. (b) Elevation map of 81P. (c) Three defined regions of the explored side: left foot & right foot (orange), Rahe (red), and Hemenway (green). (d) Depressions extracted by the edge detection technique (same as the outermost contours in Fig. 2).

In the text

Fig. 2 Depressions were detected in three regions: Left foot & right foot, Rahe, and Hemenway. The black solid lines show the results of the edge detection technique. The red contours indicate the relative elevation labeled with depth (m). The outermost contours of each region represent the borders of the detected depressions. The x- and y-axes of each region are projected distances (km). In each plot, the horizontal and vertical axes indicate arbitrary Cartesian coordinates on the triaxial ellipsoid surface.

In the text

Fig. 3 Seasonal water production rate of 81P scaled with observational data. The red, green, blue, magenta, and cyan points show the observational data from previous studies, and the black-filled circles are the average data within the 50-day bin (Table A.1). The gray and blue-filled regions represent the possible ranges of the water production rates from the entire surface (gray) and from only the explored side (blue), respectively, which is scaled by the average data. The dotted line shows the dust escape rates adopted from Pozuelos et al. (2014a). The black arrow indicates the day on which Stardust encountered 81P.

In the text

Fig. 4 Orbit-integrated water production rates of 81P. (left) Surface distribution of water production. The left and right figures preferentially show the explored northern and the unexplored southern side, respectively. (right) Latitudinal-dependent water production. The filled circles and the error bars show average sublimation rates and standard deviations of the facets inside the latitudinal bin. The dotted line shows the fraction of the explored side with respect to the entire surface inside the latitudinal bin.

In the text

	Fig. 5 Heliocentric dependencies of (left) the dust ejection rate, (center) the dust velocity of a0 = 1 cm-sized particles at a cometocentric distance of R = 20RN, and (right) the power index of the dust size distribution. The plots are adopted from Fig. 3 in Pozuelos et al. (2014a), while we excluded the data relevant to the two small outburst events observed in the 2010 apparition.
In the text

Fig. 6 Dust velocity (vd) at R = 20RN and maximum escape radius (aescmax) in our dust ejection model. The upper and lower panels represent examples at small and large heliocentric distances, respectively. The left column shows the size-dependent dust velocity at a cometocentric distance of R = 20RN . The horizontal dotted lines indicate the escape velocity (vesc) at the cometocentric distance. The vertical dotted lines indicate the maximum liftable radius (amax) at the heliocentric distance. The intersection points with horizontal dotted and colored lines correspond to the maximum escape radius at the given αv values. The right column shows the determined maximum escape radius as a function of αv. The horizontal dotted lines represent the maximum liftable radius, showing the upper limit of the maximum escape radius at given heliocentric distances.

In the text

	Fig. 7 Maximum escape radius (aescmax) at different heliocentric distances and αv values. The dotted line represents the maximum liftable radius and indicates the upper limit of the maximum escape radius; thus, the colored lines cannot exceed this limit.
In the text

	Fig. 8 Dust fallback rate (Ṁfallback) with respect to the distance from the Sun and the αv values. The black solid line represents the dust mass escape rate adopted from Pozuelos et al. (2014a).
In the text

	Fig. 9 Mass excavation rate (Ṁexcavate) from the explored side of 81P. Three cases of αv are presented. The filled range of each plot is attributed to the uncertainty in the water production rate.
In the text

	Fig. 10 Dust-to-ice mass ratio with respect to the different heliocentric distances and αv values. The dust mass is the sum of the escaping and fallback masses.
In the text

	Fig. 11 Water sublimation model from different activity models. (left & center) Orbit-integrated water sublimation of 81P’s explored side. The left shows a fully homogeneous model, and the center shows a depression-active model, respectively. (right) Time-dependent water production rates of two models.
In the text

	Fig. 12 Cumulative depression size distributions of three JFCs (9P, 67P, and 81P) and one Kuiper Belt object Arrokoth. The data used were obtained from the following studies: 9P: Thomas et al. (2013a); 67P: Ip et al. (2016); 81P: Kirk et al. (2005); and Arrokoth: Spencer et al. (2020).
In the text