© ESO 2019

1 Introduction

Comets are considered primordial objects that can reveal information about the formation and evolution of the solar system (Glassmeier et al. 2007). Their history may provide indications to physical and chemical conditions within the nebula in which our solar system formed.

The European Space Agency’s Rosetta mission was first designed to orbit and land on a comet. Launched in 2004, the Rosetta spacecraft arrived at its primary target, the Jupiter-family comet 67P/Churyumov-Gerasimenko (hereafter 67P), on 5 August 2014. The Rosetta payload included a suite of instruments that investigated the coma, the chemical interaction with radiation and the solar wind, the activity of the comet, and the nucleus and its surface morphology. The Optical, Spectroscopic and Infrared Remote Imaging System (OSIRIS, Keller et al. 2007) on board Rosetta consisted of two cameras operating from near-ultraviolet to near-infrared wavelengths. The Wide-Angle Camera (WAC) had a spatial scale of 10.1 m px⁻¹ at 100 km from the surface, and its primary objective was to image the dust and the gas directly surrounding the nucleus. The Narrow-Angle Camera (NAC) was designed to resolve the nucleus of 67P with a spatial scale of 1.86 m px⁻¹ at the same distance. The surface of 67P reveals a huge variety of terrains and geological features (Thomas et al. 2015). In particular, the comet has two striking types of landforms: smooth flat regolith plains, vertical cliffs, and talus aprons. The observations of the nucleus surface revealed layers that corresponded to terraces (planar patches of terrains arranged in staircase patterns of overimposed tabular elements) and/or as linear traces on vertical cliffs and on the wall off pits (Massironi et al. 2015; Thomas et al. 2015; Vincent et al. 2015). These surfaces, their distribution, and the magnitude of insolation (Vincent et al. 2017) seem to define the morphologies that are present on the cometary surfaces. The presence of layering can also exert an effect on the fragmentation processes at the scale of boulders, for instance by producing tabular-shaped boulders rather than spherical ones. The large blocks are one of the ubiquitous and most important geological features of 67P: these boulders can be found both isolated and in clusters, and their size distribution depends on the formation and evolution they have undergone (Pajola et al. 2015). Being exposed on the surface of the comet, these objects undergo intense thermal fatigue, and reveal imprints of erosional and geological processes (Ehlmann et al. 2008). It is therefore important to understand the origin of these objects and what this means in the framework of cometary evolution.

The goal of this study is to quantitatively describe isolated populations of boulders located in three different regions (Fig. 1): Imhotep, Hapi, and Hatmehit (Thomas et al. 2015; El-Maarry et al. 2015), using images of the OSIRIS NAC. In particular, we identified numerical parameters (circularity and solidity) with which the shape of boulders can be characterized and perform a statistical analysis of their distribution and properties in these regions.

2 Data and methods

The sublimation of ices plays a fundamental role on the nuclear surface of 67P (Auger et al. 2015). It is expected that the cometary activity increases as the comet approaches the Sun, altering the nuclear surface of 67P, and several studies have investigated the link between surface processing and activity (Hu et al. 2017; El-Maarry et al. 2017; Fornasier et al. 2019; Lai et al. 2018; Vincent et al. 2015). Fragmentation processes on a cometary body are thought to be the result of cumulated thermal fatigue (Pajola et al. 2017a), propagating fractures within the material, and the combined effect of gravitational force (Pajola et al. 2015). In this work we are interested in characterizing isolated boulder populations that are unrelated to specific niches or detachment scarps. Boulders are scattered allover the surface of the comet, and we cannot be sure about their origin. These boulders could have been the result ofpast collapses or gravitational events. However, by discarding boulder populations originating from confirmed gravitational phenomena, we can investigate other types of fragmentation, such as the thermal fatigue and the sublimation of ices. We selected three areas (Fig.1) located on the head, neck, and body of 67P, in order to study the correlation between boulder fragmentation and their location on the comet nucleus. The selected areas are located in the Imhotep, Hapi, and Hatmehit regions (Thomas et al. 2015; El-Maarry et al. 2015).

The Hatmehit region is located on the small lobe at the equator, where the sublimation of ices can be stronger than in other regions of 67P, and the erosion can be much faster (Kossacki & Czechowski 2018). The Hapi region is located in the northern hemisphere of 67P, in the neck that joins the large and small lobe of the comet. Because of its location and self-shadowing, Hapi experiences day-night cycles twice every comet rotation (Pajola et al. 2019) during the summer when the comet is at heliocentric distances >1.6 au pre-perihelion and >2.6 au post-perihelion (Keller et al. 2015). Imhotep is a complex region located close to the cometary equator. It is illuminated daily from aphelion (5.7 au) to perihelion (1.2 au), and it is relatively flat compared to the shape of the nucleus (Auger et al. 2015). Because of its location, illumination conditions, and geomorphology, this region is a good candidate for investigating sublimation erosion close to perihelion. The size of boulders located on this region has no obvious correlation with the gravitational slope of the terrain on which they stand (Auger et al. 2015), and they vary widely in size, shape, and surface texture. We emphasize that the following analysis focuses on specific and confined areas withinthese three regions of 67P, and the obtained results might not be representative of the entire regions.

