# The changing temperature of the nucleus of comet 67P induced by morphological and seasonal effects

F. Tosi, F. Capaccioni, M. T Capria, S. Mottola, A. Zinzi, A. Longobardo, M. Ciarniello, G. Filacchione, M. Hofstadter, S. Fonti, et al.

► To cite this version:

F. Tosi, F. Capaccioni, M. T Capria, S. Mottola, A. Zinzi, et al. The changing temperature of the nucleus of comet 67P induced by morphological and seasonal effects. Nature Astronomy, 2019, 3 (7), pp.649-658. 10.1038/s41550-019-0740-0. hal-02415158

HAL Id: hal-02415158

https://hal.science/hal-02415158v1

Submitted on 2 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The changing temperature of the nucleus of comet 67P induced by morphological and seasonal effects

F. Tosi<fnref n="1" />*, F. Capaccioni<fnref n="1" />, M. T. Capria<fnref n="1" />, S. Mottola², A. Zinzi³, M. Ciarniello<fnref n="1" />, G. Filacchione<fnref n="1" />, M. Hofstadter⁴, S. Fonti⁵,¹², M. Formisano<fnref n="1" />, D. Kappel²,⁶, E. Kührt², C. Leyrat⁷, J.-B. Vincent², G. Arnold², M. C. De Sanctis<fnref n="1" />, A. Longobardo<fnref n="1" />,⁸, E. Palomba<fnref n="1" />, A. Raponi<fnref n="1" />, B. Rousseau<fnref n="1" />, B. Schmitt⁹, M. A. Barucci⁷, G. Bellucci<fnref n="1" />, J. Benkhoff¹⁰, D. Bockelée-Morvan⁷, P. Cerroni<fnref n="1" />, J.-Ph. Combe¹¹, D. Despan⁷, S. Erard⁷, F. Mancarella⁵, T. B. McCord¹¹, A. Migliorini<fnref n="1" />, V. Orofino⁵ and G. Piccioni<fnref n="1" />

Knowledge of the surface temperature distribution on a comet's nucleus and its temporal evolution at different timescales is key to constraining its thermophysical properties and understanding the physical processes that take place at and below the surface. Here we report on time-resolved maps of comet 67P/Churyumov-Gerasimenko retrieved on the basis of infra-red data acquired by the Visible InfraRed and Thermal Imaging Spectrometer (VIRTIS) onboard the Rosetta orbiter in 2014, over a roughly two-month period in the pre-perihelion phase at heliocentric distances between 3.62 and 3.31 au from the Sun. We find that at a spatial resolution ≤15 m per pixel, the measured temperatures point out the major effect that self-heating, due to the complex shape of the nucleus, has on the diurnal temperature variation. The bilobate nucleus of comet 67P also induces daytime shadowing effects, which result in large thermal gradients. Over longer periods, VIRTIS-derived temperature values reveal seasonal changes driven by decreasing heliocentric distance combined with an increasing abundance of ice within the uppermost centimetre-thick layer, which implies the possibility of having a largely pristine nucleus interior already in the shallow subsurface.

Knowledge of the surface temperature distribution of a comet's nucleus proves to be of fundamental importance for a number of reasons: the evaluation of the thermophysical properties (thermal inertia and roughness, at several spatial scales), the characterization of the thermal behaviour of peculiar surface units (for example, local ice exposures, or inherently cold material units not to be confused with shaded areas), and for the understanding of the physical processes affecting the surface and the shallow subsurface layers (for example, sublimation of volatile compounds). The distribution of surface temperatures of a cometary nucleus, closely measured by a spacecraft, can be compared with theoretical models meant to predict the thermophysical properties of the nucleus at various depths and with previous thermal observations carried out by both Earth-based and space-based telescopes, with the goal of validating and/or improving those models<fnref n="1" />. Furthermore, the temporal evolution of the surface temperatures, primarily driven by insolation, is critical in triggering the activity of a comet, allowing the migration of volatile compounds from the interior of the nucleus and the production of gas and dust observable with different techniques².

In the past, thermal surveys of minor bodies were carried out using ground-based and space-based facilities. However, due to the limited spatial resolution of those observations, in most cases the derived thermal properties were only global averages, and the determined temperatures were highly model dependent. Before Rosetta, direct measurements of cometary nuclei surface temperatures were obtained during short periods for a handful of comets, namely 1P/Halley³, 19P/Borrelly⁴, 9P/Tempel 1¹ and 103P/Hartley 2⁵, with a maximum spatial resolution of approximately 30 m per pixel.

Here we study the surface temperature distribution of the nucleus of comet 67P/Churyumov-Gerasimenko as derived by the Visible InfraRed and Thermal Imaging Spectrometer, Mapping channel (VIRTIS-M, hereafter VIRTIS)⁶ in Rosetta's early global mapping phase after comet encounter. These data cover the pre-perihelion period from 1 August to 23 September 2014, when the heliocentric distance decreased from 3.62 to 3.31 au and the spacecraft was in the altitude range 61–13 km above the surface, resulting in a spatial resolution from approximately 15 to 3 m per pixel (most data showing a resolution of 13 to 15 m per pixel). In this period, the solar phase angle ranged from 17° to 93°, which

<footnote n="1">¹Istituto Nazionale di AstroFisica - Istituto di Astrofisica e Planetologia Spaziali (INAF-IAPS), Rome, Italy. ²Institute of Planetary Research, Deutsches Zentrum für Luft- und Raumfahrt (DLR), Berlin, Germany. ³Agenzia Spaziale Italiana - Space Science Data Center (ASI-SSDC), Rome, Italy. ⁴Jet Propulsion Laboratory (JPL), Pasadena, CA, USA. ⁵Dipartimento di Matematica e Fisica “Ennio De Giorgi”, Università del Salento, Lecce, Italy. ⁶Institute of Physics and Astronomy, University of Potsdam, Potsdam, Germany. ⁷Laboratoire d'Etudes Spatiales et d'Instrumentation en Astrophysique (LESIA), Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université Paris Diderot, Meudon, France. ⁸Università degli Studi di Napoli “Parthenope”, Dipartimento di Scienze e Tecnologie (DiST), Centro Direzionale, Napoli, Italy. ⁹Université Grenoble Alpes, CNRS, Institut de Planétologie et d'Astrophysique de Grenoble (IPAG), Grenoble, France. ¹⁰European Space Agency - European Space Research and Technology Centre (ESA-ESTEC), Noordwijk, Netherlands. ¹¹Bear Fight Institute, Winthrop, WA, USA. ¹²Deceased: S. Fonti. *e-mail: federico.tosi@inaf.it</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1** | VIRTIS-derived global temperature map of comet 67P. **a**, Global map of surface temperatures retrieved from VIRTIS measurements in August and September 2014, in a cylindrical projection and with a fixed angular binning of 0.5° by 0.5° in both latitude and longitude. Each bin collects the average temperature recorded during a full diurnal rotation period. We show only illuminated points that have at the same time: solar incidence angle $i < 80°$, emission angle $e < 80°$ and accuracy better than 2 K. A 'red temperature' colour code is used, such that coldest temperatures are shown in the dark red tones while warmest temperatures are whitish. The macroregions 'small lobe', 'big lobe' and 'neck' (large concavity where the two lobes join) are labelled in yellow. A white cross below the small lobe label indicates the location of the Abydos landing site of Philae. A yellow circle highlights the location of extreme temperature values. **b**, Four projections of the map shown in **a** on a digital shape model of the nucleus of comet 67P, obtained with stereo-photogrammetric methods and degraded to four million triangular plates ('facets')$^{22,47}$. We show different three-dimensional views of the nucleus, with the macroregions highlighted in yellow.</figcaption></figure>

