© ESO, 2017

1. Introduction

The brief revisit to the comet 9P/Tempel 1 by the spacecraft Stardust-NExT in 2011 led to the discovery of various changes on the nucleus surface that had occurred over a period of almost six years following the exploration by the Stardust-NExT’s predecessor Deep Impact (A’Hearn et al. 2005; Thomas et al. 2013; Veverka et al. 2013). The receding scarps, several smoothed crater rims, and possible local variations in the contrast of surface albedo, would represent the most intriguing and, yet, the only direct observations of cometary surface modifications prior to the Rosetta mission.

Since the rendezvous of Rosetta with the comet 67P/Churyumov-Gerasimenko (hereafter 67P) in mid 2014, the Optical, Spectroscopic, and Infrared Remote Imaging System (OSIRIS) onboard the spacecraft had been observing the comet at unprecedented spatial and temporal resolutions (Keller et al. 2007; Sierks et al. 2015). The observations provided critical constraints in determining the physical properties of the bi-lobed, dusty nucleus (Sierks et al. 2015; Pätzold et al. 2016), and enabled detailed morphological characterization of the surface (Thomas et al. 2015b) as well as realistic thermal analysis on the activity of 67P (Keller et al. 2015). Parts of the nucleus have been recurrently observed in detail. Various surface changes on 67P had been noted before its perihelion passage. A prominent case occurred in the smooth terrains of Imhotep¹ from late May to early July, 2015, where multiple quasi-circular, fast-expanding features were observed that probably indicate migrating scarps due to continuous collapse or removal of the top multi-meter thick surface layer (Groussin et al. 2015). Moving escarpments of possibly comparable propagating pattern and speed were detected in Hapi (Davidsson et al., in prep.), which had been among the first visibly active areas of the nucleus since mid 2014 (Lara et al. 2015; Lin et al. 2015).

In the eyes of OSIRIS, hints of surface changes on 67P were abundant and multiplying as the comet approached perihelion. A number of “honeycomb” features, distinguished by ensembles of small depressions of a few decimeters in width, were detected in late March 2015 in the Ma’at region (Shi et al. 2016a). These features resided within the formerly smooth dust deposits, and likely indicate a change in macroscopic surface roughness over time. It is possible that the honeycombs had resulted from accumulated erosion in the dust deposits due to previous activities, for example, sublimation of water ice, in the course of a few months. Given the prevalence of dust deposits in the northern hemi-nucleus of 67P, this possibility stands as a strong motivation for expanding the scrutiny over the nucleus for other similar but often subtle indications of surface changes.

This work is an exposition of the various types of surface changes associated with the dust deposits observed on 67P before perihelion. The discussions are organized as follows. In Sect. 2, we present temporal comparisons of morphological observations indicative of changes over time. In Sect. 3, we perform a photometric analysis of observations to authenticate and quantify the surface changes as a probable consequence of erosion. In Sect. 4, a thermal model is applied to analyze the role of intensifying insolation in triggering or facilitating the dust erosion via water outgassing. We discuss the global distribution of the observed surface changes in Sect. 5; in particular, we argue that the erosion and restoration of the dust deposits are a morphological expression of seasonality on 67P. The implications of the surface changes studied here are discussed in Sect. 6 in connection with some key findings from other instruments onboard Rosetta. Finally, our conclusions are summarized in Sect. 7.

2. Observation

2.1. OSIRIS observation of honeycomb features during late March flyby

The Rosetta spacecraft performed a flyby around 67P on March 28, 2015. During the close encounter, the Narrow Angle Camera (NAC) of OSIRIS was able to collect a series of images with pixel resolutions² of 0.5 m or better within 30 km from the center of nucleus. The images revealed the presence of some distinct yet unfamiliar surface features in several morphological regions in the northern hemi-nucleus³. Most of these features can be easily distinguished as macroscopically rough patches ranging in size from dozens up to hundreds of meters (Fig. 1). Each appears to be densely indented or pitted in texture. Although the indentations lack the hexagonal shape and regular arrangement, these patches are reminiscent of the cellular structure of honeycombs on decimeter to meter scales (Shi et al. 2016a). The honeycomb texture here can be loosely generalized as a resolved surface roughness or topography in contrast with the surrounding smooth texture with random and more uniform roughness.

A list of identified honeycomb features and their distribution on the nucleus is provided in Table 1. All features are located between 15◦N and 50◦N in latitude. The distribution is somewhat uneven in longitude. The gaps from 30◦E to 60◦E and from 220◦E to 300◦E partly correspond to the Hapi region deep into the concavity between the two lobes (Sierks et al. 2015; Thomas et al. 2015b). The lack of features between 120◦E and 200◦E may be in part attributed to the incomplete coverage of high-resolution observations during the period.

The regions in which honeycombs were identified include Ma’at and Serqet, on the small lobe of the nucleus, and Ash, Babi, and Seth on the large lobe. All regions display dust deposits characteristic in the northern hemi-nucleus (Thomas et al. 2015b; El-Maarry et al. 2015). The dust cover probably resulted from the deposition of dust particles ejected from the nucleus that had not escaped from the nucleus and re-accumulated on the surface (Thomas et al. 2015a,b). The thickness of the dust cover varies. Thomas et al. (2015a) estimated a thickness of between about 1 and 5 m based on the depth-to-diameter ratio of a potential impact crater in the Ash region. A similar estimate is obtained by Mottola et al. (2015) for dust deposits at the designated landing site for Philae at Agilkia (La Forgia et al. 2015). Possibly, the honeycombs occurred sporadically within smooth terrains with extant dust deposits by the time of observations. Shi et al. (2016a) infer from the pattern of surface roughness as well as the typical dimension of individual depressions that the honeycombs may be akin to the pitted dusty terrains found primarily on the small lobe of 67P (Mottola et al. 2015; El-Maarry et al. 2015; La Forgia et al. 2015).

2.2. Comparison with earlier observations of local morphologies

The plethora of honeycombs observed across the nucleus prompted the question of whether they were fresh features indicating surface changes over time. The detection of changes by comparing images of a common location is restricted by the difference in resolution, viewing geometry, and illumination condition. The first detection and the last observation before the occurrence of a honeycomb specify a time frame in which the (potential) surface change has occurred. The observations of some particular honeycomb features in different regions are presented below.

2.2.1. Ma’at

The first clear sight of honeycombs was traced back to January 2015, when MAT01 and MAT02 were already distinguishable (Fig. 2). Both were half-shielded at the foot of a steep wall or cliff. The abrupt local topography often limits the visibility of the feature locations in the images. As a result, we have found few images from before this date that show the local morphology in as much detail. An unobstructed but more slanted view of the cliff floor where MAT02 was located suggests that the feature was imperceptible on September 12, 2014 (Fig. 2a). It cannot be ruled out that, in this case, higher incidence and phase angles of Fig. 2b may have accentuated the pattern of indentations in the honeycomb. However, we note that the widths of the indentations can be roughly likened to the size of a few adjoining boulder-like features from the debris. Many indentations were in full shadow, or unilluminated at the bottom. We infer that the depths of indentations should be somewhat comparable to the widths. Therefore, when the boulders are clearly discernible, the absence of indentations is unlikely to be a visual coincidence in Fig. 2a. Instead, it may indicate the presence of smooth dust deposits at the time of observation and, hence, a change in surface roughness over time. The positive relief of MAT02 indicates that the excavation of material was not restricted to within individual indentations but occurred around the honeycombed patch (note the valley formed between MAT02 and the cliff in Fig. 2).

In Fig. 3, we adapt and present a series of three images, spanning six months from the end of September 2014, which provide arguably the clearest temporal contrast in morphologies around honeycombs MAT03-07. The features occurred within the widespread dust deposits that masked the consolidated substrate revealing rugged, fractured terrains at some large outcrops (Fig. 3a). The changes may have occurred before mid-February 2015, when the feature locations already appeared heavily etched within a few shallow depressions (Fig. 3b). We considered the comparison between Fig. 3a and b as a further indication of the textural evolution of the dust deposits, which is revealed in the latter image with lower resolution. Figure 3c shows the possibly more developed texture at the end of March 2015. It is inferred from the observations that the increase of the surface roughness was not confined to the honeycomb features but occurred over the entirety of the dust-covered areas.

2.2.2. Ash/Seth