Fig. 1 OSIRIS NAC images taken on 5 and 6 August 2014 showing the location of the Hatmehit (A, red), Hapi (B, green), and Imhotep (C, blue) regions of 67P. The first image of this set (A) was taken at 23.19.25 UT at a distance of 123.39 km and a scale of 2.30 m px−1. Image B was taken at 03.19.25 UT at a distance of 115.22 km and has a scale of 2.14 m px−1, while the last one (C) was taken at 06.19.26 UT at a distance of 109.70 km and a scale of 2.04 m px−1.

2.1 Size-frequency distribution and cumulative fractional area

Boulder counting can be used to understand and characterize surface properties (Tang et al. 2017) and the processes that result in fragmentation. Characterizing the population of boulders and the cumulative size-frequency distribution may lead to a better comprehension of the geologic history of the cometary surface (Garvin et al. 1981). Using orthorectified images from the OSIRIS-NAC, we derived the cumulative size-frequency distribution of three populations of boulders located in Hapi, Imhotep, and Hatmehit. We performed our analysis both before and after the perihelion passage (13 August 2015) to investigate the role of insolation in terms of fragmentation. Pre- and post-perihelion images have a different spatial resolution because of the different spacecraft range, which affects boulder mapping. To directly compare these two situations and to obtain a homogeneous data set, we used the lowest resolution images as reference, and degraded the highest resolution images to the same levels. The spatial reference scale for the Imhotep case study area is that obtained on 29 September 2014 at 13:29 UT and is equal to 0.35 m px⁻¹. For the Hapi area we useda spatial scale of 0.38 m px⁻¹, derived from an image dated 10 December 2014. For Hatmehit we used an image acquired on 12 December 2014, with a spatial scale of 0.36 m px⁻¹. Following the definition of “boulder” given by Pajola et al. (2015), we used the Arcgis 10.3.1 software (ESRI 2011) to manually outline all positive reliefs detectable in multiple images obtained with different observing geometries. We measured the position of boulders on the surface, the corresponding projected area, and their mean long axis. We only considered fully resolved boulders, that is, features larger than 3 pixels (Pajola et al. 2015), to minimize the possibility of misidentification (Nyquist 1928). This implies a boulder size of at least 1.08 m on Hatmehit, 1.05 m on Imhotep, and 1.14 m on Hapi. In Table 1 we present the observing log of the image dataset we used in this work to measure the size-frequency distribution of the boulders, and in Fig. 2 we show the areas of interest. As displayed in Figs. 2A and D, we decided to analyze both the entire population of the Imhotep region and a subpopulation consisting of a main boulder surrounded by its smaller fragments (hereinafter “detail”). White boxes are used to indicate this particular. The aim is to compare this situation with the entire population, assuming that the fragments derive from the parent body. As mentioned in the previous paragraph, we focused on areas that are unrelated to detachment scarps and niches. To confirm that these populations are not altered by local gravitational slopes, we report the gravitational slope maps (Fig. 3, the centrifugal force is included), that is, the angle between the local surface normal and the vector opposite to the estimated acceleration fields (Penasa et al. 2017a). All areas have gravitational slopes ranging between 0° and 20°, which is below the angle of repose of loose granular material on 67P. This angle is 30° (Vincent et al. 2016).

The finite spatial resolution of the images results in a lack of count completeness at small boulder sizes (Li et al. 2017). To compensate for this effect, we estimated the number of boulders smaller than 3 pixels by subtracting the model number (N_m) of boulders derived from the power law to the number of counted boulders (N_r), previously multiplied by the corresponding area (Li et al. 2017).

Finally, we derived the cumulative fractional area, which is the curve of area covered by boulders versus diameter represented in a log–log plot (Golombek et al. 2003). As shown by Golombek & Rapp (1997), a general cumulative area distribution of blocks can be described by the following equation: $$F(R)=C{R}^{-S},$$(1)

where F(R) is the cumulative fractional area coveredby blocks of a given diameter R or larger, C is a constant derived from the cumulative area covered by boulders greater than 3 pixels, and S is a parameter that explains how abruptly the area covered by rocks decreases with increasing diameter.

Fig. 2 OSIRIS NAC pre- and post-perihelion images used in this work (see Table 1 for image ID). (A) OSIRIS NAC image taken on 12 December 2014. The spatial scale of the image is 0.36 m px−1. This area is located in the Hatmehit region. (B) OSIRIS NAC image taken on 23 July 2016 with a spatial scale of 0.19 m px−1. This area is located in the Hatmehit region. (C) OSIRIS NAC image taken on 10 December 2014 with a spatial scale of 0.38 m px−1. This area is located in the Hapi region. (D) OSIRIS NAC image taken on 20 July 2016. The spatial scale of the image is 0.17 m px−1. This area is located in the Hapi region. (E) OSIRIS NAC image taken on 29 September 2014. The spatial scale of the image is 0.35 m px−1. This area is located in the Imhotep region. (F) OSIRIS NAC image taken on 23 July 2016 with a spatial scale of 0.28 m px−1. This area is located in the Imhotep region. In panels E and F we highlight a detail of the Imhotep case study area.