allowed variable illumination conditions to be explored (Methods and Supplementary Fig. 1).

We adopted a Bayesian approach to retrieve surface temperature from the thermal radiance measured in the infrared spectral range 4.5–5.1 µm (refs. $^{7–9}$). In all VIRTIS hyperspectral images, this retrieval is applied separately on each pixel unaffected by detector saturation in the above range. Due to its spectral range limited to $\le 5.1$ µm, combined with the very low surface albedo at visible to near-infrared wavelengths$^{10,11}$ and the high thermal emissivity, VIRTIS can sample temperatures within the uppermost radiatively active surface layer, as thick as some tens of micrometres. Within this wavelength range, the radiance is not a linear function of temperature, resulting in a preferential weighting of the hottest subpixel

features. Thus, the temperature obtained by VIRTIS is representative of the subpixel regions at higher temperatures weighted by their effective areas.

The lower limit of temperatures sensed by VIRTIS is dictated by the in-flight instrumental noise, which may vary over time depending on some instrumental parameters. In the time period studied here, this lower bound was 156 K on average (Supplementary Fig. 2), which restricts our coverage to the dayside of the nucleus, including areas that recently underwent shadowing or recently exited from shade. It should be noted that the VIRTIS-derived temperatures are not directly comparable with those measured by the Microwave Instrument for the Rosetta Orbiter (MIRO)$^{12}$, which sampled shallow subsurface temperatures (at millimetre to centimetre depths) both on the dayside and the nightside of 67P$^{13,14}$. Nonetheless, the combination of VIRTIS- and MIRO-derived temperatures can provide valuable constraints on the vertical thermal gradients.

## VIRTIS temperature maps

Due to the large obliquity of the comet$^{15}$ that causes pronounced seasonal effects, at the time of the observations northern latitudes were in the summer season, while southern latitudes below $45°$S were in the winter season and experienced polar night. Under these conditions, surface temperatures at northern latitudes reached values as high as 230 K, following the surface topography and the distribution of solar incidence angles (Fig. 1).

The surface temperatures only moderately depend on latitude. The morphological regions$^{16}$, Ash, Babi, Hapi and Seth, located in the 'neck', that is, the surface area connecting the two lobes, and the Ma'at region located in the small lobe close to the neck, are the locations where maximum surface temperatures were detected, whereas other regions at comparable latitudes attained lower temperatures. However, the surface temperature does not appear to follow black-body behaviour. In the ideal black-body approximation (black dashed line in Supplementary Fig. 3), the trend of surface temperature follows the fourth root of the cosine of the solar incidence angle, but most VIRTIS data show a deviation from this simple trend. The surface temperature for each morphological region is shown in Fig. 2 as a function of the solar incidence angle. The regions in the neck or close to it (Ash, Babi, Hapi, Ma'at and Seth) display a larger variability than expected for a black-body, which can be an indication of the larger contribution of additional thermal sources associated with macroscopic self-heating. In this case, the surface temperature is affected by the presence of large concavities, which result in local infrared thermal flux enhancement by repeated reflections from mutually facing surface areas and their thermal emission$^{17,18}$. This is also an indication that self-heating is predominant with respect to the effects of heat absorption and heat transport, which confirms a low thermal inertia for the surface layers$^{19}$. Conversely, regions located at equatorial and southern latitudes in the two lobes, which were experiencing grazing illumination conditions, such as Aker, Anubis, Atum, Imhotep and Kephry, show patterns more similar to those expected in the case of a black-body. Various regions (for example, Aker, Apis, Atum, Bastet, Khepry, Maftet, Nut and Serquet) have a bimodal pattern, with one part of the covered surface closer to the behaviour expected for a black-body and at least another part that substantially deviates from this trend. The Hatmehit region reveals three separate patterns (Fig. 2). While these multiple patterns could be ascribed to the non-uniformity of the morphology within a given region, and thus to a variable role of the mutually illuminating areas, temperature values higher than the ideal case of a black-body, observed in all of the explored morphological regions for large solar incidence angles, are evidence for small-scale roughness. Such a small-scale roughness by definition occurs beyond the spatial resolution of most optical images used to reconstruct the topography$^{17,18,20}$, and is expected based on previous thermal flux measurements of low thermal inertia, airless
<figure><img src="image_3.png" /><figcaption>**Fig. 2** | Temperature versus solar illumination in several morphological regions. Surface temperatures measured by VIRTIS in morphological regions of the nucleus¹⁶ that were experiencing sunlight in August-September 2014, as a function of the solar incidence angle. Each panel represents a different region (in alphabetical order from left to right and from top to bottom). Data points correspond to individual VIRTIS pixels. We show only illuminated points that have at the same time: solar incidence angle *i* < 80° and emission angle *e* < 80°. Colour is related to effective latitude (effective latitude and effective longitude are defined based on the orientation of the local surface normal of a point on the surface with respect to the Sun), such that highest latitudes are shown in a red colour while lowest latitudes are coloured in dark blue-violet. The black dashed line represents the ideal black-body behaviour. In the range of heliocentric distances 3.62–3.31 au, radiative equilibrium temperatures in the subsolar point (that is, the point of the surface with solar incidence angle equal to 0°) span the interval 206–215 K, having assumed an average Bond albedo of 0.0155 (ref. ²⁴). See also Supplementary Fig. 2 for a single scatterplot showing all regions together. The last panel in the lower right corner shows the placement of the different regions of the nucleus in a cylindrical projection.</figcaption></figure>

bodies⁹,²¹. Depending on illumination and observation geometry, mini-concavities can produce increased thermal fluxes due to the effects of re-radiation between wall segments⁹.

By splitting the thermal maps into different time spans, each covering two hours in true local solar time (tLST), the surface does not exhibit outstanding thermal anomalies (Supplementary Fig. 4). These maps show a fairly uniform distribution of surface temperatures, mildly dependent on the latitude, while small-scale roughness effects can largely account for residual variations. The big lobe and the neck, despite their mutual proximity, show regional differences, with the neck exhibiting a faster temperature change than

the big lobe, both in the morning and in the evening hours. Finally, Imhotep and its surrounding regions in the sunlit hemisphere show the lowest values of maximum daily temperature, driven by grazing sunlight around local noon during this season.

The 'Abydos' site where the Philae lander finally settled (8.2°S, 358.2°E in the Cheops coordinate frame²²) in this period was located in the borderland between dayside and permanent night-side. The VIRTIS measurements of a 1° by 1° area at this location (~30 m wide) showed maximum temperatures of 207 K. However, the VIRTIS pixel sampling was typically ten times as coarse as the physical size of the lander, which prevents a direct comparison with
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

**Fig. 3** | Surface temperature versus albedo. **a–c**, Scatterplots of surface temperature versus VIRTIS-derived SSA¹⁰ at three different wavelengths: 0.55 µm (**a**), 1.21 µm (**b**) and 2.00 µm (**c**). We show only illuminated points that have at the same time: solar incidence angle *i* < 20°, emission angle *e* < 80° and accuracy on temperature retrieval better than 2 K. The average temperature on the dayside of the nucleus during the time period centred on local noon is 213 ± 3 K. The main cluster of points is homogeneous and no substantial correlation emerges between temperature and albedo at the global scale. However, a secondary cluster of points can be noticed (highlighted by the red circle in **b**), with an average temperature of 219 ± 1 K, that is, at least 5 K higher than the main cluster. This group is located in the Hapi region (centre latitude 13.1° N, longitude 40.3° E). **d–f**, Density scatterplots corresponding to **a–c**, respectively, with colours accounting for the number of data points per bin of size 0.5 K × 0.005 SSA, increasing from blue to red.

