# Macro and micro structures of pebble-made cometary nuclei reconciled by seasonal evolution

Mauro Ciarniello, Marco Fulle, Andrea Raponi, Gianrico Filacchione, Fabrizio Capaccioni, Alessandra Rotundi, Giovanna Rinaldi, Michelangelo Formisano, Gianfranco Magni, Federico Tosi, et al.

► To cite this version:

Mauro Ciarniello, Marco Fulle, Andrea Raponi, Gianrico Filacchione, Fabrizio Capaccioni, et al. Macro and micro structures of pebble-made cometary nuclei reconciled by seasonal evolution. Nature Astronomy, 2022, 6, pp.546-553. 10.1038/s41550-022-01625-y. insu-03705364

HAL Id: insu-03705364

https://insu.hal.science/insu-03705364v1

Submitted on 3 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
<figure><img src="image_2.png" /></figure>

# Macro and micro structures of pebble-made cometary nuclei reconciled by seasonal evolution

Mauro Ciarniello<sup>1, ∗</sup>, Marco Fulle<sup>2</sup>, Andrea Raponi<sup>1</sup>, Gianrico Filacchione<sup>1</sup>, Fabrizio Capaccioni<sup>1</sup>, Alessandra Rotundi<sup>1,3</sup>, Giovanna Rinaldi<sup>1</sup>, Michelangelo Formisano<sup>1</sup>, Gianfranco Magni<sup>1</sup>, Federico Tosi<sup>1</sup>, Maria Cristina De Sanctis<sup>1</sup>, Maria Teresa Capria<sup>1</sup>, Andrea Longobardo<sup>1</sup>, Pierre Beck<sup>4</sup>, Sonia Fornasier<sup>5,6</sup>, David Kappel<sup>7,8</sup>, Vito Mennella<sup>9</sup>, Stefano Mottola<sup>7</sup>, Batiste Rousseau<sup>1</sup> and Gabriele Arnold<sup>7</sup>

Comets evolve due to sublimation of ices embedded inside porous dust, triggering dust emission (that is, erosion) followed by mass loss, mass redistribution and surface modifications. Surface changes were revealed by the Deep Impact and Stardust NExT missions for comet 9P/Tempel 1 (ref. <fnref n="1" />), and a full inventory of the processes modifying cometary nuclei was provided by Rosetta while it escorted comet 67P/Churyumov-Gerasimenko for approximately two years²⁻⁴. Such observations also showed puzzling water-ice-rich spots that stood out as patches optically brighter and spectrally bluer than the average cometary surface⁵⁻⁹. These are up to tens of metres large and indicate macroscopic compositional dishomogeneities apparently in contrast with the structural homogeneity above centimetre scales of pebble-made nuclei¹⁰. Here we show that the occurrence of blue patches determines the seasonal variability of the nucleus colour⁴,¹¹,¹² and gives insight into the internal structure of comets. We define a new model that links the centimetre-sized pebbles composing the nucleus¹⁰ and driving cometary activity¹³,¹⁴ to metre-sized water-ice-enriched blocks embedded in a drier matrix. The emergence of blue patches is due to the matrix erosion driven by CO₂-ice sublimation that exposes the water-ice-enriched blocks, which in turn are eroded by water-ice sublimation when exposed to sunlight. Our model explains the observed seasonal evolution of the nucleus and reconciles the available data at micro (sub-centimetre) and macro (metre) scales.

The Mapping channel of the Visible, Infrared and Thermal Imaging Spectrometer (VIRTIS-M)¹⁵ observed the surface of 67P/ Churyumov-Gerasimenko (hereafter 67P) throughout the whole Rosetta mission with its visible (VIS) channel (0.2–1 µm). Previous studies⁷,¹¹, supported by the availability of infrared observations from August 2014 to April 2015, correlated the VIS colour to superficial water-ice abundance; lower (bluer) spectral slopes indicated larger amounts of water ice. Here we exploit the whole VIRTIS-M/VIS channel dataset of nucleus observations (August 2014 to September 2016) to systematically characterize the surface water-ice seasonal variations. We produced maps of spectral slope in the 0.55–0.8 µm interval to monitor the spatial and temporal variability of 67P surface colour along the orbit, as reported in Fig. 1 (see Supplementary

Sections 1 and 2 for spectral-slope computation, maps production and data calibration refinement).

The August 2014 spectral-slope map reflects the dichotomy characterizing the 67P surface colour, with the Hapi morphological region, at northern latitudes, being distinctively bluer than the rest of the nucleus as a result of a larger content of water ice¹¹,¹⁶. Due to the inclination of the comet's rotation axis (obliquity of 52.4°) the summer of northern and southern regions occurs in opposite orbital phases¹⁷. Equatorial regions are illuminated during most of the orbit and have been systematically imaged by VIRTIS-M. These regions are characterized by a redder spectrum with respect to Hapi in August 2014, and get globally bluer approaching the Sun. The minimum spectral slope is reached around perihelion, while the original colours are gradually restored outbound from perihelion.

We selected five regions of interest (ROIs, Fig. 1) that are illuminated and observed throughout most of the mission, to monitor in more detail the spectral-slope behaviour from August 2014 to September 2016, with a temporal sampling of 30 days (Fig. 2). All the ROIs display a similar evolution of the surface colour, with maximum blueing reached around perihelion (1.23 au, August 2015) and a spectral-slope reduction of the order of 15–25% compared to August 2014 (~3.5 au, inbound). In September 2016 (~3.8 au, outbound), the spectral slope returned to approximately the same values as those measured on August 2014. These facts indicate that the progressive enrichment of water ice at smaller heliocentric distances is a global process occurring on the whole comet surface, with the exception of the Hapi region, which is not illuminated at perihelion.

Although the colour evolution implies a global water-ice enrichment approaching perihelion, water ice is characterized by a patchy distribution across the surface of comet 67P (the blue patches (BPs))⁶⁻⁸ consistent with the dust ejection rates at perihelion¹⁸ measured by GIADA (Grain Impact Analyser and Dust Accumulator)¹⁹. Cometary dust ejection is principally driven by sublimation of H₂O and CO₂ ices. Water-driven erosion ejects sub-centimetre dust¹³, whereas the dust volume distribution in 67P's coma is dominated by decimetre-sized chunks¹⁸,²⁰, which can be emitted only by CO₂-ice sublimation, occurring at greater depth with respect to water ice given the deeper sublimation front¹⁴,²¹. Activity models¹³ of water-ice sublimation and vapour diffusion inside the centimetre-sized pebbles

<footnote n="1">¹Istituto di Astrofisica e Planetologia Spaziali (IAPS), Istituto Nazionale di Astrofisica (INAF), Rome, Italy. ²Osservatorio Astronomico, INAF, Trieste, Italy.</footnote>