Fig. 3 Gravitational slopes of the selected areas. Panels A, C, and D: OSIRIS NAC pre-perihelion images used in this work (see Table 1) for image ID. Panels B, D, and F: related slopes maps are shown. The values ranges from 0 to 90degrees.

2.2 Fractal theory

Empirical data indicate that the size-frequency distribution of material expected from fractures and fragmentation follows a fractal rule (Bittelli et al. 1999). The reason for this is that fragmentation is caused by the propagation of fractures at different length scales (Li et al. 2017), and this behavior can therefore be described by connecting the number of crack initiators with the fractal dimension. In this context, a fractal fragmentation model can be used to quantify and predict the size-frequency distribution of fragments produced by a given energy input. The power-law exponent determined by the size-frequency analysis can be interpreted as a fractal dimension. This describes the exact measure for the filling capability of the fragments generator.

The size-frequency distribution of the boulders for an area combines the number of counted boulders N and their size R, with the relation (Mandelbrot 1983; Matsushita et al. 1985; Turcotte 1986) $$N(R)=k{R}^{-D},$$(2)

where N(R) is the number of boulders with a diameter larger than a fixed value R, k is a constant proportionality, and D is the power-law exponent. D characterizes the population itself and the degree of fragmentation. At the same time, Eq. (2) can be interpreted as the number of elements N at the Rth level in the hierarchy of a fractal system, k is the number of crack initiators of unit length, R represents the element size, and D is the fractal dimension (Mandelbrot 1983).

The size-frequency distributions can be interpreted as the fractal behavior of the boulder populations, produced by fragmentation processes. To test this, we used the box-count approach (Smith et al. 1989; Miloevic et al. 2013), which is a common method for analyzing complex patterns and to determine their fractal dimension in a Euclidean space. Here, the image is overlaid with a grid (Miloevic et al. 2013) and is subdivided into progressively smaller squares. Then, each box is checked to determine whether it contains any boulders. At each scale a different number of boxes N will containboulders, and that number is used to obtain the fractal dimension. If N(r) is the number of boxes of side length r required to cover the image, the fractal dimension D_f is defined as $$D_{\text{f}}={\text{lim}}_{r\to 0}\frac{\text{log}(N(r))}{\text{log}(1/r)}.$$(3)

If this limit does not exist, the analyzed image does not show fractal behavior.

2.3 Shape of boulders

The shape of boulders may provide indications on geological processes that leave morphological marks on rock surfaces at different spatial scales (Viles 2001). This is particularly true if signatures in boulder shape and surface texture can be identified and interpreted (Ehlmann et al. 2008). The aim of this section is to generate a quantitative purpose-built parameter set able to represent boulder size and shape. We consider that multiple processes may occur as timescales under consideration increase. This complicates the interpretation of erosion signatures. In the following section the most common shape parameters are listed. We provide nine metrics, some of which are dependent on each other. Then we perform some tests to verify their invariance with respect to different observation geometries.

The original shape of a particle, its lithology, its composition, and the duration of the processes responsible for the erosion are the main influence on its final shape. The lithology of a block is the primary factor that constrains whether and how a particle breaks (Yingst et al. 2007). Layered aggregates might present well-defined planes of weakness parallel to the layers, providing planar elements when subject to fragmentation. On the other hand, corners and edges are more vulnerable than faces because in situ weathering acts more effectively on the points of highest exposed surface area compared to volume (Yingst et al. 2007).

Shape factors are dimensionless quantities that are used in image analysis to quantitatively describe the shape of a particle independently of its size. They represent the degree of deviation from an ideal shape, such as circle or sphere, and are calculated from measured dimensions, such as diameter (d), perimeter (P), area (A), and major (a) and minor (b) axis. The major-to-minor axis ratio of a particle (i.e., the projection of its shape) is described by the elongation (E = b∕a). The measure of how closely a polygon shape matches that of a circle is defined by the Cox circularity (Circ = 4πA∕P²). These parameters may also indicate the presence or absence of edges and corners and may therefore describe if the polygon has a relatively simple boundary with vertices that are relatively equidistant from the centroid. Another important parameter is the solidity (S = A∕H, where H is the convex hull area of the shape; Ballan et al. 2005). It describes the extent to which the shape is convex or concave. In particular, it expresses how the border is ruffled or how many concave cavities are on the surface. The size of an object along a specified direction is quantified by the Feret diameter (d_f). In general, it can be defined as the distance between the two parallel planes that restrict the object perpendicularly to that direction.

