## Abstract

VIRTIS-M observations of the nucleus of comet 67P/Churyumov–Gerasimenko acquired from 2014 August to 2015 May have been analysed to investigate surface temporal variability at both seasonal and diurnal scales. The measured reflectance spectra are studied by means of comet spectral indicators (CSI) such as slopes in the visible and infrared ranges, and 3.2 μm band area and band centre. CSI maps derived from data acquired at different heliocentric distances (from 3.62 to 1.72 au) along the inbound leg of the comet's orbit are used to infer surface water ice abundance. We measure a global scale enrichment of water ice from 2014 August to 2015 May across the body of the comet, along with variability at small spatial scale, possibly related with the local insolation conditions. Analysis of water ice diurnal variability is performed on 2014 August observations. Water ice appears at the border of receding shadows in the neck of the comet (Hapi), sublimating in less than 1 h, after exposure to sunlight. As similar variability is not observed in other regions of the comet, we interpreted this as the expression of a diurnal cycle of sublimation and re-condensation of water ice, triggered by sudden shadowing produced on the neck by the body and the head of the nucleus.

## 1 INTRODUCTION

Launched in 2004, the Rosetta spacecraft encountered the Jupiter family comet 67P/Churyumov–Gerasimenko on 2014 August 6. The first observations from Rosetta instruments reveal an elongated nucleus, characterized by two major lobes (head and body), connected by a smoother region (neck) appearing as the main source of the comet activity (Sierks et al. 2015). After several missions that performed direct imaging of the nuclei of comets 1P/Halley, 9P/Tempel 1, 19P/Borrelly, 81P/Wild 2 and 103P/Hartley 2, Rosetta is the first to rendezvous with and escort a comet along its orbit. This allowed a continuous monitoring of 67P's activity and surface changes from a heliocentric distance of 3.6 au down to 1.24 au at perihelion, and then out again as the comet receded from the Sun.

Outgassing of volatile materials and re-condensation processes triggered by the insolation conditions (De Sanctis et al. 2015), while determining cometary activity, also continuously modify a comet's surface and its composition as well as its rotation state. Keller et al. (2015) have shown that the seasonal variability of the thermal input on 67P determines a non-homogeneous erosion of the surface, producing the observed morphological dichotomy between the north and the south hemispheres, and a variation of the rotation state. Such a dichotomy is revealed also in the surface spectral properties (Ciarniello et al. 2015; Filacchione et al. 2016a) and in the measured water and carbon dioxide distribution in the inner coma, possibly related to the different thermal history of the north and south hemispheres (Bockelée-Morvan et al. 2015; Fink et al. 2016; Migliorini et al. 2016). Groussin et al. (2015) reported morphological variability on a temporal scale of less than two months in 67P's Southern hemisphere (Imhotep region), compatible with exposure of ice rich areas. Albedo and morphological variations at small scale have been also inferred on the surface of comet Tempel 1 from the comparison of Deep Impact (Thomas et al. 2007) and Stardust images (Veverka et al. 2013). Besides temporal variability at medium and large temporal scale, results from the Rosetta mission have also evidenced diurnal effects on surface composition and activity. In particular, De Sanctis et al. (2015) measured cyclic variations of water ice content in the 67P Hapi region, driven by the light and shadow alternation during one comet rotation, while Shi et al. (2016) reported of sunset jets, sustained for about one hour after dusk, by the thermal lag in the subsurface. Moreover, Hässig et al. (2015) have shown that H2O, CO and CO2 gas emissions are characterized by diurnal variability.

The scope of this work is to investigate 67P surface composition variability both on diurnal and seasonal temporal scales, focusing on the water ice content and its distribution on the nucleus. Concerning the short-term compositional variability, we complement the work of De Sanctis et al. (2015) about the above-mentioned water ice cycle, extending the analysis to the whole observed surface of the comet, while, in the context of seasonal effects, we complete the analysis carried out in Filacchione et al. (2016a) including observations up to 2015 May.

To this aim, we used observations from the VIRTIS (Visible Infrared Thermal Imaging Spectrometer) instrument (Coradini et al. 2007) onboard the Rosetta mission. VIRTIS is an imaging spectrometer composed of two subunits, a mapping spectrometer [Visible Infrared Thermal Imaging Spectrometer (VIRTIS)-M] and a high-resolution spectrograph (VIRTIS-H). VIRTIS-M, whose data have been used in this work, is composed of two channels covering, respectively, the visible 0.24–1.04 μm (VIS) and the infrared 1–5.1 μm (IR) spectral ranges, with an angular resolution of 250 μrad × 250 μrad.

The data set used in this study is presented in Section 2, while a description of the methods adopted in this analysis is given in Section 3. Two separate sections, Sections 4 and 5, describe the results obtained from the analysis of seasonal and diurnal effects, respectively. A summary of the main findings is presented in Section 6.

## 2 DATA SET AND OBSERVATIONS

In this study, we exploit the large data set produced by the VIRTIS-M instrument onboard the Rosetta mission from 2014 August, when we performed the first disc-resolved acquisitions of the nucleus surface, until 2015 May, when the IR stopped operating due to a failure of the cryocoolers needed to operate the IR sensor. This time period corresponds to medium-term planning (MTP) observation phases going from MTP006 to MTP015 during which 977 hyperspectral images were acquired with the spatial resolution varying from 1.9 m pix−1 to 153 m pix−1 and the solar phase angle varying from 14° to 122°. In this time frame, the comet moved along its orbit, in towards the Sun, from an heliocentric distance of 3.62 to 1.72 au. Simultaneously, because of its spin axis orientation (69_{.}^{\circ}$$54 RA, 64 _{.}^{\circ}$$11 Dec, J2000 equatorial coordinates; Preusker et al. 2015), the position of the subsolar point varied from an average latitude of 24_{.}^{\circ}$$5 (MTP006, north hemisphere summer) to 3 _{.}^{\circ}$$7 (beginning of MTP015, equinox on 2015 May 5). This information is summarized in Fig. 1, where the target-spacecraft distance, phase angle, comet heliocentric distance and subsolar point latitude are shown for the different MTP mission phases.

Figure 1.

VIRTIS-M observations data set properties as a function of time. From top to bottom: solar phase angle, heliocentric distance rh, comet–spacecraft distance Δ and subsolar latitude, SS lat. Oscillations of the subsolar latitude with diurnal frequency are due to the irregular shape of the nucleus.

Figure 1.

VIRTIS-M observations data set properties as a function of time. From top to bottom: solar phase angle, heliocentric distance rh, comet–spacecraft distance Δ and subsolar latitude, SS lat. Oscillations of the subsolar latitude with diurnal frequency are due to the irregular shape of the nucleus.

## 3 METHOD

The nucleus of comet 67P, similar to other cometary nuclei (Li et al. 2007a,b, 20092013), is characterized by a very dark surface with a geometric albedo at 0.55 μm of ≈0.06 (Ciarniello et al. 2015; Fornasier et al. 2015). The average spectrum as measured by VIRTIS (Capaccioni et al. 2015; Ciarniello et al. 2015; Quirico et al. 2016), exhibits a red slope across the VIS-IR range with a prominent absorption feature centred at 3.2 μm and extending from 2.9 to 3.6 μm. Quirico et al. (2016) interpreted the red spectral slope and the low albedo as the effect of a mixture of a refractory polyaromatic carbonaceous component and opaque minerals (Fe-Ni alloys or Fe sulphides such as pyrrhotite and troilite), and suggested the COOH group in carboxylic acid or NH$$_4^+$$ ion as the best candidates for the 3.2 μm absorption. Along with these species, stretching bands from C–H in aliphatic and aromatic hydrocarbons and OH in alcohols can provide a further contribution to the wide 3.2 μm feature. The C–H absorptions are not resolved because of possible overlapping with vibration modes from other chemical groups, while alcoholic hydroxyl produces an absorption from 2.7 to 3.2 μm, which, however, is not sufficient to cover the full broadness of the 3.2 μm feature (Quirico et al. 2016).