the diurnal range 90–130 K measured in situ by Philae between 12 and 14 November 2014²³.

The temperature recorded on the dayside of the nucleus, for solar incidence angles <20° (that is, around local noon) and emission angles <80°, is 213 ± 3 K on average (Fig. 3). No correlation arises on a global scale between VIRTIS-derived temperature and single scattering albedo (SSA), both at visible and near-infrared wavelengths (Fig. 3). However, a deviation from this general behaviour may be observed at the local scale. A secondary group of points, centred in the Hapi region (13.1° N, 40.3° E), shows albedo slightly higher than average, associated with higher temperatures. This area is slightly brighter than the two main lobes¹⁰,²⁴, and is characterized by smooth terrains¹⁶. Extreme values, reaching up to 230 K, were recorded in two VIRTIS observations in the Seth region (59.5° N ± 0.3°, 204.7° E ± 0.3°), corresponding to the pit area named ‘Seth_01’²⁵ (Fig. 1, yellow circle). This location has SSA values of 0.06 at 0.55 µm, 0.10 at 1.2 µm and 0.12 at 2 µm, which are very close to the average values found at the global scale¹⁰. Therefore, these extreme temperatures are most likely due to the peculiar topography of the pit combined with the instantaneous solar illumination shortly after local noon, which enhances the self-heating effect.

## Thermophysical modelling

On 22 August 2014, VIRTIS obtained seven consecutive snapshots of the nucleus, from an altitude of 60 km above the surface (spatial resolution 15 m per pixel). These data cover large areas of the neck and the two lobes, following them during approximately 15% of the rotational period. The corresponding temperature images reflect the combination of local topography and solar illumination at the time of the observations (Fig. 4a). Temperatures predicted in this region were modelled by using a quasi-stationary approach¹⁸ described in the Methods. In Fig. 4b, we assume a 1-mm-thick dust layer and dust agglomerates of 50 µm radius over a deeper, sublimating water-ice layer (corresponding to model C in Keller et al.¹⁸), and we further add small-scale roughness in the form of 50% surface coverage of nearly hemispherical mini-concavities. This assumption provides temperature values closer to the measured data, compared with the case of exposed ‘dirty’ water ice sublimating directly from

the surface (model A in Keller et al.¹⁸) or a very thin layer of small dust agglomerates above the ice–dust mixture of the nucleus (model B in Keller et al.¹⁸). Analysis of the observed—calculated residuals demonstrates that a two-layer dust-ice model, with the addition of small-scale roughness, is also more suited than a simple grey-body stationary case with no sublimating material in reproducing the observed temperatures (Supplementary Fig. 5 and Supplementary Tables 2 and 3).

Rosetta's early images have shown that the main jets of water vapour originated from Hapi²⁶. Following the activity observed there by the Optical, Spectroscopic, and Infrared Remote Imaging System (OSIRIS)²⁷, we modelled temperature profiles for one specific location in Hapi (59.9° N, 19.1° E) using a one-dimensional model that solves the dynamic equation of heat transfer²⁸,²⁹ (Methods). Modelled results are compatible with VIRTIS-derived temperature values for that location only if subpixel surface roughness²⁰ is considered (Supplementary Fig. 6). These inferred thermal properties are consistent with spectral data of a dehydrated microporous surface³⁰⁻³² and with an overall low thermal inertia as inferred from MIRO and VIRTIS observations¹³,¹⁴,¹⁹, emphasizing the insulating role played by a mostly dusty porous surface layer.

## Thermal gradients in Hapi

Owing to the prominent concavity of the nucleus, the Hapi region is a peculiar place that undergoes large temperature gradients triggered by sudden shadowing alternately produced on the neck by the two lobes during the maximum diurnal insolation as clearly seen in the sequence shown in Fig. 4a and Fig. 5a. On airless terrestrial bodies, daytime shadowing may produce higher temporal temperature gradients than those caused by diurnal sunrise/set if the topography is optimally oriented³³. On small asteroids, an abrupt thermal gradient induced by day–night temperature variations (~2.8 K min⁻¹) was found to be an effective way to fracture rocks larger than a few centimetres more quickly than do micrometeoroid impacts³⁴.

Spatial (lateral) temperature gradients can be systematically computed as the absolute value of the difference in temperature between adjacent VIRTIS pixels. As a result, at the 15 m scale, spatial gradients of terrains experiencing fair solar illumination are on
<figure><img src="image_10.png" /><figcaption>**Fig. 4** | Measured and modelled surface temperature. **a**, Temperature image retrieved from VIRTIS hyperspectral image 11_00367329515, acquired on 22 August 2014 from an altitude of 60 km above the surface (spatial resolution 15 m per pixel). From top to bottom, seven hyperspectral snapshots were obtained, which cover approximately 15% of the rotational period. In the left edge, UTC stamps account for the start time of each of the seven sections. The coldest temperatures are shown in dark red tones while the warmest temperatures are whitish. **b**, Result of thermophysical modelling<sup>18</sup>, with addition of small-scale roughness in terms of 50% surface coverage of mini-concavities, represented on the digital shape model of comet 67P<sup>22,47</sup> for the same period of the VIRTIS observations shown in **a**. The grey (shaded) colour refers to points with modelled temperature values lower than 160 K, for which neither VIRTIS data nor the model can provide accurate values. Other spurious grey dots are due to the projection of the data onto the shape. **c**, Difference between temperature values measured by VIRTIS and represented in **a**, and theoretical temperature values modelled in **b**. In this presentation, the blue and cyan colours highlight regions of the nucleus where the modelled temperature values are significantly greater than the measured ones, and vice versa for the colours from yellow to red. The green colour represents a substantial agreement between observed and calculated data. **d**, Distribution of the solar incidence angle values on the digital shape model used to characterize the thermal data. A rainbow palette is used, such that small values are blue while large values are red. **e**, Distribution of the emission angle values on the digital shape model (same colour code as **d**). **f**, Distribution of tLST values on the digital shape model. A red-yellow-blue palette is used, such that bright tones indicate day, dark night, reddish colours are pre-noon and bluish colours are post-noon.</figcaption></figure>

average $0.8 \pm 0.5$ K per pixel or $0.05 \pm 0.03$ K m$^{-1}$ (Fig. 5b), but at the border between lit and shaded terrains in the neck region they reach up to $6.6 \pm 1.7$ K per pixel or $0.4 \pm 0.1$ K m$^{-1}$, with extreme values as high as 14 K per pixel or $0.9$ K m$^{-1}$ (Fig. 5b). However, temporal temperature gradients result from the evolution of the surface temperature at a given point of the surface at different times. From the VIRTIS data, it emerges that the dayside of the nucleus shows typical values of $0.1$ K min$^{-1}$, except in the neck area where sudden shadowing takes place. In this case, the temporal gradient increases up to $2.0$ K min$^{-1}$, that is 20 times as high as the typical value (Fig. 5c and Supplementary Fig. 7). A comparison of these data with the thermophysical modelling previously described for Hapi (Supplementary Fig. 6) indicates that temporal thermal gradients can be actually as large as $9$ K min$^{-1}$ in those regions where sudden daytime shadowing takes place, although the limitations of our instrument did not allow the minimum nighttime temperatures to be reached<sup>14</sup>.