In addition to these parameters, another important particle shape factor is its sphericity. This is the measure of how closely the shape ofan object approaches that of a mathematically perfect sphere. This is a three-dimensional parameter that can be estimated by measuring only two dimensions of a particle. This is possible because the correlation coefficient between sphericity values calculating two to three dimensions is 0.85 (Riley 1941).

Defined by Wadell (1933), the sphericity Ψ of a particle is the ratio of the surface area of a sphere, which has the same volume as the given particle, to the surface area of the particle. To address the problem of the dimensionality, particle sphericity has been defined in various ways (Zheng & Hryciw 2015). The four most commonly used formulations are:

• the area sphericity,

$${\Psi }_1=\frac{A_{\text{s}}}{A_{\text{circ}}},$$(4)

where A_s is the projected area of a particle, and A_circ is the area of the minimum circumscribing circle;

• the diameter sphericity,

$${\Psi }_2=\frac{D_{\text{c}}}{D_{\text{circ}}},$$(5)

where D_c is the diameter of a circle that has the same projected area as the particle, and D_circ is the diameter of the minimum circumscribing circle;

• the circle ratio sphericity,

$${\Psi }_3=\frac{D_{\text{ins}}}{D_{\text{circ}}},$$(6)

where D_ins is the diameter of the largest inscribing circle and D_circ is the diameter of the minimum circumscribing circle.

• The perimeter sphericity,

$${\Psi }_4=\frac{P_{\text{c}}}{P_{\text{s}}},$$(7)

where P_c is the circumference of a circle that has the same projected area as the particle and P_s is the perimeter of the particle.

The method we present here was adapted to the specific limitations of data derived from boulders in the OSIRIS NAC images. First, we retrieved these selected images from the Osiris Archive Catalogue, then we manually outlined every single boulder on the image to fix the shape and to calculate the area of all boulders, their perimeter, and the major and minor axes. After the image mapping, we analyzed every shape through a Java program called ImageJ (Schneider et al. 2012) to quantify the shape parameters.

2.4 Error sources

As described in the work of Yingst et al. (2007), one of the sources of error for particle shape analysis is poor pixel selection when the border of the boulder is sorted. Errors in outlining depend on the capability of the analyst to create an accurate contour of the particle. In the literature, the calculation of this uncertainty has been estimated by the mean values of resulting morphologic parameters derived by each one of several operators. The resulting error estimate is approximately 3% for sphericity and 20% for roundness (Yingst et al. 2007). Here, we analyzed pre- and post-perihelion images of the same regions, but because the two incidence angles vary, the inclination of the focal plane with respect to the target surface and the position of the spacecraft introduces further errors. In order to assess the reliability of the shape parameter set, we performed two tests to validate our method.

2.4.1 Two-dimensional test

We created a test image containing different shapes, ranging from the simplest, like a circle, to the most complex, like a star (see Fig. 4). We rotated and tilted the image under different angles, varying the azimuth and the elevation angles of the view in order to simulate the real conditions of the observations. This process is equivalent to applying different homographies, increasing perspective (high elevation), or changing the image plane rotation (by changing the azimuth) to the image. We summarize the test results here. We randomly rotated the image for the first test; we tilted the image by an azimuth angle of 120° and an elevation of 20° for the second test; used an azimuth of 55° and elevation of 40° for the third test; and finally, an azimuth angle of 120° and elevation of 50° for the last test. The resulting test images are shown in Fig. 4. We repeated the test several times with different values of azimuth and elevation. All the test images were then given as input to ImageJ for the shape detection and parameter measurements. Table 2 shows that the circularity and solidity are invariant shape factors. Their invariance allows us to use them to characterize the boulder populations.

Fig. 4 Variation in azimuth and elevation angle of the test image. (A) Casual rotation of the test image. (B) Azimuth = 120°, and elevation = 20°. (C) Azimuth = 55°, and elevation = 40°. (D) Azimuth = 120°, and elevation = 50°.

2.4.2 Three-dimensional test

To validate our method and to verify the invariance of the shape parameters, we created a three-dimensional test. A two-dimensional test is not sufficient to quantify the biases introduced by measuring a two-dimensional projection of a three-dimensional object. We simulated a boulder field and accordingly measured the shape parameters. The surface has a cometary texture and is composed of 10 boulders, the sizes of which vary from 1 to 5 m. Figure 5 shows an example of the synthetic boulder field. The shape of the boulders is random, and this introduces as many inconsistencies as possible. We combined different emission and incidence angles to include shadows and to measure the shape at different observational geometries. The emission values are 0°, 20°, 35°, and 60°. The incidence angles are 0°, 30°, 45°, and 60°. The zero-emission and zero-incidence cases set the shape parameters, because they are able to measure them in a nadiral way. Table 3 shows that the results are consistent with the previous test, and the solidity and circularity are invariant even in this test.

Fig. 5 Simulation of a boulder field composed of 10 boulders. (A) Emission = 0°, and incidence = 0°. (B) Emission = 0°, and incidence = 30°. (C) Emission = 0°, and incidence = 45°. (D) Emission = 0°, and incidence = 60°.