<footnote n="3">³Dipartimento di Scienze e Tecnologie, Università degli Studi di Napoli Parthenope, Naples, Italy. ⁴Université Grenoble Alpes, CNRS, IPAG, Grenoble, France. ⁵Institut Universitaire de France (IUF), Paris, France. ⁶LESIA, Observatoire de Paris, Université PSL, CNRS, Université de Paris, Sorbonne Université, Meudon, France. ⁷Institute of Planetary Research, German Aerospace Center (DLR), Berlin, Germany. ⁸Institute of Physics and Astronomy, University of Potsdam, Potsdam, Germany. ⁹Osservatorio Astronomico di Capodimonte, INAF, Naples, Italy. e-mail: mauro.ciarniello@inaf.it</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1 |** VIS spectral-slope ($S_{\text{vis}}$) maps of comet 67P/Churyumov-Gerasimenko's surface at different orbital phases. Evolution of the comet's colour along the inbound orbit, close to perihelion and during the outbound orbit (see Supplementary Table 1 for additional information). The surface colour gets progressively bluer inbound, when the comet approaches the Sun (perihelion at 1.23 au on 13 August 2015), followed by a reddening during the outbound phase. Solid coloured rectangles border five ROIs (Supplementary Table 2) for which the temporal evolution of the spectral slope is monitored in Fig. 2.</figcaption></figure>

that make up a cometary nucleus¹⁰ are consistent with the observed sub-centimetre dust fluxes¹⁸ if only up to ~6% (ref. ¹³ and Methods) of the southern hemisphere undergoes erosion as a consequence of dust ejection induced by water-ice sublimation. The water-driven activity is regulated by the refractory to water-ice mass ratio $\delta$. For values of $\delta$ lower than 5, sources in the southern hemisphere at perihelion are rich enough in water to sustain sub-centimetre dust ejection through water-ice sublimation, exposing the ice inside the pebble. Conversely, for larger values of $\delta$, water-driven erosion dies out and the surface material is dehydrated until ejected by CO₂ sublimation in the form of decimetre-sized chunks. This exposes underlying pebbles, activating their water-driven erosion or dehydration according to the corresponding $\delta$. The comet's surface can be described as a geographical mixture of bluer water-ice-rich areas (undergoing water-driven erosion) and redder water-ice-poor (dehydrating) areas, with the former being the BPs.

The observed seasonal colour evolution is then produced by the temporal evolution of the nucleus surface fraction covered by BPs. High-resolution observations of BPs on 67P's surface indicate a volumetric ice fraction as high as 50% (refs. ⁹,²²), which proves that these areas, having a low enough refractory to water-ice mass ratio, effectively sustain water-driven erosion. There is a direct correlation between surface temperature (and by consequence heliocentric distance) and the maximum value of $\delta$ still compatible

with water-driven erosion (Methods and Supplementary Fig. 2). Assuming that all the BPs are characterized by a similar characteristic water-ice content (same $\delta$), the condition for water-driven erosion in the southern hemisphere at perihelion ($\delta<5$, ref. ¹³) ensures their water-driven erosion also at larger heliocentric distances. Areas with $\delta \gg 5$ are eroded by water-ice sublimation only at large heliocentric distances, and are progressively dehydrated approaching the Sun if not refreshed by CO₂-driven erosion.

By modelling (Supplementary Section 3) the observed spectral-slope temporal evolution, we estimate the corresponding variation of the 67P surface fraction covered by BPs. We perform this computation for the ROI 1 (located in Imhotep), the region characterized by the best observation conditions and smallest uncertainties in the spectral-slope values (Fig. 3). We can compare our evaluation with the independent estimate of the fraction of areas undergoing water-driven erosion required to match the dust ejection rate ($Q_d$) measured in March 2015 by the GIADA experiment¹⁸, corresponding to approximately 0.8–1.5% of Imhotep (Supplementary Section 4). In the same period, modelling the VIRTIS results provides a BP fraction of between 0.3% and 0.9% in the ROI 1. Assuming that these values as representative of the entire Imhotep region, the estimated BP fraction is compatible with the fraction of surface ejecting sub-centimetre dust derived from the dust ejection rate. Moreover, the observed surface blueing
<figure><img src="image_4.png" /><figcaption>**Fig. 2 |** Spectral-slope ($S_{\text{VIS}}$) temporal evolution from August 2014 to September 2016, for different areas on 67P's surface. Each curve corresponds to one of the five different ROIs shown in Fig. 1, as indicated by numbers and colours. From bottom to top, the curves are offset by +5%/(100 nm). Each point and the corresponding error bars represent the spectral-slope average value and standard deviation, respectively, computed over a 30-day bin from all the observations falling in the corresponding ROI, after sigma-clipping. The Medium Term Planning (MTP) phases of the Rosetta mission and the perihelion date (dashed line) are indicated.</figcaption></figure>

is compatible with a surface fraction getting covered in BPs for a maximum of 1.4–2.2% around perihelion and decreasing at larger heliocentric distances (Fig. 3). If we extrapolate the result obtained

at perihelion for Imhotep to the southern hemisphere, the derived maximum BP areal fraction is compatible with the upper limit of areal fraction characterized by water-driven activity as estimated from observed dust ejection rates¹³,¹⁸ (~6%).

We also took into consideration an alternative scenario (Supplementary Section 5), in which the observed surface colour variability is due to the seasonal variation of the thickness of a ‘semi-transparent’ dehydrated crust on top of an active subsurface ice-rich layer (a paradigm common to most thermophysical models of comets¹⁷,²³). This would affect the intensity of the subsurface ice spectral signal reaching the observer and by consequence the surface colour. However, this approach implies that the nucleus colour should depend on its surface temperature. This result is inconsistent with VIRTIS data, which indicate that Hapi and Imhotep have similar colours at different surface temperatures and vice versa. Also, it is not consistent with a pebble-made nucleus, where the crust must be at least as thick as a pebble¹³, thereby being totally opaque to the radiation in the VIRTIS spectral range.

Explaining the colour seasonal cycle of the 67P nucleus as the effect of the occurrence of BPs implies that the dominating CO₂ erosion exposes new water-ice-enriched blocks (WEBs) present in the subsurface of the nucleus, which then become BPs at the surface (Fig. 4). The exposure of new BPs depends on the CO₂-driven erosion rate and the volumetric abundance of the WEBs in the nucleus. After exposure, BPs are eroded with a timescale primarily dependent on the water-driven erosion rate and on the typical size of the WEBs (Methods). To provide the observed blueing, the WEBs must be large enough to survive the time required to build up the observed BP areal fraction close to perihelion. Thus, the temporal evolution of a cometary nucleus colour, and in particular the observed blueing of 67P at perihelion, is due to the balance between two competing mechanisms: CO₂-driven erosion, which exposes new BPs, and water-driven erosion occurring in the BPs, which removes them (Fig. 4). This balance, in turn, depends on the surface insolation conditions, which modulate seasonally the H₂O- and CO₂-driven erosion rates. In Methods, we report a detailed modeling of the temporal evolution of the BP areal fraction on the 67P nucleus undergoing H₂O- and CO₂-driven erosions. We apply our simulations to the ROI 1 case (see Fig. 3 and complete results in Supplementary Figs. 6 and 7), assuming an average erosion in the region of ≲0.5 m during one orbit²⁴. Depending on the volumetric

