# Analysis of night-side dust activity on comet 67P observed by VIRTIS-M: a new method to constrain the thermal inertia on the surface

G. Rinaldi, M. Formisano, D. Kappel, F. Capaccioni, D. Bockelee-Morvan, Y.-C. Cheng, J.-B. Vincent, P. Deshapriya, G. Arnold, M. Capria, et al.

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

## ► To cite this version:

G. Rinaldi, M. Formisano, D. Kappel, F. Capaccioni, D. Bockelee-Morvan, et al.. Analysis of night-side dust activity on comet 67P observed by VIRTIS-M: a new method to constrain the thermal inertia on the surface. Astronomy and Astrophysics - A&A, 2019, 630, pp.A21. 10.1051/0004-6361/201834907. hal-02396688

HAL Id: hal-02396688

https://hal.science/hal-02396688

Submitted on 10 Nov 2021

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>
Rosetta mission full comet phase results

# Analysis of night-side dust activity on comet 67P observed by VIRTIS-M: a new method to constrain the thermal inertia on the surface

G. Rinaldi<fnref n="1" />, M. Formisano<fnref n="1" />, D. Kappel<fnref n="2" />,<fnref n="3" />, F. Capaccioni<fnref n="1" />, D. Bockelée-Morvan<fnref n="4" />, Y.-C. Cheng<fnref n="4" />, J.-B. Vincent<fnref n="3" />, P. Deshapriya<fnref n="4" />, G. Arnold<fnref n="3" />, M. T. Capria<fnref n="1" />, M. Ciarniello<fnref n="1" />, E. D'Aversa<fnref n="1" />, M. C. De Sanctis<fnref n="1" />, L. Doose<fnref n="5" />, S. Erard<fnref n="4" />, C. Federico<fnref n="1" />, G. Filacchione<fnref n="1" />, U. Fink<fnref n="5" />, C. Leyrat<fnref n="4" />, A. Longobardo<fnref n="1" />,<fnref n="6" />, G. Magni<fnref n="1" />, A. Migliorini<fnref n="1" />, S. Mottola<fnref n="3" />, G. Naletto<fnref n="7" />, A. Raponi<fnref n="1" />, F. Taylor<fnref n="8" />, F. Tosi<fnref n="1" />, G. P. Tozzi<fnref n="9" />, and M. Salatti<fnref n="10" />

Received 17 December 2018 / Accepted 28 March 2019

## ABSTRACT

On 2015 July 18, near perihelion at a heliocentric distance of 1.28 au, the Visible InfraRed Thermal Imaging Spectrometer (VIRTIS-M) on board the Rosetta spacecraft had the opportunity of observing dust activity in the inner coma with a view of the night side (shadowed side) of comet 67P/Churyumov-Gerasimenko. At the time of the measurements we present here, we observe a dust plume that originates on the far side of the nucleus. We are able to identify the approximate location of its source at the boundary between the Hapi and Anuket regions, and we find that it has been in darkness for some hours before the observation. Assuming that this time span is equal to the conductive time scale, we obtain a thermal inertia in the range 25–36 W K⁻<fnref n="1" /> m⁻<fnref n="2" /> s⁻<fnref n="1" />/<fnref n="2" />. These thermal inertia values can be used to verify with a 3D finite-element method (FEM) numerical code whether the surface and subsurface temperatures agree with the values found in the literature. We explored three different configurations: (1) a layer of water ice mixed with dust beneath a dust mantle of 5 mm with thermal inertia of 36 J m⁻<fnref n="2" /> K⁻<fnref n="1" /> s⁻⁰.<fnref n="5" />; (2) the same structure, but with thermal inertia of 100 J m⁻<fnref n="2" /> K⁻<fnref n="1" /> s⁻⁰.<fnref n="5" />; (3) an ice-dust mixture that is directly exposed. Of these three configurations, the first seems to be the most reasonable, both for the low thermal inertia and for the agreement with the surface and subsurface temperatures that have been found for the comet 67P/Churyumov-Gerasimenko. The spectral properties of the plume show that the visible dust color ranged from 16 ± 4.8%/100 nm to 13 ± 2.6%/100 nm, indicating that this plume has no detectable color gradient. The morphology of the plume can be classified as a narrow jet that has an estimated total ejected mass of between 6 and 19 tons when we assume size distribution indices between −2.5 and −3.

**Key words.** comets: general – comets: individual: 67P/Churyumov-Gerasimenko – infrared: planetary systems

## 1. Introduction