Such large temperature thermal gradients have strong implications for thermal stress weathering<sup>33</sup>, although the competing role of other processes, such as erosion by dust, dust deposition and water-ice sublimation<sup>35</sup>, can mitigate or intensify its contribution to the shaping of the nucleus surface.

## Thermal properties evolution

Moving away from the neck the effects of macroscopic self-heating are weaker, making it easier to highlight seasonal effects. In this respect, we considered a 10° by 10° wide area in the Imhotep region, which appears relatively smooth at the scale of optical imagery<sup>36</sup>, and we selected VIRTIS thermal data covering this area in the overall time period between August 2014 and May 2015. We limited our search to those points providing the best measurement conditions (solar incidence angle and emission angle < 80°, tLST between 11 h and 15 h, and uncertainty of retrieved temperature < 10 K). By applying these selection criteria, we found a total of 16 measurement points mostly acquired in August 2014, while the last measurements were taken approximately 8 months later in April 2015 (Supplementary Fig. 8). Among these, we selected only those points that: (1) most closely approximate a linear distribution in temperature versus time and (2) are also broadly spaced out in time. The outcomes of this selection are respectively designated 'point A' (12.98°S, 144.21°E, pixel size 21.5 m) and 'point B' (11.04°S, 141.03°E, pixel size 35.4 m) (Supplementary Fig. 8). Point A has an approximately 194 K temperature as measured by VIRTIS on 7 August 2014 at 13.1 h tLST and at 3.593 au from the Sun, whereas point B has an
<figure><img src="image_11.png" /><figcaption>**Fig. 5** | Spatial and temporal thermal gradients. **a**, Same sequence of seven thermal images of the nucleus of the 67P, acquired by VIRTIS on 22 August 2014 and already shown in Fig. 4a. **b**, Map of spatial thermal gradients, calculated on the basis of the temperature image shown in **a**. For each pixel in the image, the mean gradient is given by the absolute value of the difference between its temperature and the temperature of adjacent pixels in the horizontal and vertical directions. A rainbow palette is used, such that small values are violet while large values are red. The interface between sunlit and shaded regions in the neck area has a distinct contour, which correspondingly shows high values of the spatial thermal gradient. **c**, Map of temporal thermal gradients, obtained by using the digital shape model of the comet²². For each VIRTIS pixel as seen in **a**, we consider only the facet of the digital shape model intercepted by the pixel centre's line-of-sight extended to infinity (while each pixel covers multiple facets). When a given facet of the shape is observed multiple times by VIRTIS within a period of time ranging from minutes to tens of minutes, it is coloured other than black. The colour code is linked to the absolute value of the temporal thermal gradient, with the highest values, represented by the orange and red tones, concentrated in the area of the neck where sudden daytime shadowing occurs. Here the measured temporal gradients reach up to 2 K min⁻¹, exceeding by a factor of 20 those typical of the rest of the sunlit hemisphere. However, the actual temporal gradients in the neck can be several times greater than this value in those shaded areas where the surface temperatures are below the lower limit of VIRTIS sensitivity (see Supplementary Fig. 6).</figcaption></figure>

approximately 232 K temperature as measured by VIRTIS on 10 April 2015 at 12.2 h tLST and at 1.884 au from the Sun. Theoretical temperature profiles in these two points were modelled adopting the one-dimensional thermophysical approach used for Hapi (see details in Methods) and assuming two cases: (1) an uppermost layer composed of loose dust and (2) an uppermost layer composed by a mixture of dust and water ice. On the surface, VIRTIS point A is best modelled by assuming case 1, whereas VIRTIS point B can only be fitted by assuming case 2, which is representative of an increasing contribution of sublimation at that time (Fig. 6). We compared these profiles to an observation of the Imhotep region obtained by the submillimetre channel of MIRO on 13 August 2014 at 11.0 h tLST, yielding a shallow subsurface (down to approximately 1 cm depth) physical temperature of approximately 127 K. To reconcile the different

ground resolutions of VIRTIS and MIRO, we averaged multiple VIRTIS pixels to involve the same area observed by MIRO.

The combined VIRTIS and MIRO measurements indicates that the uppermost layers of the smooth terrain in Imhotep, in the first half of August 2014, were compatible with an assemblage of dry and loose dust grains in its uppermost layer. In the later phases, the surface temperature indicates a temporal evolution that brings a larger amount of water to the surface as indicated by point B. This is in good agreement with the results obtained by VIRTIS³⁷⁻⁴⁰ and OSIRIS²⁴,⁴¹,⁴², which showed a global enrichment in the amount of water ice in the uppermost layers throughout the pre-perihelion phase. Furthermore, a preliminary analysis based on both submillimetre- and millimetre-wavelength MIRO data confirms evidence for seasonal change in Imhotep in the two-year period 2014–2016⁴³.
<figure><img src="image_12.png" /><figcaption>**Fig. 6 |** Seasonal evolution of surface temperature in Imhotep. Theoretical temperature profiles as a function of tLST, compared with the temperatures measured by VIRTIS at 'point A' and 'point B', located in a smooth terrain within the Imhotep region. Point A was observed by VIRTIS on 7 August 2014 at 3.593 au from the Sun. Both this VIRTIS temperature determination and the MIRO measurement of the same region acquired six days later are best matched by a thermophysical model that assumes only loose dust in the uppermost cm-thick layer (red solid curve and red dashed curve, respectively). VIRTIS point B, which was observed on 10 April 2015 at 1.884 au from the Sun, is best fitted by a model that assumes a mixture of dust and ice in the uppermost layer (blue solid curve). For each measurement point, the error bar expresses the formal error related to random variations of the signal. Assuming only loose dust in the uppermost cm-thick layer in April 2015 would provide a much higher maximum temperature value (blue dash-dotted curve).</figcaption></figure>

## Implications

Our present results, combined with previous VIRTIS, MIRO and OSIRIS findings, indicate a global uniformity, both spatially and temporally, of the surface properties of the 67P nucleus. The spatial distribution of the surface temperatures, observed in roughly two months during the pre-perihelion phase and spanning approximately 0.3 au in heliocentric distance, is mainly driven by latitude and insolation, except in the areas where self-heating is important, indicating a largely homogeneous composition and physical properties of the nucleus as a whole. Over a longer timescale, the seasonal effects due to the combination of the increasing diurnal insolation and the simultaneous reduction of the heliocentric distance translate into an increasing global abundance of water ice at the surface caused either by transport mechanisms from deeper layers or by surface erosion or a combination of both.

The low thermal inertia consistent with the observed temperatures, the vertical gradients resulting from the combined VIRTIS and MIRO observations, and also the large scale uniformity, all together point to the fact that most of the active physical processes take place in the uppermost few centimetres of the nucleus body. Even the intense surface erosion caused by the removal of gas and dust due to solar heating is a process active instantaneously layer by layer, where each layer is few centimetres thick. Only the uppermost layers are then those most affected by an evolution in composition and physical properties; this implies the possibility of having a largely pristine nucleus interior already at a limited depth of few metres or less.

## Methods

**The VIRTIS instrument.** VIRTIS-M (Visible InfraRed and Thermal Imaging Spectrometer, Mapping channel) is the visible-infrared hyperspectral imaging channel aboard Rosetta¹. The instrument performs 0.25–5.1 µm spectroscopy by