3 Results

3.1 Size-frequency distribution and cumulative fractional area

The total number of identified boulders in this analysis is 11 811, of which 4581 belong to the Imhotep case study area, 641 to the Imhotep detail, 1402 to the Hapi case study area, and 5187 to the Hatmehit study area. By excluding boulders with diameters smaller than3 pixels, a power-law function can be used to fit the data from each boulder population to represent the general form of the size-frequency distribution of the boulders: $$n(R)=k{R}^{-D},$$(8)

where n(R) is the cumulative density of boulders per km² with diameter larger than R. K and D are constants. The cumulative distribution is approximately linear, and we applied a nonparametric statistical method (e.g., Kolmogorov-Smirnov method) to find the slope of the cumulative distribution (for the method, see Lamy et al. 2004; also Pajola et al. 2017b).

We derived the cumulative size-frequency distribution per km² using the equivalent area computed from the 3D shape model of 67P (Preusker et al. 2017). For each size-frequency distribution we used a constant bin-size of 0.35, 0.36, and 0.38 m for the Imhotep, Hapi, and Hatmehit study areas, respectively. This corresponds to the size of the pixel resolution of the OSIRIS NAC images.

The cumulative boulder size-frequency distribution per unit area on the Hatmehit area has a power index of −2.7 + 0.1∕ − 0.2 before perihelion and −2.8 + 0.1∕ − 0.1 after perihelion. The Imhotep population shows a power-law trend with an index of −2.4 + 0.1∕ − 0.2 before perihelion, and −2.4 + 0.1∕ − 0.1 after perihelion. The results for the detail of the Imhotep boulder population (see Fig. 2A) differ from the other populations. The power index in the pre-perihelion image is equal to −2.7 + 0.1∕ − 0.1, and −2.7 + 0.1∕ −.2 for the post-perihelion picture. The Hapi study area shows a power index of −1.7 + 0.2∕ − 0.2 before the perihelion passage, and −1.2 + 0.2∕ − 0.3 after the perihelion passage. In Table 5 we summarize the power-law indices for the areas we analyzed. In Figs. 6–9, the cumulative size-frequency distributions of boulders per km² are shown. We also included the total number of counted boulders, the area calculated trough triangulating the boulder point cloud on the 3D shape model, the corresponding fragmented area, the minimum, the maximum, and the average diameters. The squared correlation R² resulting from the least-squares fit of Eq. (8) is larger than 98% for all of the cases.

In orderto obtain complete information on the size-frequency distribution, we estimated the number of unrecognized boulders that are smaller than 3 pixels. Table 4 shows that to extrapolate the size-frequency distribution, we need to include 85 067 boulders for Hatmehit pre-perihelion, and 76 130 post-perihelion, 2794 unrecognized features for Hapi pre-perihelion, and 118 for Hapi post-perihelion. For the entire Imhotep region, we seem to have missed 35 367 boulders in the pre-perihelion image and 29 662 in the post-perihelion image. Finally, we found a number of 9380 unrecognized boulders for the Imhotep pre-perihelion detail and 4330 for the post-perihelion detail. As expected, the number of boulders that are smaller than 3 pixels is very high for the Imhotep and Hatmehit study areas because we considered all boulders whose dimensions would vary from 0.1 to 3 pixels.

Finally,we calculated the area covered by boulders with a diameter larger than R as a function of the diameter itself. Table 6 shows that the cumulative fractional area distributions for the selected areas can bebest fit by the following power-law equation: $$F(R)=C{R}^{-S(R)},$$(9)

where F(R) is the cumulative fractional area covered by boulders with diameters equal to or larger than R, C represents the total area covered by all boulders, and S(R) indicates how abruptly the area covered by boulders decreases with increasing size. As shown in Table 6, we found that the Hatmehit and Imhotep areas show a cumulative fractional area that is well fit by a power law, as expected. The Hapi case study area instead shows a low value of the square correlation. This suggests that this population does not cover the area in the same way.

Fig. 6 Cumulative size-frequency distribution of boulders per km2. Upper panel: Hatmehit pre-perihelion image. Lower panel: Hatmehit post-perihelion image. The vertical error bars indicate the root of the cumulative number of counting boulders. The fitted regression lines have a power-law index of −2.7 +0.1/−0.2 and −2.8 +0.1/−0.1.

Fig. 7 Same as Fig. 6 for Hapi. The fitted regression lines have a power-law index of −1.7 +0.2/−0.2 and −1.2 +0.2/−0.3.

Fig. 8 Same as Fig. 6 for Imhotep. The fitting regression line gives a power-law index of −2.4 +0.1/−0.2 and −2.4 +0.1/−0.1.

Fig. 9 Sameas Fig. 6 for Imhotep (detail). The fitted regression lines have a power-law index of −2.7 +0.1/−0.1 and −2.7 +0.1/−0.1.

3.2 Fractal behavior

