© ESO, 2015

1. Introduction

After a ten-year journey in the inner solar system, the European Rosetta spacecraft on 6 August 2014, entered the orbit of its final target, the Jupiter family comet 67P/Churyumov-Gerasimenko (hereafter 67P). A few months later, on 12 November 2014, Rosetta succeeded in one of the most ambitious projects ever attempted by a planetary mission: it delivered a man-made lander, named Philae (Biele et al. 2006), on the surface of this small (large lobe: 4.1 × 3.3 × 1.8 km; small lobe: 2.6 × 2.3 × 1.8 km) irregularly shaped and low-gravity (10^-3 m/s²) comet (Sierks et al. 2015).

The location of Philae’s nominal landing site on the surface of 67P was selected on the basis of a preliminary analysis of the first global mapping of the surface obtained with the Optical, Spectroscopic, and Infrared Remote Imaging System (OSIRIS; Keller et al. 2007). The area seemed favorable with good illumination (>6 h of sunlight), the majority of terrain within a square kilometer having a topographical slope <30° relative to the average local vertical and a few large hazardous boulders. The selected area appeared scientifically appealing, containing outcrops of consolidated material and offering interesting surface features to be investigated.

Philae was deployed with a velocity of 18.76 cm/s when Rosetta was at a distance of 22.5 km from the center of the comet, and the lander touched down seven hours after the separation. The ROsetta Lander Imaging System (ROLIS) onboard Philae (Mottola et al. 2007) imaged the surface of the first touchdown while Philae was between 70 and 10 m from the surface (Mottola et al. 2015b). The comparison of ROLIS images with OSIRIS images acquired during the mapping phase allowed confirming that the first touchdown (Fig. 1) occurred inside the predicted landing ellipse, possibly within 130 m from its center. Trajectory calculations and radio wave transmissions, together with data from most of the instruments have been used to reconstruct the path pursued by Philae and suggest that, after the first touchdown, Philae remained unanchored to the surface and drifted above the surface for about 1:50 h, bouncing three times and traveling for about 1 km (Ulamec et al. 2015) toward the large depression known as the Hatmehit region, perhaps sitting now just beyond its rim, possibly on the Bastet region. Throughout the paper we use the regional nomenclature outlined in Thomas et al. (2015a) and further described in El-Maarry et al. (2015).

In this framework, the OSIRIS images are of paramount importance because of their high resolution and multispectral nature, not only to constrain the lander’s position on the comet, but, more importantly, to provide a comprehensive analysis of the comet and in particular of the landing area for the unrivaled possibility of comparison with direct in-situ measurements performed by the lander instruments.

Fig. 1 OSIRIS NAC images acquired on 12 November 2014 at a) 15:18 and b) 15:43 just after the first Philae touchdown, showing evidence of the footprints of Philae’s three legs.

Before the Rosetta mission only five cometary nuclei − 1P/Halley (Keller et al. 1986), 19P/Borrelly (Soderblom et al. 2002), 81P/Wild2 (Brownlee et al. 2004), 9P/Tempel 1 (A’Hearn et al. 2005), and 103P/Hartley 2 (A’Hearn et al. 2011) − have been imaged by different ESA and NASA space missions (see Barucci et al. 2011, for a complete review). High-resolution images of these cometary nuclei revealed irregularly shaped objects showing a complexity of surface morphologies and photometric properties.

The resolution of Halley’s images did not allow a detailed nucleus morphology analysis, while each of the other comets showed its own morphological peculiarity. 19P/Borrelly, imaged at its best resolution of 47 m/pixel by the Deep Space 1 mission, showed four major morphological units: dark spots, mottled terrains, mesas, and smooth terrains, and four terrain features: ridges, troughs, pits, and hills, but did not display any sign of impact craters (Britt et al. 2004). 81P/Wild 2, imaged at the higher resolution of 14 m/pixel by the Stardust spacecraft, instead showed a variety of circular depressions interpreted as impact craters and divided into two categories: pit-halo and flat-floored features bounded by high cliffs. It also showed evidence of flat-topped mesas that have been interpreted as the results of ablation (Brownlee et al. 2004). Comet 9P/Tempel 1, imaged by the Deep Impact mission at 10 m/pixel scale, revealed two main types of terrains: pitted terrains and smoother areas, displaying various surface features, from circular craters to evidence of flows and terraced areas Thomas et al. (2013). Comet 103P/Hartley 2, observed at about the same resolution, showed a distinctive bilobed shape with knobby terrains filled by elevated rough forms on the lobes and much smoother areas on the waist (A’Hearn et al. 2011).

This morphological diversity has been interpreted both in terms of evolutionary sequence (Belton 2010) and of diverse geological outcomes, including cometary activity and eolian erosion (Cheng et al. 2013).

Despite their morphological variety, different sizes, appearance, and activity levels, all visited comets show similar bulk photometric properties, except for Borrelly (Li et al. 2013). Extensive photometry of these cometary nuclei presented in Li et al. (2007a,b, 2009, 2013) shows that they share dark surfaces with a geometric albedo ranging from 0.03 to 0.07 and featureless visible spectra with red slopes of about 7−12%/100 nm and small albedo variegation (10−20%). Comet 19P/Borrelly instead shows stronger variations in brightness of more than a factor of two and strong variations in photometric parameters and roughness (Li et al. 2007b).

The nucleus of 67P shows a specific morphology diversity presented in Thomas et al. (2015a) and its similarities with previously visited Jupiter-family comets are briefly described by El-Maarry et al. (2015). The global photometric properties of 67P, mostly consistent with that of the other visited cometary nuclei, are extensively illustrated by Fornasier et al. (2015).

We present here a geomorphological and spectrophotometric investigation focused on the landing site area on the surface of 67P. The aim is to provide a context description that is essential for interpreting the Philae measurements and for placing constraints on the investigation of the global surface of 67P and other cometary nuclei.

2. Data and methods

We performed our analysis on a dataset of OSIRIS Narrow Angle Camera (NAC) and Wide Angle Camera (WAC) images covering the nominal landing site and extending to the boundary between the Hatmehit depression and the Bastet region, which includes the possible final resting position of Philae on the head lobe of 67P, close to the equator. The overall region, referred here as Agilkia, is represented in Fig. 2 obtained as a composite of five NAC images (sets 1−5 in Table 1) acquired in the orange filter (649.2 nm). The location of the first Philae touchdown is highlighted by the white circle.

Fig. 2 Composite image of Philae landing site (Agilkia) obtained from five NAC frames in Orange filter (sets 1−5 in Table 1) acquired on 14−15 September 2014 at an average phase angle of 64.23°. The circle highlights the location of the first touchdown, while the dashed line indicates the equator.

The full dataset was acquired between 1 August 2014 and 12 November 12 2014 with a phase angle ranging from 10° to 90°. The images were acquired at distances ranging from 735 km to 10 km above the nucleus surface, resulting in a pixel scale from 13 to 0.18 m/pixel. This wide dataset was extremely useful for investigating the complexity of the region, and the varying illumination and observation conditions allowed us to study small details of the surface that would otherwise have been hidden in the shadows. The comet was at 3.6−3.0 AU from the Sun, and although 67P was found to be active already from the first unresolved studies (Tubiana et al. 2015b) and from later observations between August and September (Lin et al. 2015), the analysis of this dedicated dataset also allowed determining that no sign of activity was evident on the landing site area in this time frame (1 Aug.−12 Nov.).

We used the software ArcGis to map and identify geomorphological units and linear features.

For both the geomorphological and the spectrophotometric analysis we made use of a 3D shape model of the comet. The entire NAC dataset has provided multiple stereo coverage of the cometary surface that was used to derive accurate models of the complex shape of 67P at increasing resolutions and with different techniques such as stereo-photoclinometry (SPC; Jorda et al. 2012) and stereo-photogrammetry (SPG; Preusker et al. 2012).

