© ESO 2020

1 Introduction

The European Space Agency’s Rosetta mission was designed to orbit and land on the Jupiter-family comet 67P/Churyumov-Gerasimenko (hereafter 67P). Rosetta arrived at its primary target on 6 August, 2014. The probe was first guided into an orbit around its target to perform a first analysis and find a suitable landing site for the lander module Philae. Rosetta revealed that the nucleus of comet 67P consists of two lobes (Sierks et al. 2015) connected by a narrow neck, with a stable spin axis. The Optical, Spectroscopic, and Infrared Remote Imaging System (OSIRIS, Keller et al. 2007) was designed to study the nucleus and its dust and gas environment. The system consisted of two cameras operating from near-ultraviolet to near-infrared wavelengths. The Wide Angle Camera (WAC) imaged the dust and the gas surrounding the nucleus with a spatial scale of 10.1 m px⁻¹ at 100 km from the surface. The comet nucleus and its surface topography were investigated by the Narrow Angle Camera (NAC) with a spatial scale of 1.86 m px⁻¹ at the same distance. The OSIRIS observations revealed that the northern regions, such as Ash, Ma’at, Seth, and Hapi (Thomas et al. 2015; El-Maarry et al. 2015) are fully covered by dust. On the contrary, equatorial regions, such as Anubis, Aket, and Bastet, look different. In these regions, we observe consolidated and coarse terrain instead of dust deposits (Keller et al. 2017). The nucleus rotates with an obliquity of 52° (Keller et al. 2007). Due to the inclination of the axis of rotation, the comet experiences strong seasonal effects, resulting in significant differences in insolation between the northern and southern hemispheres. This strong dichotomy is reflected in the morphology of the two hemispheres. Southern summer coincides with the perihelion passage, hence causing the erosion in the southern hemisphere to be much stronger than in the northern regions (Jorda et al. 2016; Keller et al. 2017). The approach to perihelion causes a rise in the temperature of the nucleus, sublimating the ices. Observations of the coma revealed water to be the most abundant volatile (Gulkis et al. 2015). Keller et al. (2017) calculated the erosion due to sublimation of water ice to investigate the link between insolation, erosion, and water content of the nucleus surface. The value was found to be four times stronger on the southern hemisphere than on the northern one. The strong insolation and the water ice content in the south could erode the surface up to 20 m (Keller et al. 2017) at perihelion. Instead, the northern hemisphere, particularly Hapi, is characterized by a minimal amount of insolation, and therefore minimal erosion. These pieces of evidence confirmed that the dichotomy in appearance between the two hemispheres is linked to the dichotomy in erosion, and that the dust cover in the northern regions could be the result of transport mechanisms of particles from the southern hemisphere during the southern summer (Keller et al. 2015).

In this study, we describe a tool we developed to quantify the seasonal erosion and deposit and/or accretion in the Hapi region. We started from this region to link the time evolution of dust with the mass transfer mechanism and the erosion of the comet surface. The tool is based on the monitoring of the time evolution of boulder height, which is defined as the difference between the top of the boulders and the surrounding pebble deposit surface. This technique has led to measurements of the seasonal evolution of the deposit erosion and/or accretion of the Hapi region with a vertical accuracy of 0.2–0.9 m. The amount of erosion of the southern hemisphere, the subsequent transport of material, and then its fallout on the nucleus allow us to investigate the pristine water ice abundance in comet 67P.

2 Mass transfer on 67P

The dichotomy in appearance between the two hemispheres of comet 67P is linked to the dichotomy in erosion. The eroded southern surface is subject to a strong insolation and water sublimation (Keller et al. 2015), contributing to the release of dust particles of different sizes. OSIRIS detected both small dust particles in the size range from 3 × 10⁻³ up to 1 cm (Fulle et al. 2016) and larger particles, named chunks. These objects are defined as pieces of the nucleus of an average mass of 1 kg, in the 10 to 20 cm range (Fulle et al. 2019). The dust particles size–frequency distribution (SFD) allows us to link the ejection mechanism with the water content of 67P. Fulle et al. (2019) calculated the chunk volume ejected by 67P from 24 July, 2015, to 15 September, 2015. The lost volume is about 4 × 10⁷ m³, corresponding to an eroded southern surface of about 10 km² (Keller et al. 2015; Blum et al. 2017). This corresponds to an average erosion thickness of about 4 m (Fulle et al. 2019). Assuming that the southern erosion occurs because of the ejection of chunks implies a total erosion in average steps of about 13 cm (Fulle et al. 2019), involving a nucleus surface of 65 m² every second. The erosion model developed by Keller et al. (2015) suggests that the maximum water loss rate per unit area is 3 × 10⁻⁴ kg m⁻² s⁻¹. This implies that from a nucleus surface of 65 m², the water loss rate is at most 0.02 kg s⁻¹, which is a negligible mass fraction of the corresponding chunk loss rate (Q_v = 8.3 ± 2.1 m³ s⁻¹ ; Fulle et al. 2019; namely the mass loss rate divided by the bulk density). This result suggests that the chunk ejection from the nucleus surface is dominated by perihelion erosion of the southern hemisphere, but behaves independently of water ejection. As mentioned above, the dichotomy in appearance between the two hemispheres is linked to the dichotomy in erosion. The dust cover in the northern regions is the result of transport mechanisms of particles from the southern hemisphere during the southern summer (Keller et al. 2015). Chunks ejected at perihelion fall back over the whole nucleus, including Hapi. OSIRIS observations of the nucleus reveal a surface characterized by a varied surface granularity. Pebbles of ≈ 25 cm in size have been observed in the Sais region (Pajola et al. 2017), suggesting a deposit built up by chunks, and confirming the chunk mass distribution in the 67P coma (Fulle et al. 2016; Ott et al. 2017). The distribution of pebbles in different regions, including Hapi, can be explained with the dust fallout mechanism, which causes chunks ejected during perihelion to fall back over the whole nucleus. As the outbound equinox approaches, the southern erosion decreases, and the outgassing in the northern hemisphere self-cleans the fallout, removing the dust and leaving chunks, because the nucleus outgassing is too low to lift these objects.