As shown by Mandelbrot (1982), Perfect (1997), and Xie (1993), the size-frequency distribution of fragments resulting from fragmentation follows a fractal rule. Using the box-count method, we further investigated whether the analyzed boulder populations can be described by a fractal law. The aim of this approach is to quantify the fractal scaling, although this factor is not always known a priori. The calculation begins with an arbitrary number of boxes of size r. These elements are then applied to the dataset and counted. This optimized way of cutting an image will reveal the scaling factor.

The number N of boxes of size r needed to cover the set follows a power law: $$N=N_0{r}^{-D_{\text{f}}},$$(10)

with D_f the fractal dimension. Applying the algorithm, we determined the number of boxes N as a function of the size r, and it is represented by a solid line in Figs. 10–13. In the same figure, the scaling law n ∝r⁻² (required to completely space-fill a two-dimensional image) is represented by a dashed line. A clear separation between the two curves indicates a possible fractal behavior of the data set.

Figures 10–13 show that the analyzed populations seem to be fractal at all dimensions, without differences before and after the perihelion passage. The Hapi distribution, however, shows a possible fractal behavior for small dimensions alone. If the box dimensions increase and larger boulders are included, the fractal behavior disappears. Table 7 lists the resulting fractal dimensions D_f.

Fig. 10 Results of the box-count method for the Hatmehit region. The dashed line shows the scaling n ∝ r−2 expected for a space-filling 2D image. The solid line represents the power law N = N0r−Df. The possible fractal behavior of the dataset is indicated by the discrepancy between the two curves.

Fig. 11 Same as Fig. 10 for Hapi.

Fig. 12 Same as Fig. 10 for Imhotep.

Fig. 13 Same as Fig. 10 for Imhotep (detail).

3.3 Boulder shape

The quantitative shape factors can give meaning to physical characteristics if a reasonable comparison between populations can be made. In addition, these parameters can provide a powerful method for assessing the geological history of the surface when the data provision is limited (Yingst et al. 2007). The mean value of each boulder shape factor is considered representative of the entire population. Because we found similar average values for all boulder distributions, we can probably state that the boulder distributions are comparable in all of the inspected areas.

Table 8 shows that the average solidity of boulders ranges from 0.87 to 0.89. The concept of circularity is strictly connected to the compactness and complexity of a particle. The values range from 0.83 to 0.92, and this means that boulders at all sites are generally compact and there are small cavities at the edge of the boulders.

When we compare the mean values of the shape parameters before and after perihelion, we do not note significant differences among the analyzed regions. We graphically represent the values of circularity and solidity for all populations (Fig. 14). Finally, we report a comparison between circularity, solidity, and boulder size for each population (Figs. 15–18). The trend lines show that the circularity and solidity decrease with increasing diameter. This result would be consistent with the fragmentation process due to thermal fatigue. This process removes facets and corners and rounds smaller features more quickly because they are characterized by fewer sensing elements.

Fig. 14 Graphical representation of circularity (square) and solidity (triangle). The average values of these shape parameters suggest boulders with a moderate circularity and solidity. These values are considered as representative of the entire populations. The image in the right box (Leibrandt & Pennec 2015) has been modified to emphasize the relevant parameters.

Fig. 15 Comparison between circularity, solidity, and boulder size. Upper panel: circularity. Lower panel: solidity. The stars represent the Hatmehit pre-perihelion case (the solid line is the corresponding trend line), and the diamond represents the Hatmehit post-perihelion case. The dotted trend line refers to the post-perihelion case.

Fig. 16 Same as Fig. 15 for Hapi.

Fig. 17 Same as Fig. 15 for Imhotep.

Fig. 18 Same as Fig. 15 for Imhotep (detail).

4 Discussion

From our analysis of 11 811 boulders we found significant differences between the Imhotep, Hapi, and Hatmehit study areas in size-frequency distribution, cumulative fractional area, and fractal analysis. Surprisingly, we found no difference in size-frequency distribution and shape before and after perihelion. There may be several explanations for this: a first reason might be that the heat flux density received during perihelion passage is not enough to change the examined populations. Second, there could have been changes, but the erosion was uniform and the shape parameters can only distinguish differential erosion. This erosion occurs at irregular or varying rates caused by the differences in the resistance and hardness of surface materials. This uniformity in erosion could therefore reflect the uniformity of cometary material.

Analyzing the cumulative size-frequency distributions, we noted a substantial difference between the power index of Imhotep, Hatmehit, and Hapi populations. The first two are in the range of −2.3/−2.7. The Hapi area differs from this behavior, with a value of −1.2/−1.7 pre- and post-perihelion. As suggested by Pajola et al. (2015), these data can indicate different formation processes of the boulders. In particular, collapses and pit formation are characterized by power-law exponents of about −5 to −6.5. Indices between −3.5 and −4.5 are typical for gravitational events caused by sublimation and thermal erosion, while material formed during gravitational events and collapses that has undergone continuous sublimation typically shows power-law indices of −1 to −2 (Pajola et al. 2015). The values of −2.3 to −2.7 could therefore represent a transition between these two scenarios.