using two separate detectors: a CCD for the 0.25–1.0 µm range and a HgCdTe bidimensional array for the 1.0–5.1 µm range. The average spectral sampling is equal to 1.8 and 9.7 nm per band for these two channels, respectively. Both spectral channels share the same telescope equipped with a scan mirror to build hyperspectral images (so-called cubes). The instrumental field of view is equal to 3.67° with an instantaneous field of view of 250 µrad. The slit aperture is therefore acquired by 256 spatial samples.

**The VIRTIS-M dataset.** In the overall time period between 1 August 2014 and 3 May 2015, corresponding to Rosetta orbiter's Medium Term Plan (MTP) phases MTP006 to MTP015, VIRTIS-M acquired a total of 1,059 infrared hyperspectral images. A sudden malfunction of the VIRTIS-M infrared focal plane's cryocooler, which occurred on 3 May 2015, inhibited the use of infrared data acquired by the mapping channel from this moment on. In the remaining part of the mission, VIRTIS-M acquired data only in its visible channel 0.25–1.0 µm, while the infrared spectral range 2–5 µm was covered by the high-spectral-resolution channel VIRTIS-H, with no imaging capability.

In the aforementioned period, the heliocentric distance decreased from 3.62 to 1.72 au and the latitude of the subsolar point decreased from 44.8° to 2.7°, that is, moved from the northern summer to the equinox. Both the spatial resolution of VIRTIS-M data and the solar phase angle showed substantial variations between one MTP and another and within a given MTP, due to the changing trajectory of the Rosetta orbiter relative to the comet's nucleus. As a result, the VIRTIS-M pixel resolution varied between approximately 2 and 96 m and the solar phase angle overall ranged from approximately 9° to 120° (see Supplementary Fig. 1).

In this work, most of the data analysis relies on 266 hyperspectral images acquired by VIRTIS-M in the time period between 1 August 2014 and 23 September 2014, corresponding to MTP006 and MTP007. In MTP006, VIRTIS-M acquired data from 1 August 2014, 04:23:44 UTC to 2 September 2014, 04:45:57 UTC, while in MTP007 VIRTIS-M acquired data from 2 September 2014, 12:00:23 UTC to 23 September 2014, 05:05:59 UTC. During MTP006 and MTP007, heliocentric distances decreased from 3.62 to 3.31 au and the solar phase angle ranged from 17° to 93°. The mean spatial resolution of this VIRTIS-M subset of data is 14.8 µm per pixel, and the median spatial resolution is 13.3 µm per pixel.

Accuracy in the temperature retrieval depends on the level of the emitted thermal radiance. By considering only VIRTIS-M-retrieved temperature values above 180 K, the median of the formal errors (that is, related to random variations of the signal) associated with these values is 2.9 K, while for temperature values above 210 K the median of the formal errors decreases to 0.8 K.

VIRTIS-M scan 11_00367329515, shown in Figs. 4a and 5a, was acquired on 22 August 2014 in the time period between 11:59:45 and 13:54:52 UTC. It is a hyperspectral image of 256×658 pixels, obtained with a repetition time (that is, period between the start acquisition time of one horizontal line and the start acquisition time of the following horizontal line) of 10 s and an effective exposure time of 3 s per line. During the approximately 115 min period in which the hyperspectral image was acquired, the nucleus of the comet was observed by VIRTIS-M seven times. Supplementary Table 1 summarizes the timings for each subframe seen within the overall image.

**Thermophysical modelling of temperature images (refer to Fig. 4b and Supplementary Fig. 5).** To compute high-resolution temperature images for the nucleus of comet 67P in Fig. 4b, we adopted the simplified and thus computationally efficient thermal model described in Keller et al.¹⁸,⁴⁴⁻⁴⁶. This model assumes a two-layer medium (inert dust over sublimating ice) and a quasi-stationary state, where the equilibrium temperature is reached in a short time, compared with the typical timescales for changes in illumination. In the case of comet 67P, this condition is satisfied because of the low thermal conductivity and overall low thermal inertia. With this approach, the thermal behaviour of the comet is controlled by the instantaneous energy input (computed by considering both direct insolation and indirect visible and infrared irradiation), in addition to parameters such as the thickness of the dust layer, albedo and so on, while it does not depend on its thermal history.

In Keller et al.¹⁸, four cases are considered: (1) a simple grey-body stationary case with no sublimating material and no heat transport to the interior; (2) a stationary model with exposed 'dirty' water ice sublimating directly from the surface and no heat transport to the interior (model A); (3) a quasi-stationary two-layer model with a very thin layer (50 µm) of small dust aggregates above the ice-dust mixture of the nucleus interior, about the minimum to make the ice invisible to spectroscopy (model B); and (4) a quasi-stationary two-layer model with a 1-mm-thick dust layer with dust agglomerates of 50 µm radius, superposed to a deeper, sublimating water-ice layer (model C). Model C, which approaches the limit case of the grey-body representing a relatively inert dust surface as observed during the beginning of the Rosetta rendezvous phase with comet 67P¹⁸, was expected to produce a good match to the temperature values measured by VIRTIS. The heat transfer equation is solved in one dimension accounting for sublimation and heat transport due to conduction, radiation and gas diffusion. The solution to the heat transfer equation, consisting of the surface and ice-layer temperatures, as well as the water vapour production rate, is pre-computed once for the range of input energies occurring over the orbit. The results are stored in a look-up table,
which allows rapid access to the temperatures, without needing to re-compute the
heat transfer equation for each facet of the digital shape model.

Some parameters are not free, but they are rather derived as a function of
other quantities, in particular the size of the dust grains and the number of layers.
For example, the permeability of the gas is computed with an empirical function
that describes the permeability as a geometric quantity<sup>41,45</sup>. The actual thermal
conductivity is given by the sum of three components: (1) conductivity through
the solid medium (fixed at the value determined with laboratory experiments);
(2) conductivity through thermal emission—calculated; and (3) conductivity
through diffusion—calculated but neglected because very small compared with
the other two contributions. The resulting conductivity is a complex function that
depends on the level of activity. Upper limits under the illumination conditions of
67P are low (2 × 10<sup>-3</sup> W K<sup>-1</sup> m<sup>-1</sup> for the upper dust layer and 3 × 10<sup>-3</sup> W K<sup>-1</sup> m<sup>-1</sup> for
the dust and ice mixture below), which justifies the assumption of short thermal-
equilibrium timescales. Upper limits for thermal inertia can be
derived consequently.

To model and present synthetic temperature values, we relied on the
latest digital shape model of the nucleus of comet 67P obtained with stereo-
photogrammetric methods (SPG-SHAP7 v1.6)¹² and decimated to about 300,000
triangular plates ('facets') with an average length of approximately 3 m per side,
mostly concentrated in the region covered by the best VIRTIS-M observations
during scan I1_00367329515. In this way, we matched the resolution of the
VIRTIS-M observations, while keeping the processing time within
reasonable limits.

The thermophysical modelling carried out according to Keller et al.<sup>18</sup> produces a thermal snapshot of the nucleus at a given time. However, VIRTIS-M is a line scanner and as such it acquires the scene one line at a time, which typically requires tens of minutes for a hyperspectral image to be built. This time dependence complicates a direct comparison between observed and modelled data. To minimize the time difference between the observed and modelled temperatures, we computed thermophysical models at three epochs for each of the seven sections that make up the VIRTIS image in Fig. 4b (for a total of 21 models) and whose exact timings are reported in Supplementary Table 1. The three epochs account respectively for the exposure start time plus three minutes, the intermediate exposure time, and the stop time minus three minutes. In this way, for each VIRTIS scan line there is an output of synthetic temperature that does not deviate more than three minutes from the time at which the VIRTIS data were acquired. To project each temperature model onto the observer's plane, we used a dedicated three-dimensional data visualization tool<sup>19</sup>, which allows a proper characterization of both measured and modelled data on the bilobate nucleus of comet 67P.