<figure><img src="image_5.png" /><figcaption>**Fig. 3 |** Temporal evolution of the BP areal fraction in the ROI 1. Diamonds indicate the BP areal fraction estimated from the $S_{\text{VIS}}$ temporal profile. Error bars at each position follow from the corresponding uncertainty in $S_{\text{VIS}}$. The coloured lines are best-fit models (from Supplementary Fig. 6b,c) of the BP temporal evolution as the effect of CO₂-driven and water-driven erosion (Methods) for two different combinations of WEB size *L* and WEB volumetric fraction $f_{\text{WEB}}$ (red, $L=0.5$ m and $f_{\text{WEB}}=0.095$; blue, $L=1.0$ m and $f_{\text{WEB}}=0.055$). After December 2015, the observed BP temporal evolution is also affected by fallout of dehydrated material ejected at perihelion (Supplementary Section 6). For each model, the WEB size and the WEB volumetric abundance in the nucleus are indicated. The perihelion passage is marked by the vertical dashed line.</figcaption></figure>
<figure><img src="image_6.png" /></figure>

**Fig. 4** | 67P/Churyumov-Gerasimenko surface blueing at perihelion is due to the progressive exposure to sunlight of subsurface WEBs. The evolution of the surface and subsurface is sketched here as a function of heliocentric distance (displayed on the right). The comet nucleus is composed of two types of pebbles made of refractories and CO₂ ice plus: (1) a high content of H₂O ice (forming the WEBs); (2) a low content of H₂O ice (H₂O-ice-poor pebbles, compatible with the material observed during the occasional exposure of CO₂ ice by ref. ³⁹), comprising most of the nucleus. CO₂ ice survives under the CO₂ condensation front at depths >0.1 m (ref. ¹⁴). Its sublimation rate increases approaching perihelion, eroding the surface by chunk ejection and exposing the subsurface WEBs. After exposure, WEBs lose their CO₂ and become BPs; water-ice sublimation erodes the BPs while removing dust from their surface and preventing the formation of a dry crust¹³. The BPs survive as such until the complete sublimation of their water-ice fraction, and long enough to induce the observed surface blueing. Dehydration overcomes erosion in the H₂O-ice-poor pebbles¹³, whereas the formation of a crust is prevented by the ongoing CO₂-driven erosion. BPs and dehydrating surface eject the same water flux, because both have some water ice inside and are assumed to be at the same temperature¹³.

abundance of the WEBs in the nucleus, we find that these must be sub-metre-sized to match the observed surface colour variability up to the maximum blueing. After perihelion, the surface in the ROI 1 gets progressively redder while receding from the Sun; according to our simulations this is due to a progressive reduction of the WEBs exposure rate along the outbound orbit as the CO₂ activity decreases and due to the cumulative backfall of dehydrated material ejected at perihelion that blankets the comet surface and covers the BPs²³,²⁵,²⁶ (Supplementary Section 6).

To further constrain the typical WEB size, we computed the expected areal fraction of BPs in the active southern hemisphere at perihelion (Methods) as a function of volumetric abundance in the nucleus and size of the WEBs as obtained from the ROI 1 analysis. Assuming a total erosion from CO₂-driven chunk emission of ~4–10 m (refs. ¹³,²⁰), we find that only WEBs with a dominant size in the ~0.5–1 m range and volumetric abundance in the nucleus from ~9.5% to 5.5% are compatible with the ~6% active areas (that is, BPs) upper limit inferred from the measured dust ejection rates¹³,¹⁸. Sub-metre WEBs are compatible with CONSENT (Comet Nucleus Sounding Experiment by Radiowave Transmission) radar observations, which indicate a mostly homogeneous 67P nucleus interior structure down to a few-metres scale²⁷ and with OSIRIS Narrow Angle Camera (NAC) observations (Methods). Nonetheless, BPs up

to tens of metres have also been observed on 67P's surface⁷,⁸, sampling the large-size tail of the WEB size distribution.

As a large part of the nucleus surface displays a colour variability similar to that of the ROI 1 (Figs. 1 and 2), we infer that the mechanism described above acts on the comet surface at a global scale, with the exception of most of the northern deposits (for example, Hapi) not illuminated at perihelion and blanketed by backfall material from the southern hemisphere during the northern winter. The low spectral slope of the Hapi deposits in August 2014 is indeed compatible with an origin of this material from the southern hemisphere active at perihelion²⁵ (Supplementary Section 7).

A uniform distribution of WEBs (Fig. 4) explains the observed surface colour temporal evolution on 67P and the observed dust ejection rates. Furthermore, it reconciles the macroscopic (up to tens of metres) compositional dishomogeneities observed in comets as BPs, with the structural (although not compositional; Fig. 4) homogeneity at small (centimetre) scale characterizing a pebble-made nucleus, and with the processes at the base of water-driven activity at microscopic (sub-pebble) scales. WEBs are likely primordial, because orbital water-ice recondensation within the nucleus can operate only within very shallow depths⁴,²⁸. WEBs challenge the present comet formation models, requiring primordial macroscopic localization of water-ice-rich pebbles during the
formation of comets. Protoplanetary disk observations²⁹ suggest the presence of clouds of ice-rich dust segregated from ice-poor grains. Crystalline water ice detected in protoplanetary disks³⁰ is consistent with water-ice recondensation³¹ possibly forming pebbles enriched in water ice with respect to the average in molecular clouds³², that is, stickier than water-ice-poor pebbles. These facts, coupled with the complex evolution of the protoplanetary disk³³, may favour the formation of WEBs during the gravitational collapse triggered by streaming instabilities that formed comets and planets¹⁰.

## Methods

**Condition for water-driven activity.** According to the cometary activity model by ref. ¹³ for a nucleus made of pebbles, dust emission upon water-ice sublimation (water-driven erosion) can be described by means of two characteristic rates represented by the dust erosion rate $E(T)$ and the dehydration rate $D(T)$, both dependent on the surface temperature $T$. The former represents the thickness of the surface layer eroded into dust per unit time, whereas the latter is the thickness of the surface layer that gets dehydrated by water-ice sublimation per unit time. The condition for a given area to be active is that $E(T) \ge D(T)$. If not, in the surface upper layer dehydration occurs more rapidly than erosion, leading to a dry and insulating crust that prevents further activity if not refreshed by CO₂-driven erosion¹¹. From ref. ¹³ the dehydration rate can be expressed as

$$D(T) = (\delta + 1) \frac{Q_g(T)}{\rho_n}, \qquad (1)$$