The refractory-to-ice mass ratio

Water content and delivery to the terrestrial planets is still a subject of debate. Water is present in different bodies of the Solar System, even in the outer asteroid belt and beyond in the form of ices. Primitive meteorites such as CI-chondrites are believed to come from C-type asteroids that dominated the outer part of the asteroid belt (Burbine et al. 2002), and can have ~5–20% water by mass (O’Brien et al. 2018). Growth processes of carbonaceous CI-chondrites, characterized by a chemical composition mostly resembling the solar photosphere, allow water to be trapped in their silicates at a molecular scale (Garenne et al. 2014). Ordinary chondrites are linked to S-type asteroids, which are present in the inner asteroid belts, and contain a different amount of water. As regards comets, they are formed beyond Neptune and their water was incorporated in the form of ices embedded in the refractory matrix of their nuclei (Blum et al. 2017). The radial distribution of water and ices is recorded in the refractory-to-ice mass ratio in comets (Fulle et al. 2019). This ratio is a fundamental parameter which constrains the origin of comets and Kuiper Belt objects (KBOs). Quantifying the difference in water content between comets and other bodies allows us to distinguish between different formation processes, providing knowledge on the origin of the Solar System.

As mentioned above, chunks ejected at perihelion fall back over the whole nucleus, including Hapi. As the outbound equinox approaches, the outgassing from the southern hemisphere decreases due to seasonal changes, increasing the outgassing in the northern hemisphere where fallout occurred around perihelion (Fulle et al. 2019). Outbound, the self-cleaning in the Hapi region is negligible, and possible only from 2.5 to 4 au outbound (Fulle et al. 2019). This implies that Hapi ejects sub-centimetre dust only (Rotundi et al. 2015), acting as a chunk deposit with a thickness of metres. This can be explained as follows. Chunks ejected at perihelion have a refractory-to-ice mass ratio larger than inside the nucleus, and have an upper exposed dehydrated crust, the thickness of which increases as the refractory-to-ice mass ratio increases (Fulle et al. 2019). Fresh ice on Hapi is exposed to sunlight by water ice migration to the surface (De Sanctis et al. 2015) and by the removal of chunks of dehydrated crust. This causes outgassing coming from the interior of the chunks (Fulle et al. 2019), preventing any outgassing from below the surface.

Keller et al. (2017) calculated the water production of the nucleus along its orbit. These latter authors found that the production rates of the northern regions follow the insolation trend, and are controlled by the peculiar shape of the nucleus and the inclination of the spin axis. Hapi represents the conjunction point between the two lobes. This cavity reached insolation for short intervals of a cometary day because of the shading by the lobes (Pajola et al. 2019). This should mean that this region cannot be as active as other regions located in the south. However, Hapi appears to be the most active area during northern summer as a consequence of its water content and morphology. Because of the morphology of this region, the absorbed energy is not sufficient to produce strong outbursts, as observed elsewhere on the nucleus, but allows a modest continuous activity also far from perihelion. The activity is due to a thinner desiccated dust layer which accumulates after perihelion because of the fall out, and is eroded by ice sublimation when approaching the sun again after aphelion (Keller et al. 2017). In this context, the seasons of Hapi are fundamental to link (Fulle et al. 2019) its processed ice abundance to the pristine ice content of the metres-thick layers eroded every perihelion from the southern nucleus hemisphere.

3 Data and method