For a given illumination geometry, surface roughness significantly affects the
apparent temperature spatial distribution as seen by the observer, which depends
on the observing direction, and also affects the shape of the spectral energy
distribution (for example, the brightness temperature at different wavelengths)⁵.
Concerning airless bodies, until the mid-1990s, most models represented surface
roughness with unresolved concavities defined by their opening angle and the
fraction of the surface they occupy²⁰. Later, alternative models were proposed, such
as Gaussian random surfaces¹⁸, circular paraboloid pits¹⁹ and fractal terrains²¹. For
our study, to model subpixel surface roughness we have applied a simple approach
already used in the past to correct the thermal emission measurements of the
Martiian satellites Deimos and Phobos²⁵, which assumes a variable density of nearly
hemispherical concavities at subfacet scale, distributed randomly
across the surface.

For the purpose of comparing observed and calculated temperature images, we
have applied two different models: (1) model C of Keller et al.<sup>18</sup>, where we added
small-scale roughness, modelled by means of a 50% surface coverage of mini-
concavities, and (2) the grey-body model of Keller et al.<sup>18</sup> (which we label 'NoIce'
for simplicity), adding small-scale roughness in terms of a 10% surface coverage
of mini-concavities. The mini-concavities surface coverage percentages were
determined by fitting the respective models to the observed average temperatures.

It should be noted that, in principle, a given solution might be obtained starting
from different combinations of parameters. However, a comprehensive exploration
of the whole parameter space (model type: grey-body, A, B or C<sup>18</sup>, thickness of the
uppermost dust layer, dust grain size and roughness) would be unpractical in terms
of processing time. Therefore, we followed a pragmatic approach: we started from
two realistic scenarios, (1) model C and (2) grey-body, and changed the small-scale
surface roughness until a substantial agreement between measured and modelled
results was achieved.

The results are shown in Supplementary Fig. 5, and are summarized in
Supplementary Tables 2 and 3, respectively. All statistics are calculated only for
T > 160 K, because uncertainties on both the measured and modelled temperature
values are very large below this threshold. As expected, the 'NoIce' model
produces higher temperatures both compared with the measured data and to the
nominal model C, because it is in fact an upper limit for the thermal radiation
that can be emitted by the nucleus. In those areas where 'outliers' concentrate, the
observed – calculated difference values exceed 35 K. Such large deviations between
observed and calculated data correlate with areas experiencing large emission
angle values (Fig. 4). The statistics for the whole sequence can be computed either
including or excluding these outliers. By excluding the outliers, both the mean and

the median values of the distribution demonstrate that model C of Keller et al.<sup>18</sup>,
with the addition of small-scale roughness, provides a better approximation to
the measured data (typically by 2 K) with respect to the 'NoIce' model. However,
the percentage of outliers increases throughout the sequence going from the first
to the seventh subframe, because of the changing observation geometry due to
the rotation of the nucleus: as the shaded regions in the neck become increasingly
visible, the model hardly fits the measured data, even if the measured
values are >160 K.

Thermophysical modelling of theoretical temperature profiles (refer to Fig. 6 and Supplementary Fig. 6). Stationary and quasi-stationary thermophysical models such as Keller et al.<sup>18</sup> can not predict accurate temperature values in shaded areas, and therefore are unsuited to reproduce the thermal behaviour on the nightside of the comet's nucleus and to estimate day-night thermal gradients. To achieve this goal, one needs to use a model that solves the dynamic equation of heat transport. While a detailed thermophysical modelling of comet 67P is treated elsewhere<sup>51</sup>, we applied the 'Rome' thermophysical model<sup>28,29</sup> to compute synthetic temperature profiles and infer seasonal trends for specific locations on the nucleus. Using the latest digital shape model of comet 67P (SPG-SHAP7 v1.6)<sup>32</sup> we simulated the thermal behaviour of specific surface elements ('facets') identified by their planetocentric coordinates.

This one-dimensional model solves the heat transport equation coupled with the gas diffusion equation, expressing mass and energy conservation, respectively. The comet's nucleus is treated as a porous mixture of silicate dust grain aggregates, H₂O ice and CO₂ ice. The model accounts for the solar radiation reaching the surface, heat conduction in the interior, heat advection by gases, sublimation and recondensation of ices, amorphous-crystalline ice transition, surface erosion and gas and dust fluxes from the surface. This model may include the self-heating effect, which can not be neglected when the local topography substantially deviates from the average curvature of the body. The dust grains are distributed in different size classes, characterized by specific physical and thermal properties. They can be released by the sublimation of the ices and then undergo the drag exerted by the escaping gas, becoming part of the dust flux from the surface, or can instead redeposit on the surface. A given set of a priori parameters, describing the initial state of the body and the properties of the matter of which it is composed, defines a 'case'. In defining the boundary conditions, the body can either have a homogeneous composition throughout the entire radius at the beginning of the calculation, or separate layers with composition changing as a function of depth. In each layer, this composition can evolve over time.

For a specific surface location, SPICE-based routines<sup>52</sup> and navigational
databases, the so-called SPICE kernels, provide illumination conditions at any
time along the orbit. At the surface, the temperature is the result of the energy
balance between the insolation, the energies radiated into space, conducted into
the interior and used to sublimate ice (if present). Subpixel surface roughness
is taken into account in a very simple way through the small-scale surface
roughness parameter ξ (ref. <sup>20</sup>), which is an allowable approach as long as the
solar illumination and emission angles of the observations are small enough and
the temperature points observed at high phase angles have been filtered out and
discarded in the computation. Outputs of the model are: gas and dust fluxes and
profiles of temperature, density and composition as a function of time. The system
of equations is solved with a finite-differences method.

A direct comparison between these theoretical temperatures and VIRTIS-
measured temperatures may not be straightforward. While modelled temperatures
refer to a specific 'case', in the spectral range sampled by VIRTIS the temperature
is a nonlinear function of radiance, with a larger weight of the hottest subpixel
features (such as dust-covered, well illuminated zones) weighted by their
effective areas.

In this case the 'Rome' model was applied to three distinct points on the
nucleus, specifically one location in Hapi (59.9°N, 19.1°E) and two locations,
tagged as 'point A' and 'point B', respectively, in a smooth terrain in the Imhotep
region (12.98°S, 144.21°E; 11.04°S, 141.03°E).

In the model, the density does not necessarily coincide with the bulk density of the nucleus and may vary depending on the composition of the layer and its porosity (average diameter of the pores). An initial porosity is set at the beginning of the processing, which in turn may vary over time as a consequence of sublimation. By constraining the thermal conductivity such that it may span the interval 10<sup>-3</sup>-10<sup>-2</sup> W K<sup>-1</sup> m<sup>-1</sup>, consistent with experimental data obtained for cometary-like dust layers<sup>53</sup>, thermal inertia does not exceed 30 J K<sup>-1</sup> m<sup>-2</sup> s<sup>-0.5</sup> (thermal inertia units or TIU), which is exactly in the middle of the range of approximately 10 to 50 TIU derived by MIRO on a global scale<sup>17</sup>. In Supplementary Fig. 6, we report the theoretical temperature profiles obtained in Hapi for decreasing values of near-infrared emissivity, interpreted as due to increasing contribution of the surface roughness.