where $Q_g(T)$ is the gas emission rate and $\rho_n = 538 \text{ kg m}^{-3}$ is the nucleus density¹⁴. Unlike the dehydration rate, the erosion rate does not depend on the dust to water-ice mass ratio $\delta$, and we refer the reader to ref. ¹² for the expression of $E(T)$. Positing that $D(T) = E(T)$, we obtain from equation (1) the maximum value of $\delta$ ($\delta_{MAX}$) for which the surface can be active, which is a function of the surface temperature $T$:

$$\delta_{MAX}(T) = \rho_n \frac{E(T)}{Q_g(T)} - 1. \qquad (2)$$

By taking into account the insolation history of the ROI 1 (Supplementary Fig. 1), we estimate the temporal evolution of the corresponding maximum surface temperature as provided by the model of ref. ¹³ and obtain $\delta_{MAX}(T)$ at different times from equation (2) (Supplementary Fig. 2). In general, $\delta_{MAX}$ decreases with increasing temperatures, and for the ROI 1 it attains its minimum around perihelion. With the same approach it can be shown that in August 2014 the Hapi region is active if $\delta < 10^5$, whereas $\delta < 5$ is required for the southern hemisphere at perihelion¹³.

**Modelling of BP temporal variability.** Explaining the colour variability of a cometary nucleus as the effect of the progressive exposure of new BPs from CO₂-driven activity implies a distribution of discrete subsurface WEBs across the whole comet. The WEBs are exposed by the removal of overlying chunks lifted by CO₂ sublimation¹⁴. The exposed parts of the WEBs appear as BPs on the surface, and their abundance at a given time results from the combined effect of CO₂-driven activity, which determines the BP exposure rate, and water-ice sublimation, which erodes the BPs. We assume that all the WEBs have $\delta < 5$, which ensures that they can sustain water-driven activity up to perihelion. Both CO₂-driven activity and water-ice sublimation are modulated by the surface temperature, which varies with time, and determine the temporal evolution of the BP abundance on the surface. We note that also the sublimation of CO and O₂ might contribute to the surface erosion by chunk emission. Given the deeper sublimation front of these species with respect to CO₂, the ejected chunks would be typically larger than 0.2 m. However, the volume erosion rate for chunks with sizes >0.2 m is much lower than that for 0.1–0.2 m sizes (refs. ¹⁸,³²), so we can neglect this contribution to the total erosion.

We assume that at the onset of the CO₂ activity, when the nucleus exhibits the maximum redness, the surface of the comet is characterized by a dehydrated layer covering an internal distribution of WEBs having a constant density with depth. For simplicity, we assume that the emission of a chunk from CO₂ sublimation in the dehydrated layer directly samples depths that are already representative of the internal density of WEBs. Depending on the numerical density of WEBs and their size, it is then possible to define the fractional area of WEBs in the subsurface, $f_{WEB}$, which represents the probability of exposing a BP that results from chunk emission. $f_{WEB}$ also represents the WEB volumetric abundance in the nucleus. Let's now define $BP(t)$ as the fraction of the surface covered by BPs at the time $t$. At a given time, it is possible to express the variation of the areal fraction covered in BPs as

$$\mathrm{dBP}(t) = \mathrm{dBP}_{\mathrm{exp}}(t) - \mathrm{dBP}_{\mathrm{rem}}(t), \qquad (3)$$

where $\mathrm{dBP}_{\mathrm{exp}}(t)$ is the increase in areal fraction of exposed BPs in the interval $[t, t+dt]$ due to CO₂-driven chunk emission and $-\mathrm{dBP}_{\mathrm{rem}}(t)$ is the decrease of areal fraction due to BP removal produced by H₂O-driven and CO₂-driven erosion

taking place over the BPs. Assuming that $f_{CO_2}(t)$ is the areal erosion rate from CO₂ activity (that is, the areal fraction experiencing chunk emission), we can express $\mathrm{dBP}_{\mathrm{exp}}(t)$ as

$$\mathrm{dBP}_{\mathrm{exp}}(t) = f_{\mathrm{CO}_2}(t)(f_{\mathrm{WEB}} - \mathrm{BP}(t)) \mathrm{d}t, \qquad (4)$$

where the term $-f_{CO_2}(t) \mathrm{BP}(t)$ accounts for chunk emission taking place over already-existing BPs, which does not increase the areal fraction of exposed WEBs.

The term $\mathrm{dBP}_{\mathrm{rem}}(t)$ accounts primarily for the sublimation of the exposed BPs. A BP sublimates completely after a time $\Delta(t)$ since the start of exposure, which depends on the H₂O-driven erosion rate and the WEB size $L$. In addition, CO₂ chunks emission may in principle take place also over exposed BP, which contributes to reducing their survival time. Given this, $\mathrm{dBP}_{\mathrm{rem}}(t)$ can be expressed as