The features SET01 and 02 were situated within a shallow basin on a terrace looking down toward Hapi (Fig. 1b). The area is bordered by another higher terrace into Ash where ASH01 is located. The dust cover around the honeycombs was relatively smooth by mid-September 2014 (Fig. 4a). The surface displays corrugations on the scale of several meters that are probably either result of aeolian transport or reflect underlying subsurface topography (Thomas et al. 2015a,b). A closer but more slanting view of the same area on February 14, 2015 may hint at a change in the surface pattern (Fig. 4b). Many dim patches roughly several meters in size are visible in Fig. 4b. Some appear half-closed with shaded edges on the sunward side. It is inferred that they indicate shallow depressions with unilluminated walls shadowing the bottoms. Overall, the surface began to exhibit a pitted pattern in Fig. 4b, distinguished from the ruffled, undulating pattern in Fig. 4a. The surface may have subsequently evolved to the honeycombed texture distinct in Fig. 4c (cropped from Fig. 1b) in late March 2015. The rim of the basin occupied by SET01,02 had become sharper and more rugged at several places. A long scarp and a smaller V-shaped scarp were likely formed by displacement of materials from the inside, possibly due to slumping. The loss of dust deposits is also evident at some fresh concavities near the terrace wall toward Ash. A shallow depression emerged by SET02 after mid-February 2015 (Fig. 4b).

2.2.3. Babi

The features BAB01 and BAB02 both span roughly one hundred meters and had evolved from the formerly smooth, dusty plain fenced by extensive scarps along the boundary (Fig. D.1). A few images taken before mid-December 2014 show nearly identical surface textures, suggesting that the dust cover had probably not been altered down to the scale of 1 m (three-pixel resolution of the images) by then. The comparison of an image on December 30, 2014 with one on February 14, 2015 under similar illumination conditions revealed some localized (potential) surface changes around BAB02 (Figs. D.1c,d). For instance, a shallow concavity of a few meters in width had emerged next to an existing scarp (or, indistinguishably, a half-exposed boulder), in addition to a few other marginally resolved new depressions (Fig. D.1d). Nearby, a strip of the surface at the edge of BAB02 had coarsened perceptibly. However, there were no clear signs of significant evolution of the surface texture. The vastly different illumination conditions make it difficult to compare Figs. D.1c and d with earlier observations, such as Figs. D.1a and b. It is inferred that BAB01 and BAB02 probably did not take clear shape before mid-February 2015.

2.2.4. Serqet/Ma’at/Nut

Serqet is generally categorized as a region with consolidated materials and littered with boulder-like debris; nevertheless, it is morphologically diverse and exhibits distinct dust deposits toward its northern borders with Ma’at and Nut (El-Maarry et al. 2015). Three honeycombs were identified within the dusty portions of the region, showing the indented texture similar to other features (Fig. D.2b). More seemed to exist, for example, beyond the border with Nut, though they are not enumerated in Table 1. The proximity of these features to MAT12-13 may suggest that the surface texture of the dust deposits had been widely, if not evenly, altered across the tri-border area. It is noted, again, that the emergence of honeycombs was accompanied by numerous other potential changes in the vicinity (Fig. D.2a). For instance, the dense formation of small depressions appeared to have dispersed out from SEQ02 whose boundaries cannot be clearly distinguished. The topography nearby appeared more angular, which may reflect thinning of the dust deposits.

2.3. Other surface changes in the vicinity of honeycomb features

The emergence of honeycombs appeared to indicate an increase in the surface roughness at mid-latitudes in the northern hemi-nucleus. Numerous other surface changes have been noted that may indicate a wider expanse of surface change on 67P. We present some exemplary cases below with a crude classification of observations.

Exposure of abrupt topography.

By the end of March 2015, a few thin but sheer outcrops seemed to have been exposed from the exhausted dust deposits by the crumbled, discontinuous precipice extending toward MAT01 (Fig. 5). Observed at high incidence angles, the area around honeycombs MAT05,06 showed a rugged, furrowed topography that appeared to be buried under smooth dust deposits four months earlier (Fig. 6). As can perhaps be perceived in Fig. 3c, the ruggedness appeared pervasive and may not be easily discerned from the honeycombed textures. Numerous new depressions are visible on the neighboring terraced mound with occasionally exposed fractured edges, where the dust removal was probably non-uniform (Fig. 6b). A thin trench formed along an existing scarp suggests that the dust removal had been affected by local topography. As noted, the topography around SEQ02 and SEQ03 had apparently sharpened (Fig. D.2). The protrusion of a few boulder-like features seemed more pronounced in contrast to some more articulate depressions likely indicating retreat of the surface due to local dust removal.

Scattered depressions.

A common form of surface changes that we noted is the emergence of multiple shallow depressions. Observations of some distinct changes in Ma’at bordering Hatmehit are shown in Figs. 7, 8 and D.3, D.4 through comparison of two images taken in November 2014 and March 2015, respectively. The depressions formed more sporadically than at the honeycombs or pitted deposits. However, the size of these features, varying from (at least) decimeters to meters, is similar to the size of indentations of the honeycombs. The depressions probably resulted from the removal of dust deposits, as inferred from the excavation of materials around a few boulder-like blocks (Fig. 7). Some may have evolved from the existing concavities on the surface (Fig. 7). The fresh depressions were often accompanied by plenty of other subtle indications of changes, for example, the disappearance of dune-like structures (Fig. 7). Some changes occurred within the originally thin dust deposits that occasionally ran out to reveal rugged, fractured substrate (Fig. D.3).

Fresh scarps.

Fresh scarps had been found within the dusty deposits in several regions. The emergence of a scarp in Ma’at is easily explained by the removal of the surrounding dust deposits, as the edges of several boulders nearby had been revealed over time (Fig. 9). The rim of a trapezoidal-shaped, debris-filled depression by ASH01 had sharpened, which could have resulted from partial removal of the dust deposits, such as slumping or thinning (Fig. D.5). The deposits on the slope of the depression had been at least partially mobilized, with excavation evident around several distinguished boulders. One of the blocks appeared to have been exposed. A large irregular-shaped depression emerged near SET02 with scarp on the one side wandering into the adjacent field of talus (Fig. 10). There was a marked change in surface texture over the floor, where the previously uniform and stippled pattern of dust cover had diminished. We infer that this textural contrast may suggest a decrease in surface roughness, for example, by cleaning of large particles. Notably, a boulder at least 1 m in size had apparently been removed. Still in the Seth region but closer to Hapi, two semi-circular, likely meter-high scarps were visible in an image taken on January 23, 2015 (Fig. 11b). One feature can be confirmed to have emerged after early December 2014 (Fig. 11a). The scarp was hard to distinguish ten days later, possibly suggesting ongoing changes. During this time, the other scarp feature receded by more than 10 m. This event might be similar to the prominent changes observed in Imhotep between May and July 2015 (Groussin et al. 2015). Scarps were also observed in other areas, such as Babi and Atum (see Fig. D.6, for two new features near BAB02).

Erasure of pitted deposits.

With a similar texture, honeycombs may be akin to the pitted terrains (Shi et al. 2016a). Thomas et al. (2015a) noted that the pitted deposits are unlikely to have formed by direct dust deposition but indicate, rather, a textural evolution of the surface. The exact mechanisms of their formation are not understood. Mottola et al. (2015) proposed that the pits may have resulted from either collapse of the dust deposits ensuing the excavation of subsurface volatiles or mobilization of the deposits. At the time of writing, the pitted terrains at the border between Ma’at and Maftet had been visibly altered (Figs. 12a,b). The changes probably occurred around July 2015, that is, not long before perihelion (Figs. 12c,d,e). Heaps of pitted dust deposits were likely erased by erosion, leaving behind hollows within the residual dust cover. This serves as supporting evidence that the honeycombs indicate an evolutionary stage of the surface texture of the dust deposits. It also suggests that the erosion may both enhance and diminish the roughness of the dusty surface, whereas the mechanisms can be complicated.

2.4. Global distribution of identified surface changes

The honeycombs and the aforementioned changes of various forms are probably related (see an argument based on the timeline of the surface changes in Sect. 5). There is little question that the new outcrops were formed by the removal of the overlying dust (Figs. 5, D.5). The scattered depressions were probably akin to the indentations of honeycombs and resulted from surface erosion giving rise to the pitted texture at a later stage (Figs. 7, 8). The features SET01,02 and BAB02 appeared to have evolved as such. It can be difficult to clearly separate a honeycomb from the scattered depressions nearby (Figs. 4c, 8). Aside from an ensemble of small pits, an individual honeycomb can probably be better distinguished by the excavation of surrounding dust deposits (Figs. 1, 2). Many honeycombs, for example, MAT03-07 (Fig. 3), appeared to be isles of remnants from an expansive erosion that also affected their surroundings. We interpreted the fresh scarps, edges, and various depressions around SET02 and ASH01 as a consequence of surface erosion that had occurred beyond the honeycombs. If the honeycombs are indeed pitted remnants from erosion, they are subject to further textural variation (Fig. 12).

The distribution of the surface changes presented in this section as well as those discussed by Groussin et al. (2015) is shown in Fig. 13. The two general morphologies of smooth and rugged terrains are distinguished. Apart from the evident latitudinal concentration from about 20◦N to 40◦N, the presented changes seem to be more populated near morphological boundaries between smooth and rugged regions. It is very likely that the dust cover of 67P had been widely altered as the comet approached perihelion.