Table 5 shows that the minimum, maximum, and average values of the boulder diameters change after perihelion passage. The highest resolution of post-perihelion images allows us to outline more boulders and more details, decreasing the minimum and the average value of diameters. The maximum diameter changes because pre- and post-images have different observing conditions, and this is reflected in a diameter variation that is justified by an intrinsic error in the boulder selection method. From this, we conclude that the distribution of boulders does not show significant variations with respect to the perihelion passage, in agreement with what is found in other works (Lucchetti et al. 2017; Pajola et al. 2015).

We estimated the number of undetected boulders that are smaller than 3 pixels to obtain complete information on the size-frequency distribution, finding a large number of missing blocks for all populations, largest before the perihelion passage due to the lower image resolution. The calculated fractal dimension indicates that all study areas show a possible fractal behavior, except for the larger boulders of Hapi. We also examined the position of the boulders with respect to their size. The cumulative fractional area explains how abruptly the area covered by the boulders decreases with increasing diameters. A power law fits the Imhotep and Hatmehit populations, but the CFA of the area in Hapi cannot be described well by a power law. Here, that is, on the neck of the comet, the analyzed boulders reveal a spatial distribution that is also different from the other populations.

The results obtained so far show that the Hapi area differs from the others in size-frequency distribution, cumulative fractional area, and fractal analysis. One possible explanation could be that these boulders collapsed during gravitational events and their fragmentation in situ is controlled by thermal fatigue. This can be justified by the location of Hapi, which is a connecting region between the two lobes of comet 67P that represents a gravitational potential well. Another scenario takes into account the configuration of the Hapi region. Penasa et al. (2017b) showed that the peculiar alignment of boulders in the Hapi region is geometrically consistent with the inner layering of comet 67P (Massironi et al. 2015). The boulders appear to be aligned along a specific layer of the large lobe. The low values of the power index derived for the boulders of Hapi and the differences found in the various parameters analyzed here compared to other populations might reflect this particular relationship with a specific layer of the large lobe. This layer may have undergone fragmentation, which left behind a field of boulders that is aligned with the layering. These boulders would then represent the tops of outcrops, immersed in a deposit of back-fall material that is several tens of meters thick (Keller et al. 2017).

Finally, we defined a set of shape parameters to describe the boulder morphology. Solidity and circularity are invariant shape parameters under changes of observing geometries. We noted that boulders have a homogeneous morphology across the different populations.They consist of compact, convex, rounded, and solid boulders without many inlets. Boulders with average dimensions appear not to exhibit an enhanced tabular shape (i.e., low elongation values), which would be expected for a fragmentation that preferentially occurs along planes that are defined by a tight alternation of layers. The current rounded shape might be connected to a pattern of equally spaced planes of weakness and/or to the weathering processes that follow any kind of fragmentation (i.e., sublimation processes). This result also indicates how sublimation and thermal stress act on the surface of the comet. The near-circularity of all boulders points to a limiting case, where uniform erosion shrinks the boulders on all sides. This might also reflect that the processes involved in the fragmentation are not much influenced by pre-existing planes of weakness such as layers or fractures (and hence other factors exert larger control on the final shape), or that layering on this scale (as opposed to the global scale of the body) might be not present at all. This case would be consistent with the observations made by Belton et al. (2018), who suggested a layer thickness of about 10 meters, and with the observation that the studied boulders do not exhibit evident layering. Similarly, a tight alternation of layers characterized by a variable attitude to sublimation would have been expected to produce more inlets by favoring differential sublimation (enhancing the sublimation of the more erodible layers), suggesting that this kind of heterogeneity is not present at the boulder scale either, or that it cannot be observed due to other prevailing processes. On the other hand, larger boulders have diameters of 20–30 m, comparable to the layer thickness (or multiples of it) suggested by Belton et al. (2018). In that case, decameter-scale layers wouldprovide decameter-scale boulders with moderate to high values of circularity and solidity (Fig. 14).

In conclusion, we presented a detailed quantitative analysis of boulder populations that are located in three different regions of the cometary nucleus of 67P: Imhotep, Hapi, and Hatmehit. By using OSIRIS-NAC images, we identified 11 811 boulders, 4581 of which belong to the Imhotep population, 641 to the Imhotep detail, 5187 to the Hatmehit population, and 1402 to the Hapi region. We proposed techniques for analyzing boulder populations that can complement previous studies, analyzing the size-frequency distribution and the cumulative fractional area of these populations, and connecting these data with the fractal theory. We also outlined the boulder morphology, providing a shape parameter set that does not depend on variations of observing geometries. The analysis of the cumulative fractional area, the link with the fractal theory, and the analysis of the shapes can provide additional information about the population of boulders. As demonstrated by the Hapi case study area, anomalies can be identified, and these differences give the possibility to further investigate the origins of these boulders. The analysis was also performed to understand the role of the perihelion passage of 67P in the fragmentaryand erosive processes. We expected some variations in terms of number and shape of fragments after the perihelion passage, which is the point of the comet orbit closest to the Sun and results in the position with the highest insolation. Our results do not show variation of this type, at least within the limit of this image resolution.