In the smooth terrain within the Imhotep region we have considered two cases:
(1) one case in which the uppermost cm-thick surface layer is made up of pure
dust, thickness increasing with increasing sublimation from the underlying layers,
which provides overall higher temperatures; and (2) another case where the surface
layer is made up by a mixture of dust and water ice with a dust/ice ratio close to 1,
which yields lower surface temperatures. In both cases, the dust/ice ratio below the
surface layer is equal to 10 in mass at the beginning of the processing, and varies with time in the layers where sublimation and recondensation processes take place.

The physical temperature of 127 K, as derived by the submillimetre channel of MIRO in Imhotep in August 2014, results from an average brightness temperature of 114 K (figure 4 in Gulkis et al.¹³) vertically integrated with a weighting function extending a little above and below 1 cm depth.

## Data availability

The VIRTIS calibrated data are publicly available through the ESA's Planetary Science Archive (PSA) website (https://archives.esac.esa.int/psa/) and NASA's Planetary Data System (https://pds.nasa.gov/) in accordance with the schedule established by the Rosetta project. Other data that support the plots within this paper and other findings of this study are available from the corresponding author upon reasonable request. Readers are welcome to comment on the online version of the paper.

## Code availability

The code used to retrieve surface temperature values from VIRTIS-M infrared data is a direct implementation of a published method⁷. The code used to derive synthetic thermal images of the nucleus of comet 67P is a direct implementation of published models¹⁸,⁴⁴⁻⁴⁶,⁵¹. The code used to derive theoretical temperature profiles for specific locations of the nucleus of comet 67P is a direct implementation of a published model²⁸,²⁹.

Received: 13 April 2018; Accepted: 1 March 2019;

Published online: 22 April 2019

## References

1. Groussin, O. et al. Surface temperature of the nucleus of comet 9P/Tempel 1. *Icarus* **187**, 16–25 (2007).

2. Squyres, S. W., McKay, C. P. & Reynolds, R. T. Temperatures within comet nuclei. *J. Geophys. Res. Solid Earth* **90**, 12831–12392 (1985).

3. Emerich, C. et al. Temperature and size of the nucleus of comet P/Halley deduced from IKS infrared Vega 1 measurements. *Astron. Astrophys.* **187**, 839–842 (1987).

4. Soderblom, L. A. et al. Short-wavelength infrared (1.3–2.6 µm) observations of the nucleus of comet 19P/Borrelle. *Icarus* **167**, 100–112 (2004).

5. Groussin, O. et al. The temperature, thermal inertia, roughness and color of the nuclei of comets 103P/Hartley 2 and 9P/Tempel 1. *Icarus* **222**, 580–594 (2013).

6. Coradini, A. et al. VIRTIS: an imaging spectrometer for the Rosetta mission. *Space Sci. Rev.* **128**, 529–559 (2007).

7. Tosi, F. et al. Thermal measurements of dark and bright surface features on Vesta as derived from Dawn/VIR. *Icarus* **240**, 36–57 (2014).

8. Coradini, A. et al. The surface composition and temperature of asteroid 21 Lutetia as observed by Rosetta/VIRTIS. *Science* **334**, 492–494 (2011).

9. Keihm, S. et al. Interpretation of combined infrared, submillimeter, and millimeter thermal flux data obtained during the Rosetta fly-by of Asteroid (21) Lutetia. *Icarus* **221**, 395–404 (2012).

10. Ciarniello, C. et al. Photometric properties of comet 67P/Churyumov-Gerasimenko from VIRTIS-M onboard Rosetta. *Astron. Astrophys.* **583**, A31 (2015).

11. Longobardo, A. et al. Photometric behaviour of 67P/Churyumov-Gerasimenko and analysis of its pre-perihelion diurnal variations. *Mon. Not. R. Astron. Soc.* **469** (Suppl. 2), S346–S356 (2017).

12. Gulkis, S. et al. MIRO: Microwave Instrument for Rosetta Orbiter. *Space Sci. Rev.* **128**, 561–597 (2007).

13. Gulkis, S. et al. Subsurface properties and early activity of comet 67P/Churyumov-Gerasimenko. *Science* **347**, aaa0709 (2015).

14. Schloerb, F. P. et al. MIRO observations of subsurface temperatures of the nucleus of 67P/Churyumov-Gerasimenko. *Astron. Astrophys.* **583**, A29 (2015).

15. Jorda, L. et al. The global shape, density and rotation of comet 67P/Churyumov-Gerasimenko from preperihelion Rosetta/OSIRIS observations. *Icarus* **277**, 257–278 (2016).

16. Thomas, N. et al. The morphological diversity of comet 67P/Churyumov-Gerasimenko. *Science* **347**, aaa0440 (2015).

17. Davidsson, B. J. R. & Rickman, H. Surface roughness and three-dimensional heat conduction in thermophysical models. *Icarus* **243**, 58–77 (2014).

18. Keller, H. U. et al. Insolation, erosion, and morphology of comet 67P/Churyumov-Gerasimenko. *Astron. Astrophys.* **583**, A34 (2015).

19. Marshall, D. et al. Thermal inertia and roughness of the nucleus of comet 67P/Churyumov-Gerasimenko from MIRO and VIRTIS observations. *Astron. Astrophys.* **616**, A122 (2018).

20. Lagerros, J. S. V. Thermal physics of asteroids. III. Irregular shapes and albedo variegations. *Astron. Astrophys.* **325**, 1226–1236 (1997).

21. Davidsson, B. J. R. et al. Interpretation of thermal emission. I. The effect of roughness for spatially resolved atmosphereless bodies. *Icarus* **252**, 1–21 (2015).

22. Preusker, F. et al. Shape model, reference system definition, and cartographic mapping standards for comet 67P/Churyumov-Gerasimenko-stereo-photogrammetric analysis of Rosetta/OSIRIS image data. *Astron. Astrophys.* **583**, A33 (2015).

23. Spohn, T. et al. Thermal and mechanical properties of the near-surface layers of comet 67P/Churyumov-Gerasimenko. *Science* **349**, aab0464 (2015).

24. Fornasier, S. et al. Spectrophotometric properties of the nucleus of comet 67P/Churyumov-Gerasimenko from the OSIRIS instrument onboard the Rosetta spacecraft. *Astron. Astrophys.* **583**, A30 (2015).

25. Vincent, J.-B. et al. Large heterogeneities in comet 67P as revealed by active pits from sinkhole collapse. *Nature* **523**, 63–66 (2015).

26. Sierks, H. et al. On the nucleus structure and activity of comet 67P/Churyumov-Gerasimenko. *Science* **347**, aaa1044 (2015).

27. Keller, H. U. et al. OSIRIS — The Scientific camera system onboard Rosetta. *Space Sci. Rev.* **128**, 433–506 (2007).

28. De Sanctis, M. C., Lasue, J. & Capria, M. T. Seasonal effects on comet nuclei evolution: activity, internal structure, and dust mantle formation. *Astron. J.* **140**, 1–13 (2010).

29. Capria, M. T., Marchi, S., De Sanctis, M. C., Coradini, A. & Ammannito, E. The activity of main belt comets. *Astron. Astrophys.* **537**, A71 (2012).

30. Capaccioni, F. et al. The organic-rich surface of comet 67P/Churyumov-Gerasimenko as seen by VIRTIS/Rosetta. *Science* **347**, aaa0628 (2015).

31. Quirico, E. et al. Refractory and semi-volatile organics at the surface of comet 67P/Churyumov-Gerasimenko: insights from the VIRTIS/Rosetta imaging spectrometer. *Icarus* **272**, 32–47 (2016).

32. Rousseau, B. et al. Laboratory simulations of the Vis-NIR spectra of comet 67P using sub-µm sized cosmochemical analogues. *Icarus* **306**, 306–318 (2018).

33. Molaro, J. & Byrne, S. Rates of temperature change of airless landscapes and implications for thermal stress weathering. *J. Geophys. Res. Planets* **117**, E10011 (2012).

34. Delbo, M. et al. Thermal fatigue as the origin of regolith on small asteroids. *Nature* **508**, 233–236 (2014).

35. De Sanctis, M. C. et al. The cycle of water ice on the surface of 67P/Churyumov-Gerasimenko. *Nature* **525**, 500–503 (2015).

36. El-Maarry, M. R. et al. Regional surface morphology of comet 67P/Churyumov-Gerasimenko from Rosetta/OSIRIS images: the southern hemisphere. *Astron. Astrophys.* **593**, A110 (2016).

37. Filacchione, G. et al. Exposed water ice on the nucleus of comet 67P/Churyumov-Gerasimenko. *Nature* **529**, 368–372 (2016).

38. Filacchione, G. et al. The global surface composition of 67P/CG nucleus by Rosetta/VIRTIS. (I) Prelanding mission phase. *Icarus* **274**, 334–349 (2016).

39. Ciarniello, M. et al. The global surface composition of 67P/Churyumov-Gerasimenko nucleus by Rosetta/VIRTIS. II) Diurnal and seasonal variability. *Mon. Not. R. Astron. Soc.* **462** (Suppl. 1), S443–S458 (2016).