3. Change in surface roughness and erosion of the dust cover

Most of the identified potential changes appear to be an evolution in macroscopic roughness of the dust cover. We note that observational effects such as dissimilar viewing geometries and varied illumination conditions of images could belie or, at least, pose a challenge to the factuality of the intuited surface changes. Hence, more care is warranted in the comparison of the images as well as interpretation of any visual differences. This is the crux of our discussion in this section.

We use radiometrically-calibrated level-three images, as described by Tubiana et al. (2015). Let I_λ be the spectral radiance (in W m^{− 2} sr^{− 1} nm^{− 1}) of sunlight at wavelength, λ (in nm), scattered from the nucleus and detected by the camera. The radiance factor is defined by (Hapke 1981; Shkuratov et al. 2011), (1)where F_λ is the solar irradiance (in W m^{− 2} nm^{− 1}) at λ. The observations analyzed in this work were acquired at λ≈640 nm and 650 nm, between which F_λ varies insignificantly. Hereafter, we omit the subscript, λ, from the notation. It is illustrative to express the radiance factor as, (2)A and ℱ are the (wavelength-dependent) normal albedo and phase function, respectively. The latter is approximated here by the polynomial (Shkuratov et al. 2011; Schröder et al. 2013), (3)for some coefficients, c_k. The normal albedo, A, and phase function with c_k up to K = 3 are provided in Table 2 by fitting the Hapke model (Hapke 2002) for 67P derived by Fornasier et al. (2015). We note that c₀ = 1 by definition of Eq. (3), and that the phase angle, α, is measured in degrees. The disk function, D, is assumed to be approximated by the Lommel-Seeliger law: (4)where i and e denote the incidence and emission angles, respectively. γ is the photometric or luminace longitude (Hapke 2002; Shkuratov et al. 2011; the notation here is that of the latter). As shown in Fig. 14, γ is the angle between the camera (observer) and the component of the surface normal projected onto the “Sun-surface-observer” plane, measured positive if overlapped with α. The advantage of introducing D(α,γ) is that, when phase angle varies negligibly over the image, the disk function depends only on γ that varies with surface topography.

3.1. Macroscopic surface roughness

We applied the basic shape-from-shading technique, introduced by Rindfleisch (1966), to infer the surface roughness pattern in each image. This classical approach provides an efficient solution for topography profiling from a single image. The technique is briefly reviewed in Appendix B. Suppose the distance, r, between the (focal point of) camera and a certain reference point on the nucleus surface at P₀ in the image is known as r₀. The variation, Δr, due to the topographic variation from P₀ to any other point, P, in the direction of the Sun from the observer can be determined as (5)where ds denotes some differential step along the straight path from P₀ to P projected on the imaging plane, and where r_F is the focal length of the camera. tanγ is derived from the pixel intensities at s.

The procedure of the roughness analysis is illustrated in Fig. 15, via an image of MAT05. The Sun is projected upward onto the image plane (as can be inferred from the direction of the shadow cast downwards in Fig. 15a). The direction of integration on the right-hand side of Eq. (5) is therefore from bottom to top. The profile of topography along a single path indicated in Fig. 15a is shown in Fig. 15b. We note that the elevation varies opposite to Δr and is referred to a local plane given by the corresponding facet(s) of the existing shape model of the nucleus (Jorda et al. 2016). The constant and linear trends of the curve are removed (according to the remark C in Appendix B). Thus, the profile captures only the undulation arising from surface roughness but does not improve the measurement of absolute surface elevation from the shape model. The undulations of the curve suggest that depths of the indentations are typically a few decimeters. The depths are insignificant compared with the width of the entire honeycomb of tens of meters.

The roughness model in terms of elevation is derived by connecting profiles for all paths across the image (Fig. 15c). However, there seems to be a lack of unique method for connecting profiles. The reason is simply that, while Δr can be estimated via Eq. (5), the initial value, r₀, for each profile is arbitrary. We approximate the initial value by the distance from the camera to the local nucleus surface represented by the shape model. Errors arise from the improper initial value, discreteness of signals, and accumulation of noises along the path of integration that will sometimes cause unrealistically large undulation across profiles (Horn 1977). To alleviate this issue somewhat, we suppress any separation between adjacent profiles that can be fitted by the polynomial, Δr′ = ∑ _{k = 0}a_ks^k, up to degree five. It is not an unrealistic assumption that there should be no abrupt, large-scale topographic change, extensive cliffs for example, between two adjacent profiles.

The roughness model was used to simulate an image to be compared with the original image under the instantaneous illumination (Figs. 16a,b). Our aim is not to develop an elaborate shape model that minimizes the difference of a simulated image from the actual. Nevertheless, the similar sunlit patterns of surface in the images suggest that the derived topography is a valid characterization of the prominent surface roughness.

It is found that MAT05 was not present six months before; otherwise, it would have exhibited a resolved roughness that is not observed from the real observation which instead showed a far smoother texture (Figs. 16c,d). This analysis was repeated for many other features. Another example is provided in Fig. 17 for feature ASH01. The absence of the feature from earlier images in November 2014 cannot be attributed to visual effects, since the roughness would have been noticeable and differed from that in the real image.

We therefore confirm that the honeycombs were not enduring morphological features or surface texture, whose appearance varies significantly and can be concealed by illumination and viewing conditions. Instead, they indicate an increase in the surface roughness of the dust cover over time, on spatial scales from decimeters to meters.

3.2. Thickness of surface erosion

The surface changes in volume can be rigorously quantified by assessing the difference in shape models for the corresponding epochs of interest. It has been shown that the widths of the individual indentations are a few meters while their depths are typically a few decimeters. Thus, an accurate quantification would necessitate development and detailed comparison of multiple shape models with high resolution (Preusker et al. 2015; Jorda et al. 2016), an effort far beyond the present treatment in complexity.

We limited ourselves to the probable scenario in which the increase in surface roughness resulted from net (uneven) mass loss from the dust deposits (Schulz et al. 2015), that is, erosion caused by ice sublimation and dust ejection. This scenario is in keeping with the observed increase in activities of outgassing and dust ejection with decreasing heliocentric distance (Bockelée-Morvan et al. 2015; Della Corte et al. 2015; Fougere et al. 2016; Fulle et al. 2016c). Since the thickness of the dust cover is small relative to the width of a single honeycomb feature, it suffices to quantify the surface erosion in thickness. The topography or elevation of the roughness model is measured from an arbitrary level close to the nucleus. It is not possible to assess changes of the dust cover by directly comparing the two roughness models.

In the case of surface erosion, the thickness of the dust cover could only decrease, meaning that, the surface could only be lowered over time. Figure 18a is a comparison of the elevation profiles along the common path before and after change had occurred. The profile after the change exhibits a more abrupt and frequent variation, exemplifying the increase in surface roughness. The profile is placed below the earlier curve such that it grazes but does not cross the latter. The difference between the two profiles is 1.3 m on average. Due to multiple assumptions and many unknown sources of uncertainty, it is impractical to quantify the error of this estimate. To contain the uncertainty somewhat, we remove the quartic trend of the profiles in Fig. 18b such that the difference arises only from the higher-frequency undulations of the blue curve, that indicate an increase in surface roughness. In this case, the average difference of 0.51 m is probably a lower limit of the actual erosion.

It can be argued that 1 m is a robust order-of-magnitude estimate for the thickness of erosion. The lower limit of this estimate is strictly given by the pixel resolutions of images, which are about 0.5 m and always larger than 0.1 m. Changes below this detection limit would not have been observed (in other words, we cannot distinguish between 1 m and 0.5 m). The upper limit, on the other hand, is given by the thickness of the dust cover, estimated to be up to a few meters in general (Thomas et al. 2015a; Mottola et al. 2015). It follows that the erosion of the dust cover would not exceed a few meters. Thus, we may formally express the maximum uncertainty of this estimate as 0.5 m< Δx<5 m, while the upper bound of the dust cover is from Thomas et al. (2015a).

The estimated thickness of erosion of approximately 1 m is also reflected from other types of surface changes (Sect. 2.3). In Fig. 19, we compare the topographic profiles before and after the surface change, where the overall deviation of the curves is, again, attributed to the recess of the surface level. The potential outcropping structure shown in Fig. 5 in Ma’at could have been exposed by thinning of the dust cover by at least 1 m in six months (Fig. 19a). The fresh, scarp-like feature shown in Fig. 9b may have formed by the removal of the surrounding top few decimeters of dust deposits (Fig. 19b).

4. Accumulated insolation and modeled erosion

The surface changes could be induced by various mechanisms. As a first attempt, we have concerned ourselves with the general scenario where the loss of dust deposits was driven by the outgassing of volatiles from the nucleus.