The most abundant volatile component on the surface of 67P is represented by water ice (De Sanctis et al. 2015; Filacchione et al. 2016b). Water ice has been previously detected by Sunshine et al. (2006) on the surface of comet Tempel 1, making it the principal volatile component on comet surfaces, so far. Thermal evolution models of cometary nuclei indicate the development of a subsurface layering following the sublimation fronts of the different volatile species (Orosei et al. 1995; De Sanctis et al. 2007; Prialnik et al. 2008; Marboeuf et al. 2012). In particular, given its sublimation temperature, H2O ice is expected to be present closer to the comet surface mixed with refractory materials, while the sublimation front of higher volatility species, such as carbon dioxide and carbon monoxide, are expected to be located at greater depths (De Sanctis et al. 2007).

Water ice has been identified on 67P's surface both as exposed isolated patches (Pommerol et al. 2015; Barucci et al. 2016; Filacchione et al. 2016b) and as intimately mixed with the refractory material of the nucleus (De Sanctis et al. 2015). As evidenced by VIRTIS-M observations, the former case pertains to a few isolated small regions, while the large-scale distribution of water ice is more likely consistent with the latter, as indicated by the lack of a widespread distribution of the water ice 2.0 μm absorption feature (Filacchione et al. 2016a).

In order to investigate the compositional variability, we focus on monitoring of the water ice surface content, which is sensitive to changes of insolation conditions on diurnal and seasonal temporal scales (De Sanctis, Lasue & Capria 2010). The amount of water ice mixed with the dark refractory material representing the main component of 67P surface (dark terrain) can be inferred by means of spectral analysis. In particular, the presence of water ice in the mixture produces modifications of the spectral shape with respect to the average spectrum of the nucleus. An increasing amount of ice produces a reduction of the slope across the VIS-IR range, a deepening of the 3.2 μm absorption feature with a shift towards shorter wavelengths and an increase of the albedo. This is shown in Fig. 2(a) (see Appendix A for modelling details) where simulated spectra of mixtures of the 67P dark terrain (Ciarniello et al. 2015) and various amounts of water ice are reported.

Figure 2.

(a) Simulated SSA of 67P dark terrain intimately mixed with various amounts of water ice from 0 to 5 per cent (black to grey). Spectral models are computed by means of Hapke theory (Ciarniello et al. 2011; Hapke 2012; Raponi et al. 2016) assuming 2 μm grain size water ice particles (De Sanctis et al. 2015). SSA of 67P's dark terrain is given in Ciarniello et al. (2015), while for H2O ice it is computed from optical constants by Warren (1984), Mastrapa et al. (20082009) and Clark et al. (2012). The blue and red lines across the spectrum indicate how the spectral slope in the VIS (SVIS) and in the IR (SIR) are computed, respectively. Similarly, the green area in correspondence of the 3.2 μm band evidences the interval selected to compute the band area (BA) and the band centre (BC), which, we recall here, are computed after continuum removal. (b) CSI (comet spectral indicators) normalized at their value computed on 67P dark terrain as a function of water ice abundance: SVIS (diamonds), SIR (asterisks), BA (crosses) and SSA (triangles). BC (squares) is not normalized and its scale is reported on the right axis.

Figure 2.

(a) Simulated SSA of 67P dark terrain intimately mixed with various amounts of water ice from 0 to 5 per cent (black to grey). Spectral models are computed by means of Hapke theory (Ciarniello et al. 2011; Hapke 2012; Raponi et al. 2016) assuming 2 μm grain size water ice particles (De Sanctis et al. 2015). SSA of 67P's dark terrain is given in Ciarniello et al. (2015), while for H2O ice it is computed from optical constants by Warren (1984), Mastrapa et al. (20082009) and Clark et al. (2012). The blue and red lines across the spectrum indicate how the spectral slope in the VIS (SVIS) and in the IR (SIR) are computed, respectively. Similarly, the green area in correspondence of the 3.2 μm band evidences the interval selected to compute the band area (BA) and the band centre (BC), which, we recall here, are computed after continuum removal. (b) CSI (comet spectral indicators) normalized at their value computed on 67P dark terrain as a function of water ice abundance: SVIS (diamonds), SIR (asterisks), BA (crosses) and SSA (triangles). BC (squares) is not normalized and its scale is reported on the right axis.

Such spectral variability can be described by defining a set of comet spectral indicators (CSI), that allow us to describe the shape of the measured spectrum as well as to reduce the data set dimension. This approach, described by Filacchione et al. (2016a), is adopted in this work with a few modifications. The set of spectral indicators used in this analysis are as follows:

• SVIS: Spectral slope in the VIS range, calculated by best fitting the reflectance in the 0.55–0.8 μm range after normalization to 1 at 0.55 μm.

• SIR: Spectral slope in the IR range, calculated by best fitting the reflectance in the 1.0–2.35 μm range and after normalization to 1 at 1μm

• BA: Band area of the absorption at 3.2 μm computed after continuum removal in the 2.7–3.6 μm range.

• BC: Band centre of the absorption at 3.2 μm, computed after continuum removal as the position of the minimum of a third-degree polynomial fit in the 2.7–3.6 μm range.

In Fig. 2(b), the spectral indicators described above are computed for the mixtures of Fig. 2(a) and reported as a function of water ice abundance, after normalization to their value as measured on the average 67P dark terrain spectrum. It can be noted that the most sensitive parameter to water ice abundance is the band area (BA) at 3.2 μm, which doubles when ice reaches 3 per cent, while the VIS and IR slopes are reduced by ≈35 per cent and ≈25 per cent, respectively. For the same amount of ice, the shift of the band centre (BC) is of the order of ≈160 nm towards shorter wavelengths, corresponding to 17 VIRTIS spectral channels. These quantities are computed on VIRTIS spectra after application of the photometric correction (see Appendix A and Ciarniello et al. 2015) and removal of thermal emission (Filacchione et al. 2016a). This procedure allows us to minimize the effect of the observation geometry on the measured spectral shape, such as phase-reddening (Schröder et al. 2014; Ciarniello et al. 20152016), and to retrieve the reflected component of the measured radiance. The set of CSI defined by (Filacchione et al. 2016a) also comprises the single scattering albedo at 0.55 μm (SSA) after photometric correction and the BA of the water ice absorption feature at 2 μm. In our analysis, we do not take advantage of the absolute value of the SSA because of the potentially large systematics errors in its derivation due to unresolved shadows, variable spatial resolution between different observations and unpredictable photometric effects of subpixel topography. Such issues are minor in the case of SVIS, SIR and BA, being computed as normalized quantities (with respect to either a given wavelength or the continuum). We also do not use the 2 μm band to investigate the temporal variability of 67P nucleus properties since it is absent across the most of the surface, with the exception of few small ice-enriched regions (Filacchione et al. 2016a; Raponi et al. 2016).