Acknowledgements

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 thank Frederic Moisy for sharing the box-count work.

References

All Tables

All Figures

Fig. 1 OSIRIS NAC images taken on 5 and 6 August 2014 showing the location of the Hatmehit (A, red), Hapi (B, green), and Imhotep (C, blue) regions of 67P. The first image of this set (A) was taken at 23.19.25 UT at a distance of 123.39 km and a scale of 2.30 m px−1. Image B was taken at 03.19.25 UT at a distance of 115.22 km and has a scale of 2.14 m px−1, while the last one (C) was taken at 06.19.26 UT at a distance of 109.70 km and a scale of 2.04 m px−1.

In the text

Fig. 2 OSIRIS NAC pre- and post-perihelion images used in this work (see Table 1 for image ID). (A) OSIRIS NAC image taken on 12 December 2014. The spatial scale of the image is 0.36 m px−1. This area is located in the Hatmehit region. (B) OSIRIS NAC image taken on 23 July 2016 with a spatial scale of 0.19 m px−1. This area is located in the Hatmehit region. (C) OSIRIS NAC image taken on 10 December 2014 with a spatial scale of 0.38 m px−1. This area is located in the Hapi region. (D) OSIRIS NAC image taken on 20 July 2016. The spatial scale of the image is 0.17 m px−1. This area is located in the Hapi region. (E) OSIRIS NAC image taken on 29 September 2014. The spatial scale of the image is 0.35 m px−1. This area is located in the Imhotep region. (F) OSIRIS NAC image taken on 23 July 2016 with a spatial scale of 0.28 m px−1. This area is located in the Imhotep region. In panels E and F we highlight a detail of the Imhotep case study area.

In the text

	Fig. 3 Gravitational slopes of the selected areas. Panels A, C, and D: OSIRIS NAC pre-perihelion images used in this work (see Table 1) for image ID. Panels B, D, and F: related slopes maps are shown. The values ranges from 0 to 90degrees.
In the text

	Fig. 4 Variation in azimuth and elevation angle of the test image. (A) Casual rotation of the test image. (B) Azimuth = 120°, and elevation = 20°. (C) Azimuth = 55°, and elevation = 40°. (D) Azimuth = 120°, and elevation = 50°.
In the text

	Fig. 5 Simulation of a boulder field composed of 10 boulders. (A) Emission = 0°, and incidence = 0°. (B) Emission = 0°, and incidence = 30°. (C) Emission = 0°, and incidence = 45°. (D) Emission = 0°, and incidence = 60°.
In the text

	Fig. 6 Cumulative size-frequency distribution of boulders per km2. Upper panel: Hatmehit pre-perihelion image. Lower panel: Hatmehit post-perihelion image. The vertical error bars indicate the root of the cumulative number of counting boulders. The fitted regression lines have a power-law index of −2.7 +0.1/−0.2 and −2.8 +0.1/−0.1.
In the text

	Fig. 7 Same as Fig. 6 for Hapi. The fitted regression lines have a power-law index of −1.7 +0.2/−0.2 and −1.2 +0.2/−0.3.
In the text

	Fig. 8 Same as Fig. 6 for Imhotep. The fitting regression line gives a power-law index of −2.4 +0.1/−0.2 and −2.4 +0.1/−0.1.
In the text

	Fig. 9 Sameas Fig. 6 for Imhotep (detail). The fitted regression lines have a power-law index of −2.7 +0.1/−0.1 and −2.7 +0.1/−0.1.
In the text

	Fig. 10 Results of the box-count method for the Hatmehit region. The dashed line shows the scaling n ∝ r−2 expected for a space-filling 2D image. The solid line represents the power law N = N0r−Df. The possible fractal behavior of the dataset is indicated by the discrepancy between the two curves.
In the text

	Fig. 11 Same as Fig. 10 for Hapi.
In the text

	Fig. 12 Same as Fig. 10 for Imhotep.
In the text

	Fig. 13 Same as Fig. 10 for Imhotep (detail).
In the text

	Fig. 14 Graphical representation of circularity (square) and solidity (triangle). The average values of these shape parameters suggest boulders with a moderate circularity and solidity. These values are considered as representative of the entire populations. The image in the right box (Leibrandt & Pennec 2015) has been modified to emphasize the relevant parameters.
In the text

	Fig. 15 Comparison between circularity, solidity, and boulder size. Upper panel: circularity. Lower panel: solidity. The stars represent the Hatmehit pre-perihelion case (the solid line is the corresponding trend line), and the diamond represents the Hatmehit post-perihelion case. The dotted trend line refers to the post-perihelion case.
In the text

	Fig. 16 Same as Fig. 15 for Hapi.
In the text

	Fig. 17 Same as Fig. 15 for Imhotep.
In the text

	Fig. 18 Same as Fig. 15 for Imhotep (detail).
In the text