The sublimation of water ice, the most abundant volatile measured in the coma of 67P, is predominantly influenced by insolation and exhibits clear diurnal variation (Keller et al. 2015). The sublimation of more volatile ices, such as CO₂ and CO, is strongly dependent on the interior thermal and structural conditions of the nucleus (Yabushita 1995; Enzian et al. 1997; Belton 2010). The fact that the surface changes (in this work) all emerged between September 2014 and March 2015 suggests they were likely to have been the result of intensifying insolation due to a decreasing heliocentric distance from about 3.5 to 2 AU. In particular, the sub-solar latitude lowered from 45◦N to 15◦N during the aforementioned period, roughly the range in which the surface changes were distributed. This strongly suggests that the decrease in local solar zenith (incidence) angle had further invigorated changes of the surface.

4.1. Accumulated insolation

For a given location on the nucleus, we evaluate the accumulated energy of insolation from t₀ to t₁ by, (6)where (7)is the instant flux of absorbed energy. C_⊙ = 1361 W m^{− 2} is the solar constant and A_B = 0.05 is adopted for the surface Bond albedo. The heliocentric distance of 67P, r_⊙, and the solar incidence angle, i_⊙, are time-varying. They are calculated from the estimated position and orientation of the nucleus (approximated by a shape model) from SPICE kernels (Acton 1996). δ_⊙ = 1 if the local surface is illuminated and δ_⊙ = 0 otherwise, that is, if i_⊙ ≥ 90◦or if the local surface is shadowed by surrounding topography. r_AU denotes one AU in unites of r_⊙ (introduced for the sole purpose of nondimensionalization).

We applied local shape models truncated from the global SPG shape model with a spatial resolution of approximately 10 m (Preusker et al. 2015) to derive the accumulated insolation from early September 2014 through February 2015. We focused on the areas around honeycombs MAT02-11 and that around ASH01, SET01 and SET02, located at roughly the same latitudes on separate lobes (Fig. 13). Figures 20a and 21a show the areal pattern of accumulated insolation superposed on the respective context images. The insolation is influenced by topography. The accumulated insolation over the cliffs (as always inferred from their rugged, fractured appearance) shows clear dichotomy: the south-facing walls and the nearby shielded floors absorbed limited energy below 6 ×108 J m^{− 2} as a result of scarce sunlight during the northern summer; whereas, those looking to the north had been amply illuminated, accumulating more than 109 J m^{− 2} of energy, over six months. This contrast is also widely observable with other abrupt topographies such as outcrops and boulder-like monoliths. The honeycombs were all located within open smooth terrains with an accumulated insolation above E_⊙≈109 J m^{− 2}. More locally, some honeycombs could be distinguished by patches bounded by higher and lower accumulated insolation on opposite sides, for example, MAT03-05 (Fig. 20a). This illustrates that these honeycombs were within locally flat areas bordered by slopes or scarps well illuminated on northward side but less so on the southward side.

Higher accumulation in Fig. 20a occurred near the rim of a large pit within the dust-covered plateau to the north (see Vincent et al. 2015 and Ip et al. 2016 for a description of the pits on 67P). The enhancement toward the north was also somewhat noticeable in Fig. 21a. It can be inferred that the areas further north were more persistently, albeit softly, illuminated overhead when the comet remained further out from the Sun (before 2015).

This observation is more conspicuous on the global scale. The accumulated insolation overall increases toward higher latitudes, warmed by the trailing heat from a prolonged yet mild summer in the northern hemi-nucleus (Figs. 20b, 21b). As noted, sharp variations may occur over large-scale abrupt topographies, for example, the (north-facing) cliffs had been subject to strongest insolation. We also note that the bi-lobed shape and abrupt topography of the nucleus give rise to unintuitive lines of (equi-)latitudes and, accordingly, notable variation of insolation along the same latitudes.

4.2. Modeled erosion of water ice

In light of the surface changes as an evident consequence of intensifying insolation, it is of primary importance to assess the role of water ice in driving or facilitating the surface erosion. On the other hand, there is a lack of observational indication that the activities of supervolatiles played a vital role in triggering the surface changes (the discussion is given Sect. 6.1).

We employed a thermal model to assess the accumulated erosion of water ice in response to varying insolation. It is assumed that nucleus is covered by a dry, porous dust mantle that is insulating and (gas-)flow-resistant. The dust mantle, some millimeters up to a few centimeters thick, should be distinguished from the dust cover as an observed morphology. The subsurface below the dust mantle is a mixture of water ice and dust (Fig. C.1). Hence, it is stipulated here that the dust cover is icy (however slightly) from some shallow depths (De Sanctis et al. 2015; Spohn et al. 2015; Shi et al. 2016a). The thermal model is hereafter referred to as the “dust mantle” model.

Sublimation of water ice occurs at the ice front with water vapor diffusing through the overlying dust mantle. The mass flux per unit surface area is given by (Appendix C), (8)where is the sublimation flux of pure solid water ice beneath the dust mantle, that varies with layer thickness, X_i, size of the constituent dust particles, d_p, and temperature at the ice front, T_i (Gundlach et al. 2011). The empirical factor, 0 < f ≤ 1, loosely characterizes the “icy area fraction” and arises from the fact that the nucleus interior is composed of dusty ice (or icy dust) rather than pure ice, as clarified by Crifo (1997). For instance, a high dust-to-ice ratio gives rise to low f in the dust cover.

The topography of the nucleus is represented by a shape model. The resolution (the size of a single facet) of the shape model exceeds the diurnal thermal skin depth for 67P (likely a few centimeters at most) by at least two orders of magnitude. With the plane-parallel assumption, the profile of subsurface temperatures can be derived by solving the 1D heat equation for each facet with boundary conditions of energy balance. The reader is referred to Appendix C for details of the solution.

The thermal model is laterally homogeneous over the nucleus, that is, it does not reflect any variation of the (sub)surface structural or thermal properties with different morphologies. While it is a consensus that activity was stronger in Hapi on a global scale, the contrast is hardly striking (Fougere et al. 2016; Fink et al. 2016). For the time being, we have considered a homogeneous model viable and instructive for a preliminary study. We further assume the thickness of the dust mantle to be constant, in spite of the complex, varying vertical structures observed both experimentally and in situ (Spohn & Benkhoff 1990; De Sanctis et al. 2015; Biele et al. 2015). The erosion of the dust cover can be conceptualized as a top-down steady-state process that is always commensurate with erosion of the water ice for a given dust-to-ice ratio.

4.2.1. Choice of model parameters

The modeled temperatures and sublimation flux are sensitive to the choice of parameters, which is not unique. For both the dry dust mantle and the underlying icy dust, we adopt κ =0.002 W K^-1 m^-1, c = 1000 J K^-1 kg^-1, and ρ = 500 kg m^{− 3}. They correspond to a thermal inertia of I= (κcρ)^{1 / 2}≈30 W K^-1 m^-2 s^{1 / 2}, being the median in the range between 10 and 50 W K^-1 m^-2 s^{1 / 2} derived by MIRO (the Microwave Instrument for the Rosetta Orbiter) measurements around September 2014 (Gulkis et al. 2015). We assume that the dust mantle consists uniformly of particles 1 mm in diameter. The choice of dust mantle thickness, X_i, and the icy area fraction of the subsurface, f, will be discussed later.

4.3. Comparison with in situ measurements

The water production rate of 67P has been measured in situ by various instruments onboard Rosetta and derived from ground based observations. Figure 22 shows the data of measurements presented by Fougere et al. (2016) from a heliocentric distance of 4 AU to 1.5 AU inbound. Some estimates based on the early measurements of MIRO are appended (Gulkis et al. 2015; Biver et al. 2015). These measurements allow us to assess the feasibility of the thermal model for estimating the erosion of water ice.

We modeled the water production of 67P averaged over one comet rotation during the period of reported measurements, (9)where dS denotes differential element of the total surface area, S, of the nucleus, and where t_P ≈ 12.5 h⁴ is the rotation period of the nucleus. In this calculation, dS is given by the area of a single triangular facet of the shape model.

At first, we adopted f= 0.1 for the icy area fraction, which is roughly compatible with the ice content of about 10% relative to the refractory materials observed in Hapi (De Sanctis et al. 2015). A range of values has been tested for the thickness of the dust mantle, X_i, a parameter strongly affecting the production rate. In the first case, we assumed that the dust mantle is 5 mm thick. The model overestimates the measurements by a factor of about six. Increasing the thickness of the dust mantle reduces the production rate. For a thickness of 1 cm, the modeled production rates are best-fitting the measured. With an even thicker dust mantle of 1.5 cm, the model tends to underestimate the water production rate. Additionally, we reduced the icy fraction to f= 0.01 that corresponds to the water ice abundance of about 1% in the subsurface. It is found that the water production of the nucleus with a mantle thickness of X_i = 5 mm is similar to the results for X_i = 1 cm and f= 0.1. In all cases, the modeled production shows a temporal trend in general agreement with that of the measurements, obviously in response to the intensifying insolation as the comet approached perihelion. Hereafter, we adopt X_i = 1 cm and f= 0.1 for the thickness of the dust mantle for deriving the erosion of water ice.