The Hapi region is located in the northern hemisphere (Thomas et al. 2015), between the two lobes of comet 67P. It has been classified as a deposit of debris (Thomas et al. 2015), and is characterized by a distribution of boulders and/or outcrops (Cambianica et al. 2019) scattered all over the debris plain. Hapi is also dominated by smooth terrain (El-Maarry et al. 2015 and dune fields. Due to the presence of boulders and outcrops of tens of metres in size, El-Maarry et al. (2015) suggested a dust deposit thickness of several metres. The neck region corresponds to the gravitational minimum of the nucleus (Keller et al. 2017) and is therefore considered as the preferred location for the accumulation of the back-falling material (Keller et al. 2017). To measure the erosion and deposition of dust in this region, we monitored changes in the height of boulders, assuming that these objects are not involved in erosion processes, and that it is the surrounding terrain that gains height or is eroded.

3.1 Data selection

For the analysis of dust erosion and deposit on the Hapi region we used 19 OSIRIS NAC images (see Table 1 for imageIDs) and the photogrammetric SHAP8 V.2.1 (issued by Gaskell and Jorda in March 2018) comet shape model, which represents an advanced model based on a data set of 20679 OSIRIS-NAC images and 6072 OSIRIS-WAC images acquired between 11 July, 2014, and 30 September, 2016. The shape model can therefore be considered as a mean model of the surface, in which the evolution of the surface morphology cannot be appreciated. We divided the data set into two groups. The first set of 16 NAC images was acquired from 21 August, 2014, to 10 December, 2014, which is before the comet inbound equinox. The spatial scale ranges from 1.26 to 0.37 m px⁻¹. The second set of three NAC images was acquired from 19 June, 2016, to 30 September, 2016 (spatial scale of 0.56, 0.54, and 0.26 m px⁻¹, respectively) after the outbound equinox. We developed a MATLAB (MATLAB 2010) tool to monitor the time evolution of boulder heights from 21 August, 2014, up to the end of the mission. The height is defined as the difference between the top of the boulders and the surrounding pebble deposit surfaces. The tool is based on the measurement of the shadow length projected by the boulder on the surrounding pebbles deposit. This technique relies on the geometric Spice Kernel data (Acton 1996) and requires high-resolution images (at least 1.30 m px⁻¹). High-resolution NAC images provided global views of the Hapi region, allowing us to analyze as many boulders as possible. According to the illumination and visibility conditions, we measured the height of 22 boulders. The locations of the boulders is shown in Fig. 1. The figure reports the corresponding ID for each boulder (see Table 2 for the ID of the boulders and their latitudes and longitudes).

3.2 Surface plane definition and image alignment

To measure the height H of a boulder, we have to consider the projection of the OSIRIS NAC image on the 3D shape model of the comet. To obtain the correct projection, we defined a set of uniquely identified tie points both on the shape model and on the images. From the correspondence between these points, we derived the proper homography transformation for each projective system. This method allows us to refine the instrument attitude minimizing the residuals of the 3D points projected on the images. As a result, we can associate a corresponding 3D point on the surface photogrammetric model to each pixel of the images. Once the misalignment between thetwo is corrected, and the proper boulder is identified, it is possible to define an average surface plane Σ around the boulder in which the boulder shadow is projected (see Fig. 2). The definition of an average surface plane is fundamental to avoid the local granularity of the mesh and to smooth possible surface irregularities, since the topography of Hapi is nonhomogeneous. The definition of the average surface is defined by a grid of 3D points around the contour of the boulder and a suitable interpolation of these points. The resulting point cloud is then decomposed in eigenvectors. The eigenvector associated to the minimum eigenvalue represents the normal to the surface. Other parameters to be defined are the incidence angle i of the solar direction with respect to the normal to the plane Σ, and the length L of the boulder shadow. These three elements completely define the adopted geometry for determining the boulder height.

Fig. 1 Zoom of the image from 21 August, 2014, (NAC_2014-08-21T16.42.56.549Z_ID30_1397549300_F22). The figure shows the corresponding ID for each boulder reported in Table 2.

Fig. 2 Schematic representation of the adopted geometry and of the parameters used for determining the boulder height. Σ is the average surface plane, Π is the plane containing the normal n to Σ and the illumination vector, i is the incidence angle, Q refers to the peak of the shadow, which is defined by the highest point of the boulder P, and L is the length of the shadow. Finally, H is the heightof the boulder.

Fig. 3 Plot a: NAC view of the Hapi region (0.37 m px−1); this image was acquired in 2014. Plot b: close-up of a boulder and of its shadow. The green line represents the projection of the Sun illumination direction. Plot c: boulder section (we note the different scales on the plot axes deforming the boulder shape). The y-axis is oriented as the normal to the average plane around the peak of the shadow. The x-axis is obtained by projecting the green line in plot b on this average plane; the x-axis origin coincides with the peak of the shadow. We note that by comparing plots b and c, the shadow lengths appear different because of projection effects.

3.3 Boulder profile reconstruction and height calculation

As mentioned before, the 3D model cannot be used to directly measure the height of a given because it does not take into account variations inthe surface morphology. For this reason, we adopted a technique based on the assumption that the surface normal to Σ is considered to be locally time-invariant and based on the definition of the illumination vector, which uniquely identifies the position of the peak of the shadow Q (see Fig. 2) on the average surface Σ, and therefore on the 3D model. The peak is defined as the point of the shadow contour with the longest extension along the illumination direction, and originates from the highest point of the boulder P. Following this, it is necessary to define a plane Π perpendicular to Σ which contains the direction of the illumination direction, and passes through the point Q. This plane cuts the boulder along the direction of the solar illumination and passes through its highest point P. The bi-dimensional boulder profile obtained by sectioning the comet shape model with the plane Π finally allows us to determine the position of the peak P, which is the highest boulder point tangent to the illumination direction. The length L of the shadow on the Σ plane can then be defined as the distance between the Q point and the projection of P on Σ. Finally, the boulder height H can be calculated as follows: $$H=L\cdot \text{tan}\left(\frac{\pi }{2}-i\right),$$(1)

where L is the length of the shadow and i is the incidence angle. In Fig. 3 an example of the height measurement is shown. After the identification of the proper boulder and definition of the average surface plane, the Sun illumination direction is calculated (see the green line in Fig. 3, plot b). After the manual selection of the peak of the shadow (which corresponds to the origin of the x-axis in plot c of Fig. 3) and the projection of the illumination vector on the average plane, the height of the boulder is calculated.

Fig. 4 Gravitational slopes of the 67P nucleus surface. The centrifugal force is included. The values are restricted to the 0°–60° ranges to emphasize the slopes below the repose angle of loose materials (Groussin et al. 2015). The green ellipse defines the surface of 0.2 km2 encircling all the boulders considered in this work.

3.4 Gravitational slope

To confirm that our measurement has not been performed on a deposit the thickness of which is altered by local gravitational instability processes, for example landslides and granular flows, or on a local point of fallout accumulation on the comet surface, we investigated the surface gravitational slope. The gravitational slope is defined as the angle between the local surface normal and the vector opposite to the estimated acceleration field (Penasa et al. 2017). As shown in Fig. 4, the whole Hapi region has gravitational slopes ranging between 0° and 20°, lower than the angle of repose of loose granular materials on 67P of (45 ± 5)° (Groussin et al. 2015). This is consistent with the assumption that Hapi deposits are the result of a homogeneous deposition driven by fallout.

3.5 Method validation

To validate the method, we performed a comparison test. El-Maarry et al. (2017) measured the height of a boulder located in the Imhotep region (see Fig. 5 for the location of the region and the considered boulder) finding a value of 3.9 +0.1/−0.2 m in height. We applied our method to the same boulder and we found 4.08 ± 0.35 m, which is a value consistent with the determination of El-Maarry et al. (2017). The two methods being completely independent, this check confirmed the reliability of our method.

Fig. 5 OSIRIS-NAC image taken on 25 December, 2016 (NAC_2016-05-25T15.32.54.769Z_1397549100_F22). The white circle indicates the analyzed boulder. Bottom right panel: location of the Imhotep region on the comet nucleus.

4 Results

According to Eq. (1), we calculated the incidence angle (i) and the length of the shadow (L) to measure the height of 22 boulders. A summary of the boulder height measurements is reported in Table A.1. We provide the corresponding universal time coordinated (UTC), the calculated incidence (i) and emission (e) angles, the measured length of the shadow (L), and the height (H) of boulders with the associated average error bar (δH). As shown in Eq. (1), the height calculation does not depend on the emission angle, which is defined as the angle of camera bore-sight relative to the surface normal. However, we report both the emission and incidence angle values to show the statistical variability of the images used. To mediate the error bars, we used OSIRIS NAC high-resolution images with the best visibility conditions to compare as many heights as possible. We measured 22 boulders, but some measurements have not been included because of their large uncertainty due to the adverse illumination or visibility conditions. The reported error bars have been estimated propagating the individual errors associated to the selection of the pixel identified as the shadow peak, and the calculation of the incidence angle. The accuracy of the manual selection of the pixel depends on the ability of the operator to select the proper pixel. The accuracy of the incidence angle calculation strongly depends on the definition of the normal to the Σ plane. This is associated to the standard deviation of the distance between the surface points around the boulder and the Σ plane. To estimate the accuracy of the Σ definition, and therefore the accuracy of the incidence angle calculation, we investigated the impact of the granularity on the incidence angle. We performed a Monte-Carlo simulation which consists in defining a set of surfaces with different vertical standard deviations. These values represent the granularity of the surface. The estimation of the normal and the calculation of the incidence angle are performed by calculating the previously defined eigenvectors (see Sect. 3.2). Figure 6 shows the results of the Monte-Carlo simulation. The incidence angle being related to the definition of the above-mentioned normal, the derived standard deviations for the measured surfaces range between 0.02 and 0.1 m. As seen in Fig. 6, in this range the maximum error on the incidence angle value is lower than 0.5°. After measuring the height of each boulder, we calculated the corresponding maximum and minimum values considering the associated errors (see Table A.1). These values have been obtained as follows: $$H_{\text{min/max}}=[L\mp \delta p(\text{POG}\mp \delta l)]\text{tan}\left[\frac{\pi }{2}-(i\pm \delta i)\right],$$(2)

where L is the calculated length of the shadow, δp is the error due to the manual selection of the peak of the shadow, the POG is the pixel on ground (target–spacecraft distance multiplied by the pitch size and divided by the focal length), δl is an error related to the uncertainty of the Digital Terrain Model (δl = δh ⋅tan(e), where δh is the average uncertainty of the location of the DTM with respect to the zero level of the image, and e is the emission angle), i is the incidence angle, and the error on the incidence angle, δi, is estimated to be equal to 0.5° (see Fig. 6).

To better analyze the data, the images were divided into subsets as follows:

• Subset 1: 21–22 August, 2014;

• Subset 2: 28 August, 2014–1 September, 2014;

• Subset 3: 10–22 September, 2014;

• Subset 4: 10 December, 2014;

• Subset 5: 19 June, 2016;

• Subset 6: 30 September, 2016.

This allowed us to identify differences in height (ΔH^*) between a generic image and the height measured on 21 August, 2014. Having used 19 images, per subset there is more than one image for which the boulder heights can be measured. Therefore, we calculated the weighted average height of each boulder ⟨H⟩. The associated error δ(⟨H⟩) is calculatedas follows: $$\delta (\langle H\rangle )=\frac{1}{\sqrt{\Sigma {(1/\delta H)}^2}},$$(3)

where δH is the largest absolute value difference between the measured height and H_Min/Max.

The reported ΔH values are the differences between the < H> of two subsets, ΔH = <H_i> − <H_j>. The associated error is finally calculated as: $$\delta \Delta H=\sqrt{\delta {(<H_i>)}^2+\delta {(<H_j>)}^2}.$$(4)

In Tables 3 and 4, the time evolution of the weighted average height of each boulder and the time evolution of the boulder height difference ΔH are shown. The graphical representation of these results is shown in Fig. A.1.

The formula we used to calculate the weighted mean ΔH_weighted of the data sets ΔH_i with corresponding error bars δΔH_i shown in Table 4 is: $$\Delta H_{\text{weighted}}=\frac{{\sum }_{i=1}^n(w_i\cdot \Delta H_i)}{{\sum }_{i=1}^n{w}_i},$$(5)

where $w_i=1/{(\delta \Delta H_i)}^2$. The standard error of the weighted mean is given by $$\delta \Delta H_{\text{weighted}}={\left(\sqrt{{\sum }_{i=1}^n{w}_i}\right)}^{-1}.$$(6)

Fig. 6 Incidence angle error as a function of the surface standard deviation, which refers to the surface granularity.

5 Discussion

We present observations by the Rosetta mission of comet 67P which have allowed us to quantify the seasonal evolution and deposit and/or accretion in the Hapi region. Figure 7 shows the time evolution of the heights of boulders on the surface of the comet during the observations. Here the difference ΔH^* between the height measurement of a generic image and the height measured on 21 August, 2014, is shown as a function of time. We find a systematic increase of ΔH of about 2 m over 2014 and then a decrease by a similar amount between the equinoxes followed by a plateau up to the end of the observations. The assumption here is that any height variation is not due to boulder intrinsic changes but to erosion or accretion of the surrounding pebbles and dust deposits; due to the dominant cross-section of the fallout chunks (from 80 up to 350 cm² ; Fulle et al. 2019), only a minimal fraction of the fallout remains settled on the top of the boulders. Furthermore, following the method developed by Cambianica et al. (2019), we analyzed the location and the morphology of these boulders and outcrops, and no variation in morphology or position of the boulders was seen during the mission. The erosion of the southern hemisphere, the subsequent transport of material, and then its fallout on the nucleus (Fulle et al. 2019; Keller et al. 2017) are responsible for the time evolution of the height of the boulders, which is linked to the likely source of the deposit surrounding them. This mainly happens at the perihelion phase (Fulle et al. 2019; Keller et al. 2017), as evidenced by the similarity between the size distributions of large particles suspended in the coma (Fulle et al. 2016; Ott et al. 2017) and the cobbles in Sais and Agilkia regions (Pajola et al. 2017). The measured decrease of the boulder heights at the perihelion phase from 2014 to 2016, δΔH_2014−2016 = −1.4 ± 0.4 m, is therefore interpreted as a direct measurement of the fallout thickness. This value matches the predicted (Fulle et al. 2019) fallout thickness of 1.8 ± 0.7 m, revealing a consistency between the perihelion erosion and the total loss of nucleus mass per orbit (Pätzold et al. 2018). The height difference between the boulders and the debris field steadily increases during the inbound orbit, indicating an ongoing erosion of the deposits in the Hapi region before the comet reached its spring equinox. As shown in Fig. 7, between June and September, 2016, ΔH variations are lower than their error bars, consistent with negligible fallout and erosion after the outbound equinox. This is consistent with the process that has been suggested to be the main driver of the fallout (Fulle et al. 2019; Bertini et al. 2018), namely the expelled chunks in bound orbits are slowed by friction with the coma gas and finally collapse onto the nucleus. After the outbound equinox, the coma of 67P is too tenuous to further affect the chunk motion.

5.1 Erosion rate

The negligible erosion of the Hapi region from June to September, 2016, also suggests that Hapi maintains all its ice during the outbound orbit, up to the next inbound activity (Fulle et al. 2019; Keller et al. 2017). Any possible ice regression without dust ejection would soon build up an insulating crust of thickness of a few centimetres at most (sufficient to dump any further ice sublimation), thinner than the chunk size. Knowing the amount of eroded material (ΔH), we calculated the average erosion rates. This value can be calculated as follows: $${\text{Rate}}_{\text{Erosion}}=\frac{\Delta H\times A}{\delta t},$$(7)

where A is the area of the considered elliptical surface encircling all boulders and δt is the time interval under consideration. The average erosion rate decreases from 0.15 ± 0.07 m³ s⁻¹ in August, 2014, to 0.06 ± 0.03 m³ /s in September, 2014, and to 0.012 ± 0.010 m³ s⁻¹ in October-December, 2014. Since the adopted surface of 0.2 km² is only 10% of the surface of the Hapi region, and there is no reason for the deposits in this latter region to only be eroded around the boulders, it is reasonable to assume that these are just lower limit values of the global erosion rates.

Fig. 7 Time evolution of the boulder height differences ΔH* (as defined in text). Each time step contains information related to 22 boulders and the values correspond to the weighted average. They-axis reports the average height measurements referring to the subsets described in Sect. 4. Subset 1 refers to 21–22 August 2014. According to the illumination and visibility conditions, the time interval between Subset 1 and the perihelion is about 4 months. The time interval between the perihelion and Subset 6 is about one year. Erosion during the inbound orbit until December, 2014, nearly balances the fallout from the southern hemisphere during perihelion cometary activity. The dotted red line indicates the perihelion.

5.2 Dust fallout

The illumination and temperature of the Hapi region did not change significantly from August to October 2014 (Tosi et al. 2019), which is when we measure surface erosion of the order of metres. Therefore, a certain physical evolution of the terrain must be assumed to explain this phenomenon, and the most reasonable mechanism is the fallout self-cleaning (Fulle et al. 2019; Pajola et al. 2017). Chunks falling onto the nucleus surface are rich in ice as they have been freshly expelled from the active southern hemisphere, and, if irradiated, they outgas and self-clean. In the Hapi region, this outgassing and self-cleaning activity is absent during the fallout season, when Hapi is in winter time, and is almost negligible after outbound equinox because Hapi is in a rather shadowed location. However, outgassing and self-cleaning become more intense before reaching inbound equinox, which is exactly when we observe the erosion phenomenon. Subsequently, as the comet gets closer to the Sun and Hapi approaches its winter again, outgassing decreases in the Hapi region, is no longer able to remove the chunk dehydrated crust, and the erosion rate decreases and the process stops. At the same time, the activity of the southern hemisphere increases again, as do the total gas- and dust-loss rates (Fougere et al. 2016), and the cycle is repeated.

It is also possible to directly estimate the fallout amount on the comet nucleus associated to the activity of the neck region in August 2014. In fact, by comparing the erosion rate in the Hapi region in that month (0.15 ± 0.07 m³ s⁻¹) with the dust volume loss rate measured in the coma in the same period (0.006 m³ s⁻¹; Migliorini et al. 2016), we can see that the former is 25 times larger than the latter. Since Hapi was contributingfor the most part to the comet outgassing in August 2014, we can see that the volume loss rate measured in the coma is only 4% of the total erosion. This implies that the remaining 96% in volume of material eroded from Hapi is falling again on the comet nucleus surface. Fulle et al. (2018) estimated that the fallout at perihelion is about 80% of the total southern eroded volume. The larger fallout value measured here is probably due to the peculiar structure of the neck, which is a region surrounded by steep walls: this makes it more difficult for the particles to escape the comet because of the high probability of collision with the walls, and increases the fallout percentage (Shi et al. 2018).

5.3 Water ice fraction in the Hapi region

As mentioned before, the measurement of the deposit erosion and/or accretion in the Hapi region allows us to investigate the pristine water ice abundance in comet 67P. The power index of the differential dust size distribution at dust sizes >1 mm has been calculated to be equal to −4 (Rotundi et al. 2015), which, compared to that of the fallout from Hapi, constrains its composition in particles smaller than 1 cm, meaning that its bulk density is ${\rho }_b={800}_{-100}^{+500}$ kg m⁻³ measured insub-millimetre particles (Fulle et al. 2017).

The bulk density allows us to calculate the dust mass loss rate Q_m in August, 2014, from the considered area A, as follows: $$Q_{\text{m}}={\rho }_b\times {\text{Rate}}_{\text{erosion}}.$$(8)

The dust mass loss rate becomes 120${}_{-60}^{+160}$ kg/s. The ratio between the erosion rate and the corresponding water vapour loss rate of 1.2 kg s⁻¹ (Gulkis et al. 2015) provides Hapi’s dust-to-water mass ratio (Rotundi et al. 2015) at the erosion of 100${}_{-50}^{+140}$. Coming from the crust of the chunks in the deposit (Fulle et al. 2019), the dust ejected by Hapi is dry (Fulle et al. 2018). This is supported by the match between the water loss rates provided by local (ROSINA) and remote (MIRO) observations (Marshall et al. 2017) in August, 2014, which are <1% in mass of Hapi’s measured erosion rate. This fact allows us to infer the ice mass fraction of the Hapi region from the measured dust-to-water ratio. The inverse of the dust-to-water ratio provides Hapi’s water ice fraction of (1.2 ± 0.8)% in mass. The lower limit of the dust-to-water ratio corresponds to a water ice fraction of 2% in mass. Since the dust ejected from the surface of the Hapi region was probably greater than 0.2 km² in total volume, this value can be considered as an upper limit (a larger area means an increase of the dust volume). The fallout from Hapi is inert, but Hapi’s erosion of 1.7 ± 0.2 m during the inbound orbit is statistically diluted by a factor of about 250 (namely, the total nucleus surface divided by 0.2 km²), providing an average dry fallout of <1 cm thick over the whole nucleus. This layer is negligibly thin with respect to the total southern erosion of at least 4 m (Fulle et al. 2019), with about 97% in volume of the total ejected material in chunks of sizes >1 cm according tothe perihelion dust size distribution (Fulle et al. 2016).

6 Conclusions

In this study, we measured the seasonal evolution of the deposit erosion and/or accretion in the Hapi region of the comet 67P with a vertical accuracy of 0.2–0.9 m, quantifying the mass transfer mechanism from the southern to the northern hemisphere of the comet. To this aim, we developed a tool based on the monitoring of the time evolution of 22 boulders located in the neck of the comet, a region named Hapi. This region is located in the northern hemisphere, and represents an ideal region for applying the method. This region is considered to be the preferred location for the accumulation of material coming from the southern hemisphere (Keller et al. 2017). We find that erosion during the inbound orbit until December, 2014, nearly balances the fallout from the southern hemisphere during perihelion cometary activity. A comparison between the eroded material and the dust volume loss rate measured in the coma provided the amount of fall back material due to the morphology of the Hapi region. The fallout represents 96% of the eroded volume and is consistent with the model (Fulle et al. 2019) linking the metres-thick southern erosion of pristine nucleus material to the northern fallout. Using Eqs. (2) and (3) in Fulle et al. (2019), it is possible to estimate the pure ice and pure refractory mass ejection rates ejected in the chunks. The water ice fraction in Hapi’s deposit of (1.2 ± 0.8)% in mass provides a refractory-to-ice mass ratio ranging from 6 to 110 in the 4 × 10⁷ m³ volume of pristine nucleus material eroded at perihelion, corresponding to a pristine ice mass fraction of (8 ± 7)% in mass. The refractory-to-ice mass ratio of the eroded pristine material can be compared with the same ratio measured in CI-Chondrites and in the interstellar medium (ISM). This value is in the range of 5% (Mogi et al. 2017) to 20% (Garenne et al. 2014) measured in CI-chondrites and in the ISM (about 20%). The molecular abundance of water ice in molecular clouds (Boogert et al. 2015), and likely in the outer protoplanetary discs, H₂ O/H ≈ 10⁻⁴, and the hydrogen-to-refractory mass ratio of approximately 100 in the ISM (Spitzer 2008) imply a refractory-to-water mass ratio ≈10⁴/(18 × 100) ≈5 in the discs beyond the snow line. The differences we find in terms of water abundance of the ices incorporated into the comet 67P could be caused by one of two different scenarios. The first hypothesis is that the water abundance of the ices incorporated into comet 67P were lower than the ISM value. This is supported by observations indicating a ≤10⁻⁴ abundance of the water sublimated from ices in the hot corinos of Solar-type protostars (Ceccarelli et al. 2000; Visser et al. 2013). Another hypothesis is that comet 67P lost some water in its formation. These results imply water-trapping mechanisms that are more efficient in possible asteroidal chondritic parents than in comets (Lorek et al. 2016), a negligible water loss by the catastrophic collisions fragmenting asteroids into chondrites (unless CI-chondrites after their formation were enriched with water to values higher than the ISM average), and a relatively uniform radial distribution of water ice in the protoplanetary disc beyond the snow line, consistent with the idea that significant radial mixing of the disc explains the minerals found in comets (Fulle et al. 2017; Ogliore et al. 2009).

Acknowledgements

We thank the anonymous referee for having significantly improved the manuscript. OSIRIS was built by a consortium of the Max-Planck Institut für Sonnensystemforschüng, in Güttingen, Germany, CISAS University of Padova, Italy, the Laboratoire de Astrophysique de Marseille, France, the Instituto de Astrofísica de Andalucía, CSIC, Granada, Spain, the Research and Scientific Support Department of the European Space Agency, Noordwijk, The Netherlands, the Instituto Nacional de Tecnica Aeroespacial, Madrid, Spain, the Universidad Politechnica de Madrid, Spain, the Department of Physics and Astronomy of Uppsala University, Sweden, and the Institut für Datentechnik und Kommunikationsnetze der Technischen Universitat Braunschweig, Germany. The support of the national funding agencies of Germany (DLR), France (CNES), Italy (ASI), Spain (MEC), Sweden (SNSB), and the ESA Technical Directorate is gratefully acknowledged. We thank the ESA teams at ESAC, ESOC and ESTEC for their work in support of the Rosetta mission. We made use of Arcgis 10.3.1 software together with the Matlab, Java, and ImageJ software to perform our analysis. I.T. acknowledges the support from project GINOP-2.3.2-15-2016-00003 “Cosmic effects and hazards”.

Appendix A Boulder height measurements

Fig. A.1 Time evolution of the boulder height. The histogram shows the average height measurements for each boulder. As shown inthe legend, the results are divided into five time periods to better represent the trend.

References

All Tables

All Figures

	Fig. 1 Zoom of the image from 21 August, 2014, (NAC_2014-08-21T16.42.56.549Z_ID30_1397549300_F22). The figure shows the corresponding ID for each boulder reported in Table 2.
In the text

Fig. 2 Schematic representation of the adopted geometry and of the parameters used for determining the boulder height. Σ is the average surface plane, Π is the plane containing the normal n to Σ and the illumination vector, i is the incidence angle, Q refers to the peak of the shadow, which is defined by the highest point of the boulder P, and L is the length of the shadow. Finally, H is the heightof the boulder.

In the text

Fig. 3 Plot a: NAC view of the Hapi region (0.37 m px−1); this image was acquired in 2014. Plot b: close-up of a boulder and of its shadow. The green line represents the projection of the Sun illumination direction. Plot c: boulder section (we note the different scales on the plot axes deforming the boulder shape). The y-axis is oriented as the normal to the average plane around the peak of the shadow. The x-axis is obtained by projecting the green line in plot b on this average plane; the x-axis origin coincides with the peak of the shadow. We note that by comparing plots b and c, the shadow lengths appear different because of projection effects.

In the text

	Fig. 4 Gravitational slopes of the 67P nucleus surface. The centrifugal force is included. The values are restricted to the 0°–60° ranges to emphasize the slopes below the repose angle of loose materials (Groussin et al. 2015). The green ellipse defines the surface of 0.2 km2 encircling all the boulders considered in this work.
In the text

	Fig. 5 OSIRIS-NAC image taken on 25 December, 2016 (NAC_2016-05-25T15.32.54.769Z_1397549100_F22). The white circle indicates the analyzed boulder. Bottom right panel: location of the Imhotep region on the comet nucleus.
In the text

	Fig. 6 Incidence angle error as a function of the surface standard deviation, which refers to the surface granularity.
In the text

Fig. 7 Time evolution of the boulder height differences ΔH* (as defined in text). Each time step contains information related to 22 boulders and the values correspond to the weighted average. They-axis reports the average height measurements referring to the subsets described in Sect. 4. Subset 1 refers to 21–22 August 2014. According to the illumination and visibility conditions, the time interval between Subset 1 and the perihelion is about 4 months. The time interval between the perihelion and Subset 6 is about one year. Erosion during the inbound orbit until December, 2014, nearly balances the fallout from the southern hemisphere during perihelion cometary activity. The dotted red line indicates the perihelion.

In the text

	Fig. A.1 Time evolution of the boulder height. The histogram shows the average height measurements for each boulder. As shown inthe legend, the results are divided into five time periods to better represent the trend.
In the text