Dust activity on the night side (dark side) of comet 67P/ Churyumov-Gerasimenko (hereafter 67P) has been observed by the Optical, Spectroscopic, and Infrared Remote Imaging System (OSIRIS) on board the Rosetta mission (Knollenberg et al. 2016; Shi et al. 2016), and similar activity has been reported on other comets by previous space missions (Farnham et al. 2007; Feaga et al. 2007; Belton et al. 2008; A'Hearn et al. 2011). Some of these activities, which often occur deep in the night side, are classified as outbursts that lead to the formation of short-lived dust jets and plumes. Feaga et al. (2007) suggested that the dust activities observed near sunrise might be triggered

by sublimation of water ice and of more volatile ices (super volatiles, such as CO and CO₂) that have refrozen and accumulated on the surface during the night. Belton et al. (2008) showed that outbursts occur both in darkness or daylight, suggesting that they are not controlled by direct sunlight but rather by processes in the interior. This conclusion is strengthened by astronomical observations that show no dependence of the occurrence rate of outbursts with heliocentric distance and only a small increase in their occurrence rate post-perihelion relative to pre-perihelion. The post-perihelion activity peak, observed on 67P and other comets, suggests that the thermal environment in the immediate vicinity of the surface, in particular, the lag due to the thermal inertia of the nucleus surface layers, may play a role in

<footnote n="1">¹ Istituto di Astrofisica e Planetologia Spaziali, Istituto Nazionale di Astrofisica, Rome, Italy
e-mail: giovanna.rinaldi@inaf.it</footnote>

<footnote n="2">² Institute of Physics and Astronomy, University of Potsdam, Potsdam, Germany</footnote>

<footnote n="3">³ Institute of Planetary Research, DLR, 12489 Berlin, Germany</footnote>

<footnote n="4">⁴ LESIA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Université, Université Paris Diderot, Sorbonne Paris Cité, 5 place Jules Janssen, 92195 Meudon, France</footnote>

<footnote n="5">⁵ Lunar and Planetary Laboratory, University of Arizona, Tucson, AZ 85721, USA</footnote>

<footnote n="6">⁶ Università degli Studi di Napoli Parthenope, Dip. di Scienze e Tecnologie, CDN IC4, 80143 Naples, Italy</footnote>

<footnote n="7">⁷ Department of Physics and Astronomy, University of Padova, Via Marzolo 8, 35131 Padova, Italy</footnote>

<footnote n="8">⁸ Department of Physics, Oxford University, Oxford, UK</footnote>

<footnote n="9">⁹ INAF – Osservatorio Astrofisico of Arcetri, Firenze, Italy</footnote>

<footnote n="10">¹⁰ ASI, Via del Politecnico snc, 00133, Rome, Italy</footnote>
triggering the outburst (Gulkis et al. 2015; Spohn et al. 2015). On 2015 March 12, Knollenberg et al. (2016) observed with the OSIRIS Wide Angle Camera (WAC) one prominent event before dawn, when 67P was approaching perihelion. This event was probably due to a mini-outburst coming from the Imhotep region. To explain the mechanism that caused the outburst event, these authors suggested the build-up of high-pressure gas conditions under the devolatilized surface layers. The pressure would increase until it overcame the tensile strength of the surface, at which point the gas pocket erupted, releasing the gas and accelerating the surrounding non-volatile material. This would leave a small pit or depression that exposed fresh volatile material that continued to sublimate. The validity of such a mechanism, which is a strong candidate for explaining the observation we report here, depends on the properties of the surface and subsurface, such as its layering and the thermal inertia profile.

Thermal inertia of comet surfaces has been analyzed for several previous comet missions. Groussin et al. (2013) derived values lower than 250 W K⁻<fnref n="1" /> m⁻<fnref n="2" /> s⁻<fnref n="1" />/<fnref n="2" /> (hereafter TIU, Thermal Inertia Units) for comet 103P/Hartley 2 and lower than 45 TIU for comet 9P/Tempel 1 from the temperature rise on the morning terminator. Similar values have been found by A'Hearn et al. (2011) and Davidsson et al. (2013) for the same comets. For 67P, thermal inertia has been derived from data obtained by the Microwave Instrument on the Rosetta orbiter (MIRO; Gulkis et al. 2007) and from data by the Multipurpose Sensors for Surface and Sub-Surface Science (MUPUS) instrument package at the Philae landing site Abydos (Spohn et al. 2007). Spohn et al. (2007) derived a low value of thermal inertia, in the range between 10 to 80 TIU (Gulkis et al. 2015; Spohn et al. 2015; Schloerb et al. 2015).

The main purpose of the present work is to analyze a specific night-side dust activity observed by the Visual Infrared and Thermal Imaging Spectrometer (VIRTIS; Coradini et al. 2007) and so to obtain an independent constraint on the surface thermal inertia of comet 67P in the correspondence of the active emitting area. The mapping channel VIRTIS-M (VM hereafter) observed this night-side dust emission at a heliocentric distance of 1.28 au, one month before the perihelion passage.

Despite the large number of dust emission observations by Rosetta and other comet missions and hypotheses such as those described above, the mechanism that forms the dust is still unknown. Gas flow must be involved to accelerate the dust away from the surface, and we need a thermal analysis of the conditions to understand their role in generating the night-side activity. Our analysis of the illumination conditions of the likely emitting region provides a new method for inferring the thermal inertia of that region. For this we used the COMSOL Multiphysics 5.3a package<fnref n="1" />, in particular, the module "Surface-to-Surface Radiation" in combination with MESHLAB<fnref n="2" />, applied to the shape model SHAP5 (Jorda et al. 2016), taking into account shadowing effects. It is crucial to calculate the illumination conditions to compute the solar energy input to the surface. We investigated whether it is possible for an event to be sustained by a thermal lag in the subsurface and applied a 3D thermophysical model (using the COMSOL Multiphysics module "Heat Transfer in Solids") in order to evaluate the surface and subsurface temperature of the emitting region. Because we are only interested in the thermal history of the likely emitting region and also because this requires much computational time, the model was applied only to a specific part of the neck region of the nucleus, using the 3D

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Configuration of the nucleus of the comet and jet with respect to the VIRTIS-M slit (horizontal red line) and its direction (red arrow). The VIRTIS-M data cubes are acquired with a temporal scan in which each line corresponds to a given time (white arrows). The orange thick line shows the Sun direction. The spacecraft is approximately in a terminator orbit with a phase angle of 90°, so that one side of the comet is illuminated by the Sun and the other side is in darkness.</figcaption></figure>

radiation results as boundary condition. We analyzed configurations with and without a dust mantle over the dust-ice mixture for which we considered two different values of thermal conductivity for the dust, in order to compute the dependence of the surface and subsurface temperatures on the thermodynamic parameters such as density, specific heat, and thermal conductivity.

To distinguish between transient events and long-lasting features such as jets (which are stable for more than one comet rotation, as shown by Vincent et al. 2016), we adopted the following criteria: the transient events observed by VM were identified by their light curve (radiance at a given wavelength versus time) characterized by a sudden brightness increase in the coma that is associated with a release of gas and dust over a very short timescale, that is, 5–30 min (Bockelée-Morvan et al. 2017; Rinaldi et al. 2018). The light curve of a transient event can only be derived when the direction of the dust emission is perpendicular to the VM slit, for instance, when the scan occurs along the radial direction of the plume. The VM data are acquired with a temporal scan, as shown in Fig. 1, in which each line corresponds to a given time. This allows us to reconstruct the temporal evolution of the event. However, in the observation analyzed in this work, the dust emission direction is parallel to the VM slit, which prevents us from determining the light curve. Therefore we cannot assess whether the event is an outburst or transient, and we refer to it as a plume.

Section 2 presents a VM observation acquired on 2015 July 18. We study the physical properties of the dust ejecta (Sect. 3), including the radial profile (Sect. 3.1), color (Sect. 3.2), filling factor, and dust mass loss (Sect. 3.3). We analyze the illumination conditions of the likely emitting region and the thermal inertia of that region in Sects. 4–6. In Sect. 7 we investigate the 3D thermophysical model in order to evaluate the surface and subsurface temperature of the emitting region. Conclusions are provided in Sect. 8.

<footnote n="1">¹ www.comsol.com</footnote>

<footnote n="2">² www.meshlab.net</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Radiance at 0.55 µm on VIRTIS-M data cube V1_00395824539 acquired on 2015 July 18. *Panel A:* VM image at 0.55 µm. *Panel B:* same image, but the coma has been stretched to visualize the plume on the shadowed side (indicated by the oval). The cube details are listed in Table 1. The radiance has units of W m⁻² sr⁻¹ µm⁻¹.</figcaption></figure>

## 2. Observation and data reduction

The results presented here were derived from data obtained with the VIRTIS experiment on board the Rosetta spacecraft. VIRTIS is composed of two instruments: VIRTIS-H (VH; Drossart et al. 2000), a cross-dispersing point spectrometer providing spectra with high spectral resolution in the range 1.9–5 µm, and VIRTIS-M (VM; Coradini et al. 2007), an imaging spectrometer working in the range 0.25–5 µm by means of two channels, the visible (VIS) and the infrared (IR) channels, which cover the 0.25–1 µm and the 1–5 µm spectral intervals, respectively. Because of the failure in early May 2015 of the cryocooler, which is necessary to operate the IR channel, our analysis is restricted to VIS hyperspectral images. This dataset was continuously acquired for the entire duration of the mission. The spectra and images used for the analysis were reduced using the VIRTIS calibration pipeline (Ammannito et al. 2006; Filacchione et al. 2006), with additional corrections derived from in-flight data. We removed defective pixels and cosmic-ray strikes using a median filter despiking algorithm, which was employed only in the spatial dimensions of the data and therefore left the spectral data intact.

We here analyze a hyperspectral cube acquired by VM on 2015 July 18, which shows night-side dust activity (Figs. 1 and 2). This image cube was acquired one month before the

perihelion passage at a heliocentric distance of 1.28 au and a distance from the comet nucleus of 189 km. The corresponding Field of View (FOV) is about 12 km in the horizontal direction and 6 km in the vertical one, referred to the distance of the comet center to the spacecraft. The spacecraft was approximately on a terminator orbit so that the Sun illuminated one side and the other side was in darkness (Fig. 1).

Figure 2 displays the intensity map of the dust continuum in units of W m⁻² sr⁻¹ µm⁻¹. The map is a composite image where the comet nucleus (as an average in the wavelength range 0.45–0.55 µm) is superimposed on the image of the dust continuum averaged in a bandpass of 0.10 µm centered on 0.55 µm. The Sun is at the top of the image, and on the day side, the dust activity shows a predictable behavior that is correlated with the illumination conditions. The dust continuum image displays a narrow, collimated plume that is rather sharply bounded on the image plane. The dust plume is seen to expand as a narrow ejecta projected on the FOV plane. The maximum intensity of the dust plume is 0.098 W m⁻² sr⁻¹ µm⁻¹.

Table 1 provides the geometry information for the hyperspectral cube calculated by a routine (Acton 1996) that uses the spacecraft trajectory and orientation stored in SPICE kernels, and the 67P SHAP5 shape model for the comet nucleus (Jorda et al. 2016). The image was acquired line by line by means of a scanning mirror, taking 20 seconds per line. The final image
<table><caption>**Table 1.** VIRTIS-M observation analysed in this work.</caption><thead><tr><th>VIS File name</th><th>Cube dimension</th><th>Spatial resolution (m px<sup>-1</sup>)</th><th>Exposure time (sec)</th><th>Start time (UTC)</th><th>Duration time (sec)</th><th>Distance (km)</th><th>Phase angle (°)</th><th>Sub/Solar lat (°)</th><th>Sub/Solar long (°)</th><th>r<sub>h</sub> (au)</th></tr></thead><tbody><tr><td>V1_00395824539</td><td>256 126 432</td><td>47.25</td><td>16</td><td>07:16:57</td><td>2659</td><td>189</td><td>90</td><td>-35.22</td><td>164.04</td><td>1.28</td></tr></tbody></table>

**Note.** Column 1: observation file name. Column 2: cube size in number of samples, number of scan lines and spectral bands (432 for each channel). Column 3: pixel size at the distance of the observation. Column 4: exposure time for each line. Column 5: start time of the image cube (UTC). Column 6: total duration time for the image cube from acquisition start to stop. Column 7: distance of spacecraft from the comet center. Column 8: Observation phase angle. Column 9: subsolar latitude. Column 10: subsolar longitude. Column 11: heliocentric distance.

<table><caption>**Table 2.** Dust plume properties in the VIS channel.</caption><thead><tr><th>VIS file name</th><th>Detection time (h)</th><th>Local time (h)</th><th>Longitude range (°)</th><th>Latitude range (°)</th><th>Max radiance (W m<sup>-2</sup> sr<sup>-1</sup> µm<sup>-1</sup>)</th><th>Color %/100 nm</th></tr></thead><tbody><tr><td>V1_00395824539</td><td>[7:42, 7:48]</td><td>[20.20, 21.19]</td><td>[275, 285]</td><td>[-10, 0]</td><td>0.10</td><td>13.0 ± 2.6</td></tr></tbody></table>

**Note.** Column 1: observation file name. Column 2: start and stop detection time for the dust plume. Column 3: local time range for the source of the dust plume on the surface (see Sect. 5). Column 4: longitude range for the source of the dust plume on the surface (see Sect. 5). Column 5: latitude range for the source of the dust plume on the surface (see Sect. 5). Column 6: radiance at 0.55 µm at the emission maximum of the dust plume. Column 7: color at the radiance maximum of the dust plume (see Sect. 3.2).

is composed of a sequence of consecutive lines in the vertical direction.

In Table 2 we list all the relevant information obtained from the analysis of our observation, such as dust plume observation time, longitude and latitude of the estimated source region of the dust plume (Sect. 5), radiance level, and color at the maximum of the dust plume emission (Sect. 3.2).

The analysis of the VM dust continuum can be limited by infield stray light when the instrument slit is partially filled by the bright nucleus, and a sizable portion of the incoming photons is spread into the adjacent coma pixels. For the plume studied here, there is no stray light because the nucleus is in the dark, and its signal is therefore close to zero.

## 3. Properties of night-side dust ejecta

In this section, we analyze the physical properties of the dust plume in terms of radial profile, color, filling factor, and dust mass loss.

### 3.1. Radial profile

A particular feature of this observation is that the projected radial direction of the dust plume is along the instrument slit (horizontal red line in Fig. 1), spanning more than ten lines in the vertical direction that were acquired in about 3 min. For the first time, we therefore have the possibility to study the radial profile of a dust plume with VM data and to obtain information on the applicability of the fountain model to the cometary coma. In an ideal model, in which spherical continuity is satisfied, the radial column dust particle density profile $n(\rho)$ is expected to follow the expression (Eddington 1910)

$$n(\rho) = \frac{Q}{4\rho v(\rho)}, \qquad (1)$$

where $\rho$ is the distance from the comet center, $Q$ is the production rate of the number of particles of a given size, and $v(\rho)$ is

the particle outflow velocity. If there is no acceleration of the particles that are lifted off the surface, the observed scattered light intensity along the FOV should follow a simple $1/\rho$ profile. Deviations from a $1/\rho$ relationship reflect changes in the dust component. Several processes can modify the radiance profiles: dust acceleration or a lateral dust flow can increase the slope, while optical depth effects and extended or secondary sources act to flatten it (Gerig et al. 2018).

In Fig. 3, panel A presents the radial profile of the observed radiance at 0.55 µm, averaged over ten lines, versus $\rho$. The distance, $\rho$, is calculated using the angular distance of each pixel from the center of the nucleus and the distance of the spacecraft to the center of the comet. According to the description of the shape model by Preusker et al. (2017), this center is defined as the center of mass. The profile is characterized by a steep increase in the first 1 km from the nucleus surface, followed by a smooth decrease for 1 km after maximum was reached. Panel B shows the radiance multiplied by $\rho$ as a means for better visualizing changes in the coma. Beyond 4 km, the profile multiplied by $\rho$ has an approximately constant value, independent of $\rho$. This is consistent with a cometary coma in a steady state, with constant dust production and outflow speed together with the conservation of dust grains. The decrease observed between 2.8 and 4 km can be explained by dust acceleration. A similar behavior has been reported in other comets by Belton et al. (2013), Tubiana et al. (2015), Knollenberg et al. (2016), and Miles et al. (2016). The vanishing radiance values at distances shorter than 2 km show that the inner parts of the plume are not illuminated by sunlight because of shadowing (see Sect. 5).

### 3.2. Color

The dust emitted from a comet in a plume or jet is often characterized by the color, the normalized reflectivity gradient, or reddening, measured in percent/100 nm (Jewitt & Meech 1986), which yields information about the dust scattering properties from which factors such as particle size, shape, and composition can be inferred. The color can be calculated using the values of
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Top panel: radiance of the radial profile of the dust plume at 0.55 µm on July 18 along the distance from the nucleus center (ρ). Bottom panel: radial profile multiplied by ρ vs. ρ.</figcaption></figure>

the reflectance at two or more wavelengths. The reflectance, R, is a dimensionless quantity calculated by dividing the measured scattered light intensity I by the solar incident flux. When we take the wavelength-dependent solar flux from Kurucz (1994), the mean reflectance gradient per 100 nm, r, for a particular wavelength interval becomes

$$r = \frac{R_{\lambda_2} - R_{\lambda_1}}{\lambda_2 - \lambda_1} \times \frac{200}{R_{\lambda_2} + R_{\lambda_1}}. \quad (2)$$

The reflectances we used are averages over a narrow band-pass of 10 nm in width centered on 550 nm (λ₁) and 750 nm (λ₂). This was chosen to optimize the signal-to-noise ratio (S/N) and so minimize the internal error in the color determination. We obtained a 2D color map using the spectrum for each pixel in the image (Fig. 4). The color uncertainties were evaluated with the method used by Rinaldi et al. (2018), by propagating the formal error that is inversely proportional to the S/N. Outside the dust plume, the uncertainty is higher because the radiance and the S/N are both very low (Fig. 2). For this reason, the color fluctuations outside the dust plume region delimited by the contour lines in Fig. 2 are not realistic. Inside the plume, the radiance and the S/N give us an uncertainty of about 15–20%.

The 2D color map does not show evidence of different reddening values with respect to the surrounding coma. In Fig. 4 the black contour lines are the radiance levels at 550 nm. The red color, more than 18%/100 nm, outside the larger contour is noise. The value of the color in the upper part of the dust plume is close to 13 ± 2.6%/100 nm, which is essentially the same as the values below and above, which are around 15–16 ± 4.8%/10 nm. There is no evident color gradient in this dust plume, in contrast

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Panel A: dust continuum image of the July 18 plume at 0.55 µm. Panel B: spatial distribution of the color calculated in the yellow square of the dust image. The black contours are the radiance levels at 0.55 µm. The dust ejecta do not display an evident color gradient, as observed in other outbursts observed by VM (Rinaldi et al. 2018).</figcaption></figure>

to outburst ejectas (Bockelée-Morvan et al. 2017; Rinaldi et al. 2018).

### 3.3. Filling factor and dust mass loss

In the method applied in Rinaldi et al. (2018), the expression for computing the filling factor is

$$f_i = \frac{4\pi I_{\text{obs}}(a, \lambda)}{F_i(\lambda)A}, \quad (3)$$

where $I_{\text{obs}}(a, \lambda)$ is the observed scattered light intensity by the instrument (in W m⁻² sr⁻¹ µm⁻¹), a is the dust particle radius, $F_i(\lambda)$ is the incident solar flux, and $A = Q_{\text{sca}} \times p(g)$, where $Q_{\text{sca}}$ is the scattering efficiency and $p(g)$ is the phase function. Rinaldi et al. (2018) assumed a different albedo for the quiescent dust coma and for the outburst dust particles (Bockelée-Morvan et al. 2017; Agarwal et al. 2017) because the color map showed an evident color gradient. For this dust plume the color gradient is weak, and to compute the filling factor, we therefore assumed the quantity $wp(g)$ as albedo, which was calculated by Ciarniello et al. (2015), where $w = Q_{\text{sca}}/Q_{\text{ext}}$ is the single scattering albedo at 0.55 µm. The extinction efficiency $Q_{\text{ext}}$ is equal to 1, and so $w$ is just the scattering efficiency $Q_{\text{sca}}$. The value for the quiescent dust coma at 90° of phase angle is 0.04. The value of the filling factor was obtained at the maximum of the plume radiance, where the value is 3.0%.
<table><caption>**Table 3.** Observed total ejected dust mass.</caption><thead><tr><th rowspan="2">Bulk density (kg m<sup>-3</sup>)</th><th rowspan="2">Range of radius (μm)</th><th colspan="6">Ejected mass (tons)</th></tr><tr><th colspan="6">α</th></tr></thead><tbody><tr><td></td><td></td><td>-4.5</td><td>-4.0</td><td>-3.5</td><td>-3.0</td><td>-2.5</td><td>-2.0</td></tr><tr><td>795</td><td>0.1–1000</td><td>0.02</td><td>0.05</td><td>0.56</td><td>6.05</td><td>18.81</td><td>27.94</td></tr></tbody></table>

**Note.** Column 1: particle bulk density. Column 2: range of particle radius used for the calculations. Column 3–8: observed total ejected dust mass, computed for different differential size distribution indices $\alpha$.

To calculate the observed dust mass that was ejected, we converted the filling factor to a total ejected dust mass, assuming the differential size distribution of the dust grains ($n(a)da$), where $a$ is the grain radius, and the particle bulk density ($\rho_{bulk}$):

$$f_f = \frac{N\sigma_{geom}}{d^2}, \qquad (4)$$

where $N$ is the number of particles inside the pixel area $d^2$, and $\sigma_{geom}$ is the geometrical particle cross section. The mass of material responsible for the VIS dust emission inside the pixel area $d^2$ was computed with the equations reported by Rinaldi et al. (2018):

$$M_{pixel} = \frac{4}{3} \pi \rho_{bulk} f_f d^2 \frac{\int_{a_{min}}^{a_{max}} a^3 n(a) da}{\int_{a_{min}}^{a_{max}} \pi a^2 n(a) da}, \qquad (5)$$

where $a_{min}$ and $a_{max}$ are the minimum and maximum grain radii. The minimum and maximum sizes considered in this work are $a_{min} = 0.1$ µm and $a_{max} = 1000$ µm, which takes a wide range of particle sizes into account.

In the formulas above, the total observed dust mass that is ejected in the plume is the integral over a large area containing the dust plume. By integrating the mass for each pixel ($M_{pixel}$) over the area where signal from the dust plume is present, we can determine the observed total dust mass that is emitted during the event. In our calculations we used a constant bulk density of 795 kg m$^{-3}$ for all particles (Fulle et al. 2016).

We integrated over a rectangle extending from the closest nucleus boundary, which was chosen to avoid capturing any signal from the nucleus, to the final edge of the plume. The total mass computations are reported in Table 3.

Assuming a typical size distribution taken from the literature on 67P, with size distribution indices between -2.5 and -3 (Vincent et al. 2016), we estimate the observed total mass in the plume to be between 6 and 19 tons.

## 4. Illumination conditions

Figure 2 suggests that this dust event originates from the dark side of the neck region. In order to evaluate the illumination conditions at the observation time, we have to investigate the observation geometry that occurred before, after, and at the time of the image cube acquisition. To calculate the local illumination at any position on the surface of the comet nucleus, we need to know the orientation of the cometary rotation axis and the rotation period. As in Kömle et al. (2017), we considered the Sun as a point light source because the distance from the comet is very

large compared to the size of the Sun. The position of the Sun in the Cartesian system of the comet is

$$\left\{ \begin{aligned} x_{\text{sun}}(t) &= r_h(t) \cos \theta_{\text{sun}}(t) \cos \left(\phi - \frac{2\pi t}{P_{\text{comet}}}\right) \\ y_{\text{sun}}(t) &= r_h(t) \cos \theta_{\text{sun}}(t) \sin \left(\phi - \frac{2\pi t}{P_{\text{comet}}}\right) \\ z_{\text{sun}}(t) &= r_h(t) \sin \theta_{\text{sun}}(t), \end{aligned} \right. \qquad (6)$$

where $\theta_{sun} = -35.22^{\circ}$ is the subsolar latitude, $\phi = 164.04^{\circ}$ the subsolar longitude, and $r_h$ is the heliocentric distance at the time of observation (see Table 1). The shape model adopted in this calculation is SHAP5 (Jorda et al. 2016), adapted to a format that is suitable for the software COMSOL Multiphysics 5.3a<fnref n="3" />. In particular, we used the COMSOL module “Surface-to-Surface Radiation”.

The hemi-cube method was used to calculate the shadowing effects on the cometary surface (Formisano et al. 2018). In Figs. 5 and 6, we report the illumination conditions at 1.28 au. This corresponds to the heliocentric distance at which the dust plume observed by VIRTIS occurred on 2015 July 18. In Fig. 5 the comet is oriented as in Fig. 2, from two hours before the start time of the observations (7:16 h solar local time) to one hour after. In all the frames shown in Fig. 5, we observe that the neck, which is probably the region from which the dust plume is emitted, is entirely in shadows. For this reason it is interesting to analyze the illumination conditions on the far side of the comet, that is, the side that was not covered by the VM observation discussed here (see Fig. 6). The relevant fact is that even if the neck was not illuminated at the onset of the dust plume, it was illuminated until some hour before the observation starting time, so that the emission event could have been activated well before that time.

To validate our illumination calculation tool, we report in Appendix A the comparison between three images of Rosetta's Navigation Camera (NavCam) at different times. We also present our illumination results.

## 5. Source location of the dust plume

The calculation of the illumination conditions critically depends on the identification of the probable location of the source of the dust plume event. We selected three regions on the basis of the illumination conditions at different times (see Fig. 7): the first is in the neck at 4:35 h (panel A). The second is located in the head of the comet, at 6:15 h (panel B). Finally, the third region is at the bottom of the neck, at 6:56 h (panel C). To identify the most likely emitting region, we drew the normal (black lines in Fig. 7) to the facets covering those regions, assuming that because the dust plume is a collimated and narrow jet (Fig. 2), the ejected dust particles have been accelerated away from the surface in a direction perpendicular to the surface. By comparing the directions of the normal to the configuration of Fig. 2, it appears that the location we search for is probably that from panel A of Fig. 7, which is mainly in the Anuket region and partially in the Hapi region, as shown in Fig. 8B (red rectangle). Anuket is a consolidated region that is mostly smooth at intermediate and large scales, but with small-scale roughness that gives it a rocky appearance. Its boundary with Hapi is extremely well defined (Thomas et al. 2018), characterized by textural and topographic discontinuities in the local terrain (Vincent et al. 2016; Fornasier et al. 2019). It is tempting to conjecture that the precise location of the source of the dust plume may be at the morphological

<footnote n="3">³ www.comsol.com</footnote>
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Solar illumination on the comet oriented as in Fig. 2, at 1.28 au on 2015 July 18 from 5:16 h to 8:16 h. The color scale indicates the cosine of the solar incidence above the surface. The time step between the frames is 20 min.</figcaption></figure>

boundary. As shown in Fig. 9, such areas are characterized by discontinuities in the local terrain that are either textural or topographic, as observed in Vincent et al. (2016) and Fornasier et al. (2019).

## Geolocalization of the source of the dust plume

In order to constrain the source region better, we measured in the VM image the position and orientation of the plume, and the size of the shadow region. Scattered light from the near-nucleus part of the jet is not observed. This section of the jet indeed lies in the shadow of the large lobe and therefore is not illuminated by sunlight. These geometric properties of the jet can be used to constrain its source region.

We assumed that the dust particles within the jet escaped perpendicularly to their source region. Using the shape model version SHAP7 (500k facets, Preusker et al. 2017), we derived the plume properties geometrically for nearly 15 000 facets inside the longitude range = [260°, 300°] and latitude range = [-40°, 20°]. We also determined the three parameters mentioned above at the mid-time of the VM observation of the plume (2015-07-18T07:45:00).

The calculations were restricted to the neck regions illuminated 1.75 h before the VM observation. The reference axis was defined as the pointing vector from the nucleus center to the Sun (toward +Y in Fig. 1). The orientation of the plume was defined as the position angle between the plume and the reference axis, the shadow size was defined as the distance
<figure><img src="image_8.png" /><figcaption>Fig. 6. Solar illumination on the far side of the comet, i.e., the side that was not covered by the VM observation. The color scale indicates the cosine of the solar incidence above the surface. The time step between the frames is 20 min.</figcaption></figure>

between the intersecting point of the plume vector to the Y-axis and the maximum radiance point along the plume, and the position of the plume is the distance between this intersecting point and the nucleus center (which is on the Y-axis). Figure 10 shows scatterplots of the residuals of the three parameters (plume orientation, shadow size, and plume position, in percentage) and the combination score for each facet projected on a longitude-latitude based frame. The residual of the parameter (R) and the identifying score (S) are defined as follows:

$$R = |P_{\text{sim}} - P_{\text{obs}}|/P_{\text{obs}} \qquad (7)$$

$$S = (R_{\text{ang}}^2 + R_{\text{ss}}^2 + R_{\text{d}}^2)^{0.5}, \qquad (8)$$

where $P_{sim}$ is the derived parameter for the facet, and $P_{obs}$ is the observed value on the VM image. $R_{ang}$, $R_{ss}$, and $R_{d}$ are the residuals of the plume orientation, shadow length, and jet position in percentage. The measured position angle, shadow length, and position distance are 96.5°, 1.664 km, and 0.379 km, respectively. The possible source regions should have a score close to zero. The score panel in Fig. 10 shows that the jet source is in a very narrow region for which the latitude is between -10° and 0° and the longitude between 275° and 285° (see Table 2 and green rectangle of Fig. 8B).

The narrow source region is in the Anuket region and is characterized by a 500 m long fracture at 283–284° longitude and -7 to +18° latitude. This fracture was spotted early in the Rosetta mission, and was observed in OSIRIS images to evolve over the
<figure><img src="image_9.png" /></figure>

**Fig. 7.** Definition of the normals (black lines) to the possible emitting regions by observing the illumination conditions of Fig. 6 at different local times. The color scale refers to the cosine of the solar incidence. **Panel A:** emitting region on the neck at 4:35 h (peak illumination), see panel A of Fig. 8. **Panel B:** another possible emitting region (in the head of the comet) at 6:16 h. **Panel C:** bottom region at 6:56 h. By observing the normal directions and comparing them with the emission of Fig. 2, the emitting region of panel A seems to be a more likely candidate.

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

**Fig. 8.** **Panel A:** mean cosine of the solar incidence angle vs. rotation time of the likely emitting region (Fig. 7): the red interval ranges from about 4:35 h to 6:15 h, the time span in which the likely emitting region is illuminated. We selected 4:35 h because it corresponds to the peak of the cosine plot, while 6:15 h corresponds to the moment in which the region starts to move into shadow. The blue interval represents the observation time of the dust plume (7:42–7:48 h). We note that during the dust plume, the emitting region is completely in shadow. **Panel B:** source region of the night-side dust activity detected by VM. The red rectangle indicates the emitting region calculated with the method described in the first part of Sect. 5, with a latitude ranging from –30° to +10° and a longitude from 265 to 280°. The green rectangle indicates the source region constrained by measuring the length of the plume shadow. The region is mainly in the Anuket and Hapi regions.

two years of study. Already by 2014 December, it had extended by about 30 m, which is thought to be linked to the increasing spin rate of the comet before perihelion. Furthermore, in OSIRIS images taken in 2016 June, a new fracture that was at least 150 m long and perhaps extended up to 300 m was identified parallel to the original fracture (Thomas et al. 2018). Thomas et al. (2015) first described the neck–head boundary as probable tension fractions. The fractures and lineaments are dark-looking features with a width of a few meters and are some hundred meters long. They are characterized by sinuous outlines and variable width, suggesting that they are probably the surface expression of tension fractures and outburst sources (Skorov et al. 2016; Höfner et al. 2017).

In Sect. 3.3 we estimated the observed total mass of the dust plume to be between 6 and 19 tons, assuming a typical size distribution with indices between –2.5 and –3. The excavated volume lies in the range of 7.5–24 m³, which corresponds to an emitting area of 750–2390 m² for a vertical thickness of 10 mm (Sect. 7). When we consider the calculated emitting area as a portion of this fracture of a few meters wide (i.e., 5 m), this corresponds to the formation of a 150–478 m long fracture. When we

consider the calculated emitting area as a circular pit, we find a path with a radius of 150–270 m. Pits like this have been observed elsewhere on the comet, possibly caused by sinkhole collapse associated with the formation of previous outbursts (El-Maarry et al. 2015, 2016; Vincent et al. 2016).

## 6. Thermal inertia of 67P

Temperature measurements of the cometary surface show that the inferred thermal inertia is low, so there should be a small thermal lag that can activate a dust and gas emission a few hours after illumination (Groussin et al. 2007, 2013; Schloerb et al. 2015; Tosi et al. 2019). We here propose a new method for inferring the thermal inertia of a particular region of the surface of the comet. Previous methods are generally based on the estimating the thermal inertia from the measured temperature (e.g., Capria et al. 2014) or using daily amplitude of surface soil heat flux and temperature (Wang et al. 2010).

The key requirement is estimating the time span between the starting point, or peak, of the dust plume and the last time that
<figure><img src="image_11.png" /><figcaption>Fig. 9. OSIRIS NAC image acquired on 2016 June 8 at 15:26:38 h and projected onto the shape model of 67P. A 500 m long fracture in the neck region of comet 67P can be seen at the bottom of this frame along the +280° meridian.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 10. Simulation of the orientation, shadow length, and position of jets originating from the neck regions illuminated 1.75 h before the VIRTIS-M observations. Each point corresponds to the center of a facet in shape model version cg-dlr_spg-shap7-v1.0, with 500 k facets. In the right panel, the score (see text) is plotted. Values of the score close to zero correspond to the best fit to the observations.</figcaption></figure>

the emitting region (the border between Hapi and Anuket) is illuminated by the Sun. When we assume that this time span is equal to the conductive timescale, which is defined as

$$ \tau = \frac{\rho_m c_p L^2}{K}, \qquad (9) $$

where L is the depth of the upper boundary of the ice-dust mixture layer, K is the thermal conductivity, $\rho_m$ is the mean

density and $c_p$ is the specific heat, we can derive the thermal conductivity and consequently the thermal inertia. The thermal inertia is deduced using the following equation:

$$ I = \sqrt{K \rho_m c_p}. \qquad (10) $$

We assumed L = 5 mm (Shi et al. 2016; Hu et al. 2017), $\rho_m$ = 533 kg m⁻³ (Pätzold et al. 2016) and $c_p$ = 1000 J kg⁻¹ K⁻¹
<figure><img src="image_13.png" /><figcaption>Fig. 11. Panel A: computation domain characteristic of the neck region. The top area (surface) is 1 cm² while the height is chosen larger than the diurnal skin depth, which is on the order of one centimeter. Self-heating between the facets of the surface is taken into account. Panel B: illustration of the vertical structures assumed in this work. Left panel: dust mantle at the top (5 millimeters) and beneath a mixture of dust and water ice (blue spheres) with thermal inertia 36 TIU (model M1; derived from the method described in Sect. 6) and 100 TIU (model M2). Right panel: no dust mantle, with direct exposure of the dust-ice mixture layer (model M3).</figcaption></figure>

(Kömle et al. 2017). The time span of Eq. (9) can be constrained by observing panel A of Fig. 8, where we show the cosine of the solar incidence angle of the likely emitting region that we previously defined (Fig. 7 A). The red box marks the time span between the peak of the cosine plot (4:35 h) and the moment at which the region starts to move into the shadow (6:15 h). The blue box indicates the observation time of the dust plume (7:42–7:48 h), during which the likely emitting region is not illuminated. When we assume $τ = 187$ min (start = 4:35 h, stop = 7:42 h), we obtain a thermal conductivity of $1.2 \times 10^{-3}$ W m⁻¹ K⁻¹ and consequently a thermal inertia of 25 TIU. Conversely, when we assume $τ = 90$ min (start = 6:15 h, stop = 7:42 h), we obtain a thermal conductivity of $2.5 \times 10^{-3}$ W m⁻¹ K⁻¹ and a thermal inertia of 36 TIU. The values 25–36 TIU are compatible with the range provided by Gulkis et al. (2015) and Spohn et al. (2015) and also with the estimation provided by several previous comet missions. Our estimate is also compatible with the values found by Marshall et al. (2018), who found a thermal inertia <80 TIU for the Seth, Ash, and Aten regions. The range of values of thermal conductivity we derived (i.e., $1.2–2.5 \times 10^{-3}$ W m⁻¹ K⁻¹) is compatible with the values used in the literature in thermophysical modeling (Shi et al. 2016; Hu et al. 2017).

Observing Fig. 8A, we note that when we consider the start time as the moment in which the mean cosine of the solar incidence attains a value of 0.5 (in the increasing part of the plot), the thermal inertia estimate does not change dramatically. In this case, $τ ≈ 6$ h, and the thermal conductivity being $6 \times 10^{-4}$ W m⁻¹ K⁻¹, a thermal inertia value of 18 TIU is found. Conversely, if the start time is assumed as the moment in which the mean cosine of the solar incidence is 0.5 in the decreasing part, the time span is about 1 h, and the thermal conductivity derived is $3.7 \times 10^{-3}$ W m⁻¹ K⁻¹. Consequently, the thermal inertia is about 45 TIU. As a result of this sensitivity analysis, we can extend the thermal inertia range from 20 to 45 TIU, even if 20 TIU would correspond to a value of thermal conductivity that is very low and probably hard to justify.

## 7. Thermophysical model

We applied a thermophysical model to study the surface and subsurface temperatures of the emitting area between Hapi and

Anuket, defined in Fig. 8B, and also to test whether the thermal inertia derived in previous section is compatible with the surface and subsurface temperatures found in the literature. The model was applied only to a small portion (1 × 1 × 5 cm) of the likely emitting region. The integration domain (see Fig. 11) is characterized by a surface area of 1 cm² and a depth of 4.5 cm, which is larger than the diurnal skin depth (typically a fraction of only a few centimeters). Self-heating between the facets of the surface of the integration domain is taken into account.

### 7.1. Basic equations and boundary conditions

A 3D finite element method was applied using the software COMSOL Multiphysics<fnref n="4" /> and in particular the module “Heat transfer in Solids”, in combination with the meshing program MESHLAB. We solved the heat equation numerically, and combined it with a mass conservation equation to treat the water vapor emission, with no internal source or thermal convection, as in De Sanctis et al. (2005, 2010), Formisano et al. (2016, 2018, 2019a), Capria et al. (2017).

$$ \left( \rho c_p(T) + \frac{\Phi L(T)}{RT} \frac{\partial p_{\text{sat}}(T)}{\partial T} \right) \frac{\partial T}{\partial t} = \nabla \cdot (K(T)\nabla T) + \Phi L(T)\nabla \cdot (D(T)\nabla p_{\text{sat}}). \quad (11) $$

In this equation, $Φ$ is the porosity, $p_{\text{sat}}(T)$ is the saturation pressure, $R$ is the universal gas constant, $L(T)$ is the latent sublimation heat, and $D(T)$ is the gas diffusion coefficient. At the surface the following energy balance was imposed (Formisano et al. 2018, 2019b):

$$ S_c (1 - a) \cos(Z) + Q_{\text{SH}} = -\mathbf{n} \cdot (K(r, T)\nabla T) + f L(T) \Gamma(T) + \varepsilon \sigma T^4, \quad (12) $$

where $Q_{\text{SH}}$ is the self-heating term and takes the mutual radiative interaction between the facets of the surface into account, $a$ is the albedo, $S_c$ is the local irradiation intensity in W m⁻², $\cos(Z)$ is the cosine of the solar incidence angle, $\varepsilon$ is the emissivity (set to 0.97), $f$ is the area fraction covered by ice (if present on the surface), $\sigma$ is the Stefan-Boltzmann constant, and $\Gamma(T)$

<footnote n="4">⁴ www.comsol.com</footnote>
is the water vapor production rate in kg m⁻² s⁻¹ defined by the following formula (Delsemme & Miller 1971):

$$ \Gamma(T) = f p_{\text{sat}}(T) \sqrt{\frac{\mu}{2\pi RT}}. \qquad (13) $$

Here $\mu$ represents the water molar mass and $p_{\text{sat}}(T)$ is given by Fanale & Salvail (1984):

$$ p_{\text{sat}}(T) = A \exp(-\beta/T), \qquad (14) $$

where $A = 3.56 \times 10^{12}$ Pa and $\beta = 6141.667$ K. Local thermodynamic equilibrium was assumed. The term $fL(T)\Gamma(T)$ is included in Eq. (12) only in the case of ice on the surface. The sublimation rate calculated with Eq. (13) includes the icy area fraction factor as in Hu et al. (2017):

$$ f = \left(1 + \eta \frac{\rho_{\text{ice}}}{\rho_{\text{dust}}}\right)^{-1}, \qquad (15) $$

where $\eta$ is the dust-to-ice ratio in mass, $\rho_{\text{ice}}$ is the ice density, and $\rho_{\text{dust}}$ is the dust density. We note that $f = 1$ corresponds to a fully icy surface. At the bottom and on the vertical sides of our computational domain we imposed vanishing heat flux, and we assumed an initially constant temperature of 120 K (Kömle et al. 2017) in the whole domain, similar to the value adopted by Hu et al. (2017). The numerical grid of the computation domain shown in Fig. 11 A was chosen automatically by the COMSOL Multiphysics software, based both on the shape of the integration domain and on the physical equations involved.

## 7.2. Initial configurations

The fraction of the water ice on the surface of a comet seems to be very limited: for comet 67P, an upper limit of <1% has been estimated (Capaccioni et al. 2015). Sublimation of water ice arises from layers that are made of a mixture of dust and ice, beneath a cover of dust of a few millimeters. The dust-to-ice ratio by mass would be at least 6 (Sierks et al. 2015; Capaccioni et al. 2015; Rotundi et al. 2015; Filacchione et al. 2016).

In this paper we consider three different cases in order to evaluate the dependence of the surface and subsurface temperatures on the physical parameters adopted in this work. The first configuration (model M1) is characterized by a dust mantle of 5 mm with thermal inertia 36 TIU (or equivalently, $K = 2.5 \times 10^{-3}$ W m⁻¹ K⁻¹) above a layer made of a mixture of dust and ice (whose ratio is fixed to 6 for all the configurations). The second configuration (model M2) is similar to the first, but with a thermal inertia of 100 TIU, that is, the thermal conductivity is an order of magnitude greater ($2 \times 10^{-2}$ W m⁻¹ K⁻¹) than in the M1 case. This choice was made in order to show the effects of a high value of the thermal inertia on the surface and subsurface temperatures. The third and last configuration (model M3) represents an extreme case in which the dust-ice mixture is directly exposed without a dust layer. Model M1 uses a value of thermal inertia derived through the method we discussed in Sect. 6. Because the diurnal skin depth is very low with the typical values of density, specific heat, and thermal conductivity assumed for the comet surface, we studied only the first 4.5 cm below the surface. We recall that the diurnal skin depth is defined

(Mellon 1997) as

$$ \delta = \frac{I}{\rho c_p} \sqrt{\frac{P}{\pi}}, \qquad (16) $$

where $I$ is the thermal inertia and $P$ the rotational period. For the parameters assumed in our modeling, we obtain 0.8 cm for model M1, about 2 cm for model M2, and 0.3 cm for model M3. A schematic illustration of the vertical structures here adopted is given in Fig. 11B.

While the thermodynamic parameters ($\rho$, $c_p$ and $K$) of the dust were assumed to be independent of the temperature and the water-ice density, we assumed that the specific heat and the thermal conductivity of the water ice are functions of the temperature. In particular, the thermal conductivity of the water ice was corrected for by taking the high value of the porosity ($\Phi = 72-74\%$) into account (Pätzold et al. 2016) by applying the following equation as in Steiner & Kömle (1991):

$$ K_{\text{eff}} = (1 - \sqrt{1 - \Phi}) \Phi K_{\text{void}} + \sqrt{1 - \Phi} \left[ \lambda K_{\text{ice}} + (1 - \lambda) \left( \frac{B + 1}{B} \right) \frac{K_{\text{ice}} \times K_{\text{void}}}{K_{\text{ice}} + K_{\text{void}}} \right], \qquad (17) $$

where $K_{\text{ice}} = 567 K/T$ (Klinger 1981), $K_{\text{void}} = 4\epsilon\sigma T^3 r_p$ (Squyres et al. 1985), and the parameter $B$ is defined as (Steiner & Kömle 1991)

$$ B = 1.25 \left( \frac{1 - \Phi}{\Phi} \right)^{10/9}. \qquad (18) $$

We assumed as mean pore radius $r_p = 10^{-3}$ m (Steiner & Kömle 1991), $\Phi = 0.75$ (compatible with Pätzold et al. 2016) and $\lambda = 10^{-5}$, which represents the so-called flattening coefficient (Steiner & Kömle 1991). All physical parameters adopted in this work are listed in Table 4.

## 7.3. Temperature profiles

Figure 12 shows the calculated temperature profiles for the three different configurations we investigated. In panels A, C, and E we show the temperature at several depths as a function of time, and in panels B, D, and F, we show the temperature as a function of depth at a given time. In models M1 and M2 (panel A and C), the surface temperature is very high (up to 325 K) because the crust is covered by a dust blanket. In the case of low thermal conductivity (model M1), the layers at depths of 2.5 and 4.5 cm are not influenced by the external solar input. Conversely, in the case of model M2 (panel C), there is a general convergence of the profiles up to 0.50 cm because the thermal conductivity is higher. In this case, only the layer at 4.5 cm does not react to the solar input. In panel E the results of model M3 are shown: in this case, the exposure of the ice on the surface reduces the energy available to increase the temperature because part of it is required to sublimate the water ice. The decrease in temperature is evident: the maximum on the surface is about 220 K, in contrast to 325 K for the previous two cases.

Panels B, D, and F show the temperature profiles versus depth for several parts of a day-night cycle. All of the models show an inversion in the profile during the cometary day, linked to the formation of a recondensation front (De Sanctis et al. 2015). Model M1 exhibits a very large temperature gradient inversion (>50 K) in the first 1 cm beneath the surface; in model M2, the gradient is smaller but involves a larger depth, up to
<table><caption>**Table 4.** Physical parameters.</caption><thead><tr><th></th><th>Symbol</th><th>Value</th><th>References</th></tr></thead><tbody><tr><td colspan="4"><strong>General parameters</strong></td></tr><tr><td>Solar constant</td><td>S</td><td>1361 W m<sup>-2</sup></td><td>–</td></tr><tr><td>Initial temperature</td><td>T<sub>0</sub></td><td>120 K</td><td>Kömle et al. (2017)</td></tr><tr><td>Albedo</td><td>a</td><td>0.05</td><td>Shi et al. (2016); Hu et al. (2017)</td></tr><tr><td>Porosity</td><td>Φ</td><td>0.75</td><td>Pätzold et al. (2016)</td></tr><tr><td colspan="4"><strong>Dust parameters</strong></td></tr><tr><td>Density</td><td>ρ<sub>dust</sub></td><td>2600 kg m<sup>-3</sup></td><td>Pätzold et al. (2016)</td></tr><tr><td>Specific heat</td><td>c<sub>p,dust</sub></td><td>760 J kg<sup>-1</sup> K<sup>-1</sup></td><td>Shi et al. (2016); Kömle et al. (2017); Hu et al. (2017)</td></tr><tr><td>Thermal conductivity</td><td>K<sub>dust</sub></td><td>2 × 10<sup>-2</sup>, 2 × 10<sup>-3</sup> W m<sup>-1</sup> K<sup>-1</sup></td><td>Kömle et al. (2017); Hu et al. (2017)</td></tr><tr><td colspan="4"><strong>Ice parameters</strong></td></tr><tr><td>Ice density</td><td>ρ<sub>ice</sub></td><td>940 kg m<sup>-3</sup></td><td>Pätzold et al. (2016)</td></tr><tr><td>Ice specific heat</td><td>c<sub>p,ice</sub></td><td>7.037 T/K + 185.0 J kg<sup>-1</sup> K<sup>-1</sup></td><td>Ellsworth & Schubert (1983)</td></tr><tr><td>Thermal conductivity (*)</td><td>K<sub>ice</sub></td><td>567 K/T W m<sup>-1</sup> K<sup>-1</sup></td><td>Klinger (1981)</td></tr></tbody></table>

**Notes.** (*)The thermal conductivity of the water ice is “corrected” for the porosity according to Eq. (17).

2 cm beneath the surface. The inversion in model M3 is confined to less than 0.5 cm, and the typical temperature gradient is <20 K.

That the temperature profiles in model M1 converge from about 1.5 cm below the surface means that the deeper layers do not experience external solar energy input. This is less evident in the case of model M2, where the profiles converge from about 3 cm because the thermal conductivity is higher and a deep penetration of the thermal wave is possible. In the last case, model M3, the convergence of the thermal profiles begins from about 1 cm. The depth at which the profiles converge gives us an indication of the thermal skin depth, which is evidently not more than a few centimeters.

## 8. Conclusions

The Rosetta mission to comet 67P has given us the unprecedented opportunity of observing the regions of the coma within a few kilometers above the surface and to study dust activity on the night side (shadowed side) of the nucleus. We here analyzed a single VM hyperspectral cube. It was obtained on 2015 July 18, when the comet was close to perihelion at heliocentric distances near 1.28 au. At the time of the measurements, VM observed a dust plume that was identified as being generated on the far side of the nucleus relative to the spacecraft. Calculations of the illumination conditions at the time of the observation show that the identified spot was not directly exposed to the Sun, but that it was illuminated some hours before.

Shi et al. (2016) argued that the possible sources for these night-side dust activities are dusty terrains, and that their irregular distribution is possibly related to subsurface inhomogeneities. Our results show that this event can be sustained by a thermal lag of some hours in the shallow subsurface close to the depth of the water ice front, but it is uncertain whether it could be driven only by the sublimation of water ice.

Starting from the thermal lag, we found thermal inertia values between 25–36 J m<sup>-2</sup> K<sup>-1</sup> s<sup>-1/2</sup> for a 5 mm thick dust layer, which is compatible with the range obtained by Gulkis et al. (2015) and Spohn et al. (2015) and with the temperature measurements of the cometary surface performed by Groussin et al. (2007, 2013) and Tosi et al. (2019). In order to evaluate if these thermal inertia values are compatible with the estimation of

the surface and subsurface temperatures found in the literature, we applied a 3D Finite-Element Method (FEM) thermophysical model to a portion of the likely emitting region.

We explored three different configurations: (1) a layer of water ice mixed with dust beneath a dust mantle that is 5 mm thick and has thermal inertia of 36 J m<sup>-2</sup> K<sup>-1</sup> s<sup>-0.5</sup>; (2) the same structure, but with thermal inertia of 100 J m<sup>-2</sup> K<sup>-1</sup> s<sup>-0.5</sup>; and (3) an ice-dust mixture that is directly exposed. Temperature variations of about 50 K in the first 0.5 cm are predicted for the case of low thermal conductivity and the presence of a dust mantle, while layers at a depth greater than 2 cm are not influenced by the solar input. Increasing the thermal conductivity by a factor 10 gives temperatures of the layers in the first 0.5 cm that are very similar, but now the deeper layers (down to 2.5 cm) are affected by the solar input. The surface temperatures are in agreement with the results of Kömle et al. (2017).

An inversion in the temperature profile is present, in agreement with the results obtained by De Sanctis et al. (2015). In the M1 model, the inversion is present in the first 1 cm with the consequent formation of an “accumulation water ice front”. In case of an exposed dust-ice mixture and no dust coverage (M3 model), temperatures are greatly reduced (by up to 100 K) as a result of the energy that is used to sublimate the directly exposed water ice. This result also agrees with Kömle et al. (2017).

Of the three models we studied, M1 could represent the most likely configuration, both due to the low thermal inertia and because of a dust blanket that covers the ice-dust mixture. The model with high thermal inertia (M2), even if the surface temperatures are quite similar to those of M1 case, is probably less “realistic” than M1, in agreement with the results of Marshall et al. (2018). A configuration as in M3, where the ice-dust mixture is not covered by a dust mantle, could be possible, for example, if the mechanism of dust removal was more efficient than ice sublimation. This could be possible for short periods when the comet approaches perihelion. These numerical estimates could offer a more complete physical picture of the emitting source region.

The spectral properties of the dust in the dust plume show no detectable color gradient either in the plume itself or relative to the other parts of the coma that are visible in the observations
<figure><img src="image_14.png" /></figure>

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

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

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

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

<figure><img src="image_19.png" /><figcaption>**Fig. 12.** Temperature profiles vs. time and depth for the three models. Panels A and B: M1. Panels C and D: M2. Panels E and F: M3. The times selected in panels B, D, and F refer to the last rotation. Only for the main case, M1, do we report an inset in panel B with the surface temperature profile vs. time.</figcaption></figure>

analyzed here. The dust plume morphology can be classified as a narrow ejecta. The observed total ejected mass is estimated to be between 6 and 19 tons, assuming size distribution indices between −2.5 and −3. This would correspond to the formation of a fracture of a few meters wide and 150–478 m long, or to the formation of a circular pit on the surface with a radius of 150–270 m; pits like this have been observed elsewhere on the comet and might be due to sinkhole collapses associated with the formation of previous outbursts.

**Acknowledgements.** The authors wish to thank Norbert Kömle for his helpful suggestions and questions during the review process. We thank the following institutions and agencies for support of this work: Italian Space Agency (ASI, Italy) contract number I/024/12/1, Centre National d'Études Spatiales (CNES, France), DLR (Germany), NASA (USA) Rosetta Program, and Science and Technology Facilities Council (UK). VIRTIS was built by a consortium, which includes Italy, France, and Germany, under the scientific responsibility of the Istituto di Astrofisica e Planetologia Spaziali of INAF, Italy, which also guides the scientific operations. The VIRTIS instrument development, led by the prime contractor Leonardo-Finmeccanica (Florence, Italy), has been funded and managed by ASI, with contributions from Observatoire de Meudon financed by CNES,
and from DLR. We thank the Rosetta Science Ground Segment and the Rosetta Mission Operations Centre for their support throughout all the phases of the mission. The VIRTIS calibrated data will be available through the ESA's Planetary Science Archive Website (www.rssd.esa.int) and are available upon request until posted to the archive. This research has made use of the scientific software shapeViewer (www.comet-toolbox.com).

## References

Acton, C. H. 1996, Planet. Space Sci., 44, 65

Agarwal, J., Della Corte, V., Feldman, P. D., et al. 2017, MNRAS, 469, s606

A'Hearn, M. F., Belton, M. J. S., Delamere, W. A., et al. 2011, Science, 332, 1396

Ammannito, E., Filacchione, G., Coradini, A., et al. 2006, Rev. Sci. Instrum., 77, 093109

Belton, M. J. S., Feldman, P. D., A'Hearn, M. F., & Carcich, B. 2008, Icarus, 198, 189

Belton, M. J. S., Thomas, P., Carcich, B., et al. 2013, Icarus, 222, 477

Bockelée-Morvan, D., Rinaldi, G., Erard, S., et al. 2017, MNRAS, 469, S443

Capaccioni, F., Coradini, A., Filacchione, G., et al. 2015, Science, 347, aaa0628

Capria, M. T., Tosi, F., De Sanctis, M. C., et al. 2014, Geophys. Res. Lett., 41, 1438

Capria, M. T., Capaccioni, F., Filacchione, G., et al. 2017, MNRAS, 469, S685

Ciarniello, M., Capaccioni, F., Filacchione, G., et al. 2015, A&A, 583, A31

Coradini, A., Capaccioni, F., Drossart, P., et al. 2007, Space Sci. Rev., 128, 529

Davidsson, B. J. R., Gutiérrez, P. J., Groussin, O., et al. 2013, Icarus, 224, 154

De Sanctis, M. C., Capria, M. T., & Coradini, A. 2005, A&A, 444, 605

De Sanctis, M., Lasue, J., Capria, M., et al. 2010, Icarus, 207, 341

De Sanctis, M. C., Capaccioni, F., Ciarniello, M., et al. 2015, Nature, 525, 500

Delsemme, A., & Miller, D. 1971, Planet. Space Sci., 19, 1229

Drossart, P., Semery, A., Bouye, M., et al. 2000, in Infrared Spaceborne Remote Sensing VIII, eds. M. Strojnik, & B. F. Andresen, Proc. SPIE, 4131, 78

Eddington, A. S. 1910, MNRAS, 70, 442

El-Maarry, M. R., Thomas, N., Gracia-Berná, A., et al. 2015, Geophys. Res. Lett., 42, 5170

El-Maarry, M. R., Thomas, N., Gracia-Berná, A., et al. 2016, A&A, 593, A110

Ellsworth, K., & Schubert, G. 1983, Icarus, 54, 490

Fanale, F. P., & Salvail, J. R. 1984, Icarus, 60, 476

Farnham, T. L., Wellnitz, D. D., Hampton, D. L., et al. 2007, Icarus, 187, 26

Feaga, L. M., A'Hearn, M. F., Sunshine, J. M., Groussin, O., & Farnham, T. L. 2007, Icarus, 190, 345

Filacchione, G., Ammannito, E., Coradini, A., et al. 2006, Rev. Sci. Instrum., 77, 103106

Filacchione, G., de Sanctis, M. C., Capaccioni, F., et al. 2016, Nature, 529, 368

Formisano, M., De Sanctis, M. C., Magni, G., Federico, C., & Capria, M. T. 2016, MNRAS, 455, 1892

Formisano, M., Federico, C., De Sanctis, M. C., et al. 2018, J. Geophys. Res. Planets, 123, 2445

Formisano, M., Sanctis, M. D., Angelis, S. D., Carpenter, J., & Sefton-Nash, E. 2019a, Planet. Space Sci., 169, 8

Formisano, M., Federico, C., Magni, G., et al. 2019b, J. Geophys. Res. Planets, 124, 2

Fornasier, S., Hoang, V. H., Hasselmann, P. H., et al. 2019, A&A, 630, A7 (Rosetta 2 SI)

Fulle, M., Della Corte, V., Rotundi, A., et al. 2016, MNRAS, 462, S132

Gerig, S.-B., Marschall, R., Thomas, N., et al. 2018, Icarus, 311, 1

Groussin, O., A'Hearn, M. F., Li, J.-Y., et al. 2007, Icarus, 187, 16

Groussin, O., Sunshine, J. M., Feaga, L. M., et al. 2013, Icarus, 222, 580

Gulkis, S., Frerking, M., Crovisier, J., et al. 2007, Space Sci. Rev., 128, 561

Gulkis, S., Allen, M., von Allmen, P., et al. 2015, Science, 347, aaa0709

Höfner, S., Vincent, J.-B., Blum, J., et al. 2017, A&A, 608, A121

Hu, X., Shi, X., Sierks, H., et al. 2017, MNRAS, 469, S295

Jewitt, D., & Meech, K. J. 1986, ApJ, 310, 937

Jorda, L., Gaskell, R., Capanna, C., et al. 2016, Icarus, 277, 257

Klinger, J. 1981, Icarus, 47, 320

Knollenberg, J., Lin, Z. Y., Hviid, S. F., et al. 2016, A&A, 596, A89

Kömle, N. I., Macher, W., Tiefenbacher, P., et al. 2017, MNRAS, 469, S2

Kurucz, R. L. 1994, in Infrared Solar Physics, eds. D. M. Rabin, J. T. Jefferies, & C Lindsey, IAU Symp., 154, 523

Marshall, D., Groussin, O., Vincent, J.-B., et al. 2018, A&A, 616, A122

Mellon, M. 1997, J. Geophys. Res. Planets, 102, 25617

Miles, R., Faillace, G. A., Mottola, S., et al. 2016, Icarus, 272, 327

Pätzold, M., Andert, T., Hahn, M., et al. 2016, Nature, 530, 63

Preusker, F., Scholten, F., Matz, K.-D., et al. 2017, A&A, 607, L1

Rinaldi, G., Bockelée-Morvan, D., Ciarniello, M., et al. 2018, MNRAS, 481, 1235

Rotundi, A., Sierks, H., Della Corte, V., et al. 2015, Science, 347, aaa3905

Schloerb, F. P., Keihm, S., von Allmen, P., et al. 2015, A&A, 583, A29

Shi, X., Hu, X., Sierks, H., et al. 2016, A&A, 586, A7

Sierks, H., Barbieri, C., Lamy, P. L., et al. 2015, Science, 347, aaa1044

Skorov, Y. V., Rezac, L., Hartogh, P., Bazilevsky, A. T., & Keller, H. U. 2016, A&A, 593, A76

Spohn, T., Seiferlin, K., Hagermann, A., et al. 2007, Space Sci. Rev., 128, 339

Spohn, T., Knollenberg, J., Ball, A. J., et al. 2015, Science, 349, 6247

Squyres, S. W., P., M. C., & T., R. R. 1985, J. Geophys. Res. Solid Earth, 90, 12381

Steiner, G., & Kömle, N. I. 1991, Planet. Space Sci., 39, 507

Thomas, N., Sierks, H., Barbieri, C., et al. 2015, Science, 347, aaa0440

Thomas, N., El Maarry, M., Theologou, P., et al. 2018, Planet. Space Sci., 164, 19

Tosi, F., Capaccioni, F., Capria, M., et al. 2019, Nat. Astron., in press

Tubiana, C., Snodgrass, C., Bertini, I., et al. 2015, A&A, 573, A62

Vincent, J.-B., A'Hearn, M. F., Lin, Z.-Y., et al. 2016, MNRAS, 462, S184

Wang, J., Bras, R. L., Sivandran, G., & Knox, R. G. 2010, Geophys. Res. Lett., 37, L05404
# Appendix A: Validation of the illumination calculations

<table><caption>**Table A.1.** NavCam observations.</caption><thead><tr><th>NavCam file name</th><th>Start time (s)</th><th>Exposure time (s)</th><th>Distance (km)</th><th>Phase angle (°)</th><th>r<sub>h</sub> (au)</th></tr></thead><tbody><tr><td>ROSCAM1_20150718T060217</td><td>06:02:17</td><td>0.01</td><td>189</td><td>90</td><td>1.28</td></tr><tr><td>ROS CAM1_20150718T064703</td><td>06:47:03</td><td>0.01</td><td>189</td><td>90</td><td>1.28</td></tr><tr><td>ROSCAM1_20150718T090509</td><td>09:05:09</td><td>0.01</td><td>189</td><td>90</td><td>1.28</td></tr></tbody></table>

**Note.** Column 1: observation file name. Column 2: start time of the image cube (UT). Column 3: exposure time for each image. Column 4: distance from the comet center. Column 5: observation phase angle. Column 6: Heliocentric distance.

On 2015 July 18 the Rosetta Navigation Camera (NavCam) performed two observations at 6:02 and at 6:47 UTC, before the outburst ejection, and one after, at 9:05 UTC (Table A.1).

In order to validate our illumination calculation tool, we have compared the three NavCam images with our numerical results. Figure A.1 shows that a good agreement is met.

<figure><img src="image_20.png" /><figcaption>**Fig. A.1.** Left column: illumination of the cometary surface oriented as the NavCam observations shown on the right. The color scale indicates the cosine of the solar incidence. Right column: three observations acquired by the Rosetta NavCam on 2015 July 18 at 6:02, 6:47 and 9:05 UTC (Table A.1).</figcaption></figure>