4.4. Results

The total erosion of water ice from time t₀ to t₁ is calculated as, (10)where the instant outgassing flux, Z, is given by Eq. (8).

Figures 23 and 24 show the modeled water ice erosion in the areas around MAT02-11, and near ASH01, SET01 and SET02, respectively, accumulated in six months from early September 2014. The graphs are superposed on the same context images as in Figs. 20 and 21, respectively. The pattern of Δm_H2O resembles that of insolation, E_⊙. This is expected, since the sublimation of water ice is essentially insolation driven. The north-facing cliffs were subject to the strongest erosion of nearly 10 kg m^-2, while the minimum occurred at the least insolated locations, for example, southward walls, the nearby floors as well as the bottom of pits. The open areas where the honeycombs emerged were subject to a moderate amount of erosion, mostly in excess of 5 kg m^-2 (induced by E_⊙>109 J m^-2; see Figs. 20 and 21). The total erosion overall increases further north, as a direct result of higher accumulation of insolation.

It can be shown that 5 kg m^-2 is a reasonable model estimate for the lower limit of water ice erosion at any locations where the surface changes were detected. The derived pattern of Δm_H2O, which conforms to that of accumulated insolation, supports the interpretation of the observed surface changes that the erosion had probably occurred widely over the northern hemi-nucleus of 67P.

5. Seasonal erosion and restoration of the dust cover

Most of the surface changes were distributed from about 20◦N to 40◦N (Table 1, Fig. 13). The fact that the changes occurred⁵ concurrently with the movement of the sub-solar point across the same latitudes suggests the change over the dust cover is a seasonal phenomenon. We recall a timeline of detections of certain cases, which, however, may not specify their temporal occurrence:

• 1

  MAT01, MAT02 were detected in January while MAT03-08were detected in mid February 2015 (though MAT06 and MAT08were less distinguishable, Fig. 3b).

• 2

  Changes around ASH01, SET01 and SET02 by the Ash-Seth border were perceptible in February 2015 (Fig. 4b); the moving scarp in Seth occurred at the end of January (Figs. 11b and c).

• 3

  There was a lack of significant changes around BAB01, BAB02 before the features were detected at the end of March 2015. However, hints of changes around BAB02 date back to February (Figs. D.1c and d).

We also note:

• 4

  The prominent changes of moving scarps in Imhotep occurredbetween late May and early July in 2015 (Groussinet al. 2015). The removal of thepitted dust deposits at the border between Ma’at and Maftetprobably occurred around the same time.

Meanwhile, it is shown that (Table 1; Fig. 13):

• i

  MAT01, MAT02 are at around 33◦N; MAT03-08 are furthersouth at 25±5◦N.

• ii

  ASH01, SET01 and SET02 are located at about 25◦N.

• iii

  BAB01 and BAB02 are at about 21◦N and 15◦N, respectively, among other southernmost features (MAT11 and SEQ03, for example).

• iv

  Changes in Imhotep are further south than all surface changes analyzed in this work (Groussin et al. 2015). The erased pitted terrains were located near the equator.

Inspecting the above points by the number (numeral), we infer that the features far north, MAT01 and MAT02, appeared first, followed by ASH01, SET01, SET02, and MAT03-08, slightly lower in latitude. BAB01 and BAB02 were possibly the last honeycombs to emerge with local insolation to peak still later. All these changes predated the changes in Imhotep as well as the erasure of the pitted terrains in Ma’at and Maftet. Most of the other honeycombs were detected at the end of March 2015 and scattered in latitudes. The timing of their occurrences proved more difficult, owing to the fewer images found and often long gaps in between. The inferred temporal occurrence of the changes supports the scenario that they were directly responding to the seasonal migration of the sub-solar point and the decreasing heliocentric distance.

5.1. Transition strip of global morphology

More insight into this seasonal phenomenon is gained by noting that the latitudinal strip between about 20◦N to 40◦N where most of the surface changes took place marks a global transition in surface roughness between the general morphologies in the northern and southern hemi-nuclei (Fig. 25). Smooth dust deposits are prominent and likely perennial at higher latitudes on each lobe; whereas, they were diminished further south revealing the prevalent, more rugged terrains. The transition strip is where the dust cover overall tapers off toward the south. However, due to the sheer topography, local variations are expected to be common.

5.2. Thinning of the dust cover pre-perihelion

By the northern autumnal equinox, the surface roughness had increased along the transition strip, as evidenced by the emergence of honeycombs and other changes in plenty, marginally (un)covering the subsurface and showing a moderate surface roughness. Our thermal modeling analysis suggests that the areas further north had been subject to greater accumulated insolation and hence erosion. If the erosion had occurred uniformly as modeled, that is, due to water ice sublimation driven purely by insolation and below a dust mantle of constant thickness, it is then likely that the initial thickness of the dust cover was not uniform. This means that the stronger erosion had not exhausted the dust cover further north, while the thinner deposits between about 20◦N and 40◦N had been rapidly depleted.

Another possibility should not be dismissed, however. The prominent change along the transition strip might be attributable to the non-uniform erosion over different latitudes. In this case, the erosion would increase more rapidly or steeply than induced only by decreasing heliocentric distance. One possible indication is the late emergence of all changes from early 2015. Recall, for instance, that the earliest, subtle textual development with BAB01, ASH01, SET01 and SET02, did not appear until February 2015. Gundlach et al. (2015) showed that 2.5AU (in January 2015) may mark the onset of dust ejection efficiently driven by water outgassing. The ejected particles were likely sized around 1 dm. The significant ejection of large particles could have enhanced the erosion of the dust cover from early 2015. The present thermal model with fixed parametrization did not reflect this “acceleration”, as it neglected the micro-physics of ejection of the dust particles by outgassing and, not least of all, assumed constant thickness of overlying dry dust mantle for simplification. We note that the two possibilities, namely, the non-uniform dust cover and accelerated erosions, are not mutually exclusive. However, it is beyond the scope of the present analysis to disentangle the two scenarios.

5.3. Restoration of the dust cover post-perihelion

The dust cover probably consists of large particles (Mottola et al. 2015). Thomas et al. (2015a) demonstrated that the local and global redistribution of dust particles over the nucleus is plausible and probably a general phenomenon. Particles not accelerated to escape speed by gas drag tend to fall back and accumulate on the nucleus. Large and heavy particles are likely the main source of contribution⁶. We also refer the reader to the dedicated analysis on dust transport on 67P by Lai et al. (2016).

Keller et al. (2015) envisioned a prominent dust deposition process in the northern hemi-nucleus during the upcoming perihelion passage of 67P. The particles would originate dominantly from the south, which is subject to the strongest insolation and, hence, surface erosion during perihelion. This is substantiated by the more recent observations several months after perihelion, where the pitted texture of honeycombs was no longer visible. The disappearance of some pronounced features is shown in Fig. 26. The cluster of honeycombs in Ma’at and three features at the Ash/Seth border had vanished. The scarp in Seth had been erased. The sheer edges exposed by thinning of the overlying dust deposits before perihelion had been smoothed out. The surface texture of the feature locations had reverted to that of smooth dust deposits with more random and uniform roughness. Here, we have omitted the roughness analysis as discussed in Sect. 3.1. It suffices to note that the disappearance of features likely suggests that the lower limit of dust deposition during the perihelion is comparable to that of erosion before perihelion, which is approximately 1m.

5.4. Seasonal retreat and expansion of the dust cover

The pre-perihelion erosion of the dust cover in the northern hemi-nucleus of 67P was compensated by a prominent dust deposition at perihelion. These notable changes are a manifestation of the seasonality in the surface texture as well as in the activity of 67P.

It is along the transition strip that the surface changes would appear most prominent. As noted, this prominence could result from the thinner dust cover than higher latitudes; additionally, it may result from the more significant dust erosion over mid latitudes, for example, due to removal of sizable particles. Both cases accommodate the significant disturbance and even depletion of the dust cover along the transition strip with the deposits further north remaining thick and smooth (Fig. 25). Similarly, the dust deposition around perihelion altered more significantly the surface texture along the transition strip, when the reinstated dust cover masked the honeycombs distinguishable after the overlying deposits had thinned out before perihelion. The textural change of dust deposits in higher latitudes would probably be less notable, since the further enrichment to the (already) eminent dust cover would not alter its surface roughness. In effect, the edges of the dust cover retreat and advance on each lobe of the nucleus, as the dust deposits are removed and restored along the transition strip. We may expect this morphological phenomenon to repeat in each orbit (assuming the current orbital and rotational states to remain).

6. Discussion

6.1. Plausibility of surface erosion caused by sublimation of water ice