40. Raponi, A. et al. The temporal evolution of exposed water ice-rich areas on the surface of 67P/Churyumov-Gerasimenko: spectral analysis. *Mon. Not. R. Astron. Soc.* **462** (Suppl. 1), S476–S490 (2017).

41. Barucci, M. A. et al. Detection of exposed H₂O ice on the nucleus of comet 67P/Churyumov-Gerasimenko as observed by Rosetta OSIRIS and VIRTIS instruments. *Astron. Astrophys.* **595**, A102 (2016).

42. Fornasier, S. et al. Rosetta's comet 67P/Churyumov-Gerasimenko sheds its dusty mantle to reveal its icy nature. *Science* **354**, 1566–1570 (2016).

43. Lethuillier, A. et al. Analysis of observations of the Imhotep region of 67P/C-G performed by MIRO/Rosetta in 2014 and 2016 and derived constraints on the close subsurface properties. In *American Astronomical Society, DPS Meeting 49* 415.03 (2017).

44. Skorov, Yu. & Blum, J. Dust release and tensile strength of the non-volatile layer of cometary nuclei. *Icarus* **221**, 1–11 (2012).

45. Skorov, Yu. V., Lieshout, R. V., Blum, J. & Keller, H. U. Activity of comets: gas transport in the near-surface porous layers of a cometary nucleus. *Icarus* **212**, 867–876 (2011).

46. Gundlach, B., Skorov, Yu. V. & Blum, J. Outgassing of icy bodies in the Solar System—I. The sublimation of hexagonal water ice through dust layers. *Icarus* **213**, 710–719 (2011).

47. Zinzi, A. et al. The SSDC contribution to the improvement of knowledge by means of 3D data projections of minor bodies. *Adv. Space Res.* **62**, 2306–2316 (2018).

48. Lagerros, J. S. V. Thermal physics of asteroids. IV. Thermal infrared beaming. *Astron. Astrophys.* **332**, 1123–1132 (1998).

49. Davidsson, B. J. R., Gutiérrez, P. J. & Rickman, H. Physical properties of morphological units on comet 9P/Tempel 1 derived from near-IR Deep Impact spectra. *Icarus* **201**, 335–357 (2009).

50. Giese, B. & Kührt, E. Theoretical interpretation of infrared measurements at Deimos in the framework of crater radiation. *Icarus* **88**, 372–379 (1990).

51. Capria, M. T. et al. How pristine is the interior of the comet 67P/Churyumov-Gerasimenko? *Mon. Not. R. Astron. Soc.* **469** (Suppl. 2), S685–S694 (2017).
52. Acton, C. H. Ancillary data services of NASA's navigation and ancillary information facility. *Planet. Space Sci.* **44**, 65–70 (1996).

53. Grün, E. et al. Development of a dust mantle on the surface of an insulated ice-dust mixture: results from the KOSI-9 experiment. *J. Geophys. Res.* **98**, 15091–15104 (1993).

## Acknowledgements

The authors thank the following institutions and agencies, which supported this work: Italian Space Agency (ASI-Italy), Centre National d'Etudes Spatiales (CNES-France), Deutsches Zentrum für Luft- und Raumfahrt (DLR-Germany), National Aeronautic and Space Administration (NASA-USA). VIRTIS was built by a consortium from Italy, France and Germany, under the scientific responsibility of IAPS, Istituto di Astrofisica e Planetologia Spaziali of INAF, Rome, Italy, which led also the scientific operations. The VIRTIS instrument development for ESA has been funded and managed by ASI, with contributions from Observatoire de Meudon financed by CNES and from DLR. The VIRTIS instrument industrial prime contractor was former Officine Galileo, now Leonardo SpA in Campi Bisenzio, Florence, Italy. The authors thank the Rosetta Liaison Scientists, the Rosetta Science Ground Segment and the Rosetta Mission Operations Centre for their support in planning the VIRTIS observations. We also thank the MIRO science and MIRO archiving teams for making MIRO data available to us before their public release. This research has made use of NASA's Astrophysics Data System. D.K. acknowledges DFG-grant KA 3757/2-1. This work is dedicated to Angioletta Coradini (1946–2011), conceiver of the VIRTIS instrument, and to Sergio Fonti (1945–2018), co-author of this Article and active contributor in the development of VIRTIS. The first author dedicates this work also to Luca Malagutti (1965–2017), who was a brilliant researcher at the University of Milan.

## Author contributions

F.T. carried out the surface temperature retrieval from VIRTIS-M infrared data, derived geometric information for those data, and led the analysis of VIRTIS-

measured temperature data, writing major sections of the main text and Methods. F.C. is the principal investigator of the VIRTIS instrument; he designed the overall study and wrote part of the main text. S.M., M.T.C. and M.F. carried out the thermophysical modelling and wrote part of the Methods. M.C. led the spectrophotometric analysis and derived the single scattering albedo values. G.F. was responsible for the VIRTIS calibration pipeline and flagged the saturated spectral pixels. M.H. is the principal investigator of the MIRO instrument; he granted MIRO data calibrated with the latest responsivity function, and provided key information for their proper interpretation. F.C., G.F., S.E., D.B.-M. and C.L. planned VIRTIS observations. S.E. and G.A. are respectively the French and German group leaders within the VIRTIS Science Team. The other authors are all VIRTIS co-investigators and associates who participated in the study and/or reviewed the manuscript, providing edits, comments and suggestions that led to substantial improvement of the paper.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** is available for this paper at https://doi.org/10.1038/s41550-019-0740-0.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Correspondence and requests for materials** should be addressed to F.T.

**Journal peer review information:** *Nature Astronomy* thanks Ben Rozitis and the other anonymous reviewer(s) for their contribution to the peer review of this work.

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© The Author(s), under exclusive licence to Springer Nature Limited 2019