The SPC model, downsized to 1.5 × 10⁵ points, was used to derive the gravitational potential of this highly irregular comet, assuming a homogeneous nucleus with a density of 470 kg/m³ (Sierks et al. 2015) following the theory of Werner & Scheeres (1997) and accounting for the centrifugal force due to the comet rotation (Rossi & Fulchignoni 1999). To describe the gravitational framework of Agilkia for geomorphological interpretation, we analyzed the gravitational slopes of the terrain shown in Fig. 3. The gravitational slope is defined as the angle between the local surface normal pointing inside the nucleus and the gravity vector. The average gravitational slope across the area is 15.6° with wide areas showing slopes <10° and less than 10% of the terrain showing slopes >30°.

Fig. 3 Gravitational slope map of Agilkia obtained using the SPC shape model downsized to 1.5 × 105 points and the gravitational potential model derived assuming a homogeneous nucleus and accounting for the centrifugal force due to the nucleus rotation.

For the spectrophotometric application we instead preferred the SPG model (Preusker et al. 2015) because it makes no assumptions on the albedo. As a result of the importance of the lander operations success, the landing site has been studied in detail, and a dedicated shape model subset with 1 m nominal lateral spacing and an effective lateral resolution of about 2 m has been constructed. It consists of a total of about 2.5 million vertices and 5 million facets. Apart from local exceptions, this subset comprises a description of the surface with a typical vertical accuracy at the submeter level. The SPG shape model subset was used, together with the most recent NAIF-SPICE kernels and its IDL code toolkit (Acton Jr 1996), to correct the surface for topographic effects. A detailed description of the photometric modeling used is given in Sect. 4.

For the spectrophotometric analysis we made use of observations performed with different narrowband filters with central wavelengths ranging from 250 nm to 1000 nm (see Keller et al. (2007) for filter details) to generate color images. The images in different filters are not acquired simultaneously but sequentially, and because of the rotation of the comet and the spacecraft movements, they need to be coregistered for a correct spectral analysis. Because this task is difficult and given the large amount of images, we restricted our dataset to the spectrophotometric analysis and concentrated on nine NAC color sequences that are summarized in Table 1. They were acquired between 14 and 20 September with five filters: near-UV (360.0 nm), blue (480.7 nm), green (535.7 nm), orange (649.2 nm), and near-IR (882.1 nm).

All OSIRIS images are calibrated in radiance (expressed in W m^-2 nm^-1 sr^-1) using the OSIRIS standard calibration pipeline described in Tubiana et al. (2015a) (see also Küppers et al. 2007; Magrin et al. 2015) and were then converted to reflectance (also known as radiance factor) using $\mathit{I}\mathit{/}{\mathit{F}}_{\mathrm{meas}}\mathrm{=}\frac{\mathit{\pi }{\mathit{r}}_{\mathrm{h}}^{\mathrm{2}}{\mathit{I}}_{\mathrm{meas}}}{{\mathit{F}}_{\mathrm{\odot }}}\mathit{,}$(1)where I_meas is the measured radiance, F_⊙/π is the solar irradiance at 1 AU measured at the central wavelength for each filter, and r_h is the heliocentric distance of the comet in AU.

The Integrated Software for Imagers and Spectrometers (ISIS3; Anderson et al. 2004) from the United States Geological Survey (USGS) was used to coregister the color images using a maximum correlation-matching algorithm with a goodness of fit greater than 0.9 and a third-degree polynomial in the linear regression fit of warping. The window sizes of search and pattern were varied in different data sequences to have the best solution for each data sequence. Once aligned, each color sequence was collected in a multispectral image, hereafter called spectral cube, and was used for spectrophotometric studies.

In the following section we describe the identification of different geological units and linear features on the landing site, we investigate the possible origins of these features, and provide a geomorphological map of the area. In Sect. 4 we introduce the photometric modeling performed to correct for topographic and illumination conditions, and we evaluate the improvements of the correction. In Sect. 5 we discuss the photometric variations of the terrain, and in Sect. 6 we investigate the spectrophotometric properties and their relation with the morphology. A discussion of the results and conclusions follow in Sect. 7.

3. Geomorphological analysis

Fig. 4 Geological map of 67P observed from the top of the head lobe and nomenclature of the geological regions (modified from Fig. 2 El-Maarry et al. 2015).

Figures 2 and 4 show that the nominal landing site, and therefore the location of the first touchdown, is located on the Ma’at region, whereas the final position of Philae is presumably located on the other side of Hatmehit, in proximity of the boundary with Bastet (see El-Maarry et al. 2015; and Giacomini et al., in prep., for a description of the regional morphology).

As a result of the uncertainties in the final position of Philae and to give an extensive description of the landing site, we performed a detailed geomorphological mapping of all the area referred to as Agilkia (Fig. 5). This area includes part of the Ma-at, Maftet, Nut, and Bastet regions and the entire Hatmehit region. The main geological units were identified on the basis of their morphological properties.

The parts of Maftet, Ma’at, and Nut that are enclosed in Agilkia show a similar geomorphological character, even thought they are separate physiographic regions (El-Maarry et al. 2015). In contrast, Hatmehit and Bastet show different surface properties that mean that they were very likely subjected to a different formation process. For this reason we describe them separately.

Fig. 5 Geomorphological map of Agilkia identifying the main geological units on the basis of their morphological properties, and various linear features. The map is produced on the composite image in Fig. 2.

3.1. Maftet, Ma’at, and Nut regions

The areas belonging to Maftet, Ma’at, and Nut are mainly covered by smooth deposits of fine material that very likely are the result of airfall due to activity producing low-velocity particles that are unable to escape (Thomas et al. 2015a,b). The size of the grains constituting the surface layer are currently constrained by the resolution limit of OSIRIS images (18 cm), although the ROLIS image revealed that at the location of the first touchdown the surface is covered by fine-grained material on a scale of 9 mm/pix (Mottola et al. 2015b). The thickness of these deposits is unknown, although an impact crater of 35 m diameter detected in the Ash region allowed estimating that the smooth deposit thickness in that region might reach 5 m. Nevertheless, this superficial layer may be highly variable around the comet surface (Thomas et al. 2015a) and seems to be very thin when observed from a lateral view over walls of more consolidated material (see, for example, Figs. 2, 3 in Thomas et al. 2015b).

Fig. 6 Examples of longitudinal dunes observed in Ma’at and Hatmehit. NAC image acquired on 2014-09-14 at a) 17:47:55 and b) 16:56:07 in orange filter.

Fig. 7 Examples of aligned pit-like features visible on smooth deposits in the Ma’at region. NAC images acquired on 2014-09-30 at 09:26:53 (left panel) and on 2014-10-19 at 12:22:15 (right panel) in orange filter. Enclosed by the green dashed line, we show the region for which the local horizon has been extracted, and the red dashed line denotes the considered orientation of the features for comparison with the solar orientation (see text for details).

In the smooth areas of Maftet, Ma’at, and Nut several longitudinal dune-like features are detected (some examples are shown in Fig. 6). These features show a general uniform trend of their longest axis throughout the region (see map in Fig. 5). Since they are located on nearly flat areas in correspondence of <10° gravity slopes (see Fig. 3), a down slope origin can be excluded. The dune-like structures on the surface of the nucleus were already mentioned by Thomas et al. (2015a) from very early observations and are further discussed by Thomas et al. (2015b). The mechanism responsible for the formation of these features is still unknown, although surface dust transport, correlated with comet gas activity, is one of the most plausible candidates. This is supported by the evidence of airfall deposits and other features connected to dust transport phenomena, such as ventifacts, aeolian ripples, moats, and wind tails (see Thomas et al. 2015b). Characteristic features were also observed by the ROLIS instrument at the landing site and are interpreted as wind tail structures (Mottola et al. 2015a). An alternative but less likely hypothesis is that these formations might be related to other kinds of thermal induced processes such as sublimation. However, there is no observational evidence that supports this origin. Therefore, assuming that these dune-like features are the result of wind-driven transport, on the basis of their form and orientation it is possible to derive which wind regime could have formed them. Indeed, longitudinal dunes occur when winds blow from more than one direction, although one of these must predominate. The dominant wind determines the orientation of the dunes, but a crosswind from another direction is needed to form this type of dunes (Pye & Tsoar 2009; Melosh 2011, and references therein). The uniform orientation of these features throughout the Maftet, Ma’at, and Nut regions suggests that the dominant wind direction is constant throughout the whole area. Fine dynamical modeling of the gas flows at the location of the dunes combined with the nucleus rotation are needed to further prove this hypothesis.