In our thermal analysis it was assumed that the sublimation of water ice was the main driver of the observed dust erosion (Sect. 4). The validity of this assumption lies with the source of observed dust activity on 67P.

Thermal modeling by Keller et al. (2015) illustrated that the dust activity driven by water outgassing follows insolation over the nucleus. Vincent et al. (2016) showed that the visual pattern of jet collimation is indeed repetitive on a diurnal basis. Kramer & Noack (2016) demonstrated that the pattern of coma brightness conforms to that of dust emission from a homogeneous nucleus, echoing the earlier conclusion by Crifo et al. (2002) that the topography per se could already produce complicated near-nucleus dust coma structures, such as collimated jets. Through comprehensive modeling of the gas coma, Fougere et al. (2016) found that a homogeneous outgassing model could be improved by taking into account an enhancement from Hapi to better fit the water column densities observed by VIRTIS (the Visible and Infrared Thermal Imaging Spectrometer).

That the dust activity, in all probability, occurs extensively (albeit non-uniformly) indicates that the dust cover as a prevalent morphology was active. It may also partly explain the challenge and limitation with determining local active areas or sources of dust activities from the nucleus (Bruck Syal et al. 2013; Farnham et al. 2013; Lara et al. 2015; Lin et al. 2015; Shi et al. 2016a; Vincent et al. 2016). The dust emissions from the nucleus are usually distinguished against the dark background of rarified dust coma. In some cases, the high-resolution OSIRIS observations reveal faint localized structures of enhanced brightness against the shadowed nucleus. Such brightness, sometimes with distinct spiked pattern fading out immediately from the terminator, likely results from the sunlit dust particles aloft in the foreground (Thomas et al. 2015a; Vincent et al. 2016; Shi et al. 2016a,b).

The dust activities shown in Fig. D.8 occurred near the instantaneous dusk terminator. The area around honeycomb MAT02 as well as the nearby smooth dust deposits were probably the sources of dust emission. We note that the “sunset jets”, observed in the vicinity as early as February 2015, were probably continuation of day-side activities that would eventually cease following the subsidence of water ice activity from the shallow subsurface (Shi et al. 2016a). We infer that the observed dust activity in Fig. D.8 near the dusk terminator were also induced by water ice sublimation.

Near-surface activities were also observable in other areas. The entire tri-border area between Ma’at, Nut, and Serqet was visibly active in an image taken at the end of March 2015 (Fig. D.9). The dust emissions were contrasted against the shadows cast by adjacent topographies such as cliffs and monoliths. Similar observations were found in earlier images, as well. For instance, Fig. D.10 shows two cases of dust emission against the shadows, likely from the bright foreground close to the terminator. This observation was taken one hour after Fig. 11b showing the retreating scarp in the same area in the Seth region. We note that the activities shown in Figs. D.9 and D.10 occurred by the dawn terminator (that is, shadows were receding locally). The activity was likely induced by the sublimation of water ice that had recondensed close to or even at the surface overnight (De Sanctis et al. 2015). Dust emissions around the dawn terminator have been widely observed by OSIRIS, likely suggesting the prevalence of the water cycle in the shallow subsurface and, thus, its fundamental role in triggering dust activity (Shi et al. 2016b).

Secondary role of super-volatiles.

While it is difficult to exclude definitively any possibility, we consider it unlikely for super-volatiles to be the main driver of the observed changes. The observational evidence thus far indicates higher abundance of CO₂ from the unilluminated southern hemi-nucleus (Hässig et al. 2015) and, in particular, its depletion in Ma’at, Ash, Seth(-Hapi), and (Aten-)Babi regions in the northern hemi-nucleus where surface changes were detected (Bockelée-Morvan et al. 2015; Fink et al. 2016; Migliorini et al. 2016). Podolak et al. (2016) suggested that the measured CO₂ production rate, a few percent of water production, is insufficient to account for the typical speed of dust grains at meters per second.

The thickness and timescale of the observed surface changes also shed some light on the role of super-volatiles. The timescale of heat conduction is given by (Huebner et al. 2006), τ ~ (xI /κ)². The depth of x≈0.5 m, which is less than the inferred thickness of erosion, can be roughly conceptualized as a “seasonal” level where τ ~ 1 yr. At x<0.5 m, the activity should respond to increasing insolation within a few months, potentially compatible with the timescale of surface changes. In this case, however, super-volatile ices would have been present in the dust cover, either preserved in the dust particles ejected from the south around the previous perihelion or having condensed out of gases from the interior or coma. Both mechanisms seem elusive and have yet to be substantiated by observations or modeling. Heat propagation to x>0.5 m would take longer than the time frame of surface changes, in which case, the observed changes would imply a (spatial) distribution of super-volatile ices in deeper subsurface. This scenario is contested by observations, if not excluded (Bockelée-Morvan et al. 2015; Fink et al. 2016; Migliorini et al. 2016).

6.2. Mass ratio of dust and water ice in the dust cover