## 4 SPECTRAL INDICATORS MAPPING AND SEASONAL VARIABILITY

Spectral indicators as computed in Section 3 are mapped in a cylindrical projection with a latitude–longitude resolution of 0_{.}^{\circ}$$5 × 0 _{.}^{\circ}$$5, as was done in Filacchione et al. (2016a). VIRTIS pixels are projected according to their IFOV (instantaneous field of view) and the corresponding latitude and longitude values are computed by means of spice-based routines (Acton 1996) using the digital shape model (SHAP5) of 67P's nucleus (Jorda et al. 2016). Due to the high redundancy of the VIRTIS data set (see Fig. 3), many points on the surface were observed more than once: in these cases, we report on the map the median value derived from all available observations. Pixels in shadow or characterized by extremely large incidence angle or emission angle (i, e, >80°) are rejected.

Figure 3.

Redundancy (number of available VIRTIS observations for a given position) maps for the different observation phases. From top to bottom: MTP006, MTP007-MTP008, MTP009-MTP010, MTP011-MTP012 and MTP014-MTP015. See Figs 4 and 5 for 67P's morphological regions.

Figure 3.

Redundancy (number of available VIRTIS observations for a given position) maps for the different observation phases. From top to bottom: MTP006, MTP007-MTP008, MTP009-MTP010, MTP011-MTP012 and MTP014-MTP015. See Figs 4 and 5 for 67P's morphological regions.

Spectral indicators maps enable the investigation of surface spatial variability, and through temporal filtering of the data set give the possibility to follow surface evolution both on local and global scales. With the aim to investigate the long-term variability of 67P surface properties, spectral indicators have been mapped for the different phases of the mission, in order to describe their evolution in the period from 2014 August to 2015 May.

In order to provide an as homogeneous as possible spatial coverage of the surface and to improve the statistics, different MTPs have been merged, with the exception of MTP006, which in itself provides a good spatial coverage for the corresponding time frame. Because of the relatively small number of nucleus observations available during MTP013 and the corresponding unfavourable observation geometries (large emission angles), this sequence has not been included in the following analysis. The final arrangement of the observation sequences is reported in Table 1 and the corresponding redundancy maps are shown in Fig. 3.

Table 1.

Observation sequences arrangement for the spectral parameters maps: MTPs, the corresponding time frame and heliocentric distance rh.

MTPs Time frame rh (au)
MTP006 2014 Aug 8–2014 Sept 2 3.62–3.44
MTP007-MP008 2014 Sept 2–2014 Oct 24 3.44–3.12
MTP009-MTP010 2014 Oct 24–2014 Dec 19 3.11–2.75
MTP011-MTP012 2014 Dec 19–2015 Feb 10 2.73–2.35
MTP014-MTP015 2015 Mar 10–2015 May 5 2.13–1.72
MTPs Time frame rh (au)
MTP006 2014 Aug 8–2014 Sept 2 3.62–3.44
MTP007-MP008 2014 Sept 2–2014 Oct 24 3.44–3.12
MTP009-MTP010 2014 Oct 24–2014 Dec 19 3.11–2.75
MTP011-MTP012 2014 Dec 19–2015 Feb 10 2.73–2.35
MTP014-MTP015 2015 Mar 10–2015 May 5 2.13–1.72

### 4.1 SVIS and SIR maps