Fig. 8 Solar orientation with respect to the aligned pit features. Top panels show the elevation of the Sun over the local horizon, while the bottom panels show the angular separation between the orientation of the features (red dashed line in Fig. 7) and the projected solar orientation on the local horizon. Both angles are given as a function of the heliocentric distance (left panels) and time (right panels) for a period covering one full comet orbit with steps of 30 days for revolution and 30 min for rotation (see text for details).

Fig. 9 Examples of outcrops of consolidated material. NAC images acquired in orange filter on a) 2014-09-14 at 17:31:24, b) 2014-09-14 at 04:28:01, c) 2014-11-12 at 17:08:52, and d) 2014-11-12 at 16:53:51.

Chains of pit-like features generally aligned with the longest axis of the dunes have been observed (Fig. 7) in some regions of Ma’at and Maftet (see map in Fig. 5). They have dimensions from smaller than a meter to 4−5 m in diameter. The mechanism forming these particular features is unknown, although their appearance suggests a post-depositional modification of the surface (Thomas et al. 2015b). The morphology of these forms is comparable to wind-erosion forms such as flutes and pits observed on Earth. Flutes are usually observed on rock surface and originate at pits that become shallower on the downwind side. These pits result from the erosion of wind vortical movements (Richardson 1968; Whitney 1979), and their formation seems to be controlled by the inclination of slopes with respect to the wind direction (Greeley & Iversen 1986). The morphological similarity with the aligned pits observed on Agilkia suggests a similar wind-erosion origin for these features.

However, other processes might have been responsible. Several bright spots were detected on smooth deposits throughout the surface of 67P and have been interpreted as icy chunks (Pommerol et al. 2015; Thomas et al. 2015b). Some of these volatile-rich blocks could have been buried by airfall and/or aeolian transport and been subsequently sublimated. This might have caused the collapse or ejection of the fine non-volatile material with the consequent formation of the pits (El-Maarry et al. 2015; Thomas et al. 2015b). Laboratory experiments on water ices mixed with mineral contaminants under simulated comet surface conditions (see Pommerol et al. 2015, and references therein) indeed revealed the ejection of large to cm-sized grains with the consequent formation of pits in the superficial layer (see Fig. 10 Pommerol et al. 2015). To investigate this possibility, we studied the orientation of the Sun with respect to these features throughout the orbital period of 67P. Figure 8 shows the elevation of the Sun (top panels) over the average local horizon extracted using 180 facets of the SPG shape model across the region enclosed by the green dashed line in Fig. 7. The data were retrieved for a full rotation of the comet with steps of 0.5 h and a full orbital period, including the last perihelion (March 2009), with steps of 30 days. The elevation is shown as function of heliocentric distance (left panel) and time (right panel). The red lines in Fig. 8 indicate the solar maximum elevation (local culmination) for each rotation. Even though Agilkia is close to the equator, the 52° obliquity of the axis of 67P (Sierks et al. 2015) makes this region quite unfavorable to sunlight at perihelion. The highest elevation of the Sun is <10° at the perihelion. This might explain why the pits are not normal to the surface, but seem more orientated along the surface declivity. We also investigated the angular separation between the orientation of the observed feature (identified by the red dashed line in Fig. 7) and the projection of the direction of the Sun on the local horizon (bottom panels in Fig. 8), both as a function of heliocentric distance (left panel) and time (right panel). The average angular separation at the culmination is low (<30°) and comparable with the average orientation variations of these features, suggesting the existence of a correlation between the solar direction and the alignment of the pits. This is an indication that the sublimation might be responsible for the formation of the aligned pits even though a more detailed analysis of the surface slopes and the thermal conductivity of the material is necessary to confidently conclude on in this origin.

Smooth deposits do not cover all of Maftet and Ma’at regions. Outcrops of terrain are widespread (examples in Fig. 9) in correspondence of cuestas (see Fig. 10). Cuestas are landforms showing a steep wall on one side and a gentle dip or back slope on the other. These features are the expression of differential erosion of gently dipping strata on Earth, where hard layers alternate with weaker layers (Fig. 10a) (Simonett 1968; Easterbrook 1999). The continuous stratification on the cometary surface has been established by various observational and geological evidence throughout the nucleus (Massironi et al. 2015). Various parallel linear features (strata heads) have also been established on Agilkia, extending for several meters on the steep walls. The cuestas features on the comet therefore most likely arise from differential erosion by sublimation of strata with different content in volatile material that are then differently erodible. A similar process of differential rates of sublimation erosion was previously invoked to explain the formation and morphology of geological units on the surface of comet Borrelly by Britt et al. (2004). We here suggest that the diverse rate of sublimation is mainly driven by different volatiles content among strata. The frequent overhangs of layers overlying more erodible ones (Groussin et al. 2015) are furthermore in favor of the differential sublimation origin of cuestas landforms.

Fractures have been detected on the steep walls, suggesting a brittle nature of the outcropping terrain. These fractures often cross the strata heads (see Fig. 11) and are randomly oriented (see map in Fig. 5). Insolation weathering related to either thermal fatigue or thermal shock seems to be the responsible for most of the fracturing in 67P (Thomas et al. 2015b).

Fig. 10 Cuesta morphology a) sketch of cuesta landform. Cuestas occur in areas of tilted strata and are caused by the differential erosion of the harder top layer and the underlying layer that is weaker and therefore more erodible (figure after Castiglioni 2000). b)−c) Examples of cuestas on the Agilkia site (indicated by black arrows). The hard layer probably consists of material poor in volatiles, whereas the soft underlying layer is richer in volatile material and therefore more erodible by sublimation. NAC image acquired onb) 2014-09-14 at 17:32:15 in near-IR filter and on c) 2014-11-12 at 16:08:51 in orange filter.

Fig. 11 Examples of fractures (in cyan) and strata heads (in green) detected on the Agilkia site. Usually fractures cut strata heads, revealing that they form after the layering. NAC image acquired on 2014-11-12 at 16:23:51.

Steep walls are found in correspondence with high gravitational slopes (see Fig. 3) of 20° to 40° and usually show unconsolidated material at their foot, which therefore have been interpreted as the results of collapses due to erosion by sublimation. These deposits consist of elements of different dimensions, from boulders of 30 m (Pajola et al. 2015) to very fine particles. We classified these different collapsing deposits (see Fig. 5) on the basis of their distribution, texture, and source area.

Talus and gravitational deposits (examples in Fig. 12) have both a gravitational source and are interpreted as rockfall material deposits resulting from the erosion of the adjacent cliffs. Talus deposits generally consist of fine material uniformly located all along the foot of the walls. Gravitational accumulation deposits are instead composed by particles of a wide range of dimensions. They are not well sorted and include large boulders that are supported by fine material-matrix. They generally form fans located at the foot of cliffs. Their source area is clearly defined and associated with rough semicircular niches formed after the material detached from the steep walls. Finally, diamicton deposits involve poorly sorted material that has no unequivocal source. These deposits have been detected in Nut and fill part of the depressed region.

Fig. 12 Different types of deposits. a) Example of talus deposit. b) Example of gravitational accumulation deposit. NAC image acquired on 2014-09-14 at 05:05:28 in orange filter.