The mass loss of surface erosion per unit surface area is ΔM = Δm_dust + Δm_H2O (kg m^-2, where Δm_dust and Δm_H2O are the mass components of dust and water ice, respectively. Let the dust-to-(water-)ice mass ratio be denoted by μ = Δm_dust/ Δm_H2O, such that the total erosion is expressible as ΔM = (1 + μ)Δm_H2O. The dust-to-ice ratio can be derived as (11)where ρ is the density of the dust cover (as a dust-ice mixture) and Δx the thickness of erosion. We adopt Δx = 1 m as an order-of-magnitude approximation. It is assumed that the density of the dust cover conforms to the bulk density of the nucleus, ρ = 530 kg m^-3 (Pätzold et al. 2016). Therefore, the total mass of erosion is approximated as ΔM = 530 kg m^-2.

We discuss two instructive cases of μ. We note that μ is always referred only to the dust cover. That is, we always distinguish μ from the dust-to-ice ratio computed as the ratio between the total dust and gas loss rates provided by observations performed in the coma (Rotundi et al. 2015; Fulle et al. 2016c).

6.2.1. Case A – Lower limit: μ = 4

On a local scale, the lower limit of the dust-to-ice ratio is given by the maximum erosion of water ice. To derive a reasonable estimate, an alternative thermal model to the dust mantle model in Sect. 4 is applied, which is termed the dusty ice model as a distinction hereafter. It is assumed that ice is exposed and sublimates on the nucleus surface. While similar to Models A and D in Keller et al. (2015), the dusty ice model here accounts specifically for the influence of the factor, f, for the icy area fraction (Crifo 1997), a parameter also incorporated in the dust mantle model. The dusty ice model coincides with Model D in Keller et al. (2015) for f= 1.

We use the honeycomb feature MAT05 as an example, representative of changes along the same latitudes around 30◦N (for example, ASH01, SET01 and 02, and etc.). This is reasonable because all changes discussed in this work were located in well illuminated areas. In the extreme case of f= 1, that is, if the comet nucleus consists only of water ice, the erosion accumulated in six months from September 1, 2014 is Δm_H2O = 175 kg m^-2. Although this estimate seems to yield according to Eq. (11), a contradiction arises from the assumption Δm_dust = 0 that would lead to μ = 0. Therefore, μ = 2 underestimates the lower limit (in the presence of dust).

We performed calculations for a range of discrete values of f_(k) for k = 1,2,..., to search for a consistent μ. Again, referring to Eq. (C.6), each f_(k) corresponds to a tentative dust-to-ice ratio of . Then, Δm_H2O was derived which yields via Eq. (11). A valid solution ought to satisfy . The comparison of the curve for the analytic expression μ = f^-1 − 1 and the resulting discrete are shown in Fig. 27. The two curves intersect at around f= 0.2, which corresponds to the lower limit of the dust-to-ice ratio of μ = 3.8.

This lower limit is consistent with other findings from Rosetta instruments. Firstly, the water ice content over the nucleus surface of 67P is up to 1% on average according to VIRTIS measurements (Capaccioni et al. 2015). Thus, the icy area fraction of 20% cannot exist everywhere on the nucleus surface, and a dust-to-ice ratio of μ = 4 would correspond to a strong localized water-ice enrichment. We note that the dusty ice model assuming f= 0.2 over-estimates the global water production by a factor of about 20. Secondly, in the most water-rich region of Hapi in the north, transient water ice up to 10-15% relative to the non-volatile materials had been observed (De Sanctis et al. 2015), that does not exceed the inferred icy area fraction of 20% on the surface. Thirdly, water ice was exposed at scattered locations on the nucleus, which measured up to 4% in abundance (Filacchione et al. 2016; Barucci et al. 2016). The detected water ice consists of large, millimeter-sized grains that were likely sourced from sintered, ice-rich layers (Filacchione et al. 2016), suggesting that f= 0.2 may be an overestimate in the subsurface, as well.

Therefore, it can be justified that μ ≈ 4 represents a lower limit for the dust-to-ice ratio at any location of surface changes described in Sect. 2. It should also be emphasized that the observational evidence of exposure of water ice at the specific locations of detected surface changes is absent thus far. Hence, μ ≈ 4 is likely still an underestimate for the lower limit. For instance, if a dry dust mantle overlies the dust-ice mixture, the curve for in Fig. 27 would shift upward and, thus, increase the lower limit of dust-to-ice ratio.

6.2.2. Case B – Homogeneous dust cover: 50 ≤μ≤ 100

Alternatively, we assume that the dust cover is homogeneous at spatial scales larger than 10 m, which is the resolution of the shape model for the thermal analysis (Sect. 4). This assumption is inferred from the observation that the dust cover had been globally renewed after perihelion, where the re-deposition of dust particles ejected from the south was probably undifferentiated in terms of ice abundance (note another justification of the assumption of homogeneity in Sect. 6.1). The widely-distributed surface changes then enabled a “sampling” of the constant dust-to-ice ratio across the dust cover. We adopt 5 ≤ Δm_H2O≤10 kg m^-2 as derived via the dust mantle model (Fig. 23), representative of the erosion at the locations of most changes. The dust-to-ice ratio in the homogeneous dust cover falls in the range of 50 ≤ μ ≤ 100. This suggests that the dust cover is overall water-depleted at least in the first meter below the surface (Capaccioni et al. 2015).

6.2.3. Chessboard scenario of locally varying dust-to-ice ratio

The cases A and B above distinguish two possibilities of extreme inhomogeneity and perfect homogeneity of the dust cover in terms of water ice abundance. We note that, while the lower limit of dust-to-ice ratio is given by μ = 4 (case A), μ ∈ [ 50,100 ] in the latter case should not be misinterpreted as the upper limit of dust-to-ice ratio. The dust-to-ice ratio is always unbounded from the upper end, that is, μ = ∞ if the dust cover is completed desiccated.

The global distribution of the surface changes suggests that the dust cover had been widely eroded. On the other hand, the erosion may be non-uniform on a local scale, as indicated by the increased surface roughness over time for example. The non-uniform changes may have resulted from inhomogeneities in composition, for example, the abundance of water and other volatile ices as well as difference in material strengths. As noted in Sect. 5.2, since all changes analyzed in this paper were associated with the dust cover, non-uniformity may simply indicate uneven thickness of the deposits over more rugged substrate. As pointed out, we are unable to resolve inhomogeneities in dust-to-ice ratio on scales below 10 m.

The possibility of a hybrid scenario certainly exists, such that higher erosion occurred at scattered positions, for example, small individual pits, with locally higher ice abundance and, thus, concentration of water activity. In effect, this corresponds to a “chessboard” pattern in terms of ice abundance across the dust cover. The dust-to-ice ratio is higher on average and more uniform over larger scales.

6.3. Non-escaping dust

The in situ measurements of dust and water productions revealed a dust-to-ice ratio of 4( ± 2) averaged over the sunlit surface of 67P up to perihelion (Rotundi et al. 2015; Fulle et al. 2016c). The ratio refers to mass loss of the comet, that is, the mass of dust and water escaping from the nucleus. Not all dust removed in the surface erosion will escape. We introduce the non-escaping percentage, η, as the mass fraction of dust ejected from yet still bound by the nucleus, (12)where [ Δm_dust ] = ^∫_SΔmdS kg refers to the mass loss over the surface area of S. refers to the mass loss of the escaping dust. The non-escaping mass of η [ Δm_dust ] may enrich the dust coma, or/and fall back to the nucleus (Agarwal et al. 2016).

Again, we assume that the nucleus is homogeneous at scales greater than 10 m, and let for simplicity⁷. Taking the average dust-to-ice ratio, μ, over the nucleus surface as a variable, the non-escaping portion is given by η = 1 − 4μ^-1 (shown in Fig. 28).

6.3.1. Case A – Dust ejected from exposed ice: full-escaping

The limiting case of η = 0 (negative η is forbidden) where all dust escapes the nucleus is μ = 4. Referring to case A in Sect. 6.2.1, we see that this corresponds to the scenario that dust production was concentrated from enduringly exposed water ice on the nucleus surface (Filacchione et al. 2016; Barucci et al. 2016).

6.3.2. Case B – Dust ejected from the south: half-escaping

The bulk dust-to-ice ratio of the entire nucleus is found to be 8.5 (Fulle et al. 2016b). The ratio is probably representative of dust ejected from the southern hemi-nucleus, where the pristine constituents were exposed with the removal of top layers by intense activity of water and super-volatiles around perihelion. Then, it can be inferred that slightly more than 50% of the dust ejected from the south was still nucleus-bound, part of which had subsequently accumulated in the north replenishing the dust cover observed after perihelion (Fig. 26).

6.3.3. Case C – Dust ejected from the northern dust cover: limited-escaping

For a homogeneous dust cover where 50 ≤ μ ≤ 100, less than 10% of the ejected dust would escape from the nucleus. Excluding Hapi and Imhotep, the dust cover spans roughly 8 km² (Thomas et al. 2015b). For instance, if the entire deposits had been eroded by m on average, the mass loss would be in six months from September 2014. If = 1 m (Fig. 19), a total mass ejection of 4 × 10⁹kg would have resulted from an average ejection rate in excess of 250 kg s^-1. Based on GIADA (Grain Impact Analyser and Dust Accumulator) and OSIRIS observations, Fulle et al. (2016c) found that dust particles up to 0.1 kg totaled a mass loss rate of 70 kg s^-1 at the end of March 2015. The time-average of 35 kg s^-1 in the previous six months amounts to no more than 15% of the dust ejection rate from the dust cover. In the extreme case, if = 0.1 m which corresponds to μ = 8.5 (Fulle et al. 2016b)⁸, the average dust ejection rate is 25 kg s^-1. Since the dust cover analyzed in this paper accounts for less than half of surface area in the northern hemi-nucleus, less than half of the ejected dust from the dust cover may leave the comet (which is similar to the dust ejected from the south).

6.3.4. Global re-distribution of dust

It is inferred that the ejected particles from the dust cover may not have been a significant source of mass loss from the comet. Large particles were not energized to escape in bulk or even reach high altitudes, but were accumulated in the (inner-)coma or redistributed back onto the nucleus surface. For instance, particles drifting lower than spacecraft around the nucleus were not collected by the onboard instruments and, thus, filtered out in the derived dust production. The dust-to-water ratio of 4( ± 2) suggested by Rotundi et al. (2015) indicates a lower limit (Fulle et al. 2016b). An implication is that the dust-to-water ratio provided by coma observations could be non-uniform. A higher value may result when measured in the inner coma than at distances of Rosetta from the nucleus.

We suspect that the redistribution of dust is a general mechanism (Thomas et al. 2015a; Kramer & Noack 2015; Lai et al. 2016). The south-to-north transport prevailed during perihelion (Keller et al. 2015), where about half of the mass did not escape from the comet. The non-escaping portion could be even larger for dust ejected from the deposits in the northern hemi-nucleus before perihelion. The dust could be bound for other parts of the nucleus surface, more likely over the unlit southern hemi-nucleus and Hapi (Lai et al. 2016). Assuming that all non-escaping dust must fall back to the nucleus surface, the thickness of redeposition is a linear function of average thickness of erosion over the dust cover, that is, , which is proportional to the average dust-to-ice ratio (Fig. 28). If = 1 m, the ejecta would have been spread over the southern hemi-nucleus and Hapi by about 0.3 m in thickness. We note that the exact average thickness of erosion cannot be determined here, which needs to be constrained via a dedicated analysis in the future.

7. Conclusion

The main conclusions of this study are summarized as follows:

• 1

  Numerous changes had been observed by OSIRIS, thatoccurred within the dust deposits between September 2014through March 2015. The surface changes signify a widespreaderosion over the dust cover in the northern hemi-nucleus as 67Papproached perihelion, as predicted by Schulzet al. (2015). The surface erosion atmost locations where changes had been detected may have wellreached 1 m and, at least, several decimeters.

• 2

  The pre-perihelion surface changes analyzed here were concentrated in latitudes between about 20◦N and 40◦N, where the dust deposits thin out and transition into the rugged terrains prevalent toward the south. The occurrence of changes followed temporally the southward migration of the sub-solar point across the same latitudes.

• 3

  The dust cover was active before perihelion, at least since early 2015. The distribution and timeline of changes suggest that the erosion was more likely driven by the sublimation of water ice than that of more volatile species. The exact mechanisms of dust ejection and surface erosion are unknown.

• 4

  We applied a dust-mantle thermal model assuming the nucleus is a mixture of dust and ice overlain uniformly by a dry mantle. The mantle thickness and the icy area fraction of the subsurface, which is assumed to be governed by the dust-to-ice ratio, are adjusted to yield the water production rates which could fit the Rosetta measurements before perihelion. The locations of surface changes were subject to ample insolation and an accumulated erosion of water ice over 5 kg m^-2 in six months from September 2014.

• 5

  If the dust cover was homogeneous, with an average bulk density of about 500 kg m^-3, the dust-to-ice ratio ranges between 50 and 100, notably higher than 4( ± 2) as reported by Rotundi et al. (2015) derived from coma observations and indicating a lower limit (Fulle et al. 2016b). It cannot be ruled out that the ice abundance is inhomogeneous on spatial scales below a few tens of meters across the dust cover and the local dust-to-ice ratio is lower. Our study does not constrain the dust-to-ice ratio in a distinctly inhomogeneous dust cover that consists of significant inactive areas. For instance, a low dust-to-ice ratio (4, for example) and the infinite one are indistinguishable beneath a thick dust mantle that effectively quenches any water outgassing.

• 6

  The dust cover had been restored around perihelion, as predicted by Keller et al. (2015). Our results suggest that the dust cover contained a few percent of water ice in mass on average.

---

Acknowledgments

X. Hu acknowledges the financial support from the International Max Planck Research School for Solar System Science in Göttingen. We thank the referee, the Language Editor, and the Associate Editor for their comments and suggestions that helped improve the quality of the paper. OSIRIS was built by a consortium led by the Max-Planck- Institut für Sonnensystemforschung, Göttingen, Germany, in collaboration with CISAS, University of Padova, Italy, the Laboratoire d’Astrophysique de Marseille, France, the Instituto de Astrofísica de Andalucia, CSIC, Granada, Spain, the Scientific Support Office of the European Space Agency, Noordwijk, The Netherlands, the Instituto Nacional de Técnica Aeroespacial, Madrid, Spain, the Universidad Politéchnica de Madrid, Spain, the Department of Physics and Astronomy of Uppsala University, Sweden, and the Institut für Datentechnik und Kommunikationsnetze der Technischen Universität Braunschweig, Germany. The support of the national funding agencies of Germany (DLR), France (CNES), Italy (ASI), Spain (MEC), Sweden (SNSB), and the ESA Technical Directorate is gratefully acknowledged.

References

Appendix A: Catalog and reference of images

Information of all the OSIRIS images used in this work is listed in Table A.1.

Appendix B: Method of shape-from-shading using a single image

Denote the vector from the focal point of camera to the certain feature on the nucleus surface by , where represents a unit directional vector. The distance, r (in m), between the camera and the nucleus surface varies by, (B.1)for some differential step, ds, away from the “zero-phase” point within the image plane. The zero-phase point is defined as the pixel at which the phase angle would be zero, that is, where , with being the unit vector from the camera to the Sun. is the unit vector along the bore-sight of the camera. denotes the unit normal vector to the “Sun-surface-observer” plane. The zero-phase point lies along the intersection of a Sun-surface-observer plane with the (extended) image plane. r_F is the focal length of the camera. The photometric angle, γ, is as defined in Eq. (3), and tanγ is measured as (Rindfleisch 1966), (B.2)where cosi/ cose is, according to Eq. (2), derived from the radiance factor, ℜ, (B.3)We note that the sign of γ is opposite to that defined by Rindfleisch (1966; which explains the omission of a minus sign on the right-hand side of Eq. (B.1) compared to the original expression, Eq. (32), in the reference).

If the distance from the camera to a certain surface point, say P₀, is known as r₀, Eq. (B.1) can be integrated to find the distance to any other point, P, such as, . The paths of integration are straight lines in the image radiating from the zero-phase point along which is constant (Rindfleisch 1966).

Alternatively, Eq. (B.1) can be integrated directly to find, (B.4)such that r = r₀ + Δr, in case dr (or Δr) is small relative to r. In the present analysis, r₀ is roughly the distance between the spacecraft and the nucleus surface and is always greater than 10 km; whereas, Δr, as inferred from the shaded pixels, is on the order of decimeters and rarely exceeds a few meters for a typical indentation of honeycomb. Therefore, the use of Eq. (B.4) is legitimate. We note that the path of integration for (B.4) also originate from the zero-phase point in the image.

The analysis is further simplified by the following considerations:

• A.

  r₀ is regarded as constant within small surface areas (for example, of less than 100 m in width).

• B.

  With a detailed (but imperfect) shape model of the nucleus and with camera position and orientation extracted from SPICE, r₀ can be approximated as with small errors δr. The approximation of r₀ by in the integration of Eq. (B.4) is justified as long as δr is small compared with r₀.

• C.

  The apparent depths (in pixels) of indentations in honeycombs are insignificant compared with the expanse of an entire honeycomb feature. Therefore, for the sake of simplicity, we confine ourselves to the roughness pattern of the honeycombs in terms of Δr with respect to the plane parallel to the local nucleus surface.

We note that Δr is referred to the level indicated by that is arbitrarily close to the true but unknown level given by r₀ due to errors.

Appendix C: Thermal modeling

The 1D heat equation applies to a certain location on the nucleus, represented by a facet of the shape model, (C.1)where κ, c, and ρ denote the thermal conductivity, heat capacity, and density of the subsurface materials, respectively. At the surface, the energy input on each facet is solar irradiation, Q_⊙, and (unless otherwise noted) “self-heating”, Q_⊕, which accounts for the direct thermal radiation from nearby, unobstructed facets (Ivanova & Shulman 2006; Davidsson & Rickman 2014; Keller et al. 2015), (C.2)where σ is the Stefan-Boltzmann constant and ϵ ≈ 1 the emissivity of the surface. κ_d is referred to the dust mantle (Fig. C.1). At the bottom of the dust mantle, x = X_i, part of the energy or heat flux is consumed by sublimation of water ice (Kührt & Keller 1994; Skorov et al. 2011; Keller et al. 2015), (C.3)κ_i is the thermal conductivity of the icy dust, and ℓ the latent heat of water ice. The mass flux of sublimation of solid pure water ice below the dust mantle is given by (Gundlach et al. 2011), (C.4)where , or the Hertz-Knudsen formula, applies to the case of ice sublimating from the surface into vacuum. k_B is the Boltzmann constant and m_M the molecular weight of water. The saturation vapor pressure, P_V, and the sublimation coefficient, α, are temperature dependent and calculated using formulas determined by Gundlach et al. (2011). The dust mantle suppresses the sublimation flux by a factor of Ψ = (1 + 0.14X_i/d_P)^-1, where d_P is the diameter of the uniform-sized constituting dust paritcles (Gundlach et al. 2011). We introduce a factor, f, to account for the obvious fact that sublimation flux from the ice-dust mixture will be reduced with respect to pure ice, that is, (C.5)This factor is formally defined by Crifo (1997) as, (C.6)where ϱ_i, ϱ_d denote the densities of the (solid) ice and of the dust, respectively, and μ is the dust-to-ice mass ratio in the mixture. The gas diffusion through the pores of dust-ice mixture yields a contribution to the total water production (Davidsson & Skorov 2002), which is neglected here. It is readily evident that the water production of a dusty subsurface with μ ≫ 1 is reduced compared to an icy subsurface with μ = 0.

Numerical treatment: integration and interpolation.

We apply the Crank-Nicolson method for solving the temperature profile and the sublimation flux for each facet. Given certain boundary conditions, the diurnal or orbital variations of temperatures are solved iteratively until results converge, when the temperature profile coincides with that before (or after) one period of nucleus rotation or orbital motion. The time frame in which honeycombs appeared is narrow compared with the orbital period of more than six years for 67P. It is thus justified to derive the accumulated erosion by propagating or integrating diurnal solutions over a few months.

For diurnal solutions, we impose the isothermal condition, ∂T/∂t = 0, at the lower boundary, X_T. We note that we do not employ the adiabatic condition of vanishing temperature gradient here, as there may exist (nearly) constant heat flux (inwards) at X_T. A full orbital solution is performed to find an approximate value for X_T. We refer the local shape model to a global shape model of reduced resolution (1000 facets) and select a corresponding facet of the latter model. The comet is assumed to have a homogeneous initial temperature of 15 K at aphelion. In this case, the adiabatic condition, ∂T/∂x = 0, is imposed at a depth of 20 m (though actually fulfilled at shallower depths), and self-heating is neglected for the upper boundary at the surface. X_T is solved iteratively as the depth at which T varies by less than 1 K over the entire orbit. Fewer than 15 iterations are needed for orbital solutions to converge. We used T = 100 K at X_T = 1 m for the subsequent diurnal solutions. X_T is greater than the orbital skin depth, approximated by κI^-1(P/π)^{1 / 2} with P being the orbital period of 67P (Huebner et al. 2006), that characterizes the attenuation depth of temperature variation by e^-1 with respect to the surface.

To reduce calculation time, we estimated diurnal erosions every ten (terrestrial) days and interpolate the estimation in between.

Appendix D: Supplementary images

All Tables

All Figures