In Figs 6(a) and (b), the maps of the visible and infrared slopes SVIS and SIR are presented. In both cases, the mapped surface shows a marked dichotomy between the Hapi neck region (see El-Maarry et al. 2016 and Figs 4 and 5 for classification of 67P's morphological regions) characterized by smaller values of the spectral slopes and the rest of the comet. Most of the spatial and temporal variability is contained within 1.35–2.55 μm−1 and 0.23–0.32 μm−1 for SVIS and SIR, respectively. A spatial and temporal correlation at global scale emerges between the two indicators. As an example, in the Imhotep region, a reduction of both slopes can be noticed between MTP006 and MTP014-MTP015 observation phases, with average values going from 2.13 to 1.97 μm− 1 in the VIS and from 0.28 to 0.25 μm− 1 in the IR. As a general trend, both parameters undergo a general reduction with time across the body of the comet, indicating a global flattening of the spectral slope, as 67P was moving along the inbound orbit towards perihelion.

Figure 4.

Morphological regions for the nucleus of 67P adapted from El-Maarry et al. (2016). White areas are the subregions investigated in Section 4.4 along with Hapi: roundish features (R.F.), mid-Imhotep (M.I.), north Ash (N.A.) and smooth Seth (S.S.). The reference system is the same as used for the maps in Fig. 3.

Figure 4.

Morphological regions for the nucleus of 67P adapted from El-Maarry et al. (2016). White areas are the subregions investigated in Section 4.4 along with Hapi: roundish features (R.F.), mid-Imhotep (M.I.), north Ash (N.A.) and smooth Seth (S.S.). The reference system is the same as used for the maps in Fig. 3.

Figure 5.

Two different views of the morphological regions of 67P projected on the nucleus shape model (a and c) and of the subregions investigated in Section 4.3 (white areas in b and d panels) along with Hapi: roundish features (R.F.), mid-Imhotep (M.I.), north Ash (N.A.) and smooth Seth (S.S.).

Figure 5.

Two different views of the morphological regions of 67P projected on the nucleus shape model (a and c) and of the subregions investigated in Section 4.3 (white areas in b and d panels) along with Hapi: roundish features (R.F.), mid-Imhotep (M.I.), north Ash (N.A.) and smooth Seth (S.S.).

Figure 6.

SVIS (left-hand column, a) and SIR (right-hand column, b) maps. From top to bottom: MTP006, MTP007-MTP008, MTP009-MTP010, MTP011-MTP012 and MTP014-MTP015.

Figure 6.

SVIS (left-hand column, a) and SIR (right-hand column, b) maps. From top to bottom: MTP006, MTP007-MTP008, MTP009-MTP010, MTP011-MTP012 and MTP014-MTP015.

Figure 7.

BA (left-hand column, a) and BC (right-hand column, b) maps. From top to bottom: MTP006, MTP007-MTP008, MTP009-MTP010, MTP011-MTP012 and MTP014-MTP015.

Figure 7.

BA (left-hand column, a) and BC (right-hand column, b) maps. From top to bottom: MTP006, MTP007-MTP008, MTP009-MTP010, MTP011-MTP012 and MTP014-MTP015.

Because of the reduction of the subsolar latitude with time (see Fig. 1), passing from an average value of ≈24° in 2014 August to ≈8° in 2015 April, the neck region is barely illuminated during the MTP014-MTP015 period since it enters the winter season. Conversely, a larger portion of the Southern hemisphere encompassing the Bes, Geb, Khonsu and Anhur regions becomes illuminated and observable.

We note that the maps relative to the intermediate sequences MTP007-MTP008, MTP009-MTP010 and MTP011-MTP012 are much more contrasted with respect to the ones relative to the MTP006 and the MTP014-MTP015 sequences. This is primarily due to a larger contribution of shadows in the pixels’ field of view since the former has been acquired at higher solar phase angles (Fig. 1). In addition, the instrument response in proximity of highly contrasted light-shadow transitions can produce spurious effects. These are more important for the IR channel, which is characterized by a wider instrumental point spread function compared to the VIS (Filacchione 2006). Although this effect can introduce small-scale artefacts in correspondence of regions characterized by complex morphology, the large-scale spatial distribution of the spectral indicators remains meaningful.

### 4.2 BA and BC maps

In Figs  7(a) and (b), the maps of the BA and BC are shown. The dichotomy between the neck and the rest of the nucleus observed in the spectral slope maps is still evident in both BC and BA maps, with Hapi showing a BC located at shorter wavelengths and exhibiting a larger BA. As for SVIS and SIR, also BC and BA appear to be well correlated, both spatially and temporally across the whole surface. In particular, the increase of the BA with time in the Imhotep region is followed by a shift of the BC towards shorter wavelengths. Similarly, the reduction of BA in the eastern part of Hapi is followed by a shift of BC towards longer wavelengths. As discussed above, in correspondence of structures with complex morphology as the southernmost part of Ash or across the head (centred on Hatmehit), spurious signatures may arise, for both parameters. However, from the obtained spatial distributions it can be noted that both BA and BC are less affected by this issue than SIR. Along with this, residuals from the photometric correction and thermal removal introduce a further source of systematic variability.

Comparing the maps of the spectral slopes and of the 3.2 μm feature indicators, we find an overall correlation of the four parameters with slopes decreasing for larger BAs and smaller wavelength BCs as was also shown by Filacchione et al. (2016a) for the mission phase before the Philae landing. Referring to Fig. 2(b), this concomitant behaviour of the CSI is in agreement with a variation of the water ice abundance.

### 4.3 Water ice abundance maps

From the reported results of 67P surface spectral modelling (De Sanctis et al. 2015; Barucci et al. 2016; Filacchione et al. 2016b; Raponi et al. 2016), the lack of a widespread water ice 2 μm feature across the nucleus, with the exception of exposed water-ice-enriched area, indicates that ice, if present, is intimately mixed with the nucleus average material (dark terrain). Typical water ice abundances are of the order of few per cent, with maximum values of ≈10 per cent reached in the Hapi region, in proximity of receding shadows (De Sanctis et al. 2015). Assuming this compositional paradigm, it is then possible to convert spectral indicators maps into water ice abundance maps. For this purpose, we applied Hapke's spectral modelling (Ciarniello et al. 2011; De Sanctis et al. 2015; Raponi et al. 2016) to compute simulated spectral indicators as a function of the abundance of water ice intimately mixed with the dark terrain material.

This is achieved by deriving synthetic spectra as a linear combination of the SSA of the end-members (see Appendix A). The SSA spectrum of the ice fraction was computed as a function of the grain size and its optical constants, while for the dark terrain, the spectrum derived in Ciarniello et al. (2015) is used. The spectral SSA of the mixture can be directly compared to VIRTIS spectra after the application of the photometric correction derived by Ciarniello et al. (2015), that converts the measured reflectance into SSA (see Appendix A). For a given set of spectral indicators corresponding to a certain position in the spectral indicators maps, is then possible to determine the proper combination of water ice abundance and grain size in the simulated spectrum that provides the best match for that set. This can be done by minimizing a chi-square variable χ2, which is built as $$\chi ^2=(\frac{S_{{\rm VIS}}^m-S_{{\rm VIS}}}{\sigma _{S_{{\rm VIS}}}})^2+(\frac{S_{{\rm IR}}^m-S_{{\rm IR}}}{\sigma _{S_{{\rm IR}}}})^2+(\frac{{\rm BA}^m-{\rm BA}}{\sigma _{{\rm BA}}})^2$$, where the superscript m indicates modelled quantities and $$\sigma _{S_{{\rm VIS}}}$$, $$\sigma _{S_{{\rm IR}}}$$ and σBA are the standard deviations of the mean value of each parameter. The BC is not included in this analysis because its absolute position can be potentially affected by residual spectral artefacts and a systematic shift introduced by the thermal removal procedure. Conversely, these issues do not affect SVIS and SIR and have a negligible effect on BA.

Although theoretically possible, the simultaneous determination of water ice abundance and grain size with this method has proven impractical, because of the small number (three) of fitted values. This is shown in Fig. 8, where simulated CSI are reported for different grain sizes and water ice abundances. It can be noted that SVIS and SIR are weakly sensitive to grain size in the investigated range whereas they are affected by minimal amounts of water ice while BA shows a certain dependence with d only at very small sizes <2 μm. For this reason, we assume a fixed water ice grain size d = 2 μm in our analysis, compatible with the result derived by (De Sanctis et al. 2015). It is worth mentioning that since the spectral modelling is performed using the average 67P dark terrain as an end-member of the mixture, the derived water ice abundances must be considered as an additional contribution to the average ice content of the dark terrain, which at this stage is unknown, although extremely low, given the lack of evident H2O features.

Figure 8.

Plots of lookup tables depicting the simulated CSI as functions of grain size d (1–5 μm) and water ice abundance (0–4  per cent): SVIS (a), SIR (b) and BA (c). CSI are computed from synthetic spectra of dark terrain–water ice intimate mixtures as derived from Hapke's theory.

Figure 8.

Plots of lookup tables depicting the simulated CSI as functions of grain size d (1–5 μm) and water ice abundance (0–4  per cent): SVIS (a), SIR (b) and BA (c). CSI are computed from synthetic spectra of dark terrain–water ice intimate mixtures as derived from Hapke's theory.

Water ice abundance maps as derived with the method described above are shown in Fig. 9(a). As explained in Appendix A, the water ice abundance we report must be considered as the ice cross-section fraction seen by the photons within the medium. The typical error on the derived abundance is computed by the χ2 distribution and it is below 0.1 per cent. Here, we give a global description of the maps for each investigated time period.

• MTP006 (2014 Aug 8–2014 Sept 2): the observed dichotomy of CSI values between Hapi region and the rest of the nucleus in the early mission phases is reproduced in the water ice content. In particular, during the MTP006 sequence, an average water ice content of ≈0.7 per cent is measured on Hapi, reducing to 0.4 per cent in the northern part of Seth (smooth Seth, Fig. 4) while most of the rest of the nucleus is consistent with the absence of ice on the surface. Isolated ice-enriched areas can be recognized in the Imhotep region, at lon. 117°–lat. 13°, lon. 181°–lat. −5° and lon. 189°–lat. −5°, described as exposed ice patches in Filacchione et al. (2016a) and Barucci et al. (2016), respectively. In these cases, the spectral modelling paradigm adopted in this work cannot be rigorously applied, since Filacchione et al. (2016b) have shown that in these areas, ice is present both as an intimate and areal mixture with the refractory material. This implies that, in these cases, the derived water ice abundance cannot be directly compared to Filacchione et al. (2016b) results and is possibly biased by the different adopted modelling.

An excess of water ice with respect to the surrounding terrains is measured in an extended portion of the surface located approximately at lon. 150°–180°, lat. −30° to 0°, which corresponds to the roundish features described by Auger et al. (2015) with abundances of the order of ≈0.1 per cent. Auger et al. (2015) correlate the morphology of this structure to remnants of volatiles degassing conducts, which could be compatible with the measured water ice enrichment. Outside of Hapi and Imhotep, areas with increased water ice content are distributed at the boundaries between Seth, Atum, Anubis and Ash, characterized by a complex morphology with average values of ≈0.4 per cent and between Ash and Aten (≈0.6 per cent) at lon. 115_{.}^{\circ}$$5–lat. 28°. • MTP007-MTP008 (2014 Sept 2–2014 Oct 24): the global distribution of water ice remains relatively unchanged with respect to the previous MTP006 data. None the less, the average content of ice is slightly increased, in particular in the neck region. A similar increase is also measured in the Imhotep region and on the northern part of Ash. • MTP009-MTP010 (2014 Oct 24–2014 Dec 19): the trend of increasing abundance of water ice on the surface is confirmed also in this time frame. The spatial coverage of this observational phase allows one to perform a direct comparison with previous observation phases for Hapi, the northern part of Seth and Imhotep, which all show an increase of water ice abundance. Within the Imhotep region, an enrichment of ice is also measured in correspondence of the roundish features. Also the northern part of Ash is relatively well mapped; however, in this region, we observe a reduction of the water ice abundance. • MTP011-MTP012 (2014 Dec 19–2015 Feb 10): the spatial coverage of this observational phase is the worst of the whole data set. None the less, it is possible to notice a reduction of the water ice abundance in the Hapi region, while the scarce data available in Imhotep suggest a further increase of the water ice content, which is shown also in the portion of the roundish features accessed by the observations. • MTP014-MTP015 (2015 Mar 10–2015 May 5): during this observation phase the Southern hemisphere of the nucleus became observable, while the northern one enters the winter season. Distribution of water ice in Imhotep is broader with respect to the previous phases, leading to global larger abundances. The region of the roundish features is well covered by this set of observations and represents the area with the largest amount of water ice in Imhotep. An increase of ice is also measured in the head of the comet, mainly in Wosret and Ma'at. The observed Southern hemisphere, represented by the Bes and Anhur regions, shows large abundances of ice up to 2 per cent. However, these values must be considered with caution because the observing geometry (large emission and/or incidence angle) in this area was unfavourable. Figure 9. Water ice abundance maps (a) and solar flux maps (b). Figure 9. Water ice abundance maps (a) and solar flux maps (b). ### 4.4 Water ice abundance and solar flux In Fig. 9(b), we show maps of the median solar flux (SF) at 67P's surface. This quantity can be expressed as$$\frac{\cos (i)}{r_{\rm h}^2}F_{\bigoplus }$$, where i is the incidence angle, rh is the comet's heliocentric distance expressed in au and$$F_{\bigoplus }$$is the solar constant. SF maps are computed according to the VIRTIS pixel's geometry and time of acquisition and are thus related to the average instantaneous solar flux at the observation conditions of the sampled region. It can be noted that the overall SF increases from MTP006 to MTP015, due to the reduction of the heliocentric distance. This increase is modulated by the seasonal effect due to the shift of the subsolar point towards smaller latitudes, which produces a reduction of the SF in Hapi and an increase in Imhotep. In order to explore a possible correlation between insolation at the surface and derived water ice abundance, five areas with good spatial and temporal coverage and relatively uniform illumination conditions during the different observation sequences have been selected, and their average water ice content has been compared to the corresponding average SF (Fig. 10a). The selected regions (Figs 4, 5b,d)) correspond to the centre of the Imhotep region (mid-Imhotep), the northern portion of Ash (north Ash), the northern and smoother part of Seth (smooth Seth), a portion of the roundish features and the whole Hapi region and are discussed below. • Hapi: In this area, the abundance of water ice is the largest among the selected regions. The SF value increases from 2014 August (90 W m−2) to November (125 W m−2), with a reduction in 2015 January–February (120 W m−2) before, entering the winter season due to the shift of the subsolar point to smaller latitudes that compensates, in terms of solar flux, the increase provided by the reduction of the heliocentric distance. The water ice abundance increases from ≈0.75 per cent in 2014 August up to ≈0.95 per cent in 2014 November, and decreases to 0.85 per cent in 2015 January, following the trend of SF. Given the very low amount of observations during the MTP014-MTP015 phases, this data set has not been included in the analysis. Figure 10. (a) Water ice abundance (diamonds and solid lines) and SF (triangles and dashed lines) as functions of time for the five selected regions. (b) ISFrot (solid lines) and daylight hours (dashed lines) as derived from spice simulations as functions of time for the five selected regions. Figure 10. (a) Water ice abundance (diamonds and solid lines) and SF (triangles and dashed lines) as functions of time for the five selected regions. (b) ISFrot (solid lines) and daylight hours (dashed lines) as derived from spice simulations as functions of time for the five selected regions. • Smooth Seth: Water ice abundance varies in the ≈0.4–0.5 per cent range from 2014 August to 2015 January–February, with a significant increase during 2015 April–May at ≈0.6 per cent, corresponding to the maximum value of the SF (225 W m−2). • Mid-Imhotep: This area has the lowest abundance of water ice among the five investigated across the whole analysed time interval, being null in 2014 August and growing to ≈0.15 per cent in 2015 February–January until 2015 April–May, while showing a positive trend in SF (from 45 to 270 W m−2). Along with the local enrichment in water ice in the selected zone, the water ice distribution in Imhotep broadens with time, thus resulting in an increase of the total ice content. • Roundish features: This area is not evenly covered during all the observation sequences, none the less a significant increase of the water ice abundance can be measured from 2014 August to 2015 January–February, with average contents going from ≈0.1 to 0.9 per cent, showing the largest increase rate with respect to the other selected regions. In the 2015 April–May observations, a reduction of the average water ice content is measured down to 0.5 per cent. Conversely, SF increases monotonically within the whole investigated time frame across this area. We want to mention that because of the limited spatial coverage and redundancy of the observations taken through 2014 November–December and 2015 January–February, the derived values of the water ice abundance must be considered with care and could possibly be overestimated. Conversely, the good observation conditions in 2015 April–May provide a more robust determination of the real water ice content and indicate also an enlargement of the ice-enriched area with respect to the previous observations during MTP006 and MTP007-MTP008, similar to what is observed in the rest of the Imhotep region. • North Ash: The selected area within Ash corresponds to the portion of this region that, because of its orientation, receives the lowest amount of instantaneous solar flux across all the analysed observation phases (Fig. 9b). After the mid-Imhotep, area it is the least enriched in water ice. Its abundance increases from ≈0.1 per cent in 2014 August to 0.3 per cent in 2015 November–December, and then decreases to ≈0.2 per cent in 2015 January–February and to ≈0.15 per cent 2015 April–May. SF at the time of the observations doubles from MTP006 (80 W m−2) to MTP014-MTP015 (165 W m−2) with an increment that is the lowest among the different selected regions, except Hapi. To give a more exhaustive description of the insolation conditions of these five areas, we report in Fig. 10(b), the integrated solar flux per comet rotation ISFrot and the total daylight hours as derived from spice simulations for one selected position within each of the investigated regions, assuming SHAP5 nucleus’ digital terrain model (Jorda et al. 2016). It can be noted that both roundish features area and mid-Imhotep have a very similar thermal history, because of their close position in the central part of the Imhotep region, which is characterized by an almost constant amount of daylight hours (≈5.5 h) across the whole 2014 August–2015 May period and an increasing ISFrot. The north Ash region is permanently illuminated from 2014 August to 2015 March, thus receiving a larger amount of solar flux with the comet approaching the Sun. Conversely, daylight hours in Hapi progressively decrease up to 2015 May, while ISFrot is nearly constant because of the simultaneous reduction of heliocentric distance. A similar trend for the daylight hours is observed in smooth Seth, where the ISFrot increases with time up to the beginning of 2015 April and then rapidly decreases towards the equinox. ### 4.5 Discussion A first comparison of the temporal water ice evolution (Fig. 9a) between mid-Imhotep and the roundish features terrains shows a completely different behaviour despite the almost identical illumination history. The most significant difference is seen in the absolute value of H2O being distinctively larger in the roundish features since 2014 August, indicating a surface and subsurface lateral heterogeneity that is maximally expressed in the exposition of water ice rich patches (Barucci et al. 2016; Filacchione et al. 2016b; Raponi et al. 2016). None the less, both mid-Imhotep and roundish terrains areas show an increase of the water ice abundance with time, which has to be related to the increase of insolation and, simultaneously, of the sublimation rate. The observed trend can be understood when assuming that an enhanced activity causes removal of dust from the top surface, exposing layers that are richer in ice. However, this is a complex process since only a part of the solar energy penetrates the dust layer and triggers sublimation. The available amount of energy for sublimation depends on the heliocentric distance and is related to the complex thermal balance and the varying dust mantle thickness and is not proportional to the incoming radiation (Marboeuf & Schmitt 2014). Moreover, the mechanism of lifting off a dust layer by subliming gas molecules depends on many parameters, such as the surface gas velocity and dust mass loading that are not fully understood yet while the amount of ice below the dust may be limited. Between the MTP011-MTP012 and MTP014-MTP015 phases, both SF and ISFrot increase by a factor of 3 while the water ice abundance in Imhotep remains constant, possibly indicating that the sublimation front reached a depth with a vertically constant dust–ice mixing ratio. Conversely, the roundish terrain area shows a reduction of the water ice content during MTP014-MTP015. This can be interpreted as being due to the depletion of a former localized water ice reservoir, similar to what was found by Raponi et al. (2016) for isolated exposed water ice patches. However, given the uneven spatial coverage of observations across the roundish terrain region this argument must be considered with caution. As discussed above, the north Ash area was permanently illuminated from 2014 August up to the beginning of 2015 March. This indicates that during a single rotation, the insolation at the surface does not change strongly, thus the abundance of water ice should not exhibit large variability related to the instantaneous illumination conditions. Conversely, a correlation with the total solar flux received during one rotation is more likely. This hypothesis is supported by the fact that the north Ash area shows an increase of water ice from 2014 August up to the second half of 2014 November and then a reduction, qualitatively similar to the ISFrot trend, while SF, which is related to the average instantaneous illumination conditions, grows monotonically with time. Moving to smooth Seth, the amount of ice on the surface is fairly constant with time, showing an increase towards 2015 April that corresponds to the maximum of ISFrot and of SF. This makes it difficult to disentangle the effects of instantaneous illumination condition and integrated flux. The last investigated region is Hapi, corresponding to the neck of the comet. The water ice abundance trend follows the behaviour of both SF and ISFrot, with a maximum during MTP009-MTP010 followed by a reduction during MTP011-MTP012. Keller et al. (2015) estimated that self-heating due to the concave shape of the neck could contribute up to 50 per cent of the thermal flux the surface, thus the calculated SF and ISFrot must be considered as lower limits of the corresponding energetic inputs to the surface. Hapi region is on average the richest in ice among the selected areas. This does not seem to be related to the local insolation conditions, since, as reported by Sierks et al. (2015), ISFrot is 15 per cent lower than the most illuminated regions and the total integrated flux per orbit is the lowest across 67P's surface. Conversely, it points towards an intrinsic compositional difference in Hapi. A possible explanation for this can be given assuming it is a re-deposition region of dust mixed with icy grains, similar to what has been suggested by A'Hearn et al. (2011) for the ‘waist’ region of comet Hartley 2. In the case of 67P, the material could have been emitted from the poorly illuminated comet's Southern hemisphere, possibly driven by CO2 activity during southern winter (Bockelée-Morvan et al. 2016; Fink et al. 2016; Migliorini et al. 2016). Such a hypothesis is compatible with the fact that Hapi corresponds to a minimum of the comet's gravitational potential and exhibits a very smooth surface texture (Sierks et al. 2015). ## 5 DIURNAL VARIABILITY The variability of the surface water ice content on the comet's rotation period time-scale is investigated, similar to Section 4, by means of spectral parameters mapping. In this case, observations from MTP006 are used, providing the best spatial and local time coverage of the surface. Maps are produced for 12 different rotational phases, which approximately correspond to 1 h steps. This procedure allows us to monitor the variability of surface spectral properties versus the instantaneous illumination conditions and the light-shadow cycle at each position. The selection of data from a single observation sequence (MTP006, i.e. one month) minimizes the seasonal effects of surface variability although it also limits the signal to noise ratio of the measurements. For this reason, only the BA parameter is mapped, being the one most sensitive to the water ice abundance. Three different regions have been investigated separately: the neck, the head and the body of the comet. • Neck: BA maps are shown in Fig. 11 with an orthographic projection centred at the comet's North pole. Given the bi-lobed shape of 67P, a complex pattern of shadows is cast on to the neck region. Three fronts with a larger value of BA can be noticed as following receding shadows. This behaviour confirms results by De Sanctis et al. (2015) that measured an enrichment of water ice in Hapi for selected areas coming out from shadows. It is important to mention that fronts with increased BA are not detected at the border of advancing shadows. This is explained by assuming that the observed ice is produced by re-condensation of water vapour sublimating from the subsurface, when the corresponding area stops being illuminated. Such a mechanism is effective in the neck where the surface undergoes sudden shadowing, produced by the two lobes of 67P, which determines a positive temperature gradient with depth (De Sanctis et al. 2015). When the surface is exposed to sunlight again, ice from the top layers sublimates on a time-scale smaller than 1 h, which represents the temporal resolution of this study. An analogous mechanism is also responsible for the formation of sunset jets, as shown by (Shi et al. 2016), for the shadowed areas in Ma'at region when the comet was at 1.8 au (late 2015 April). In this case, the observed activity is directly sustained by the vertical temperature gradient due to the subsurface thermal lag. Figure 11. BA orthographic maps centred at the North pole of 67P. The rotational phase increases by ≈ 1 h/step from top to bottom and left to right. Grey areas represent missing data and regions not mapped, while black represent shadows as derived from spice calculations (Acton 1996). Shadows move counterclockwise because of the comet's rotation. The red, green and yellow arrows indicate three fronts showing an increased BA at places exposed by receding shadows, which is interpreted as an increased water ice abundance. Figure 11. BA orthographic maps centred at the North pole of 67P. The rotational phase increases by ≈ 1 h/step from top to bottom and left to right. Grey areas represent missing data and regions not mapped, while black represent shadows as derived from spice calculations (Acton 1996). Shadows move counterclockwise because of the comet's rotation. The red, green and yellow arrows indicate three fronts showing an increased BA at places exposed by receding shadows, which is interpreted as an increased water ice abundance. • Body: BA maps are shown in Fig. 12 with a cylindrical projection centred in Imhotep at 0° lat–180° lon. Differently from the neck area, this region undergoes one day-night cycle per rotation, without the occurrence of sudden shadowing. The overall value of BA is lower when compared to the neck and does not show significant variations with time related to the shadow pattern. One exception is represented by an extended increase of BA on Aten(≈15° lat to ≈90° lon), when exiting from shadow, followed by a sudden drop of BA. Notably, the roundish features area, which exhibits larger values of BA with respect to the Imhotep average (interpreted as a larger content of water ice), does not show any relevant temporal variability. Figure 12. BA cylindrical maps centred in Imhotep at 0° lat–180° lon. The rotational phase increases by ≈1 h/step from top to bottom and left to right. Grey areas represent missing data and regions not mapped, while black represents shadows as derived from spice calculations (Acton 1996). The shadow pattern proceeds from left to right. The green arrow indicates the Aten region showing an increase of BA. Figure 12. BA cylindrical maps centred in Imhotep at 0° lat–180° lon. The rotational phase increases by ≈1 h/step from top to bottom and left to right. Grey areas represent missing data and regions not mapped, while black represents shadows as derived from spice calculations (Acton 1996). The shadow pattern proceeds from left to right. The green arrow indicates the Aten region showing an increase of BA. • Head: BA maps are shown in Fig. 13 with a cylindrical projection centred in Hatmehit at 15° lat–345° lon. Similar to the body area, this region undergoes one day-night cycle per rotation and does not experience sudden shadowing. The overall value of BA is similar to the body and does not show significant variations with time related to the shadow pattern. None the less minor variations of BA can be observed in the northern portion of the head, possibly related to effects produced by extreme observational geometries and unresolved shadows. Figure 13. BA cylindrical maps centred in Hatmehit region at 15° lat–345° lon. The rotational phase increases by ≈1h/step from top to bottom and left to right. Grey areas represent missing data and regions not mapped, while black represents shadows as derived from spice calculations (Acton 1996). The shadow front moves from positive to negative latitudes and backward. Figure 13. BA cylindrical maps centred in Hatmehit region at 15° lat–345° lon. The rotational phase increases by ≈1h/step from top to bottom and left to right. Grey areas represent missing data and regions not mapped, while black represents shadows as derived from spice calculations (Acton 1996). The shadow front moves from positive to negative latitudes and backward. ### 5.1 Discussion The diurnal variability of water ice across the 67P's surface reflects the compositional dichotomy between the neck and the rest of the nucleus. Only the Hapi region shows evidence of a water ice cycle during MTP006. Such a process was identified for the first time by De Sanctis et al. (2015) in a limited area of the neck, where an increased amount of water ice was detected along shadows, with maximum values comprised between 5 and 14 per cent. This was connected to the effect of sudden shadowing, able to produce a vertical temperature gradient that sustains sublimation from subsurface and re-condensation on the top layers. Here, we show that the same mechanism is active across the whole Hapi region, with the formations of multiple water-ice-enriched fronts which follows receding shadows. The ubiquity of this process across the neck region reinforces the role of the light-shadow cycle as the main driver of the water ice abundance variability on the comet's surface. This mechanism is enhanced by the favourable illumination condition during MTP006 that produced in Hapi, a relatively large instantaneous solar flux with respect to the Southern hemisphere. Assuming the mechanism described above, a water ice cycle is prevented in areas that are permanently illuminated during MTP006 such as Ash and Ma'at. Lack of such an effect in Imhotep and in the part of the head that undergoes a day-night cycle can be instead related to the absence of sudden shadowing. In these cases, the solar incidence angle progressively increases towards the terminator and the consequent gradual reduction of the energy input at the surface prevents the formation of a strong vertical thermal gradient. The same mechanism acting in Hapi can explain the increase of BA observed in Aten when coming out of shadow. This region, because of its morphology and position is possibly subject to sudden shadowing. The smaller magnitude of the effect in this case can be related to lower instantaneous solar flux received during the northern summer with respect to Hapi and the intrinsically lower content of water ice. ## 6 SUMMARY AND CONCLUSIONS We have investigated the variability of the 67P nucleus surface properties at seasonal and diurnal temporal scales by means of VIRTIS-M data acquired from 2014 August (MTP006) up to the beginning of 2015 May (MTP015). To this aim, a set of four CSI (Filacchione et al. 2016a) has been adopted to characterize the shape of the measured spectrum and reduce the spectral dimension of the data set: SVIS (VIS spectral slope), SIR (IR spectral slope), BA (band area of the 3.2 μm absorption feature) and BC (band centre of the 3.2 μm absorption feature). The seasonal variability of the surface properties has been studied by producing global maps of the CSI for the different mission phases, corresponding to the inbound leg of the comet's orbit (3.62–1.72 au) (see Figs 6a, b, 7a and b). For each selected time frame, a general spatial correlation of the CSI has been recognized, with smaller SVIS and SIR corresponding to larger values of BA and smaller BC. This is interpreted as an effect of the variability of the amount of water ice intimately mixed with the refractory material composing the surface of 67P. Hapke spectral modelling of the CSI has been performed to derive the amount of water ice at the surface, indicating Hapi as the region richest in H2O ice (0.75–0.95 per cent, see Fig. 9a). Another area with a relatively large amount of water ice is represented by the roundish features region (Auger et al. 2015), with abundances ranging in the 0.1–0.9 per cent interval. As a general trend water ice abundances increase across the body and the head of the nucleus from MTP006 to MTP015, with the comet approaching the Sun (from 0 to 0.15 per cent in mid-Imhotep). None the less, an analysis of selected areas reveals the complexity of the processes acting at smaller spatial scales, where the seasonal effects due to the shift of the subsolar point towards smaller latitudes (change of the instantaneous solar flux to the surface and variation of the amount of daylight hours) and the reduction of the heliocentric distance both drive the water ice sublimation rate and are coupled to the intrinsic non-homogeneity of the water ice distribution. In this context, the roundish features area shows an increment of the water ice content (up to 2015 January) followed by a reduction (2015 April), despite an increase of energy input, that is possibly related to the depletion of a subsurface water ice reservoir. A similar trend is observed in Hapi and north Ash, which in contrast seems to be related to the seasonal reduction of the integrated solar flux per comet rotation. In contrast, smooth Seth and mid-Imhotep show an increase of the water ice content, in agreement with the increase of both instantaneous and diurnal solar flux. In order to investigate the diurnal variability of the water ice content during one comet rotation, BA maps at 12 different rotational phases (≈1 h/step) have been produced for three separate regions: neck, body, head. It emerges that only the neck area exhibits an evident diurnal cycle of water ice at global scale, with H2O-enriched fronts following receding shadows and sublimating in less than 1 h when illuminated. The occurrence of such a diurnal cycle is related, as indicated in De Sanctis et al. (2015), to the effect of the sudden shadowing alternatively produced by each lobe of the nucleus in the neck region, allowing the formation of a strong positive temperature gradient from the top colder surface to the deeper warmer layers, where water ice sublimates to re-condense at the surface. ## Acknowledgments We thank the following institutions and agencies for support of this work: Italian Space Agency (ASI, Italy, contract number I/024/12/1), Centre National d’Études Spatiales (CNES, France), DLR (Germany), NASA (USA) Rosetta Program and Science and Technology Facilities Council (UK). VIRTIS was built by a consortium that includes Italy, France and Germany, under the scientific responsibility of the Istituto di Astrofisica e Planetologia Spaziali of INAF, Italy, which also guides the scientific operations. The VIRTIS instrument development led by the prime contractor Leonardo-Finmeccanica (Florence, Italy) has been funded and managed by ASI, with contributions from Observatoire de Meudon financed by CNES, and from DLR. We thank the Rosetta Science Ground Segment and the Rosetta Mission Operations Centre for their support throughout all the phases of the mission. The VIRTIS calibrated data will be available through the ESA’s Planetary Science Archive website (http://www.rssd.esa.int/index.php?project=PSA&page=index) and is available upon request until posted to the archive. ## REFERENCES A'Hearn M. F. et al. , 2011 , Science , 332 , 1396 Acton C. , 1996 , Planet. Space Sci. , 44 , 65 Auger A.-T. et al. , 2015 , A&A , 583 , A35 Barucci M. A. et al. , 2016 , A&A , 595 , A102 Bockelée-Morvan D. et al. , 2015 , A&A , 583 , A6 Bockelée-Morvan D. et al. , 2016 , MNRAS , 462 , 170 Capaccioni F. et al. , 2015 , Science , 347 , 628 Ciarniello M. et al. , 2011 , Icarus , 214 , 541 Ciarniello M. et al. , 2015 , A&A , 583 , A31 Ciarniello M. et al. , 2016 , A&A , 462 , S476 Clark R. N. et al. , 2012 , Icarus , 218 , 831 Coradini A. et al. , 2007 , Space Sci. Rev. , 128 , 529 De Sanctis M. C. , Capria M. T. , Coradini A. , Ammannito E. , 2007 , AJ , 133 , 1836 De Sanctis M. C. , Lasue J. , Capria M. T. , 2010 , AJ , 140 , 1 De Sanctis M. C. et al. , 2015 , Nature , 525 , 500 El-Maarry M. R. et al. , 2016 , A&A , 593 , A110 Filacchione G. , 2006 , PhD thesis , INAF-IASF, Rome, Italy Filacchione G. et al. , 2016a , Icarus , 274 , 334 Filacchione G. et al. , 2016b , Nature , 529 , 368 Fink U. et al. , 2016 , Icarus , 277 , 78 Fornasier S. et al. , 2015 , A&A , 583 , A30 Groussin O. et al. , 2015 , A&A , 583 , A36 Hapke B. , 2012 , Theory of Reflectance and Emittance Spectroscopy Cambridge Univ. Press , Cambridge Hässig M. et al. , 2015 , Science , 347 , 276 Henyey L. G. , Greenstein J. L. , 1941 , ApJ , 93 , 70 Jorda L. et al. , 2016 , Icarus , 277 , 257 Keller H. U. et al. , 2015 , A&A , 583 , A34 Li J.-Y. et al. , 2007a , Icarus , 187 , 41 Li J.-Y. , A'Hearn M. F. , McFadden L. A. , Belton M. J. S. , 2007b , Icarus , 188 , 195 Li J.-Y. , A'Hearn M. F. , Farnham T. L. , McFadden L. A. , 2009 , Icarus , 204 , 209 Li J.-Y. et al. , 2013 , Icarus , 222 , 559 Marboeuf U. , Schmitt B. , 2014 , Icarus , 242 , 225 Marboeuf U. , Schmitt B. , Petit J.-M. , Mousis O. , Fray N. , 2012 , A&A , 542 , A82 Mastrapa R. , Bernstein M. , Sandford S. , Roush T. , Cruikshank D. , Ore C. D. , 2008 , Icarus , 197 , 307 Mastrapa R. , Sandford S. A. , Roush T. L. , Cruikshank D. P. , Ore C. M. D. , 2009 , ApJ , 701 , 1347 Migliorini A. et al. , 2016 , A&A , 589 , A45 Orosei R. , Capaccioni F. , Capria M. T. , Coradini A. , Espinasse S. , Federico C. , Salomone M. , Schwehm G. H. , 1995 , A&A , 301 , 613 Pommerol A. et al. , 2015 , A&A , 583 , A25 Preusker F. et al. , 2015 , A&A , 583 , A33 Prialnik D. , Sarid G. , Rosenberg E. D. , Merk R. , 2008 , Space Sci. Rev. , 138 , 147 Quirico E. et al. , 2016 , Icarus , 272 , 32 Raponi A. et al. , 2016 , MNRAS , 462 , S476 Schröder S. E. , Grynko Y. , Pommerol A. , Keller H. U. , Thomas N. , Roush T. L. , 2014 , Icarus , 239 , 201 Shi X. et al. , 2016 , A&A , 586 , A7 Sierks H. et al. , 2015 , Science , 347 , 1044 Sunshine J. M. et al. , 2006 , Science , 311 , 1453 Thomas P. C. et al. , 2007 , Icarus , 187 , 4 Veverka J. et al. , 2013 , Icarus , 222 , 424 Warren S. G. , 1984 , Appl. Opt. , 23 , 1206 ### APPENDIX A: PHOTOMETRIC CORRECTION AND SSA SPECTRAL MODELLING As described by Ciarniello et al. (2015), the radiance factor I/F of 67P/Churyumov–Gerasimenko comet can be expressed as a function of the observation geometry using a simplified Hapke model (Hapke 2012) with equation: $$\frac{I}{F}=\frac{w}{4}\frac{\mu _{0{\rm eff}}}{\mu _{0{\rm eff}}+\mu _{{\rm eff}}}p(\alpha )S(i,e,\alpha ,\bar{\theta }),$$ (A1) where • i, e and α are the incidence, emission and phase angle, respectively; • w is the single scattering albedo (SSA); • p(α) is the single scattering phase function (SPPF); • μ0eff and μeff are the effective cosines of the incidence and emission angle, respectively; •$$S(i,e,\alpha ,\bar{\theta })$$is the shadowing function; •$$\bar{\theta }$$is the average surface slope. The SPPF is modelled according to Henyey & Greenstein (1941) expression and depends on the asymmetry parameter b. In equation (A1), the multiple scattering terms is neglected, given the low albedo of the surface, while the opposition effect functions are posed to 1, because the observations were performed at relatively large phase angles (α > 27 _{.}^{\circ}$$2). The spectral behaviour of w and b, and the slope parameter $$\bar{\theta }$$ were derived by Ciarniello et al. (2015) from global analysis of VIRTIS observations, and are assumed as representative of the average properties of the surface of the comet. In particular, the two photometric parameters b and $$\bar{\theta }$$ can be used to invert equation (A1) and derive the SSA from the I/F measured by VIRTIS for each pixel, allowing us to determine the intrinsic variability of the surface and minimizing observation geometry effects (photometric correction).