A large quantity of boulders are visible in all types of deposits. They have been identified and analyzed originally for the lander safety requirements and subsequently for geomorphological purposes. We detected the boulders on the Ma’at, Maftet, and Nut parts on Agilkia in the NAC image that was acquired on 2014-09-14 at 17:48:08 with a scale of 0.53 cm/px. The boulders were identified with the same technique as in Pajola et al. (2015), who performed global and localized statistics of boulders ≥7 m throughout the comet. We considered as statistically meaningful boulders larger than 1.5 m in diameter, that is, at least three pixel sampling, which minimizes the likelihood of misidentifications (Nyquist 1928). We identified a total of 2065 boulders in this area, the largest with a dimension of 24.33 m and a spatial distribution as presented in Fig. 13. We then derived the cumulative size-frequency distribution by using bins of 0.5 m (Fig. 14). The resulting power-law fit returns a value of −2.9 +0.1/−0.3. This agrees well with analyses that used a similar resolution and were performed on images obtained with the lander ROLIS camera while it was approaching the surface (Mottola et al. 2015b). However, when this index is compared to the results of the global and localized statistics presented in Pajola et al. (2015), it is possible to see that it is substantially lower than the boulder power-law index derived from pure gravitational events triggered by localized sublimation. The main explanation of this behavior is that the considered area is characterized by different types of deposit, such as smooth fine deposits and gravitational accumulation deposits located at the foot of terrace margins (see, e.g., Fig 12). Therefore, the morphological diversity between the several considered units mixes lower power-law indices (−1 to −2) expected for type 3 boulder clusters, such as more evolved material deposits, and higher power-law indices (−3.5 to −4) expected for type 2 boulder clusters, such as gravitational falls triggered by sublimation and/or thermal fracturing (see Pajola et al. 2015, for further details). The resulting cumulative distribution of the entire area is then consistent with a globally and mostly flat terrain covered with fine deposits on layered material that does not present a heavily fractured texture. In addition, the absence of activity detected in this area throughout the OSIRIS dataset (see also Lin et al. 2015) together with the absence of recent processes affecting the superficial and subsurface layers, can explain the lower value of the cumulative boulder power-law index, when compared with areas that are more active, fractured, and rich in meter-sized boulders.

Fig. 13 Spatial distribution of boulders on the Ma’at, Maftet, and Nut parts of Agilkia. NAC image acquired on 2014-09-14 at 17:48:08 in orange filter.

Fig. 14 Cumulative size-frequency distribution for the total 2065 boulders identified on Maftet, Ma’at, and Nut inside Agilkia.

3.2. Hatmehit region

Fig. 15 Hatmehit region observed from the top and front. NAC image acquired from on 2014-09-02 at 07:42:53 (left panel) and on 2014-09-01 at 08:42:52 (right panel) in orange filter.

The boundaries of the Hatmehit region roughly correspond to the edges of a 0.49 km² circular depression of about 900 m in diameter and 150 m in depth (Fig. 15) located just across the comet’s equator. It is the largest irregularly shaped depression observed on the comet, apart from Imhotep on the larger lobe of the comet (see Auger et al. 2015 for a detailed study of the Imhotep region). This feature is characterized by a rather flat floor bounded by steep walls.

Fig. 16 Examples of terraces detected on the Hatmehit depression rim (indicated by white arrows in panels a) and b)). NAC image acquired on 2014-09-14 at 16.56.07.

Fig. 17 Diamicton deposit in Hatmehit. NAC image acquired on 2014-09-14 at 16:56:07 in orange filter.

Similarly as for the Ma’at and Maftet regions, steep walls show outcropping-layered terrains (Fig. 9a) often affected by fractures that cross the strata heads (Fig. 11). In some areas, steep walls present terraces (Fig. 16) that probably are the results of a differential erosion of nearly horizontal strata. Locally, walls collapsed, forming gravitational accumulation deposits and leaving semicircular niches that are clearly visible all around the rim of the depression (see map in Fig. 5).

The floor of the circular depression of Hatmehit is covered throughout by fine material that forms dune-like features whose axes show a uniform orientation similar to those observed on Maftet and Ma’at. However, the greater part of the floor is covered by a divers deposit of fine material mixed with boulders of various sizes. The source area of this deposit is not clearly defined from a geomorphological analysis and was therefore classified as diamicton (Fig. 17). It is reasonable to assume that it has been generated by diverse concurrent processes such as gravitational accumulation, regressive erosion, and sublimation variously combined. The size-frequency distribution of the boulders inside the depression (Pajola et al. 2015) show a bimodal behavior, presenting two power indicess, −3.4 and −1.0. This is indeed indicative of a mixing deposit due to the concurrent origin from collapse, gravity, and sublimation.

The floor of Hatmehit is crossed by a narrow, linear feature 300 m long that extends through the smooth deposits and probably below the diamicton deposit. The origin of this formation is not known, although its asymmetric transverse profile suggests that it is probably a buried terrace (see also Massironi et al. 2015).

The nature of the Hatmehit depression is still unknown. The crater hypothesis has been examined by Pajola et al. (2015) by relating the putative crater diameter with the largest boulder fragment found in the crater, and they found a possible association. However, the morphology of the depression, in particular the extensive evidence of coherent centripetally arranged terraces and strata heads all around the depression disagree with the expected appearance of an impact crater on a layered terrain, which would otherwise show overturned flaps due to the impact (Melosh 1989; Kenkmann et al. 2014). Moreover, the linear feature that is interpreted as a buried terrace and is not visible outside the depression, likewise contradicts the crater hypothesis. If the central linear feature is a buried terrace margin, it is unlikely that it survived the impact, while if it were a different feature instead, for instance, a ridge produced after the impact, it should be visible across the depression rims and outside. The irregular shape, the depth-to-diameter ratio, and the very low estimated cometary collisional rate (Belton et al. 2013; Vincent et al. 2015) also highly discourage the impact crater hypothesis.

On the basis of the boulder distribution, an alternative hypothesis of the formation of Hatmehit can be proposed (see also Pajola et al. 2015). It might be due to the sublimation of underlying volatile-rich strata covered by volatile-poor layers that cause the collapse of the uppermost strata. After the putative collapse, the whole area could have undergone areal sublimation and regressive erosion at the bordering cliffs. The final result could have been a dismantling and elimination of smaller boulders in the central part of the depression and a progressive erosion and retreat of the surrounding walls with production of debris continuously renovated by blocks and grains.

3.3. Bastet region

The portion of Bastet visible in Agilkia appears to be almost free of fine deposits and consists of consolidated material. This seems to also be supported by the instruments on Philae. Although the power of the hammer of MUPUS instrument was gradually increased, it was apparently unable to deeply penetrate the surface, revealing a solid hardness (DLR 2014). A more detailed analysis of these data is needed, however, to conclude on the strength of the comet surface.

The consolidated terrain appears to be heavily affected by long fractures that are parallel to each other and show a preferred orientation. Finally, the analyzed area lacks large boulders, talus, or gravitational accumulation deposits.

4. Photometric modeling

The photometric properties of airless bodies of the solar system are determined both by the physical and chemical characteristics of the material composing the surface (mineralogical composition, size distribution, and shapes of particles) and by geometrical parameters such as solar phase angle, viewing, and illumination angles. The photometric modeling is used to distinguish the physical characteristics of the material from the mere geometric conditions.