$$\mathrm{dBP}_{\mathrm{rem}}(t) = \mathrm{dBP}_{\mathrm{exp}}(t') = f_{\mathrm{CO}_2}(t')(f_{\mathrm{WEB}} - \mathrm{BP}(t')) \mathrm{d}t', \qquad (5)$$

where $\mathrm{dBP}_{\mathrm{exp}}(t')$ is the areal fractional of BPs exposed in the interval $[t', t' + dt']$, with $t' = t - \Delta(t)$ and with $\Delta(t)$ being the time required to remove at time $t$ a BP first exposed at time $t'$, due to erosion from water-ice sublimation and CO₂ activity. The emission of chunks from CO₂ activity, which decreases the depth $L$ of a given BP by the typical chunk size $L_c$, reduces also its sublimation time. Taking this into account, we define an effective depth $L_{\mathrm{eff}}(t)$ of the BPs as follows:

$$L_{\mathrm{eff}}(t) = L - L_c \times \int_{t'}^{t} f_{\mathrm{CO}_2}(t^*) \mathrm{d}t^*, \qquad (6)$$

which represents the residual average depth of BPs exposed at time $t'$, and where the term $L_c \times \int_{t'}^{t} f_{\mathrm{CO}_2}(t^*) \mathrm{d}t^*$ is the average depth eroded by CO₂-driven activity. Finally, $\Delta(t)$ is given by the implicit relation

$$L_{\mathrm{eff}}(t) = \int_{t-\Delta(t)}^{t} E(t^*) \mathrm{d}t^*, \qquad (7)$$

where $E(t)$ is the erosion rate from water-ice activity in a BP.

For given values of $L$, $f_{WEB}$ and $L_c$, the set of equations (3)–(7) can be solved numerically, once $f_{CO_2}(t)$ and $E(t)$ are provided. For the ROI 1, $E(t)$ can be determined from ref. ¹³ by taking into account the insolation history of the region as reported in Supplementary Fig. 1. We compute the average erosion rate with a 10-day running box and obtain the result of Supplementary Fig. 3.

The quantity $f_{CO_2}(t)$ depends on the vertical erosion rate from CO₂-driven activity. To our knowledge, this has been modelled so far only for the southern hemisphere at perihelion¹⁴. Given this, for the purpose of this work we test different formulations of $f_{CO_2}(t)$ computed starting from the H₂O-driven erosion rate $E(t)$, with the parametric expression $f_{CO_2}(t) = K \frac{E(t)^\alpha}{\mathrm{MAX}(E(t))^\alpha}$, where $K$ is the areal erosion rate at maximum activity. This appears to be a reasonable approximation, as $f_{CO_2}(t)$, like $E(t)$, is expected to be driven by the insolation conditions, reaching a maximum around perihelion and decreasing monotonically with heliocentric distance. We report typical behaviours of the temporal evolution of the BP areal fraction $BP(t)$ for different values of $\alpha$ in Supplementary Fig. 4. It can be noted that the $\alpha=1$ case, after an initial transient phase, reaches a stationary regime. The transient phase ends with a local maximum well before perihelion, which is determined by the BP areal fraction build-up due to the exposure of the first WEBs in November 2014 up to February–early March 2015, when the first BPs were completely eroded. This build-up behaviour can be noted also for the $\alpha=2$ case, whereas for $\alpha=3$ it is less evident as the end of the transient regime occurs at low BP areal fractions. As the end of the transient phase depends on the sublimation time of the WEBs, it occurs earlier for smaller $L$ (see the set of simulations for different $L$ values shown in Supplementary Fig. 6, where the end of the transient phase appears to shift in time with the WEB size).

The occurrence of the stationary regime can be explained with the following reasoning. In the stationary condition we have $\mathrm{dBP}_{\mathrm{exp}}(t) - \mathrm{dBP}_{\mathrm{rem}}(t) = 0$. Assuming that the areal coverage of BPs is far from saturation (that is, the areal fraction of exposed BPs is much smaller than the areal fraction $f_{WEB}$ in the interior of the comet, $BP(t) \ll f_{WEB}$), the stationary condition implies that $f_{CO_2}(t) f_{WEB}(t) = f_{CO_2}(t') f_{WEB}(t')$. In the time interval $[t, t+dt]$, BPs exposed in the interval $[t', t' + dt]$ finally sublimate. Then, in a regime of BP erosion dominated by water ice, we have

$$L = \int_{t'}^{t} E(t^*) \mathrm{d}t^* = \int_{t' + dt'}^{t + dt} E(t^*) \mathrm{d}t^*, \qquad (8)$$

which implies that

$$\int_{t'}^{t' + dt'} E(t^*) \mathrm{d}t^* = \int_{t}^{t + dt} E(t^*) \mathrm{d}t^*, \qquad (9)$$
and finally $E(t') dt' = E(t) dt$. Given this, the stationary condition becomes
$$f_{CO_2}(t)f_{WEB}(t) dt = \frac{f_{CO_2}(t')}{E(t')} f_{WEB}(t) dt. \text{ If } f_{CO_2}(t) \propto E(t)^a, \text{ we obtain } \frac{E(t')^a}{E(t')} = \frac{E(t')^a}{E(t')}$$
which is trivially valid for every $t$ and $t'$ if $\alpha=1$.

Conversely, if $\alpha \neq 1$, the stationary points of BP(t) must satisfy
$$E(t)^{\alpha-1} = E(t')^{\alpha-1}.$$

Position of the maximum of BP(t). If $E(t)$ is characterized by a local maximum at $t_{MAX}$, as in the investigated case, and $\alpha \neq 1$, the condition $E(t)^{\alpha-1} = E(t')^{\alpha-1}$ can be met for a certain combination of $t$ and $t'$, with $t > t_{MAX}$ and $t' < t_{MAX}$. Assuming that in the neighborhood of $t_{MAX}$ the BP sublimation time is approximately constant, that is, $\Delta(t) \approx \Delta(t_{MAX})$, we have $dt' \approx dt$. Given this, and expanding both $E(t)^{\alpha-1}$ and $E(t')^{\alpha-1}$ around $t_{MAX}$, we obtain

$$E(t)^{\alpha-1} = E(t_{MAX})^{\alpha-1} + \frac{1}{2} \frac{d^2(E(t_{MAX})^{\alpha-1})}{dt^2} \times (t - t_{MAX})^2 \quad (10)$$

$$E(t')^{\alpha-1} = E(t_{MAX})^{\alpha-1} + \frac{1}{2} \frac{d^2(E(t_{MAX})^{\alpha-1})}{dt^2} \times (t' - t_{MAX})^2 \quad (11)$$

as $\frac{d(E(t_{MAX})^{\alpha-1})}{dt} = 0$ and where it can be shown that

$$\frac{d^2(E(t_{MAX})^{\alpha-1})}{dL^2} = (\alpha - 1)E(t_{MAX})^{\alpha-2} \frac{d^2(E(t_{MAX})^{\alpha-1})}{dL^2}.$$

$$(t - t_{MAX})^2 = (t' - t_{MAX})^2, \quad (12)$$

which is the case

$$t = t_{MAX} + \Delta/2 \approx t_{MAX} + \frac{L}{2E(t_{MAX})}. \quad (13)$$

Interestingly, as $t_{MAX}$ is the position of a local maximum, we have $\frac{d^2E(t_{MAX})}{dt^2} < 0$ that, for $\alpha > 1$, provides the stationary point $t = t_{MAX} + \Delta/2$ as the position of a local maximum of BP(t), whereas for $\alpha < 1$ it identifies a local minimum. Given this, for the study of 67P's BP temporal evolution, only cases with $\alpha > 1$ are of interest. The discussion above indicates that the maximum of BP(t) is delayed with respect to the maximum of CO₂ activity by an amount that depends on the sublimation time of the BPs and that is related to the size L of the WEBs. Thus, the quantification of such delay from the BP(t) curve inferred from VIRTIS data would in principle provide an estimate of the WEB size. Assuming that $\alpha=2$, we illustrate in Supplementary Fig. 5 the effect on the final shape of BP(t) for different sizes of the WEB by integrating numerically our set of equations. It is worth noting that if the condition BP(t) $\ll f_{WEB}$ is not met, that is, the areal fraction of BP's is close to saturation, the considerations made above are not necessarily valid, and the position of the maximum of BP(t) can be shifted at times $\le t_{MAX} + \Delta/2$, preventing the determination of the WEB size from the position of the BP(t) maximum (see also below).

Matching the BP's temporal evolution for the ROI 1. The set of equations described above is applied to model the temporal evolution of the BP abundance in the ROI 1. For simplicity, we assume that $\alpha=2$ (as shown above, only cases with $\alpha > 1$ are of interest for 67P's BP temporal evolution) and investigate four different WEB sizes, $L=0.1$ m, 0.5 m, 1 m and 2 m. As described above, along with $L$ and $\alpha$ the model depends on the three additional parameters $L_c$ and $f_{WEB}$. The chunk size $L_c$, determined by the depth of the CO₂ sublimation front, may vary with the heliocentric distance. However, here we adopt a constant value of $L_c=0.1$ m, which is comparable to values for the southern hemisphere at perihelion¹⁴. In this respect, a time-dependent formulation of $L_c$ does not appear as strictly necessary, given the modelling uncertainty in $f_{CO_2}(t)$. Moreover, the adopted value of $L_c$ is representative for the perihelion phase, where most of the CO₂-driven erosion takes place. The maximum areal erosion rate from CO₂ activity K can be fixed by taking into account the total vertical erosion measured in the ROI 1. For the Imhotep region, the upper limit on the total erosion is 0.5 m, as determined by analysis of OSIRIS images²⁴. According to this, we set our simulations to provide an erosion of ~0.4 m through CO₂ erosion, which means requiring that $L_c \times \int_{Aug 2014}^{Sep 2016} f_{CO_2}(t) dt \approx 0.4$ m. We verified that such a value of CO₂-driven erosion is compatible with a total erosion of ~0.5 m when accounting for the water-ice-driven erosion in the BPs. For $\alpha=2$, this condition fixes $K=0.0262$ per day in our model, and by consequence the BP(t) curve can be matched by varying the free parameter $f_{WEB}$. Such a value of K results in a maximum CO₂-driven erosion rate at perihelion of $K \times L_c = 0.262$ cm per day, which is much smaller than the maximum water-driven erosion rate of 5.5 cm per day in the BPs. In Supplementary Fig. 6 we report fits of the BP(t) curve for the ROI 1 from VIRTIS data with the results of our simulations.

As post-perihelion fallout of dehydrated material emitted around perihelion²⁶ may cover exposed BPs, stop water-driven erosion and reduce the resulting abundance of BPs on the surface (Supplementary Section 6), the fit is performed by matching only the BP(t) curve up to one month after perihelion, where the maximum blueing is measured. A good match of the BP curve can be achieved for $\alpha=2$ with basically all the explored values of L. In particular,

the simulations indicate that the smaller the size of the WEBs, the larger the $f_{WEB}$ required to match the BP(t) curve. This is explained by the fact that smaller WEBs sublimate faster, and then a larger abundance in the subsurface is required to sustain the increase of exposed BPs at perihelion. All the simulations peak close to perihelion, whereas the maximum on the BP(t) curve for the ROI 1 appears shifted to one month after perihelion. However, given the error bars on the BP(t) curve, and the coarse temporal resolution of the measured BP(t) curve (30 days), we cannot consider such a one-month shift as diagnostic of the WEB size. In this respect, we note that the spectral-slope curve of the ROI 3 shows the maximum blueing slightly before perihelion, whereas ROI 4 displays a broader spectral-slope minimum centred at perihelion (Fig. 2). This is indicative of the uncertainty associated with the position at which maximum blueing occurs. It is worth noting that also the simulations performed with a large WEB size ($L=2$ m, Supplementary Fig. 6d), for which a noticeable shift of the BP maximum towards larger times after perihelion can be expected, still peak close to perihelion. This is explained by the fact that the required values of $f_{WEB}$ are close to the maximum value of the BP curve, thus indicating that the BP abundance on the surface was approaching a saturation regime, in which case equation (13) cannot be applied directly. Conversely, simulations run with $L=0.1$ m (Supplementary Fig. 6a) provide $f_{WEB} \approx 0.4$, which is much larger than the maximum value registered in BP(t), which is of the order of ~2%. As $f_{WEB}$ increases with decreasing WEBs size, it is then possible to define a limit size of the WEBs for which $f_{WEB}=1$ is required to match the BP temporal evolution in the ROI 1. This would represent a lower limit of the typical size of the WEBs for comet 67P. In this regime, where $L_c > L$, the emission of a chunk removes the BPs present in its area and simultaneously exposes new WEBs below. In this case, $f_{WEB}$ represents the average areal fraction of BPs at the scale of the area exposed by the chunk emission $\approx L_c^2$. To treat this case we simplify our model slightly by neglecting the term BP(t) $\ll f_{WEB}$ in equations (4) and (5) and by setting $L_{e1}(t)=L$, because in this regime exposed BPs can be eroded by water-ice sublimation only. In Supplementary Fig. 7 we show the result of a simulation with $\alpha=2$ and the smallest values of $L$ ($L_{min}$) for which the observed trend is still matched. We obtain that the simulated temporal evolution of BPs on the surface is still compatible with the observed one if $L_{min}=0.032$ m. Smaller values of L would not be able to match the temporal evolution in Imhotep, and as a consequence $L_{min}$ provides an estimate of the lower limit of the WEB size. We also stress that in the ideal limit of $f_{WEB}=1$, WEBs would no longer be separate units with a finite size. Thus, the provided value of $L_{min}$ represents the asymptotic WEB size in this regime and underestimates the physical WEB size lower limit.

An upper limit for the WEB size can also be provided. In May 2015, when VIRTIS already measured a detectable blueing in the ROI 1, OSIRIS²⁶ NAC observations with a nominal spatial resolution of 2.4 m per pixel show blueing over an extended area in the central part of the Imhotep region¹⁴. Considering a reference point spread function full width at half maximum for the OSIRIS NAC of ~1.8 pixels²⁷,²⁸, we can conclude that BPs are almost evenly distributed over a spatial scale of $2.4 \times 1.8 \text{ m} = 4.32 \text{ m}$ in this region. This implies that each surface element corresponding to the equivalent OSIRIS pixel, having an area $A_{eq}=(4.32 \text{ m})^2=18.7 \text{ m}^2$, approximately hosts the BP areal fraction observed by VIRTIS. This also indicates that roughly each surface element hosts at least one WEB partially or fully exposed, depending on the size. As $f_{WEB}$ represents the areal fraction of WEBs at depth, and assuming a uniform distribution of the WEBs in the subsurface, with a given size L we have on average one WEB for each surface element of area $A=L^2/f_{WEB}$. To be compatible with the May 2015 OSIRIS NAC observations of the central part of the Imhotep region, A must be smaller than $A_{eq}$, thus granting that a large part of the equivalent pixels are populated by WEBs. For $L=1$ m and the corresponding $f_{WEB}=5.5\%$, as obtained by our simulations, we have $A \approx 18-19 \text{ m}^2$. Given that $A_{eq}$ represents a lower limit of the real projected area of the equivalent pixel (nadir observation geometry) and that OSIRIS spectral-slope images of this region are not fully uniform and show residual pixel-by-pixel variability, we can conclude that a WEB size of $L=1$ m is compatible with OSIRIS observations. Larger WEB sizes require smaller values of $f_{WEB}$, resulting in a larger A. For example, with $L=2$ m and $f_{WEB}=3.5\%$ we obtain $A \approx 114 \text{ m}^2 > A_{eq}$. This indicates that $L \approx 1$ m is also a good estimate of the maximum WEB size. After May 2015, OSIRIS NAC observations have a nominal spatial resolution worse than 2 m per pixel roughly until mid-December 2015, when the BP abundance is characterized by values similar to those in May 2015. In later periods, the nominal spatial resolution of OSIRIS NAC observations is typically better than 2 m per pixel. However, the BP areal fraction rapidly reduces to ~0 in early February 2016. This prevents us from further constraining the upper limit of the WEB size by the approach described above, given the limited occurrence of BPs when the OSIRIS NAC observes at better spatial resolution.

Southern hemisphere erosion and expected blueing at perihelion. In the period from 24 July 2015 to 15 September 2015, the active part of the southern hemisphere of comet 67P underwent an average erosion by chunks ejection of ~4–10 m (refs. 13,20). Assuming such a value for erosion, we can estimate the expected BP fraction on the southern hemisphere surface by means of our modelling (see above). For simplicity, we assume constant isolation conditions in the southern hemisphere for this period, which is polar summer¹⁷, with a
water-driven erosion of E=13 cm per day (ref. ¹³). In stationary conditions, the increase of BP in the time interval [t, t+dt] can be expressed as

$$ \frac{dB\text{P}_{\text{South}}}{dt} = f_{\text{CO}_2} (f_{\text{WEB}} - \text{BP}_{\text{South}}) - \frac{\text{BP}_{\text{South}}}{\Delta}, \quad (14) $$

where BP<sub>South</sub> is the surface fraction of the active southern hemisphere covered in BP<sub>South</sub> in the stationary regime, Δ is the time of erosion of a BP due to the combined effect of H₂O, and CO₂ erosion and the term BP<sub>South</sub>/Δ is the areal fraction of BP removed per unit time. The latter stems from the fact that all the BP<sub>South</sub> on the surface at a given time t-Δ are completely eroded by the time t, and thus the BP removal rate is BP<sub>South</sub>/Δ. Assuming a constant CO₂ areal fraction erosion rate f<sub>CO₂</sub> = K, we have an average vertical CO₂ erosion rate given by K×L<sub>c</sub>, where L<sub>c</sub> is the emitted chunk size. K is defined from the total erosion E<sub>tot</sub> in chunks of size L<sub>c</sub> in the 53-day period P from 24 July 2015 to 15 September 2015, which is expressed by E<sub>tot</sub> = K×P×L<sub>c</sub>, and gives K = E<sub>tot</sub>/P×L<sub>c</sub>. Finally, Δ depends on the WEB size L and the H₂O and CO₂ erosion rates as Δ = L/(K×L<sub>c</sub>+E). In the stationary regime we have dB<sub>P<sub>South</sub></sub>/dt=0, and equation (14) can be solved straightforwardly, providing

$$ \text{BP}_{\text{South}} = \frac{K}{1/\Delta + K} f_{\text{WEB}}. \quad (15) $$

We have defined combinations of L and f<sub>WEB</sub> that provide a BP temporal evolution compatible with the one observed in Imhotep. Our results indicate that WEB sizes smaller than ~1 m are required, with f<sub>WEB</sub> decreasing with increasing L. Assuming the same WEB physical properties across the whole nucleus, we can use the values of L and f<sub>WEB</sub> inferred from the ROI 1 analysis to compute the expected surface fraction of BP in the active southern hemisphere from equation (15). We investigated four possible combinations accounting for the two extreme erosion cases, E<sub>tot</sub>=4 m (ref. ²⁰) and E<sub>tot</sub>=10 m (ref. ¹³), and two slightly different values of the chunk size, L<sub>c</sub>=0.1 m (as for the ROI 1 models) and L<sub>c</sub>=0.13 m (more representative of the southern hemisphere at perihelion¹⁰). In Supplementary Fig. 8, BP<sub>South</sub> is reported against L and is compared with the upper limit of the water-driven eroded areas (δ<5). This upper limit was originally estimated to ~4% (ref. ¹³), but we adopt a more conservative value of 6%, including uncertainties affecting the dust and water loss rates¹³, BP<sub>South</sub> decreases with increasing L. These results indicate that only L≥0.5 m provides values of BP<sub>South</sub> compatible with the upper limit on water-active areas. As a matter of fact, this constrains the WEB size to ~0.5–1 m, as WEBs larger than 1 m are not compatible with the upper limit on the WEB size, from the ROI 1 analysis. For L≈0.5–1 m we obtain BP<sub>South</sub>≈4–6%, which is relatively close to the value measured in the ROI 1 in the same period (~2%) and in qualitative agreement with the observed southern hemisphere colour at perihelion being similar to that of the rest of the nucleus¹. If only half of the southern hemisphere were active at perihelion¹³, and active and non-active areas were mixed at a few-metre scale, then the final amount of BP<sub>South</sub> would be half of BP<sub>South</sub>, thus even closer to the values in Imhotep. However, there is no clear observational evidence in this respect. WEBs much smaller than 0.5 m, along with exceedance of the upper limit of water-driven eroded areas, would result in a southern hemisphere colour at perihelion much bluer than the rest of the surface. The possibility of describing with a single set of WEB properties the occurrence of BP<sub>South</sub> (that is, water-driven eroded areas) across different regions of the comet (Imhotep and the southern hemisphere) supports the scenario of a uniform distribution of WEBs in the 67P nucleus, as per our proposed cometary nucleus model.

## Data availability

The VIRTIS calibrated data are publicly available through the European Space Agency's Planetary Science Archive website (https://archives.esac.esa.int/psa/). Source data are provided with this paper.

## Code availability

The computer code used to produce VIS spectral-slope maps of 67P is a direct implementation of a published method¹¹. The computer codes used to perform spectral modelling and simulations of the BP fraction temporal evolution are direct implementations of the models described in the present paper.

Received: 18 May 2021; Accepted: 2 February 2022;
Published online: 17 March 2022

## References

1. Veverka, J. et al. Return to comet Tempel 1: overview of stardust-NEXT results. *Icarus* **222**, 424–435 (2013).

2. Groussin, O. et al. Temporal morphological changes in the Imhotep region of comet 67P/Churyumov-Gerasimenko. *Astron. Astrophys.* **583**, A36 (2015).

3. El-Maarry, M. R. et al. Surface changes on comet 67P/Churyumov-Gerasimenko suggest a more active past. *Science* **355**, 1392–1395 (2017).

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

5. Sunshine, J. M. et al. Exposed water ice deposits on the surface of comet 9P/ Tempel 1. *Science* **311**, 1453–1455 (2006).

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

7. 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**, S476–S490 (2016).

8. 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).

9. Oklay, N. et al. Long-term survival of surface water ice on comet 67P. *Mon. Not. R. Astron. Soc.* **469**, S582–S597 (2017).

10. Blum, J. et al. Evidence for the formation of comet 67P/ Churyumov-Gerasimenko through gravitational collapse of a bound clump of pebbles. *Mon. Not. R. Astron. Soc.* **469**, S755–S773 (2017).

11. 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**, S443–S458 (2016).

12. Filaccione, G. et al. An orbital water-ice cycle on comet 67P from colour changes. *Nature* **578**, 49–52 (2020).

13. Fulle, M. et al. How comets work: nucleus erosion versus dehydration. *Mon. Not. R. Astron. Soc.* **493**, 4039–4044 (2020).

14. Gundlach, B., Fulle, M. & Blum, J. On the activity of comets: understanding the gas and dust emission from comet 67P/Churyumov-Gerasimenko's south-pole region during perihelion. *Mon. Not. R. Astron. Soc.* **493**, 3690–3715 (2020).

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

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

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

18. Fulle, M. et al. Evolution of the dust size distribution of comet 67P/ Churyumov-Gerasimenko from 2.2 au to perihelion. *Astrophys. J.* **821**, 19 (2016).

19. Colangeli, L. et al. The Grain Impact Analyser and Dust Accumulator (GIADA) experiment for the Rosetta mission: design, performances and first results. *Space Sci. Rev.* **128**, 803–821 (2007).

20. Fulle, M. et al. The refractory-to-ice mass ratio in comets. *Mon. Not. R. Astron. Soc.* **482**, 3326–3340 (2019).

21. Capria, M. T. et al. How pristine is the interior of the comet 67P/ Churyumov-Gerasimenko? *Mon. Not. R. Astron. Soc.* **469**, S685–S694 (2017).

22. O'Rourke, L. et al. The Philae lander reveals low-strength primitive ice inside cometary boulders. *Nature* **586**, 697–701 (2020).

23. Hu, X. et al. Seasonal erosion and restoration of the dust cover on comet 67P/Churyumov-Gerasimenko as observed by OSIRIS onboard Rosetta. *Astron. Astrophys.* **604**, A114 (2017).

24. Cambianica, P. et al. Long-term measurements of the erosion and accretion of dust deposits on comet 67P/Churyumov-Gerasimenko with the OSIRIS instrument. *Mon. Not. R. Astron. Soc.* **504**, 2895–2910 (2021).

25. Keller, H. U. et al. Seasonal mass transfer on the nucleus of comet 67P/ Churyumov-Gerasimenko. *Mon. Not. R. Astron. Soc.* **469**, S357–S371 (2017).

26. Bertini, I. et al. The backscattering ratio of comet 67P/Churyumov-Gerasimenko dust coma as seen by OSIRIS onboard Rosetta. *Mon. Not. R. Astron. Soc.* **482**, 2924–2933 (2019).

27. Ciarletti, V. et al. CONSERT constrains the internal structure of 67P at a few metres size scale. *Mon. Not. R. Astron. Soc.* **469**, S805–S817 (2017).

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

29. Terada, H. et al. Detection of water ice in edge-on protoplanetary disks: HK Tauri B and HV Tauri C. *Astrophys. J.* **667**, 303–307 (2007).

30. Min, M. et al. The abundance and thermal history of water ice in the disk surrounding HD142527 from the DIGIT Herschel Key Program. *Astron. Astrophys.* **593**, A11 (2016).

31. Schoonenberg, D. & Ormel, C. W. Planetesimal formation near the snowline: in or out? *Astron. Astrophys.* **602**, A21 (2017).

32. Cambianica, P. et al. Time evolution of dust deposits in the Hapi region of comet 67P/Churyumov-Gerasimenko. *Astron. Astrophys.* **636**, A91 (2020).

33. Brownlee, D. et al. Comet 81P/Wild 2 under a microscope. *Science* **314**, 1711–1716 (2006).

34. Pätzold, M. et al. The nucleus of comet 67P/Churyumov-Gerasimenko – Part I: the global view – nucleus mass, mass-loss, porosity, and implications. *Mon. Not. R. Astron. Soc.* **483**, 2337–2346 (2019).

35. Ott, T. et al. Dust mass distribution around comet 67P/Churyumov-Gerasimenko determined via parallax measurements using Rosetta's OSIRIS cameras. *Mon. Not. R. Astron. Soc.* **469**, S276–S284 (2017).
36. Keller, H. U. et al. OSIRIS – the scientific camera system onboard Rosetta. *Space Sci. Rev.* **128**, 433–506 (2007).

37. Bertini, I. et al. Search for satellites near comet 67P/Churyumov-Gerasimenko using Rosetta/OSIRIS images. *Astron. Astrophys.* **583**, A19 (2015).

38. Magrin, S. et al. Pre-hibernation performances of the OSIRIS cameras onboard the Rosetta spacecraft. *Astron. Astrophys.* **574**, A123 (2015).

39. Filacchione, G. et al. Seasonal exposure of carbon dioxide ice on the nucleus of comet 67P/Churyumov-Gerasimenko. *Science* **354**, 1563–1566 (2016).

## Acknowledgements

We thank the Italian Space Agency (ASI, Italy; ASI-INAF agreements I/032/05/0 and I/024/12/0), Centre National d'Etudes Spatiales (CNES, France) and Deutsches Zentrum für Luft und Raumfahrt (DLR, Germany) for supporting this work. VIRTIS was built by a consortium from Italy, France and Germany, under the scientific responsibility of Istituto di Astrofisica e Planetologia Spaziali (IAPS) of INAF, Rome, which also led the scientific operations. The VIRTIS instrument development for ESA has been funded and managed by ASI (Italy), with contributions from Observatoire de Meudon (France) financed by CNES and from DLR (Germany). The VIRTIS instrument industrial prime contractor was former Officine Galileo, now the Leonardo Company, in Campi Bisenzio, Florence, Italy. Part of this research was supported by the ESA Express Procurement (EXPRO) request for proposal for IPL-PSS/JD/190.2016. D.K. acknowledges DFG grant no. KA 3757/2-1. This work was supported by the International Space Science Institute (ISSI) through the ISSI International Team 'Characterization of cometary activity of 67P/ Churyumov-Gerasimenko comet'. This research has made use of NASA's Astrophysics Data System.

## Author contributions

M.C. wrote the manuscript, performed data analysis, modelling and interpretation, and contributed to VIRTIS data calibration. M. Fulle contributed to model conceptualization, data interpretation and manuscript drafting. A. Raponi, G.F., F.C., A. Rotundi and G.R. contributed to data interpretation. A. Raponi supported spectral modelling and contributed to VIRTIS data calibration. G.F. provided VIRTIS data calibration. F.C. managed the VIRTIS experiment. A. Rotundi managed the GIADA experiment. F.T. provided geometric files for VIRTIS nucleus observations. All authors, including M. Formisano, G.M., M.C.D.S., M.T.C., A.L., P.B., S.F., D.K., V.M., S.M., B.R. and G.A., contributed to the discussion of the results and helped with the manuscript preparation.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41550-022-01625-y.

**Correspondence and requests for materials** should be addressed to Mauro Ciarniello.

**Peer review information** *Nature Astronomy* thanks the anonymous reviewers for their contribution to the peer review of this work.

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

**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 2022