The spectrum of the SSA of the nucleus was modelled by applying Hapke's theory (Hapke 2012) and assuming an intimate mixture of the average composition of the comet (dark terrain), as derived by Ciarniello et al. (2015) and water ice. The SSA of water ice (wI) in the VIRTIS wavelength range was simulated as a function of the grain size using optical constants from Warren (1984), Mastrapa et al. (20082009) and Clark et al. (2012). In the case of an intimate mixture, particles of different end-members are in contact with each other and the modelled SSA of the surface wm is given by
$$w_m=pw_{I}+(1-p)w,$$
(A2)
where p and 1 − p represent the cross-section fraction, seen by the photons propagating within the medium, of water ice and dark terrain, respectively. It can be shown that
\begin{eqnarray} {p = \frac{N_I\sigma _I}{N_I\sigma _I+N\sigma }} \nonumber \\ 1-p = \frac{N\sigma }{N_I\sigma _I+N\sigma }, \end{eqnarray}
(A3)
where NI, N are the number of particles per unit volume and σI, σ the particle cross-section for water ice and dark terrain, respectively. The grain size of the dark material is unknown, thus preventing us to determine the volumetric abundance of water ice with respect to the dark terrain. However, assuming the same grain size for the two end-members, the p and 1 − p provide a direct estimate of volumetric abundances.