Widely used photometric models (Hapke 1981, 1984, 2002) describe these dependencies with a complex theory of bidirectional reflectance by particulate surfaces including several free parameters. The small phase angle coverage of our dataset (4°) does not allow accurately modeling the observations with the Hapke multiparametric function. However, Fornasier et al. (2015) performed this analysis using a dataset with lower resolution, but wider phase angle coverage on the overall cometary body and derived the Hapke parameters for the orange filter (reported here in Table 2). Therefore, assuming that these parameters are also representative of the examined area, we calculated the photometric corrected reflectance of Agilkia using the Hapke model in the orange filter. We applied for each pixel $\mathit{I}\mathit{/}{\mathit{F}}_{\mathrm{Hcorr}}\mathrm{\left(}\mathit{i,e,\alpha }\mathrm{\right)}\mathrm{=}\mathit{I}\mathit{/}{\mathit{F}}_{\mathrm{meas}}\frac{\mathit{I}\mathit{/}{\mathit{F}}_{\mathrm{Hapke}}\mathrm{\left(}{\mathrm{0}}^{\mathrm{◦}}\mathit{,}{\mathrm{0}}^{\mathrm{◦}}\mathit{,}{\mathrm{0}}^{\mathrm{◦}}\mathrm{\right)}}{\mathit{I}\mathit{/}{\mathit{F}}_{\mathrm{Hapke}}\mathrm{\left(}\mathit{i,e,\alpha }\mathrm{\right)}}\mathit{,}$(2)as in Fornasier et al. (2015), where i, e, and α are the observed incidence, emergence, and phase angles for each pixel, respectively. The resulting I/F_Hcorr represents the normal albedo, while I/F_Hapke is the modeled reflectance as from Hapke (2002), including the macroscopic roughness correction described in Hapke (1984), and is given by $\begin{array}{ccc}\mathit{I}\mathit{/}{\mathit{F}}_{\mathrm{Hcorr}}\mathrm{\left(}\mathit{i,e,\alpha }\mathrm{\right)}\mathrm{=}\frac{\mathit{w}}{\mathrm{4}\mathit{\pi }}\frac{{\mathit{\mu }}_{\mathrm{0}}}{{\mathit{\mu }}_{\mathrm{0}}\mathrm{+}\mathit{\mu }}\mathrm{·}[{\mathit{P}}_{\mathrm{HG}}\mathrm{\left(}\mathit{g}\mathrm{\right)}{\mathit{B}}_{\mathrm{SH}}\mathrm{\left(}\mathit{B}{\mathrm{0}}_{\mathrm{SH}}\mathit{,}{\mathit{h}}_{\mathrm{SH}}\mathrm{\right)}\mathrm{+}\mathit{M}]\mathrm{·}& & \\ {\mathit{B}}_{\mathrm{CB}}\mathrm{\left(}\mathit{B}{\mathrm{0}}_{\mathrm{CB}}\mathit{,}{\mathit{h}}_{\mathrm{CB}}\mathrm{\right)}\mathit{S}\mathrm{\left(}\mathit{\theta ̅}\mathrm{\right)}& & \end{array}$(3)where μ and μ₀ are the effective cosine of the emergence and incidence angle, respectively (Hapke 1984), P_HG(g) is the single-term Henyey-Greenstein function, B_SH(B0_SH,h_SH) is the shadow hiding opposition effect term, B_CB(B0_CB,h_CB) is the coherent backscattering opposition effect term, and M is the multiple scattering term, all described in Hapke (2002), while is the roughness correction term Hapke (1984). This is a function depending on seven parameters: the single-scattering albedo (w), the phase function cosine asymmetric factor (g), the shadow-hiding amplitude and angular width (B0_SH and h_SH), the macroscopic roughness (), and the coherent backscattering amplitude and angular width (B0_CB and h_CB). They have been fitted globally by Fornasier et al. (2015) for 67P and are reported in Table 2.

The computation of i, e, and α for each pixel required a high-resolution shape model of the surface to correctly take the local topography into account. For this application we used the SPG model (Preusker et al. 2015) and in particular the landing-site-dedicated shape model subset described in Sect. 2.

The SPG shape model subset was then used together with the most recent NAIF-SPICE kernels and its IDL code toolkit (Acton Jr 1996) to produce high-resolution synthetic images and maps of incidence, emergence, and phase angles. The synthetic images showed an offset of a few pixels with respect to original images, probably due to small errors in the pointing of the spacecraft and/or small differences between the geodetic datum of the shape model reference frame compared to that of the recent SPICE kernels we used. The perfect alignment between original and synthetic images was achieved with a cross-correlation semi-automatic method looking for χ² minima. We excluded all the non-illuminated (i> 90°) and the non-observable (e> 90°) regions from the analysis.

The resulting Hapke-corrected albedo map in orange filter is shown in Fig. 18 for a composite of the same images as were used in Fig. 2. From the comparison with Fig. 2 (uncorrected reflectance map in the same filter), the effect of photometric correction is evident. The topographical features are clearly reduced and the surface reflectance appears more uniform, even though the boundaries of shadowed regions remain visible. This is expected because we modeled a real surface, which might deviate from a theoretical behavior, also locally.

The histogram of the albedo distribution was best-fit by a Gaussian function with a peak at 0.067 and a standard deviation of 0.01. Thus the Agilkia surface displays an average albedo of 6.7% with variations of ~15%. The average albedo is slightly higher than the overall nucleus derived by Fornasier et al. (2015) (6.4%), while the variations are smaller than the global nucleus albedo variations. The slightly higher albedo with respect to the average nucleus is interpreted as mainly due to grain size effects. As described in Sect. 3, this area is indeed mostly smooth and covered by fine-grained particles. Small particles have an extended scattering surface, and the resulting visible surface brightness is higher than that of a surface made of a mixture of smooth terrains, outcrops, and consolidated materials as the global the nucleus measured by Fornasier et al. (2015).

Fig. 18 Composite albedo map of Agilkia in orange filter (649.2 nm) obtained with the Hapke photometric correction.

Since Hapke parameters were not available for the other filters and the phase coverage of our dataset did not allow an independent Hapke modeling, we were forced to used other photometric models for spectrophotometric investigations.

Other widely used photometric models approximate the reflectance as the product of a term depending on the illumination and viewing conditions (the disk function), and a term only depending on the phase angle (the phase function; Kaasalainen et al. 2001; Shkuratov et al. 2011).

Following this method, the radiance factor can be described by $\mathit{I}\mathit{/}\mathit{F}\mathrm{\left(}\mathit{\alpha ,i,e,\lambda }\mathrm{\right)}\mathrm{=}{\mathit{A}}_{\mathrm{eq}}\mathrm{\left(}\mathit{\alpha ,\lambda }\mathrm{\right)}\mathit{D}\mathrm{\left(}\mathit{i,e,\alpha ,\lambda }\mathrm{\right)}$(4)where A_eq is the equigonal albedo (Shkuratov et al. 2011), defined as the radiance factor for i = α/ 2 and e = α/ 2, while D is the disk function, which describes how the reflectance varies over the surface at a constant phase angle.

Given the small phase angle coverage of our data, we decided to disregard the phase function effect and only considered the disk function. The simplest disk function is the Lambert law for isotropically scattering surfaces, but this is not well suited for airless bodies (Schröder et al. 2013). Widely used functions are the Lommel-Seeliger law, adopted by Fornasier et al. (2015) for the global photometric correction of 67P, its variant with a Lambert term, and the Minnaert law (Minnaert 1941), which introduce one free parameter (see Shkuratov et al. (2011) for a complete review of the common disk functions). The Akimov disk function (Akimov 1976, 1988; Shkuratov et al. 1994; Shkuratov et al. 2003) and in particular its parameter-free version was derived theoretically for an extremely rough surface slightly randomly undulated (Shkuratov et al. 2011). All of the cited disk functions have been tested for the photometric correction of surface of asteroid Vesta (Schröder et al. 2013), and it was found that the parameterless Akimov model performs almost as well as the parametrized functions with the advantage that it has no parameters.

We tested the Lommel-Seeliger law and the Akimov parameterless law on the orange filter and compared the resulting photometric corrected maps with the Hapke-corrected albedo map to determine the best approximation to the Hapke model. Since the disk functions only correct for incidence and emergence angle, but not for the phase angle, we used as a reference for the Hapke model I/F_Hapke(0°,0°,α), where α is the observed average phase angle of the images, 64.23°. The resulting histograms of the reflectance distribution for the three models are reported in Fig. 19 together with the histogram of the measured reflectance. The figure shows that while the Lommel-Seeliger histogram is slightly shifted with respect to the Hapke model and has a higher peak brightness, the Akimov parameterless function yields the best approximation to the Hapke model: both have a peak intensity difference lower than 0.7% and a standard deviation difference lower than 3%. The relative residuals of the Hapke and Akimov corrections are also displayed in Fig. 20 and show that only 10% of the pixels have relative residual higher than 0.12 and that the average residual is lower than 0.05.

Fig. 19 Histograms of the reflectance distribution for various photometric corrections applied.

Fig. 20 Plot of residuals of the Hapke (0°,0°,α) and Akimov corrections. 90% of the pixels, indicated by the dashed line, have a relative residual lower than 0.12, the average residual over all the images, indicated by the red solid line, is lower than 0.05.

On the basis of this result and since it has the advantage of having no parameters, we adopted the parameterless Akimov function to derive the photometry-corrected reflectance in all the available filters. The disk function model we used is therefore given by $\mathit{D}\mathrm{\left(}\mathit{\alpha ,\beta ,\gamma }\mathrm{\right)}\mathrm{=}\cos \frac{\mathit{\alpha }}{\mathrm{2}}\cos \left(\frac{\mathit{\pi }}{\mathit{\pi }\mathrm{-}\mathit{\alpha }}\left(\mathit{\gamma }\mathrm{-}\frac{\mathit{\alpha }}{\mathrm{2}}\right)\right)\frac{\mathrm{(}\cos \mathit{\beta }{\mathrm{)}}^{\mathit{\alpha }\mathit{/}\mathrm{(}\mathit{\pi }\mathrm{-}\mathit{\alpha }\mathrm{)}}}{\cos \mathit{\gamma }}\mathit{,}$(5)where β and γ are the photometric latitude and longitude, respectively, which in turn depend on the incidence, emergence, and phase angles as follows: $\begin{array}{ccc}\cos \mathit{i}& \mathrm{=}& \cos \mathit{\beta }\cos \mathrm{(}\mathit{\alpha }\mathrm{-}\mathit{\gamma }\mathrm{)}\\ \cos \mathit{e}& \mathrm{=}& \cos \mathit{\beta }\cos \mathit{\gamma }\mathit{.}\end{array}$The reflectance was then photometrically corrected for a pixel-by-pixel Akimov disk function $\mathit{I}\mathit{/}{\mathit{F}}_{\mathrm{Acorr}\mathit{,\lambda }}\mathrm{=}\frac{\mathit{I}\mathit{/}{\mathit{F}}_{\mathrm{meas}\mathit{,\lambda }}}{\mathit{D}\mathrm{\left(}\mathit{\alpha ,i,e,\lambda }\mathrm{\right)}}\mathit{,}$(8)to eliminate the effects of topography, illumination, and observing conditions.

Similarly to Walter et al. (2012), we investigated the standard deviation to assess the improvement of the photometric correction. The average uncorrected reflectance in orange filter over the illuminated region is I/F_meas = 7.6 × 10^-3 with a standard deviation of 1.8 × 10^-3, while the average corrected reflectance is I/F_Acorr = 9.6 × 10^-3 with a standard deviation of 0.8 × 10^-3. The standard deviation is reduced by more than a factor of 2 by the Akimov photometric correction, while the average corrected reflectance is slightly increased. This increase is due to the apparent effect of flattening the high-slope terrains, performed by the photometric correction, that otherwise reflect a much lower fraction of the incident light. A more detailed evaluation of the photometric correction improvement was performed by investigating the relation between the reflectance and the incidence and emergence angles (see Fig. 21). The measured (top panels) and corrected (bottom panels) reflectance of each pixel of the NAC image in set 4 of Table 1 is plotted as an example (the analysis was performed for all cubes) as a function of its emergence (left panels) and incidence angles (right panel). The measured reflectance shows a linear dependence on both angles, as expected, since high incidence and emergence angles cause the surface to reflect only a fraction of the light that would reflect if illuminated and observed at opposition (i = e = 0). The Akimov-corrected reflectance instead shows a much smaller scatter around its average value independent of the specific illumination and observing conditions.

Fig. 21 Evaluation of the photometric correction improvement. In the top panels the I/Fmeas of the NAC image in orange filter of set 4 in Table 1 is plotted as an example as a function of emergence (left) and incidence (right) angles in a density slice plot. I/Fmeas shows a linear dependence on both angles. The linear fit is overplotted. In the bottom panels the Akimov-photometrically corrected I/FAcorr of the same image is plotted as a function of the same angles. The scatter is reduced now to only the intrinsic reflectance variations, and the linear dependence is corrected, being I/FAcorr equally distributed at all incidence and emergence angles. The average value and ± σ scatter width is indicated by the solid and dashed lines, respectively. Colored contours indicate the density slices for values from 10 to 90 bins with intervals of 10. Angles bins are of 1° , while reflectance bins are of 0.001.

5. Photometric variegation

The variegation observed across a few pixels was ignored because it is affected by small topographical features that at such high observing phase angle might cause small artifacts. In particular on rough materials that are crowded with boundaries and small-scale features, the smallest uncertainty in the shape model and/or in the alignment may be relevant. We therefore instead concentrated on large-scale variations, which are accurately corrected for by the shape model. A large-scale albedo variation of about 15% is indeed detectable across the surface. This was also compared with the incidence and emergence angle maps, and the absence of clear correlation precludes that it is caused by geometrical conditions. It is therefore interpreted to be due to an actual photometric variation of the material composing the surface. The sole photometric variation might be due to both composition variegation and/or texture and grain size effects.

A few very bright spots or patches with a reflectance of up to 2% were identified in this area as well as throughout the comet and were interpreted as boulders probably consisting of ice (Pommerol et al. 2015). We could not confirm the presence of bright spots on Agilkia from the photometry-corrected reflectance map since the features, which are only a few pixels wide, might suffer much larger uncertainties. For these features the measured reflectance is probably the most reliable measure, at least until a much more detailed DTM of the area will be released at much higher resolution. There are a few bright spots mostly over the consolidated material around the rims of Hatmehit, in the Bastet region, and on Maftet (see examples in Fig. 22).

The comparison of the reflectance map with the geomorphological map in Fig. 5 allowed retrieving a general correlation between the reflectance variation and the morphological units and textures defined in Sect 3. The fine-particle deposits show the highest reflectance with little variations across the region. The other deposits from finer to coarser materials, that is, from taluses to gravitational accumulations and diamicton deposits, show a gradually decreasing brightness. Finally, the outcrops and consolidated materials are associated with the lowest reflectance. The geomorphology-photometric correlation suggests a dependence on the texture of the material. Assuming a homogeneous composition, the fine particles reflect the light much more than grainy materials and rocks because they have a greater reflecting surface, as confirmed for example by experiments on carbonaceous chondrites of different granulometry by Johnson & Fanale (1973) and Clark (1999). It is therefore not surprising that smooth deposits are brighter than outcrops. However, a composition variegation cannot be excluded from the reflectance analysis alone.

Fig. 22 Bright patches and boulders detected in a) the Maftet region, b) over the rims of Hatmehit, and c) in the Bastet region. NAC images acquired on a) 2014-10-19 at 13:19:53, b) 2014-09-14 at 17.31:36, and c) 2014-09-15 at 05:42:24.

6. Spectrophotometric analysis

We investigated the surface color using the spectral cubes in Table 1. We analyzed several colors as filter ratios for the different combination of filters in the dataset. However, since the average spectrum of the region is featureless with a constant red slope at the observed wavelengths (see next paragraph), the filter ratios are only diagnostic of the spectral slope. The spectral slope map between near-IR (882.1 nm) and blue (480.7 nm) was used as representative of the spectral behavior. The spectral slope, expressed in %/100 nm, computed as (Jewitt & Meech 1986; A’Hearn et al. 1984) $\mathit{S}\mathrm{=}\frac{\mathit{I}\mathit{/}{\mathit{F}}_{\mathrm{NearIR}}\mathrm{-}\mathit{I}\mathit{/}{\mathit{F}}_{\mathrm{Blue}}}{{\mathit{\lambda }}_{\mathrm{NearIR}}\mathrm{-}{\mathit{\lambda }}_{\mathrm{Blue}}}\mathrm{·}\frac{\mathrm{20}\hspace{0.17em}\mathrm{000}}{\mathit{I}\mathit{/}{\mathit{F}}_{\mathrm{NearIR}}\mathrm{+}\mathit{I}\mathit{/}{\mathit{F}}_{\mathrm{Blue}}}\mathit{,}$(9)where I/F_NearIR and I/F_Blue are intended as the Akimov-corrected reflectance in the two filters, is shown in Fig. 23. The color bar was chosen so that regions with a lower slope than average are colored in blue, and regions where it is higher than average are red.

Fig. 23 Composite spectral slope map computed with the formula in Eq. (9). Unobservable (e> 90°) and non-illuminated regions (i> 90°) are removed from the data and are represented in white in the map. Regions where near-IR/blue is lower than average are colored in blue, and regions where the ratio is higher than average are red.

The average spectral slope across the area is 15.2%/100 nm at the observed phase angle of about 64°, consistent with the average spectral slope of the nucleus of 67P at similar phase angles measured by Fornasier et al. (2015). The entire region shows a significantly steep spectral slope, but the eastern rim of Hatmehit is the reddest regions on Agilkia, while the smoother area on Ma’at toward the north pole shows lower slopes. The spectral slope map agrees with spectral maps obtained by Fornasier et al. (2015) and by the VIRTIS instrument (Capaccioni et al. 2015; Filacchione et al. 2015). Based on their results and ours, we might speculate that this type of neutral-to-red slope is compatible with the typical slope of primitive organic materials (e.g., Moroz et al. 2003 and references therein).

A significant spectral slope variation of more than 5% is visible across the surface.

The comparison of the spectral slope map with the geomorphological units allowed us to identify an existing correlation of the morphology also with colors in addition to reflectance. The fine material deposits show the bluest colors within the area. From fine to rough material, the colors of the map gradually change to red, indicating material with a steeper spectrum. The outcrops and consolidated terrains show the reddest colors. However, some color differences are additionally observed within the same geological unit. The fine material deposits in particular show the highest variability, changing from 13% to about 14.5% going from the northern boundary of Ma’at to the equator.

To better investigate the correlation between geomorphological and photometric properties of the surface area, we distinguished the different morphological terrains, in particular, the fine dust deposits, the taluses, the gravitational accumulations, the diamictons, and the outcropping competent material (as in the map in Fig. 5), on the spectral cubes described in Table 1. For each type of terrain, a wide set of 20 to 30 low-resolution spectra was obtained, spread across the unit regions of about 12 square meters, to have a good spatial sampling and avoid small-scale topographic features.

For each morphological terrain we applied a resistant mean, excluding values exceeding ±3σ with respect to the average, producing a of low-resolution spectrum representative of each terrain type. The resulting spectra for all terrains are shown in Fig. 24, together with the average spectrum of the entire region.

Fig. 24 Absolute (top panel) and normalized (bottom panels) low-resolution spectra obtained with the Akimov-corrected reflectance in all available filters on different types of terrains, together with the overall average spectrum. Error bars in the normalized spectra are inside the symbol dimension.

The area globally shows a featureless spectrum (except for a small inflection at 535.7 nm discussed below), without evident absorption bands and a neutral-to-red spectral slope. The absence of distinctive features makes any compositional interpretation of the visible spectrum uncertain, because it might be compatible with a wide range of materials. However, the VIRTIS instrument allowed the detection of a wide infrared band in the 2.9−3.6 μm range across the entire illuminated surface (Capaccioni et al. 2015). Capaccioni and collaborators inferred a composition of refractory organic materials, such as polyaromatic organic solids mixed with sulfides and Fe-Ni alloys (Quirico et al. 2015). Our visible featureless red spectra agree with the VIRTIS conclusions.

Figure 24 shows the absolute spectra of individual terrains (top panel) and the normalized spectra at 535.7 nm (green filter). The spectra confirm the observed correlation with the morphological units, showing a slightly diverse slope trend. The fine deposits are the highest reflectance terrains in all filters and show the more shallow, that is, bluest, spectrum, immediately followed by taluses and gravitational accumulations, which are slightly darker and redder. The diamictons have the widest spread in spectral behavior, as indicated by the largest error bars, which is expected considering the wide variety of textures and source origins. However, they perfectly agree with the average reflectance of the area and are placed just between the finer deposits and the outcrops. The latter are the darkest and reddest material in the area. The strongest variation in slope is about 5%, but the gradual relation with morphological terrains validates this variegation. As pointed out previously, almost all the spectra in Fig. 24 show a slight inflection in the green filter at about 535 nm with a depth of about 2−3%. A possible interpretation could invoke an absorption feature due to charge transfer of Fe3+ ion in silicates (Lazzarin et al. 2009). However, the small amplitude of this inflection is not sufficient to confidently assess that the spectrum shows an absorption. This conclusion is also supported by the fact that none of the overall nucleus spectra studied in (Fornasier et al. 2015) show any indication of such an absorption.

7. Discussion and conclusions

We identified different typologies of terrains on the landing site of Philae on comet 67P. Part of the region is covered by finely grained deposits, while outcrops of more consolidated materials are visible that are cuestas and steep walls, which often exhibit terraces and strata heads, representative of underlying layered terrains (Massironi et al. 2015). At the foot of the walls, various deposits are visible of both fine and coarser textures, which were identified as taluses and gravitational accumulations on the basis of the local gravitational slopes and their textures. Poorly sorted deposits are also present across the area, showing fine material mixed with large boulders of unclear origin. The smooth deposits show dune-like features with a uniform trend of their longest axis all around Agilkia, suggesting a constant dominant wind regime if a gas-flow transport is invoked as the cause for the origin of these features. Chains of pit-like features, generally aligned with the dune-like features, were observed, and a sublimation origin was suggested on the basis of the solar elevation and angular separation with respect to these features along the orbital period of 67P and in particular at perihelion. The boulder size-frequency distribution analysis across the area agrees with a mostly flat terrain, suggesting a low activity and small fracturing of the region with respect to other areas on the nucleus. NAC observations in the period 1 Aug.−12 Nov. do not show any evidence of activity, which supports this result.

The average normal albedo of the site is 6.7%, consistent with the general low albedo of cometary nuclei, but slightly higher than the albedo of the overall nucleus found by Fornasier et al. (2015). This small difference was interpreted in terms of grain size effects. The average spectrum of the area shows a featureless behavior with a significant red slope of about 15.2%, in agreement with the high phase reddening observed by Fornasier et al. (2015). The featureless visible spectrum is compatible with an organic-rich composition as hypothesized from VIRTIS instrument data that display a deep 3 μm band that was interpreted as due to organics (Capaccioni et al. 2015; Filacchione et al. 2015). We observed a significant variation in albedo of about 15% and a weak variation in color and spectral slope slightly above the uncertainties. However, the authenticity of photometric differentiation across the area is supported by the very good correlation with the morphological map and by the careful statistical analysis performed on the spectra. The finest deposits show the highest reflectance and the bluest spectra with respect to the darker and redder consolidated material. The geological units with intermediate textures are placed between these two extremes.

Refractories usually show a significant reddening at decreasing particle size (e.g., Johnson & Fanale 1973; Clark 1999). The color investigation performed on the Baetica region of asteroid Lutetia, for example (Magrin et al. 2012), shows an association between fine materials and higher spectral slopes. This is valid in particular for grain sizes comparable with the observed wavelength range. ROLIS images (Mottola et al. 2015b) revealed, however, that particles constituting the fine deposits at the location of the first touchdown on Ma’at are mostly granular, with a millimeter-size distribution. In this particle-size regime it is not possible to spectroscopically distinguish smooth deposits from outcrops of the same composing material. We therefore interpreted the observed spectral slope variation as an indication of a chemical composition variegation. The lower steepness of the spectrum in the blue regions in Fig. 23 might be tentatively interpreted as due to a small fraction of neutral-to-bluer spectral components, which are able to flatten the visible spectral slope. The general bluer color of the neck-Hapi region has similarly been interpreted by Fornasier et al. (2015), indicating the presence of water ice that has been recognized in transient icy patches through a VIRTIS infrared spectrum of small regions close to the shadows (De Sanctis et al. 2015).

Another possibility is that fine particles and outcrops have been exposed to different amounts of space weathering, which altered their spectral appearance. Laboratory measurements by (Moroz et al. 2004) showed that in certain conditions, organic material may show a decrease in the steepness if exposed to ion bombardment that simulates space weathering. Because of the relatively fast resurfacing to which comets are subjected during perihelion passages and considering the short orbital period of 67P, this hypothesis is relatively unlikely.

If the spectral variegation is therefore actually due to a chemical compositional variegation, this in turn suggests that the source of the finest material is not only the sublimation and consequent fragmentation of local consolidated competent material. It seems rather that the fine grains may have a different source region. The efficiency in dust-transport processes acting on the comet suggests that a continuous mixing of material occurs across the surface of the comet due to activity and jets. Even though Agilkia does not show any sign of activity up to date, strong jets have instead been observed (Lin et al. 2015) coming from the neck-Hapi region on the comet, which shows a generally bluer color than the rest of the nucleus (Fornasier et al. 2015). Simulations in Thomas et al. (2015b) showed that this activity would be extensively able to transport the necessary amount of dust to explain all the dust-covered regions on the nucleus of 67P. The smooth deposits visible on Agilkia may therefore easily be the result of airfall of material produced by jets and originate from various regions on the nucleus, including the bluer Hapi-neck region. If the fine material comes from higher activity regions and it is enriched in neutral-to-blue slope component, the redder slope of the steep walls and margins that is usually connected with higher gravitational slopes, such as Hatmehit rim, might therefore be interpreted as due to the incapacity of the dust particles to rest on these sharp terrains, while flatter regions are more suited to detain fine deposit accumulation.

The upcoming observations that Rosetta will perform in the next months, as the comet approaches the Sun and the solar heat gradually changes its appearance, will allow a deeper investigation of the discussed hypotheses and will furthermore reveal how the activity influences the changing surface. More highly resolved observations with high temporal resolution will uncover the main processes that shape and alter the surface, and they will be used in future analyses.

Acknowledgments

OSIRIS was built by a consortium of the Max-Planck-Institut für Sonnensystemforschung, in Göettingen, Germany, the CISAS-University of Padova, Italy, the Laboratoire d’Astrophysique de Marseille, France, the Instituto de Astrofísica de Andalucia, CSIC, Granada, Spain, the Research and Scientific Support Department of the European Space Agency, Noordwijk, The Netherlands, the Instituto Nacional de Técnica Aeroespacial, Madrid, Spain, the Universidad Politécnica 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. We thank the ESA teams at ESAC, ESOC and ESTEC for their work in support of the Rosetta mission. This research has made use of the USGS Integrated Software for Imagers and Spectrometers (ISIS). We gratefully acknowledge the developers of SPICE and NAIF/PDS resources.

References

All Tables

All Figures

	Fig. 1 OSIRIS NAC images acquired on 12 November 2014 at a) 15:18 and b) 15:43 just after the first Philae touchdown, showing evidence of the footprints of Philae’s three legs.
In the text

	Fig. 2 Composite image of Philae landing site (Agilkia) obtained from five NAC frames in Orange filter (sets 1−5 in Table 1) acquired on 14−15 September 2014 at an average phase angle of 64.23°. The circle highlights the location of the first touchdown, while the dashed line indicates the equator.
In the text

	Fig. 3 Gravitational slope map of Agilkia obtained using the SPC shape model downsized to 1.5 × 105 points and the gravitational potential model derived assuming a homogeneous nucleus and accounting for the centrifugal force due to the nucleus rotation.
In the text

	Fig. 4 Geological map of 67P observed from the top of the head lobe and nomenclature of the geological regions (modified from Fig. 2 El-Maarry et al. 2015).
In the text

	Fig. 5 Geomorphological map of Agilkia identifying the main geological units on the basis of their morphological properties, and various linear features. The map is produced on the composite image in Fig. 2.
In the text

	Fig. 6 Examples of longitudinal dunes observed in Ma’at and Hatmehit. NAC image acquired on 2014-09-14 at a) 17:47:55 and b) 16:56:07 in orange filter.
In the text

Fig. 7 Examples of aligned pit-like features visible on smooth deposits in the Ma’at region. NAC images acquired on 2014-09-30 at 09:26:53 (left panel) and on 2014-10-19 at 12:22:15 (right panel) in orange filter. Enclosed by the green dashed line, we show the region for which the local horizon has been extracted, and the red dashed line denotes the considered orientation of the features for comparison with the solar orientation (see text for details).

In the text

Fig. 8 Solar orientation with respect to the aligned pit features. Top panels show the elevation of the Sun over the local horizon, while the bottom panels show the angular separation between the orientation of the features (red dashed line in Fig. 7) and the projected solar orientation on the local horizon. Both angles are given as a function of the heliocentric distance (left panels) and time (right panels) for a period covering one full comet orbit with steps of 30 days for revolution and 30 min for rotation (see text for details).

In the text

	Fig. 9 Examples of outcrops of consolidated material. NAC images acquired in orange filter on a) 2014-09-14 at 17:31:24, b) 2014-09-14 at 04:28:01, c) 2014-11-12 at 17:08:52, and d) 2014-11-12 at 16:53:51.
In the text

Fig. 10 Cuesta morphology a) sketch of cuesta landform. Cuestas occur in areas of tilted strata and are caused by the differential erosion of the harder top layer and the underlying layer that is weaker and therefore more erodible (figure after Castiglioni 2000). b)−c) Examples of cuestas on the Agilkia site (indicated by black arrows). The hard layer probably consists of material poor in volatiles, whereas the soft underlying layer is richer in volatile material and therefore more erodible by sublimation. NAC image acquired onb) 2014-09-14 at 17:32:15 in near-IR filter and on c) 2014-11-12 at 16:08:51 in orange filter.

In the text

	Fig. 11 Examples of fractures (in cyan) and strata heads (in green) detected on the Agilkia site. Usually fractures cut strata heads, revealing that they form after the layering. NAC image acquired on 2014-11-12 at 16:23:51.
In the text

	Fig. 12 Different types of deposits. a) Example of talus deposit. b) Example of gravitational accumulation deposit. NAC image acquired on 2014-09-14 at 05:05:28 in orange filter.
In the text

	Fig. 13 Spatial distribution of boulders on the Ma’at, Maftet, and Nut parts of Agilkia. NAC image acquired on 2014-09-14 at 17:48:08 in orange filter.
In the text

	Fig. 14 Cumulative size-frequency distribution for the total 2065 boulders identified on Maftet, Ma’at, and Nut inside Agilkia.
In the text

	Fig. 15 Hatmehit region observed from the top and front. NAC image acquired from on 2014-09-02 at 07:42:53 (left panel) and on 2014-09-01 at 08:42:52 (right panel) in orange filter.
In the text

	Fig. 16 Examples of terraces detected on the Hatmehit depression rim (indicated by white arrows in panels a) and b)). NAC image acquired on 2014-09-14 at 16.56.07.
In the text

	Fig. 17 Diamicton deposit in Hatmehit. NAC image acquired on 2014-09-14 at 16:56:07 in orange filter.
In the text

	Fig. 18 Composite albedo map of Agilkia in orange filter (649.2 nm) obtained with the Hapke photometric correction.
In the text

	Fig. 19 Histograms of the reflectance distribution for various photometric corrections applied.
In the text

	Fig. 20 Plot of residuals of the Hapke (0°,0°,α) and Akimov corrections. 90% of the pixels, indicated by the dashed line, have a relative residual lower than 0.12, the average residual over all the images, indicated by the red solid line, is lower than 0.05.
In the text

Fig. 21 Evaluation of the photometric correction improvement. In the top panels the I/Fmeas of the NAC image in orange filter of set 4 in Table 1 is plotted as an example as a function of emergence (left) and incidence (right) angles in a density slice plot. I/Fmeas shows a linear dependence on both angles. The linear fit is overplotted. In the bottom panels the Akimov-photometrically corrected I/FAcorr of the same image is plotted as a function of the same angles. The scatter is reduced now to only the intrinsic reflectance variations, and the linear dependence is corrected, being I/FAcorr equally distributed at all incidence and emergence angles. The average value and ± σ scatter width is indicated by the solid and dashed lines, respectively. Colored contours indicate the density slices for values from 10 to 90 bins with intervals of 10. Angles bins are of 1° , while reflectance bins are of 0.001.

In the text

	Fig. 22 Bright patches and boulders detected in a) the Maftet region, b) over the rims of Hatmehit, and c) in the Bastet region. NAC images acquired on a) 2014-10-19 at 13:19:53, b) 2014-09-14 at 17.31:36, and c) 2014-09-15 at 05:42:24.
In the text

	Fig. 23 Composite spectral slope map computed with the formula in Eq. (9). Unobservable (e> 90°) and non-illuminated regions (i> 90°) are removed from the data and are represented in white in the map. Regions where near-IR/blue is lower than average are colored in blue, and regions where the ratio is higher than average are red.
In the text

	Fig. 24 Absolute (top panel) and normalized (bottom panels) low-resolution spectra obtained with the Akimov-corrected reflectance in all available filters on different types of terrains, together with the overall average spectrum. Error bars in the normalized spectra are inside the symbol dimension.
In the text