# Land of gas and dust – exploring bursting cavities on comet 67P

Daniel R. Müller<sup>1,★</sup> Kathrin Altwegg<sup>1</sup>, Jean-Jacques Berthelier<sup>1,2</sup> Robin F. Bonny<sup>1</sup>

Michael R. Combi<sup>3</sup>, Johan De Keyser<sup>3,4</sup> Antea C. Doriot<sup>1</sup>, Stephen A. Fuselier<sup>5,6</sup> Nora Hänni<sup>1</sup>

Martin Rubin<sup>1</sup>, Susanne F. Wampfler<sup>7</sup> and Peter Wurz<sup>1,7</sup>

Accepted 2025 January 27. in original form 2025 January 6

## ABSTRACT

Gas and dust outbursts are recurring phenomena on comets, offering critical insights into their subsurface activities. On comet 67P/Churyumov–Gerasimenko, two distinct outburst types have been identified: CO₂-dominated ‘summer fireworks’ near perihelion and water-driven events often linked to cliff collapses outside the perihelion period. While CO₂-dominated outbursts are thought to originate from subsurface gas cavities, the properties of these cavities remain poorly understood. In this study, we modelled the outgassing dynamics and dust velocities of outburst events using Rosetta/ROSINA data to estimate the characteristics of subsurface gas cavities and their impact on ejected particle dynamics. Our results indicate that CO₂-dominated events involve subsurface cavities with radii ranging from 15 to 62 m for an equivalent half-sphere geometry, depending on gas distribution assumptions. Conversely, water-driven outbursts would require subsurface temperatures far above equilibrium, supporting the hypothesis of mechanical processes like cliff collapses exposing ices to sublimation. Dust velocities in CO₂-dominated events – while aligning with results from other Rosetta instruments – were notably higher across all grain sizes compared to water-driven events, reflecting distinct dynamics in dust ejection. These findings highlight the critical role of subsurface gas reservoirs in driving explosive outbursts and suggest a strong connection between cometary activity, volatile distribution, and structural conditions. This study emphasizes the need for high-resolution data on subsurface volatiles from future missions and more refined modelling and experiments to further elucidate these mechanisms, with potential broader implications for our understanding of cometary activity.

**Key words:** instrumentation: detectors—methods: data analysis—comets: general—comets: individual: 67P/Churyumov–Gerasimenko.

## 1 INTRODUCTION

Cometary outbursts are transient and dynamic phenomena observed across various comets, marked by sudden releases of gas and dust. These events offer valuable insights into the physical properties and internal structure of cometary nuclei, as well as the processes that drive their activity (Hughes 1991). Despite extensive research, the exact mechanisms that trigger these outbursts remain elusive. Proposed explanations range from subsurface pressure build-up due to volatile sublimation to structural changes like cliff collapse and fracturing of the nucleus (e.g. Hughes 1975; Prialnik, A’Hearn & Meech 2008; Vincent et al. 2016; Bockelée-Morvan et al. 2022). Another potential mechanism, proposed for comet 1P/Halley, involves the exothermic amorphous-to-crystalline ice phase transition, which can induce runaway internal heating, leading to explosive activity

(Prialnik & Bar-Nun 1992). Additionally, clathrate destabilization, associated with amorphous ice crystallisation, has been modelled as a possible driver of pit formation through sinkholes or violent outbursts (Mousis et al. 2015).

Comet 67P/Churyumov–Gerasimenko (hereafter 67P), studied extensively during the Rosetta mission, exhibited numerous outbursts. The mission provided a unique opportunity to explore drivers of outbursts in better detail. Previous studies have attributed outbursts on 67P to a variety of factors, such as mechanical stresses or the release of volatile gases from subsurface cavities (e.g. Skorov et al. 2016; Pajola et al. 2017). The continuous monitoring by Rosetta’s suite of instruments has been crucial in capturing both the dust and gas dynamics during these events (e.g. Grün et al. 2016; Agarwal et al. 2017).

In a recent study, Müller et al. (2024) analysed gas emissions during outbursts on comet 67P using data from the ROSINA Double Focusing Mass Spectrometer (ROSINA/DFMS). Their findings revealed key patterns in the comet’s gas composition, attributing some

<footnote n="1"><sup>1</sup>Space Research and Planetary Sciences, Physics Institute, University of Bern, Sidlerstrasse 5, CH-3012 Bern, Switzerland</footnote>

<footnote n="2"><sup>2</sup>Laboratoire Atmosphères, Milieux, Observations Spatiales (LATMOS), 4 Avenue de Neptune, F-94100 Saint-Maur, France</footnote>

<footnote n="3"><sup>3</sup>Department of Climate and Space Sciences and Engineering, University of Michigan, 2455 Hayward, Ann Arbor, MI 48109, USA</footnote>

<footnote n="4"><sup>4</sup>Royal Belgian Institute for Space Aeronomy, BIRA-IASB, Ringlaan 3, B-1180 Brussels, Belgium</footnote>

<footnote n="5"><sup>5</sup>Space Science Division, Southwest Research Institute, 6220 Culebra Road, San Antonio, TX 78228, USA</footnote>

<footnote n="6"><sup>6</sup>Department of Physics and Astronomy, The University of Texas at San Antonio, San Antonio, TX 78249, USA</footnote>

<footnote n="7"><sup>7</sup>Center for Space and Habitability, University of Bern, Gesellschaftsstrasse 6, CH-3012 Bern, Switzerland</footnote>

<footnote>*E-mail: daniel.mueller@unibe.ch</footnote>
outbursts to subsurface gas reservoirs while linking others to surface processes such as cliff collapses. However, the properties of these subsurface gas reservoirs remain poorly understood.

In this study, we build on these findings by investigating gas production rates and outburst dynamics from ROSINA/DFMS data, focusing exclusively on the subsurface pressure build-up mechanism. Our analysis centres on CO₂-dominated outbursts observed during the comet’s perihelion period in the summer of 2015, using a model to simulate gas outflow driven by subsurface pressure build-up. This approach allows us to estimate the sizes and characteristics of these cavities and their role in driving explosive events. Additionally, we constrain the velocities of dust particles ejected during these outbursts, shedding light on the underlying mechanisms. The following sections outline the instrumentation, data processing, and modelling techniques used in this study, and present new findings that enhance our understanding of cometary outburst dynamics.

## 2 METHODS

### 2.1 ROSINA/DFMS instrumentation and data treatment

The ROSINA/DFMS instrument aboard the Rosetta spacecraft is a high-resolution mass spectrometer designed to measure the composition of cometary gases. The instrument uses a Nier-Johnson double-focusing ion-optical configuration, including a toroidal electrostatic analyser and a permanent magnet to separate ions by their mass-to-charge (m/z) ratio. The instrument achieves a mass resolution of 3000, measured at 1 per cent of the peak amplitude for m/z 28, as detailed by Balsiger et al. (2007).

DFMS scans through a range of m/z values using specific voltages applied to its ion optics. Ions are detected via a multichannel plate (MCP) detector with detection along the mass-dispersive direction. The detector gain is adjusted to obtain a high dynamic range. A typical scan, covering an m/z range from 13 to 100, takes about 45 min. Data are processed following established methods for calibration, species identification, and signal integration, and are normalized to neutral gas densities measured by ROSINA’s COmet Pressure Sensor (COPS) (Le Roy et al. 2015; Calmonte et al. 2016; De Keyser et al. 2019; Rubin et al. 2019). Detailed analysis methods, especially for outburst events, are disclosed in Müller et al. (2024).

### 2.2 Data analysis

To study cometary outbursts, we selected DFMS data from 67P’s perihelion period during the summer of 2015, when the so-called summer fireworks occurred (Vincent et al. 2016), and periods when other similar outbursts have been reported (Feldman et al. 2016; Grün et al. 2016; Agarwal et al. 2017; Rinaldi et al. 2018; Noonan et al. 2021). We focused on times when Rosetta was positioned within ±25° of an outburst source’s latitude and longitude. This angular window was selected to take into account both the physical processes driving the expanding gas dynamics (e.g. Combi et al. 2012; Marschall et al. 2020) and the DFMS observation timing which controls in particular the spatial resolution of measurements. Yet, the actual complexity of the nucleus shape is not taken into account as we consider an overall spherical geometry for simplicity reasons. Data normalization by spacecraft-source distance, corrected for varied viewing geometries, required careful data evaluation to link ROSINA measurements to specific events. For further details on the specific selection of events and data analysis protocols, the reader is referred to Müller et al. (2024).

<figure><img src="image_1.png" /><figcaption>**Figure 1.** Schematics of the spherical cone cap distribution. The outgassing behaviour of the outburst manifests as a conical expansion with a half angle α up to the points at which the gas flow rate is measured at the spacecraft, situated at a distance $r_{CGSC}$ from the comet nucleus. Using this outgassing model, the outbursting gas is distributed homogeneously across the surface of the grey cone cap.</figcaption></figure>

### 2.3 Modelling approach

In Müller et al. (2024), we proposed that subsurface gas reservoirs could explain the elevated gas production rates observed during outburst events. A similar concept had previously been discussed by Agarwal et al. (2017), who developed a model to explain such reservoirs and estimate their sizes. Following these approaches, we implemented a similar model in this work. According to Sharipov & Kozak (2009), the gas mass flow rate, or gas production rate, of the gas escaping from a pressurized container at pressure $p_{in}$ and temperature $T_{in}$ through an aperture can be expressed analytically. These authors consider a slit of width a and length l, but since the precise geometry of the opening is unknown, we represent it with an effective radius $r_{crack}$ in this work. The modelled gas mass flow rate is then given by

$$ \dot{M} = \frac{W\sqrt{m}}{\sqrt{2\pi k_B}} \frac{p_{in}r_{crack}^2\pi}{\sqrt{T_{in}}} = \frac{W\sqrt{\pi m}}{\sqrt{2k_B}} \frac{p_{in}r_{crack}^2}{\sqrt{T_{in}}}. \quad (1) $$

Here, m is the molecular mass of the gas, W is a dimensionless parameter that characterizes the flow regime (in this case, we assume $W = 1.5$ for viscous flow), and $k_B$ is Boltzmann’s constant. The pressure inside the container, $p_{in}$, is approximated by the sublimation pressure $P_{subl}(T)$, following the model by Fray & Schmitt (2009). The model from equation (1) assumes an expansion of the gas pocket against vacuum. In general, there might be a collisional thermalized region close to the comet surface, especially for high activity comets, which would lead to a reduced outward mass flux. However, assuming expansion against vacuum can be justified due to the low gas production of 67P and the strong outgassing during the outbursts.

For each outburst event, we computed the gas production rate using gas density measurements from ROSINA/DFMS and ROSINA/COPS. The temperature $T_{in}$ and corresponding internal pressure $p_{in}$ were then calculated to satisfy equation (1).

The peak gas flow rate was determined by using the peak gas density observed at the nearest measurement point relative to the event source and outburst time. Müller et al. (2024) previously defined a field of view (FOV) to evaluate whether ROSINA measured gases originating from the outburst or the quiescent coma, and the same FOV is used here. Therefore, the gas density measured by ROSINA/DFMS is assumed to be homogeneously distributed within the 25° half-angle cone (grey area in Fig. 1) and the gas flow rate at the spacecraft distance, $Q_{sc}$, is computed as follows:

$$ Q_{sc} = v_{sc}n_{sc}A_{sc} = v_{sc}n_{sc}2\pi r_{cgsc}^2(1 - \cos\alpha), \quad (2) $$
where $v_{sc}$ and $n_{sc}$ are the gas velocity and density at the spacecraft's location, respectively, and $A_{sc}$ represents the calculated area of the gas plume determined from the spacecraft's distance from the comet $r_{cgsc}$ and the cone opening angle $\alpha$.

The gas flow rate from the crack on the comet's surface is:

$$Q_{\text{crack}} = v_{\text{crack}}n_{\text{crack}}A_{\text{crack}} = v_{\text{crack}}n_{\text{crack}}\pi r_{\text{crack}}^2, \quad (3)$$

where $v_{crack}$ and $n_{crack}$ are the gas velocity and density at the crack, respectively, and $A_{crack}$ is the area of the crack opening, determined by the radius of the crack $r_{crack}$. With the assumption of thermalization inside the pocket, the gas velocity at the crack is derived from the gas effusion equation, given by $v_{crack} = \sqrt{\pi k_B T_{in}/2m}$.

To conserve the mass of the outgassing material, the two flow rates $Q_{sc}$ and $Q_{crack}$ must be equal, which leads to the following expression for the gas density at the crack, $n_{crack}$:

$$n_{\text{crack}} = \frac{v_{\text{sc}}}{v_{\text{crack}}} \frac{2r_{\text{cgsc}}^2(1 - \cos \alpha)}{r_{\text{crack}}^2} n_{\text{sc}}, \quad (4)$$

Using this equation, the mass production rate at the crack can be expressed as

$$\begin{aligned} \dot{M}_{\text{crack}} &= Q_{\text{crack}}m = v_{\text{crack}}n_{\text{crack}}m \pi r_{\text{crack}}^2 \\ &= v_{\text{sc}}n_{\text{sc}}m 2\pi r_{\text{cgsc}}^2(1 - \cos \alpha). \end{aligned} \quad (5)$$

The mass production rate at the crack is thus independent of surface parameters or other model assumptions, except for the cone angle $\alpha$, as defined by Müller et al. (2024). The gas velocity at the spacecraft, $v_{sc}$, is assumed to be 600 m s$^{-1}$, consistent with values during outburst events reported in the literature (e.g. Grün et al. 2016; Marschall et al. 2016).

Since only the outburst-related overproduction is considered, the mass production rate due to the outburst is calculated as follows: $\dot{M}_{out} = \dot{M}_{peak} - \dot{M}_{bg}$. The background gas production follows the same approach as in equations (2) to (5). The only difference to equation (2) is that the background gas flow is not limited to the conical outgassing, but is distributed over the whole hemisphere because of its origin from the whole nucleus. The cavity temperature $T_{in}$ and pressure $p_{in}$ are computed by fitting equation (1) to the gas production rate.

As the exact size of the cavity opening is unknown, we vary its parameter to identify the value that ensures that the cavity pressure remains below the surface tensile strength of comet 67P, estimated between 3 to 150 Pa (Groussin et al. 2015; Vincent et al. 2015; Basilevsky et al. 2016), and the temperature remains below the ice sublimation equilibrium temperature.

To estimate the ice sublimation equilibrium temperature, we use the model from Keller et al. (2015). This model assumes that the energy absorbed by the comet's surface is partially used to sublimate subsurface volatiles, with the surface temperature derived by balancing incident solar irradiation, thermal re-radiation, sublimation losses, and heat conduction. The simplest assumption is that exposed water ice sublimates directly from a flat surface which is the actual geometry of the observations with Rosetta looking towards a sunlit nucleus. Neglecting heat conduction (justified by the low heat conduction observed; Gulkis et al. 2015), the energy balance is expressed as (Keller et al. 2015):

$$(1 - A_v) I(t) = \epsilon \sigma T^4 + Z(T) L_{\text{ice}}, \quad (6)$$

where $I$ is the Solar irradiation intensity, $A_v$ is the Bond albedo (assumed to be 0.01; Keller et al. 2015), $\epsilon$ is the emissivity (0.9; Keller et al. 2015), $\sigma$ is the Stefan-Boltzmann constant, $T$ is the surface temperature, $Z$ is the sublimation rate, and $L_{ice}$ is the latent

heat of sublimation. We assume $L_{ice}$ to be constant, neglecting its slight temperature dependence. The sublimation rate is calculated using the Hertz-Knudsen formula (Keller et al. 2015), with the water vapour pressure derived from Fray & Schmitt (2009). The corresponding equilibrium temperatures for all studied events are provided in Table A5.

Assuming that the reservoir is not replenished during the outburst and that the cavity temperature remains constant, the pressure inside the cavity and the mass flux through the crack decay exponentially. The e-folding decay time, $\tau$, is proportional to the cavity volume and is given by

$$\tau = \sqrt{\frac{2m}{k_B T_{in}\pi} \frac{V}{W r_{\text{crack}}^2}}. \quad (7)$$

Müller et al. (2024) observed that gas enhancements lasted longer than the dust outbursts captured by Rosetta's cameras, likely due to continued sublimation from exposed ice surfaces on the order of tens of meters (e.g. Agarwal et al. 2017) following the removal of the dust layer. In our model, however, we assume that the duration of the gas outburst matches the dust event observed by the cameras. This approach focuses on the bursting phase of the event rather than the prolonged sublimation of newly exposed ice which has been used in earlier publications (e.g. Grün et al. 2016; Agarwal et al. 2017). The summer fireworks events require a $\tau$ value between 40 and 260 s to reduce the gas production rate to 1/1000 of its initial value (Agarwal et al. 2017) within 5 to 30 min (Vincent et al. 2016). Events outside the perihelion summer fireworks period necessitate different $\tau$ values, depending on the specific observation times of the individual dust events. All $\tau$ values are listed in Table A6.

From equation (7), the gas cavity volume is determined. Assuming the cavity to be a half-sphere with radius $r_{cavity}$, the equivalent half-sphere radius is calculated as

$$r_{\text{cavity}} = \left( 3\tau W r_{\text{crack}}^2 \sqrt{\frac{k_B T_{\text{in}}}{8\pi m}} \right)^{1/3}. \quad (8)$$

The values of $\tau$ and $r_{crack}$ are varied, and the results are presented in Section 3. The shape of the cavity is assumed to be a half-sphere because of previous publications using this shape as an equivalent reference shape and to facilitate the calculations over more complex internal structures such as branching structures beneath the surface (Agarwal et al. 2017).

Uncertainties have been calculated using Gaussian error propagation applied to all relevant variables of the formulas above. Most importantly, the uncertainty of the ROSINA/DFMS measurements (≈20 per cent) and the uncertainty in the initial thermal velocity (≈30 per cent) of the gas molecules had the largest influence on the error propagation.

# 3 RESULTS

## 3.1 Subsurface gas cavities

This subsection explores the characteristics of subsurface gas cavities and shows the differences between the perihelion summer fireworks and other outbursts, encompassing events outside this period. Through detailed analysis of gas production rates and cavity dimensions as introduced in Section 2.3, we seek a better understanding of the role of subsurface cavities in driving cometary outbursts and the implications these cavities have for cometary evolution.
<table><caption>**Table 1.** Gas production rates for the studied events during the summer fireworks described by Vincent et al. (2016).</caption><thead><tr><th>Gas</th><th>Min. (kg s<sup>-1</sup>)</th><th>Max. (kg s<sup>-1</sup>)</th><th>Mean ± SD (kg s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>H<sub>2</sub>O</td><td>3.8</td><td>10.2</td><td>6.7 ± 2.0</td></tr><tr><td>CO<sub>2</sub></td><td>0.9</td><td>7.9</td><td>3.2 ± 2.8</td></tr><tr><td>CO</td><td>0.1</td><td>1.1</td><td>0.4 ± 0.3</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>**Figure 2.** Gas production rate for H₂O and CO₂ for events where Rosetta flew over the corresponding event source region shortly after the event was detected. The horizontal lines and the shaded areas represent the mean value and the standard deviation for both gases, respectively. The points have been connected by dash-dotted lines for improved readability. The event IDs correspond to the values from Vincent et al. (2016).</figcaption></figure>

### 3.1.1 Summer fireworks

Our investigation concentrated on the immediate effects on the gas coma following dust outburst events, rather than on the long-term monitoring of the source regions. Therefore, we selected eight specific outburst events (see Table A1) in the data set of Vincent et al. (2016) and Müller et al. (2024), where Rosetta observed the source regions within a day of the event. Notably, six of these eight events were observed within a single full rotation period of comet 67P.

Our analysis of these eight events revealed that the H₂O gas production rates varied between 3.8 and 10.2 kg s⁻¹, with an average of 6.7 kg s⁻¹. The production of CO₂ ranged from 0.9 to 7.9 kg s⁻¹, with a mean value of 3.2 kg s⁻¹. In contrast, CO production was negligible, typically accounting for 1–7 per cent of the total gas production of the three gases only Table A1). A summary of these gas production rates is presented in Table 1, while the specific production rates for H₂O and CO₂ for each event are shown in Fig. 2.

These results indicate that H₂O production generally exceeded CO₂ production, suggesting a stronger water contribution. However, compared to the background production values, CO₂ showed a significant increase by a factor of 1.22 ± 0.11, while H₂O increased by only a factor of 1.05 ± 0.03, and CO rose by a factor of 1.07 ± 0.03 (Fig. 3). As already stated, these findings suggest that the outburst events were predominantly driven by CO₂, consistent with the conclusions of Müller et al. (2024).

If these summer fireworks events had been triggered predominantly by water, the required effective cavity opening radius would exceed 70 metres to reach the equilibrium temperature of approximately 202 K (Table A5). Such a large opening would likely have been detected by Rosetta's cameras after the event. In contrast, for CO₂ sublimation at temperatures below 200 K, an effective cavity opening radius of just 0.05 metres would suffice (Fig. 4). However, this would produce a sublimation pressure of approximately 130 MPa – far surpassing the tensile strength of comet 67P's surface layers,

<figure><img src="image_3.png" /><figcaption>**Figure 3.** Gas production rate at the crack during an outburst (Ṁc,peak) compared to background gas production rate (Ṁc,bg) for H₂O, CO₂, and CO for events where Rosetta flew over the corresponding event source region shortly after detection. The horizontal lines indicate the mean value of the outburst-to-background ratio for each gas. The uncertainties of the individual points and the average values have been omitted and the points have been connected by dash-dotted lines for improved readability. The event IDs correspond to the values from Vincent et al. (2016).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 4.** Cavity temperature as a function of the effective cavity opening radius, calculated using equation (1) for H₂O and CO₂. The uncertainties represent the standard deviation for the summer fireworks outburst events.</figcaption></figure>

which ranges from 3 to 150 Pa (Groussin et al. 2015; Vincent et al. 2015; Basilevsky et al. 2016). To remain within this pressure range, an average cavity opening radius of at least 1.4 metres would be required, resulting in an average pressure of 142 Pa and an average temperature of 138 K (Fig. 5, bottom). In the case of H₂O, using the same cavity radius would lead to pressures approximately four to five times higher (Fig. 5, top).

We also calculated the gas cavity radius, assuming a half-spherical pocket, based on equation (8) and varying the effective cavity opening radius and time constant τ in accordance with the model described in Section 2.3. With an effective cavity opening radius of 1.4 m – consistent with the pressure and temperature constraints – the resulting cavity radius ranges between 15 and 30 metres for an equivalent half-sphere cavity (Fig. 6).

Although outbursts are localized and focused phenomena, some studies adopt a hemispherical gas distribution (Grün et al. 2016; Agarwal et al. 2017; Noonan et al. 2021) rather than the spherical cone cap distribution we used here. The hemispherical distribution has the advantage that it is straightforward to define but ignores the focused outgassing behaviour of outbursts as seen on the Rosetta images. To compare with other studies, we recalculated the gas production assuming a hemispherical distribution by setting the cone angle α in equation (5) to 90° (see Table A2 for the results).
<figure><img src="image_5.png" /><figcaption>**Figure 5.** Pressure in the cavity as function of the effective cavity opening radius according to equation (1) and the corresponding sublimation pressure from Fray & Schmitt (2009) for H₂O (top) and CO₂ (bottom) for the summer fireworks outburst events. The mean pressure and standard deviation for these events are shown with the black curve and shaded area. The maximum tensile strength, R<sub>m</sub>, is shown by the horizontal line (dark red). The plots only show effective cavity opening radii starting at 1 m for improved visibility as smaller radii lead to much larger cavity pressures. The uncertainties of the individual points have been omitted for improved readability.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 6.** Cavity radius for an equivalent half-sphere pocket as a function of the effective cavity opening radius and the exponential decay constant τ for the summer fireworks events. The cavity radii vary minimally between events (Table A1), so only the mean value for each τ is provided here. Uncertainties of the individual points have been omitted for improved readability.</figcaption></figure>

In this case, H₂O production ranged between 41 and 110 kg s⁻¹, with an average of 71 kg s⁻¹, while CO₂ production varied between 9 and 84 kg s⁻¹, with an average of 35 kg s⁻¹. These values are approximately an order of magnitude higher than those obtained with the cone cap distribution. Under these conditions, an effective cavity radius of 4.5 metres would be required to match the pressure and temperature constraints – three times larger than with the cone cap distribution. Consequently, the equivalent half-sphere cavity radius

would be between 33 and 62 metres, twice as large as in the cone cap scenario.

Thus, although different distribution models yield significantly different gas production rates, the impact on the calculated cavity and opening size remains relatively modest – a factor of two to three – considering the many approximations involved in the model.

Vincent et al. (2016) estimated an ejected dust mass flux of 60 kg s⁻¹ for particles in the size range of 1–10 µm or 260 kg s⁻¹ for particles between 1–50 µm during the event observed on 2015 July 29. This estimation was based on the polygonal area defined by the ejected dust visible in OSIRIS images and the calculated dust mass flux. By assuming the outburst duration based on the cadence of OSIRIS images, they determined that the total dust mass ejected in this event was approximately 20–80 tons. Unfortunately, this event is not included in this study, as the first ROSINA measurement over the event’s source region occurred more than a day after the outburst observed by OSIRIS. Consequently, no direct comparison is possible; however, the estimated dust mass flux can be contextualized alongside other summer fireworks outbursts. Using the hemispherical gas distribution approach, we estimated a gas mass flux of 50–205 kg s⁻¹ resulting from H₂O, CO₂, and CO, while the cone cap distribution yielded a lower estimate of 5–20 kg s⁻¹. The gas mass flux derived from the hemispherical distribution aligns reasonably well with the dust mass flux reported by Vincent et al. (2016). Note, however, that the dust-to-gas mass ratio remains a topic of active debate (Choukroun et al. 2020) and is likely to change between the case of a fully active nucleus surface and the case of a much smaller crack area.

### 3.1.2 Other outbursts

In addition to the outbursts described by Vincent et al. (2016), several other outbursts were detected on comet 67P during the Rosetta mission. The detection of these events by the ROSINA instruments is discussed in Müller et al. (2024). These outbursts are also analysed here to assess potential subsurface gas cavities, similar to those inferred in the summer fireworks events (Section 3.1.1). Among these events, the outburst from 2015 November 7 (Noonan et al. 2021) stands out as the only CO₂-dominated event, with a significant increase in CO₂ production relative to H₂O (Table A3).

Assuming CO₂ as the primary driving force and using a cone cap outgassing distribution, we determined that a potential subsurface cavity would require a temperature of 143 K and an effective cavity opening radius of 0.75 m to match the surface’s tensile strength threshold. This corresponds to an equivalent half-spherical cavity with a radius between 15 and 24 metres. These values align with those derived for the summer fireworks events, although the effective cavity opening radius for this CO₂-dominated event is approximately half that of the summer fireworks cavities (Section 3.1.1 and Table A1). If a hemispherical outgassing distribution is used instead, an effective cavity opening radius of 2.5 m would be necessary, corresponding to an equivalent hemispherical cavity with a radius between 34 and 53 m (Table A4). These cavity dimensions are comparable to those associated with the summer fireworks events, with the effective cavity opening differing by a factor of 2.

The other outbursts that were outside the perihelion period, were primarily water-dominated and likely due to cliff collapses (Grün et al. 2016; Rinaldi et al. 2018; Müller et al. 2024). If modelled similarly to the summer fireworks events but for H₂O pockets instead of CO₂ pockets, temperatures of approximately 260 K would be required, which is significantly above the equilibrium surface
temperature (Table A5). Except for one event in 2016 July (Agarwal et al. 2017), the necessary cavity sizes for these events would need to be 1.6 to 2.3 times larger for both outgassing models (Tables A3 and A4). However, the high temperatures required in these pockets (≈260 K), well above the surface equilibrium temperature, suggest that H₂O pockets alone are an unlikely source of water-dominated outbursts, reinforcing prior propositions that these events were more likely due to cliff collapses. For completeness, if these outbursts were driven by CO₂ cavities instead, the associated temperatures and cavity sizes would be comparable to those in the summer fireworks evens. None the less, the strong H₂O outgassing observed in these events makes any interpretation of CO₂ as the dominant driver unlikely.

## 3.2 Dust velocity

The motion of dust particles in the comet's vicinity is primarily influenced by two forces: the gas drag accelerating the dust away from the comet's surface, and the gravitational pull of the comet's nucleus acting as an opposing force. The effect of solar radiation pressure is considered negligible within several hundred kilometres of the nucleus (Tenishev, Combi & Rubin 2011). For simplicity, dust grains are assumed to be spherical in shape (e.g. Tenishev et al. 2011; Marschall et al. 2016), although there is strong evidence, particularly for larger grains, that many are porous fluffy aggregates (Fulle et al. 2015; Schulz et al. 2015) and likely aspherical (Ivanovski et al. 2017). However, to first order, modelling the dust as spherical particles is justified since the key factors influencing dust trajectories are mass and cross-section and owing to the simplicity of our global model with, for example, a spherically shaped comet nucleus. As a result, the values of these parameters represent effective spheres. The motion of an individual dust grain, starting with zero initial velocity, is described by the following equation of motion (Gombosi, Nagy & Cravens 1986; Skorov & Rickman 1999; Molina, Moreno & Jiménez-Fernández 2008):

$$ \frac{4}{3}\pi a^3 \rho_d \frac{dv_d}{dt} = \pi a^2 \frac{C_D}{2} \rho_g(v_g - v_d) |v_g - v_d| - \frac{4}{3}\pi a^3 \rho_d \frac{GM_n}{r^2} \frac{r}{r}, \quad (9) $$

where G is the gravitational constant, M_n is the mass of the comet's nucleus, r is the dust grain's position relative to the nucleus, C_D is the drag coefficient, v_g is the bulk velocity of the gas in the coma, v_d is the dust grain velocity, a is the grain radius, and ρ_g and ρ_d are the densities of the surrounding gas and dust, respectively. Thus, the particles are simulated to run through the time-dependent cone-cap distribution to calculate ρ_g at each point of the simulation. For typical conditions in the coma, the drag coefficient is approximated as C_D = 2 (Gombosi et al. 1986; Grün et al. 1989). The uncertainties presented in Fig. 7 and Table 2 are calculated incrementally for each time-step of the ejection simulation. This is achieved using Gaussian error propagation applied to all relevant variables (v_g, v_d, r, ρ_g), incorporating a 30 per cent uncertainty in the initial thermal velocity of the gas molecules. This initial uncertainty arises from the uncertainty in the subsurface temperature. At each time-step, the propagated uncertainty of each variable is updated and subsequently used as an input for the next step. This iterative process is continued until the final dust grain velocity is determined.

### 3.2.1 Summer fireworks

For the dust involved in the summer fireworks outbursts, a density of 440 kg m⁻³ was assumed for all particle sizes. Fulle et al. (2015) estimated that very fluffy particles, ranging from 0.2 to 2.5 mm in

<figure><img src="image_7.png" /><figcaption>Figure 7. Average dust velocities as a function of radial distance from the comet centre. The dust grains are assumed to be spherical, with radii ranging from 10 nm to 1 mm as indicated in the legend, and a bulk dust density of 440 kg m⁻³. The shaded uncertainty regions represent the standard deviation of the summer fireworks outburst events plus the uncertainty from the error propagation as explained in Section 3.2.</figcaption></figure>

size, detected by the Grain Impact Analyser and Dust Accumulator (GIADA), have densities below 1 kg m⁻³.

At distances of around 10 nucleus radii, dust particles become fully decoupled from the gas phase in the coma (Tenishev et al. 2011). The terminal speed of dust grains varies significantly across different particle sizes.

Fig. 7 and Table 2 present the average dust velocity distributions for the summer fireworks outburst events, covering particle sizes from 10⁻⁸ to 10⁻³ m. The smallest grains, with radii of 10 nm, can reach velocities of 331 ± 14 m s⁻¹, while larger particles, with radii of 0.1 mm, achieve velocities of 12 ± 3 m s⁻¹.

Vincent et al. (2016) estimated a minimum dust ejection velocity of 13 m s⁻¹, based on the plume shape observed in successive images. Lin et al. (2016) suggested that larger grains (0.1 to 1 mm in diameter) likely dominate the curved jets, while smaller particles may populate straight jets. Rotundi et al. (2015) indicated that the optical scattering of dust is primarily driven by particles between 100 µm and several millimetres in size, though this was observed at a heliocentric distance of about 3.5 au and was expected to change closer to perihelion (Fulle et al. 2015). Later, Agarwal et al. (2017) found particles ranging from 20 to 200 µm in a 2016 July outburst event. In line with Vincent et al. (2016), our calculations suggest that the dust velocity for grains around 100 µm reaches 12 ± 3 m s⁻¹, with smaller grains achieving even higher velocities.

### 3.2.2 Other outbursts

Several instruments observed the event on 2016 July 3 (Agarwal et al. 2017), including the UV spectrometer ALICE, the Cometary Secondary Ion Mass Analyser (COSIMA), GIADA, and the Star Tracker B (STR-B) of Rosetta's attitude control system. These instruments provided different values for dust grain velocities, densities, and particle sizes. According to ALICE and STR-B, the fastest particles reached 25 ± 10 m s⁻¹, while the slowest particles travelled at 0.41 ± 0.05 m s⁻¹. ALICE data indicated the presence of sub-micron-sized water ice particles in the dust outburst. However, GIADA detected larger particles, likely in the range of hundreds of microns, with speeds below 3 m s⁻¹, suggesting compact particles with densities around 800 kg m⁻³ (Fulle et al. 2016). One of the particles could be more specifically determined to be 312 µm in size, 800 kg m⁻³ in density, and reached a velocity of 1.37 ± 0.08 m s⁻¹. COSIMA data suggested that particle tensile strengths were in the range of several
<table><caption>**Table 2.** Terminal dust velocity for all events considered in this work across different dust grain sizes having a bulk dust mass density of 440 kg m⁻³.</caption><thead><tr><th rowspan="2">Grain radius<br>(µm)</th><th rowspan="2">Summer fireworks</th><th colspan="7">Velocity (m s<sup>-1</sup>)</th></tr><tr><th>2015-05-23</th><th>2015-09-13</th><th>2015-09-14</th><th>2015-11-07</th><th>2016-01-06</th><th>2016-02-19</th><th>2016-07-03</th></tr></thead><tbody><tr><td>0.01</td><td>331 ± 14</td><td>203 ± 55</td><td>209 ± 61</td><td>214 ± 65</td><td>211 ± 62</td><td>172 ± 30</td><td>162 ± 25</td><td>79 ± 5</td></tr><tr><td>0.10</td><td>217 ± 24</td><td>130 ± 20</td><td>141 ± 23</td><td>150 ± 27</td><td>143 ± 24</td><td>90 ± 9</td><td>81 ± 8</td><td>30 ± 2</td></tr><tr><td>1.00</td><td>98 ± 17</td><td>58 ± 9</td><td>66 ± 11</td><td>73 ± 13</td><td>68 ± 11</td><td>36 ± 4</td><td>31 ± 4</td><td>10 ± 1</td></tr><tr><td>10.00</td><td>36 ± 7</td><td>21 ± 3</td><td>25 ± 4</td><td>28 ± 5</td><td>25 ± 4</td><td>12 ± 2</td><td>11 ± 2</td><td>3 ± 1</td></tr><tr><td>100.00</td><td>12 ± 3</td><td>7 ± 1</td><td>8 ± 1</td><td>9 ± 2</td><td>9 ± 2</td><td>4 ± 1</td><td>3 ± 1</td><td>1 ± 1</td></tr><tr><td>300.00</td><td>7 ± 2</td><td>4 ± 1</td><td>5 ± 1</td><td>5 ± 1</td><td>5 ± 1</td><td>2 ± 1</td><td>2 ± 1</td><td>0.6 ± 0.4</td></tr></tbody></table>

hundred Pa, leading to velocities between 2 and 5 m s⁻¹ for particle sizes slightly below 100 µm and particle densities of 250 kg m⁻³ (Hornung et al. 2016).

Our calculated dust grain velocities (Fig. 8, Panel a and Table 2) are consistent with those reported by Agarwal et al. (2017). For grains with densities of 800 kg m⁻³ and radii of 300 µm, we obtain velocities of 0.44 ± 0.04 m s⁻¹, about one third of the GIADA measurement. Grains with densities of 250 kg m⁻³ and radii slightly smaller than 100 µm, reach velocities of 1.61 ± 0.04 m s⁻¹, consistent with the COSIMA observations.

For the event on 2016 February 19, Grün et al. (2016) reported that the fastest particles on the order of 10 µm travelled at speeds of at least 25 m s⁻¹, while GIADA detected 100 µm particles moving at around 6 m s⁻¹. OSIRIS also observed even larger grains, up to several centimetres in size. Our estimates (Fig. 8, Panel b and Table 2) indicate that 100 µm particles reached speeds of 2.5 ± 0.3 m s⁻¹, about 2 to 3 times slower than the values reported by Grün et al. (2016). Similarly, smaller particles, approximately 10 µm in size, reached speeds over 8 m s⁻¹, also slower by a factor of 3.

Bockelée-Morvan et al. (2017) and Rinaldi et al. (2018) detected two significant outbursts on 2015 September 13 and 14 and retrieved dust velocities from spectral data collected by the Rosetta Visible InfraRed Thermal Imaging Spectrometer (VIRTIS). They reported velocities between 20 and 32 m s⁻¹ for outburst B (Rinaldi et al. 2018) and similar values for event F. From spectral analysis (Bockelée-Morvan et al. 2017), these values are likely for particles with sizes less than 0.3 µm. In addition, Rinaldi et al. (2018) computed particle speeds for dust sizes ranging from 0.1 to 10 µm, with velocities between 110 and 12 m s⁻¹, respectively (Fig. 9, black line). They explained the difference between these computed speeds and the velocities derived from spectral data by noting that velocities inferred from light curves represent the projected particle velocity, whereas the dust model calculates the total particle velocity, which can be considered an upper limit. Additionally, different particle shapes can lead to variations in velocity (Ivanovski et al. 2017; Rinaldi et al. 2018).

Our estimation of dust grain velocities (Fig. 8, Panel c and Table 2) closely aligns with the particle speeds simulated by Rinaldi et al. (2018). For dust particles with radii ranging from 0.1 to 10 µm, we determined velocities between 118 ± 16 m s⁻¹ and 18 ± 3 m s⁻¹, which correspond well to their reported values. As shown in Fig. 9, the comparison reveals a reasonably good agreement, with our results exceeding those of Rinaldi et al. (2018) by no more than a factor of 1.5. However, for particles with a size of 0.3 µm, we get 88 ± 9 m s⁻¹. This is about a factor 3 times more than the 20 and 32 m s⁻¹ Rinaldi et al. (2018) retrieved from the VIRTIS measurements. Therefore, in line with the statement by Rinaldi et al. (2018), the dust model calculates the total particle velocity, and our values may also be upper limits for these events as our data matches well with their dust model but overestimates the VIRTIS measurements.

<figure><img src="image_8.png" /><figcaption>**Figure 8.** Dust velocities as a function of radial distance from the comet centre for selected events outside the perihelion period and literature values for comparison. Panel (a) Event observed on 2016 July 3 (Agarwal et al. 2017). Panel (b) Event observed on 2016 February 19 (Grün et al. 2016). Panel (c) Event observed on 2015 September 13–14 (Rinaldi et al. 2018). For panel (a), two different dust bulk densities were used: 250 kg m⁻³ (dashed lines; Hornung et al. 2016) and 800 kg m⁻³ (solid lines; Fulle et al. 2016). In panels (b) and (c), a dust bulk density of 800 kg m⁻³ was assumed to compare the results to literature values. The literature values have all been determined where the particles reached their terminal velocity. They are, however, plotted over the whole distance for visibility considerations.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 9.** Speeds of spherically shaped particles as function of the particle radius estimated from the forces exerted by the gas flow measured with ROSINA/DFMS in comparison to the computed values by Rinaldi et al. (2018) for particles with initial temperature of 250 K. A bulk density of 800 kg m⁻³ has been used here.</figcaption></figure>

When compared with other studies, the particle speeds derived from ROSINA gas data are generally in good agreement with the velocities reported by Agarwal et al. (2017). However, our estimated dust velocities are lower by a factor of 2 to 3 for the event described by Grün et al. (2016), while exceeding the values for the grains analysed by Rinaldi et al. (2018) by about 50 per cent. Despite these discrepancies, the overall results for dust grain velocities across different outburst events are consistent with data obtained from other instruments. These findings validate the robustness of our approach, which incorporates the calculation of subsurface cavity gas densities, cavity opening radii, and gas pocket sizes.

# 4 DISCUSSION

This study provides a detailed analysis of gas-driven dust ejec-
tions observed on comet 67P/Churyumov–Gerasimenko during the
Rosetta mission, focusing on two distinct outburst categories: CO₂-
dominated perihelion events, such as the summer fireworks and
preferably arising from subsurface gas cavities, and water-dominated
outbursts outside the perihelion period. By modelling the outgassing
dynamics and dust velocities associated with these events, we
estimated the characteristics of subsurface gas cavities and the
velocity distribution of ejected particles.

## 4.1 Comparison between summer fireworks and other events

Our findings reveal significant differences between perihelion-period and off-perihelion outbursts. The summer fireworks events appear to be driven by CO₂ sublimation, with subsurface cavities estimated to have equivalent half-sphere radii between 15 and 62 m, depending on the assumed outgassing distribution model. Subsurface temperatures at cavity rupture were approximately 137 K, consistent with CO₂ sublimation thresholds. The broad range of cavity radii highlights the sensitivity of results to modelling assumptions, with hemispherical outgassing distributions likely providing upper bounds due to the localized nature of these events, as indicated by OSIRIS images. However, the precise timing of peak gas production remains uncertain and measurements by ROSINA/DFMS on Rosetta were taken after some level of gas dissipation and lateral diffusion. Therefore, the actual cavity radius likely falls within this range, justifying the rough approximations made in modelling and accounting for the variability in radius estimates.

Dust velocities during the summer fireworks events also varied widely, with smaller grains (<10 nm) reaching high velocities of up to 330 m s⁻¹, while larger particles (100 µm to 1 mm) exhibited velocities closer to 12 m s⁻¹. These values align with OSIRIS observations (Vincent et al. 2016), supporting the robustness of our modelling framework.

In contrast, water-driven outbursts outside the perihelion period were likely mechanical in origin, possibly resulting from cliff collapse or other collapsing structures within the nucleus. The higher temperatures required for H₂O sublimation in a subsurface cavity – around 260 K and well above the equilibrium surface temperature – suggest that these events involve different mechanisms compared to CO₂-driven outbursts. Dust velocities for these events were generally lower, ranging from a factor of 0.3 to 0.1 below those observed during summer fireworks. These results support the hypothesis of distinct physical processes underlying the two types of outbursts.

## 4.2 Implications of subsurface cavity formation

The observed cavity sizes and dust ejection velocities provide critical insights into subsurface processes on 67P. CO₂-dominated cavities are consistent with models of gas reservoirs accumulating beneath an impermeable surface layer or pocket. Similar processes may occur on other comets with significant volatile content, such as 29P/Schwassmann-Wachmann (Miles et al. 2016; Lisse et al. 2022) or 17P/Holmes (Lin et al. 2009; Gronkowski & Sacharczuk 2010), both of which show massive outbursts associated with a strong increase of CO and other highly volatile components. These composition signatures from other comets suggest that volatile-driven outbursts may be a common phenomenon.

Morphological heterogeneity on 67P further emphasizes the dynamical nature of cometary surfaces. Transitions between smooth and rough terrains reflect substantial dust transport (Vincent et al. 2021), with volatile sublimation driving intra- and inter-regional redistribution of sediment (Barrington et al. 2023). Localized zones of erosion and accumulation, as documented by Jindal et al. (2024), highlight the uneven impact of these processes, even within small areas.

Several depressions identified in smooth and harder terrains by Vincent et al. (2021) provide intriguing clues. These depressions, typically 10 m in diameter and 1–2 m in depth, are smaller than the modelled subsurface cavities but are comparable to the estimated cavity openings. In the Imhotep region, depressions found in a basin associated with outbursts suggest the presence of volatile reservoirs just below the surface. While depressions in other areas lack direct links to outbursts, some events may have been missed due to insufficient observations (Vincent et al. 2021).

Lamy et al. (2024) further investigated potential icy cavities using OSIRIS image anaglyphs, identifying bright patches 15–30 m across with high reflectance and visible spectral slopes indicative of subsurface water ice. Infrared absorption features supported this interpretation, with the ice potentially being exposed due to outbursts that removed the overlying dust layers. Although the nature of the detected ice – whether water or more volatile ices – remains debatable, Lamy et al. (2024) proposed that these patches represent pristine mixtures of water ice and refractory grains at cavity bottoms. High-resolution photogrammetry revealed cavity depths of 20–47 m and suggested lifetimes of up to two years, with potential links to jets or outbursts. These cavity dimensions align well with the modelled subsurface cavities for CO₂-driven events, implying that they may be remnants of outburst-related features.
## 4.3 Model refinement and future outlook

While the current model offers valuable insights into gas-driven outbursts, several refinements could enhance its accuracy. Incorporating irregular particle geometries and tensile strengths, as well as refining assumptions about the comet's shape and gravity field, would improve dust velocity predictions. Higher-resolution data on subsurface volatiles, obtained through future missions using advanced radar, microwave, or thermal imaging, could further constrain subsurface cavity characteristics. Notably, exploring the THz frequency region, offering deeper penetration capabilities (Stöckli et al. 2025), will enable direct measurements of the volatile content and the structure of the gas pockets. Such advancements could improve our understanding of subsurface processes that drive these dynamic outbursts.

Fundamental questions persist regarding the formation and evolution of subsurface cavities. Investigating how these cavities form – whether through gas build-up, thermal gradients, or structural weaknesses – would refine our understanding of cometary activity and enhance predictive models. Environmental influences, such as Solar radiation and perihelion heating, are particularly relevant to cavity development. Current hypotheses on cavity formation may involve repetitive CO₂ frost cycles, which have been observed on comet 67P and shown to operate on seasonal time-scales (Filacchione et al. 2016; Rubin et al. 2023). This process could explain why areas exhibiting outbursts are often active in CO₂ well in advance (Läuter et al. 2020). Frost cycles could lead to the closure of cavities, while removal of the dust cover during an outburst may account for their long-term gas release after the dust plume is gone (Läuter et al. 2020).

Another potential mechanism, which might even be connected to the frost cycle, involves the formation of a CO₂ ice layer coating the porous internal structure. CO₂ ice layers are observed in laboratory experiments under Martian conditions (Portyankina et al. 2019). Although the conditions on comets differ significantly – particularly in terms of pressure, temperature, porosity, composition, and scale – a similar process might occur. In this scenario, CO₂ forms an impermeable layer over porous material, creating a structure with considerable tensile strength capable of confining a high-pressure gas pocket. The interplay between CO₂ frost cycles and ice layer formation – or a combination of the two – should be further explored to validate these potential mechanisms of subsurface cavity formation.

## 5 CONCLUSION

This study advances our understanding of gas-driven outbursts on comet 67P/Churyumov-Gerasimenko, strengthening our previous conclusions about distinct mechanisms for perihelion-period and off-perihelion events. CO₂-dominated outbursts, such as the summer fireworks, are driven by sublimation of subsurface volatiles, resulting in the formation and rupture of cavities with equivalent half-sphere radii between 15 and 62 m. In contrast, water-dominated events outside perihelion are likely mechanical in origin, involving processes such as cliff collapse or localized heating.

The differences in dust ejection velocities and cavity sizes between the two event types reflect the diverse physical processes shaping cometary nuclei. CO₂-driven outbursts provide insights into the role of volatile reservoirs, while water-driven events highlight the importance of structural integrity and mechanical collapse. These findings underscore the complex interplay of thermal, chemical, and mechanical processes governing cometary activity.

Future research should prioritize refining models of subsurface cavity formation, leveraging data from advanced instruments and

laboratory studies. By integrating high-resolution observations from future missions with comparative studies across comet populations, we can further unravel the mechanisms driving cometary outbursts, surface evolution, and the diversity of cometary activity across the Solar System.

## ACKNOWLEDGEMENTS

We gratefully acknowledge the work of the many engineers, technicians, and scientists involved in the Rosetta mission and in the ROSINA (Rosetta Orbiter Spectrometer for Ion and Neutral Analysis) instrument in particular. Without their contributions, ROSINA would not have produced such outstanding results. Rosetta is a European Space Agency (ESA) mission with contributions from its member states and the National Aeronautics and Space Administration (NASA). Work at the University of Bern was funded by the State of Bern and the Swiss National Science Foundation (200020_207312). SFW acknowledges the financial support of the SNSF Eccellenza Professorial Fellowship (PCEFP2_181150).

## DATA AVAILABILITY

All Rosetta/ROSINA data are available through the NASA Planetary Data System (PDS) and ESA's Planetary Science Archive (PSA).

## REFERENCES

* Agarwal J. et al., 2017, MNRAS, 469, s606

* Balsiger H. et al., 2007, Space Sci. Rev., 128, 745

* Barrington M. N., Birch S. P. D., Jindal A., Hayes A. G., Corlies P., Vincent J. B., 2023, J. Geophys. Res. (Planets), 128, e2022JE007723

* Basilevsky A. T., Krasil'nikov S. S., Shiryaev A. A., Mall U., Keller H. U., Skorov Y. V., Mottola S., Hviid S. F., 2016, Sol. Syst. Res., 50, 225

* Bockelée-Morvan D. et al., 2017, MNRAS, 469, S443

* Bockelée-Morvan D. et al., 2022, A&A, 664, A95

* Calmonte U. et al., 2016, MNRAS, 462, S253

* Choukroun M. et al., 2020, Space Sci. Rev., 216, 44

* Combi M. R., Tenishev V. M., Rubin M., Fougere N., Gombosi T. I., 2012, ApJ, 749, 29

* De Keyser J. et al., 2019, Int. J. Mass Spect., 446, 116232

* Feldman P. D. et al., 2016, ApJ, 825, L8

* Filacchione G. et al., 2016, Science, 354, 1563

* Fray N., Schmitt B., 2009, Planet. Space Sci., 57, 2053

* Fulle M. et al., 2015, ApJ, 802, L12

* Fulle M., Altobelli N., Buratti B., Choukroun M., Fulchignoni M., Grün E., Taylor M. G. G. T., Weissman P., 2016, MNRAS, 462, S2

* Gombosi T. I., Nagy A. F., Cravens T. E., 1986, Rev. Geophys., 24, 667

* Gronkowski P., Sacharczuk Z., 2010, MNRAS, 408, 1207

* Groussin O. et al., 2015, A&A, 583, A32

* Grün E. et al., 1989, Adv. Space Res., 9, 133

* Grün E. et al., 2016, MNRAS, 462, S220

* Gulkis S. et al., 2015, Science, 347, aaa0709

* Hornung K. et al., 2016, Planet. Space Sci., 133, 63

* Hughes D. W., 1975, Q. J. R. Astron. Soc., 16, 410

* Hughes D. W., 1991, in Newburn R. L. J., Neugebauer M., Rahe J., eds, Astrophysics and Space Science Library, Vol. 167, IAU Colloq. 116: Comets in the post-Halley era. Kluwer Academic Publishers, Dordrecht, The Netherlands, p. 825

* Ivanovski S. L. et al., 2017, MNRAS, 469, S774

* Jindal A. S. et al., 2024, J. Geophys. Res. (Planets), 129, e2023JE008089

* Keller H. U. et al., 2015, A&A, 583, A34

* Lamy P., Faury G., Romeuf D., Groussin O., 2024, MNRAS, 531, 2494

* Läuter M., Kramer T., Rubin M., Altwegg K., 2020, MNRAS, 498, 3995

* Le Roy L. et al., 2015, A&A, 583, A1

* Lin Z.-Y., Lin C.-S., Ip W.-H., Lara L. M., 2009, AJ, 138, 625
Lin Z. Y. et al., 2016, A&A, 588, L3

Lisse C. M. et al., 2022, Planet. Sci. J., 3, 251

Marschall R. et al., 2016, A&A, 589, A90

Marschall R., Liao Y., Thomas N., Wu J.-S., 2020, Icarus, 346, 113742

Miles R., Faillace G. A., Mottola S., Raab H., Roche P., Soulier J.-F., Watkins A., 2016, Icarus, 272, 327

Molina A., Moreno F., Jiménez-Fernández F. J., 2008, Earth Moon and Planets, 102, 521

Mousis O. et al., 2015, ApJ, 814, L5

Müller D. R. et al., 2024, MNRAS, 529, 2763

Noonan J. W. et al., 2021, AJ, 162, 4

Pajola M. et al., 2017, Nat. Astron., 1, 0092

Portyankina G., Merrison J., Iversen J. J., Yoldi Z., Hansen C. J., Aye K. M., Pommerol A., Thomas N., 2019, Icarus, 322, 210

Prialnik D., Bar-Nun A., 1992, A&A, 258, L9

Prialnik D., A'Hearn M. F., Meech K. J., 2008, MNRAS, 388, L20

Rinaldi G. et al., 2018, MNRAS, 481, 1235

Rotundi A. et al., 2015, Science, 347, aaa3905

Rubin M. et al., 2019, MNRAS, 489, 594

Rubin M. et al., 2023, MNRAS, 526, 4209

Schulz R. et al., 2015, Nature, 518, 216

Sharipov F., Kozak D. V., 2009, J. Vac. Sci. Technol. A, 27, 479

Skorov Y. V., Rickman H., 1999, Planet. Space Sci., 47, 935

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

Stöckli L. L., Brändli M., Piazza D., Ottersberg R., Pommerol A., Murk A., Thomas N., 2025, Rev. Sci. Instrum., (under review)

Tenishev V., Combi M. R., Rubin M., 2011, ApJ, 732, 104

Vincent J.-B. et al., 2015, Nature, 523, 63

Vincent J. B. et al., 2016, MNRAS, 462, S184

Vincent J.-B., Kruk S., Fanara L., Birch S., Jindal A., 2021, European Planetary Science Congress 2021. EPSC2021–525

APPENDIX A: EXTENDED DATA

<table><caption>**Table A1.** Excess gas production rates due to outburst events for H₂O and CO₂, cavity temperature, cavity pressure, and equivalent half-sphere cavity radius for CO₂ for each event of the summer fireworks (Vincent et al. 2016) considered in this work. The event IDs correspond to the values from Vincent et al. (2016). The excess gas production rates due to the outburst events are calculated by assuming a hemispherical cone cap distribution with a half-angle of 25°. The temperature, pressure, and equivalent half-sphere cavity radius are given for an effective cavity opening radius of 1.4 m matching best with the tensile strength of the surface layers of comet 67P. For the equivalent half-sphere cavity radius the two extreme τ-values are used in the calculations.</caption><thead><tr><th>ID</th><th>Event date (UTC)</th><th>r<sub>cgsc</sub><br>(km)</th><th>Lat<br>(°)</th><th>Lon<br>(°)</th><th>Gas prod. rate<br>(kg s<sup>-1</sup>)</th><th>H<sub>2</sub>O<br>CO<sub>2</sub></th><th>CO</th><th>T<sub>in</sub>, CO<sub>2</sub> (K)<br>r<sub>crack</sub> = 1.4 m</th><th>p<sub>in</sub>, CO<sub>2</sub> (Pa)<br>r<sub>crack</sub> = 1.4 m</th><th>r<sub>cavity</sub>, CO<sub>2</sub> (m)<br>τ = 40 s</th><th>r<sub>cavity</sub>, CO<sub>2</sub> (m)<br>τ = 260 s</th></tr></thead><tbody><tr><td>3</td><td>2015-07-26T20:22:42</td><td>168</td><td>-36</td><td>75</td><td>6.8</td><td>0.9</td><td>0.1</td><td>131</td><td>37</td><td>15.2</td><td>28.4</td></tr><tr><td>4</td><td>2015-07-27T00:14:29</td><td>168</td><td>-31</td><td>333</td><td>5.4</td><td>1.2</td><td>0.2</td><td>132</td><td>50</td><td>15.2</td><td>28.4</td></tr><tr><td>5</td><td>2015-07-28T05:23:43</td><td>181</td><td>-4</td><td>264</td><td>5.1</td><td>2.2</td><td>0.3</td><td>136</td><td>97</td><td>15.3</td><td>28.6</td></tr><tr><td>7</td><td>2015-08-01T10:53:15</td><td>214</td><td>-12</td><td>196</td><td>10.2</td><td>2.1</td><td>0.3</td><td>136</td><td>90</td><td>15.3</td><td>28.6</td></tr><tr><td>8</td><td>2015-08-01T15:44:50</td><td>211</td><td>-28</td><td>34</td><td>3.8</td><td>2.4</td><td>0.2</td><td>137</td><td>105</td><td>15.3</td><td>28.6</td></tr><tr><td>17</td><td>2015-08-22T23:46:21</td><td>334</td><td>-25</td><td>316</td><td>7.9</td><td>1.8</td><td>0.3</td><td>135</td><td>76</td><td>15.3</td><td>28.5</td></tr><tr><td>20</td><td>2015-08-26T07:51:04</td><td>417</td><td>-41</td><td>42</td><td>7.3</td><td>7.9</td><td>0.6</td><td>144</td><td>352</td><td>15.5</td><td>28.9</td></tr><tr><td>23</td><td>2015-08-28T10:10:57</td><td>410</td><td>-31</td><td>229</td><td>6.8</td><td>7.6</td><td>1.1</td><td>144</td><td>341</td><td>15.5</td><td>28.8</td></tr></tbody></table>

<table><caption>**Table A2.** The same as in Table A1 but the excess gas production rates are calculated by assuming a hemispherical distribution. Therefore, the temperature, pressure, and equivalent half-sphere cavity radius are given for an effective cavity opening radius of 4.5 m matching best with the tensile strength of the surface layers of comet 67P.</caption><thead><tr><th>ID</th><th>Event date (UTC)</th><th>r<sub>cgsc</sub><br>(km)</th><th>Lat<br>(°)</th><th>Lon<br>(°)</th><th>Gas prod. rate<br>(kg s<sup>-1</sup>)</th><th>H<sub>2</sub>O<br>CO<sub>2</sub></th><th>CO</th><th>T<sub>in</sub>, CO<sub>2</sub> (K)<br>r<sub>crack</sub> = 4.5 m</th><th>p<sub>in</sub>, CO<sub>2</sub> (Pa)<br>r<sub>crack</sub> = 4.5 m</th><th>r<sub>cavity</sub>, CO<sub>2</sub> (m)<br>τ = 40 s</th><th>r<sub>cavity</sub>, CO<sub>2</sub> (m)<br>τ = 260 s</th></tr></thead><tbody><tr><td>3</td><td>2015-07-26T20:22:42</td><td>168</td><td>-36</td><td>75</td><td>72.2</td><td>9.2</td><td>0.8</td><td>131</td><td>38</td><td>33.1</td><td>61.8</td></tr><tr><td>4</td><td>2015-07-27T00:14:29</td><td>168</td><td>-31</td><td>333</td><td>57.4</td><td>12.4</td><td>1.8</td><td>133</td><td>51</td><td>33.2</td><td>62.0</td></tr><tr><td>5</td><td>2015-07-28T05:23:43</td><td>181</td><td>-4</td><td>264</td><td>54.3</td><td>23.7</td><td>2.7</td><td>136</td><td>100</td><td>33.4</td><td>62.3</td></tr><tr><td>7</td><td>2015-08-01T10:53:15</td><td>214</td><td>-12</td><td>196</td><td>109.2</td><td>22.2</td><td>2.7</td><td>136</td><td>93</td><td>33.3</td><td>62.2</td></tr><tr><td>8</td><td>2015-08-01T15:44:50</td><td>211</td><td>-28</td><td>34</td><td>41.0</td><td>25.7</td><td>1.8</td><td>137</td><td>108</td><td>33.4</td><td>62.3</td></tr><tr><td>17</td><td>2015-08-22T23:46:21</td><td>334</td><td>-25</td><td>316</td><td>83.9</td><td>18.8</td><td>3.5</td><td>135</td><td>79</td><td>33.3</td><td>62.2</td></tr><tr><td>20</td><td>2015-08-26T07:51:04</td><td>417</td><td>-41</td><td>42</td><td>78.4</td><td>83.8</td><td>6.5</td><td>144</td><td>363</td><td>33.7</td><td>62.9</td></tr><tr><td>23</td><td>2015-08-28T10:10:57</td><td>410</td><td>-31</td><td>229</td><td>72.1</td><td>81.3</td><td>11.6</td><td>144</td><td>352</td><td>33.7</td><td>62.8</td></tr></tbody></table>
<table><caption>**Table A3.** The same as in Table A1 but for the events considered in this work that are not discussed in Vincent et al. (2016). The excess gas production rates are calculated assuming a hemispherical cone cap distribution with a half-angle of 25°. The temperature, effective cavity opening radius, and equivalent half-sphere cavity radius are shown for the most enhanced gas and have been calculated such that they are matching best with the tensile strength of the surface layers of comet 67P.</caption><thead><tr><th rowspan="2">Event date (UTC)</th><th rowspan="2">Reference</th><th rowspan="2">r<sub>cgsc</sub><br>(km)</th><th colspan="3">Gas prod. rate (kg s<sup>-1</sup>)</th><th rowspan="2">Most enhanced gas</th><th rowspan="2">T<sub>in</sub><br>(K)</th><th rowspan="2">r<sub>crack</sub><br>(m)</th><th colspan="2">r<sub>cavity</sub> (m)</th></tr><tr><th>H<sub>2</sub>O</th><th>CO<sub>2</sub></th><th>CO</th><th>τ<sub>min</sub></th><th>τ<sub>max</sub></th></tr></thead><tbody><tr><td>2015-05-23T12:42</td><td>Feldman et al. (2016)</td><td>143</td><td>3.8</td><td>0.2</td><td>0.1</td><td>H<sub>2</sub>O</td><td>259</td><td>2.0</td><td>74.6</td><td>84.9</td></tr><tr><td>2015-09-13T13:36</td><td>Rinaldi et al. (2018)</td><td>314</td><td>4.3</td><td>0.2</td><td>0.1</td><td>H<sub>2</sub>O</td><td>260</td><td>2.0</td><td>43.0</td><td>59.1</td></tr><tr><td>2015-09-14T18:47</td><td>Rinaldi et al. (2018)</td><td>316</td><td>6.0</td><td>0.3</td><td>0.1</td><td>H<sub>2</sub>O</td><td>264</td><td>2.0</td><td>34.7</td><td>53.8</td></tr><tr><td>2015-11-07T17:32</td><td>Noonan et al. (2021)</td><td>229</td><td>5.1</td><td>1.9</td><td>0.1</td><td>CO<sub>2</sub></td><td>143</td><td>0.75</td><td>15.1</td><td>23.9</td></tr><tr><td>2016-01-06T12:00</td><td>Agarwal et al. (2017)</td><td>83.9</td><td>1.0</td><td>0.6</td><td>0.02</td><td>H<sub>2</sub>O</td><td>259</td><td>1.0</td><td>22.8</td><td>38.8</td></tr><tr><td>2016-02-19T09:40</td><td>Grün et al. (2016)</td><td>34.5</td><td>0.8</td><td>0.2</td><td>0.03</td><td>H<sub>2</sub>O</td><td>263</td><td>0.75</td><td>35.5</td><td>41.8</td></tr><tr><td>2016-07-03T07:30</td><td>Agarwal et al. (2017)</td><td>8.5</td><td>0.02</td><td>0.15</td><td>0.01</td><td>H<sub>2</sub>O</td><td>260</td><td>0.15</td><td>6.4</td><td>11.0</td></tr></tbody></table>

<table><caption>**Table A4.** The same as in Table A3 but the excess gas production rates are calculated by assuming a hemispherical distribution. In addition, the estimated dust mass-loss from the corresponding literature reference is given. It is given here as comparison to the calculation assuming a hemispherical gas distribution because all authors used a hemispherical gas and dust distribution except for Grün et al. (2016) who used a solid angle of $\Omega_d = 1$ sr. Here, their value has been normalized to the solid angle of $\Omega_d = 2\pi$ sr which is used by all other reference publications.</caption><thead><tr><th rowspan="2">Event date (UTC)</th><th rowspan="2">Reference</th><th rowspan="2">r<sub>cgsc</sub><br>(km)</th><th rowspan="2">Dust mass-loss<br>(tons)</th><th colspan="3">Gas prod. rate (kg s<sup>-1</sup>)</th><th rowspan="2">Most enhanced gas</th><th rowspan="2">T<sub>in</sub><br>(K)</th><th rowspan="2">r<sub>crack</sub><br>(m)</th><th colspan="2">r<sub>cavity</sub> (m)</th></tr><tr><th>H<sub>2</sub>O</th><th>CO<sub>2</sub></th><th>CO</th><th>τ<sub>min</sub></th><th>τ<sub>max</sub></th></tr></thead><tbody><tr><td>2015-05-23T12:42</td><td>Feldman et al. (2016)</td><td>143</td><td>–</td><td>40.2</td><td>1.6</td><td>0.6</td><td>H<sub>2</sub>O</td><td>265</td><td>5.0</td><td>137.9</td><td>157.0</td></tr><tr><td>2015-09-13T13:36</td><td>Rinaldi et al. (2018)</td><td>314</td><td>70–230</td><td>46.3</td><td>2.1</td><td>0.9</td><td>H<sub>2</sub>O</td><td>264</td><td>5.5</td><td>84.5</td><td>116.2</td></tr><tr><td>2015-09-14T18:47</td><td>Rinaldi et al. (2018)</td><td>316</td><td>120–380<sup>a</sup></td><td>63.8</td><td>3.6</td><td>0.7</td><td>H<sub>2</sub>O</td><td>264</td><td>6.5</td><td>76.2</td><td>118.0</td></tr><tr><td>2015-11-07T17:32</td><td>Noonan et al. (2021)</td><td>229</td><td>–</td><td>54.1</td><td>20.3</td><td>1.1</td><td>CO<sub>2</sub></td><td>143</td><td>2.5</td><td>33.7</td><td>53.4</td></tr><tr><td>2016-01-06T12:00</td><td>Agarwal et al. (2017)</td><td>83.9</td><td>–</td><td>10.1</td><td>5.9</td><td>0.3</td><td>H<sub>2</sub>O</td><td>260</td><td>3.0</td><td>47.5</td><td>80.7</td></tr><tr><td>2016-02-19T09:40</td><td>Grün et al. (2016)</td><td>34.5</td><td>10<sup>b</sup></td><td>8.4</td><td>2.3</td><td>0.3</td><td>H<sub>2</sub>O</td><td>263</td><td>2.5</td><td>79.2</td><td>93.2</td></tr><tr><td>2016-07-03T07:30</td><td>Agarwal et al. (2017)</td><td>8.5</td><td>6.5–118<sup>c</sup></td><td>0.3</td><td>1.6</td><td>0.1</td><td>H<sub>2</sub>O</td><td>260</td><td>0.5</td><td>14.4</td><td>24.4</td></tr></tbody></table>

<table><caption>**Table A5.** Equilibrium temperatures for water ice sublimation on a flat surface for all considered outburst events according to equation (6).</caption><thead><tr><th>Event date (UTC)</th><th>Reference</th><th>Hel. Distance (au)</th><th>T<sub>eq</sub>, H<sub>2</sub>O (K)</th></tr></thead><tbody><tr><td>2015-07-26T20:22</td><td>Vincent et al. (2016)</td><td>1.26</td><td>202</td></tr><tr><td>2015-07-27T00:14</td><td>Vincent et al. (2016)</td><td>1.26</td><td>202</td></tr><tr><td>2015-07-28T05:23</td><td>Vincent et al. (2016)</td><td>1.26</td><td>202</td></tr><tr><td>2015-08-01T10:53</td><td>Vincent et al. (2016)</td><td>1.25</td><td>202</td></tr><tr><td>2015-08-01T15:44</td><td>Vincent et al. (2016)</td><td>1.25</td><td>202</td></tr><tr><td>2015-08-22T23:46</td><td>Vincent et al. (2016)</td><td>1.25</td><td>202</td></tr><tr><td>2015-08-26T07:51</td><td>Vincent et al. (2016)</td><td>1.25</td><td>202</td></tr><tr><td>2015-08-28T10:10</td><td>Vincent et al. (2016)</td><td>1.26</td><td>202</td></tr><tr><td>2015-05-23T12:42</td><td>Feldman et al. (2016)</td><td>1.58</td><td>198</td></tr><tr><td>2015-09-13T13:36</td><td>Rinaldi et al. (2018)</td><td>1.30</td><td>201</td></tr><tr><td>2015-09-14T18:47</td><td>Rinaldi et al. (2018)</td><td>1.31</td><td>201</td></tr><tr><td>2015-11-07T17:32</td><td>Noonan et al. (2021)</td><td>1.61</td><td>198</td></tr><tr><td>2016-01-06T12:00</td><td>Agarwal et al. (2017)</td><td>2.06</td><td>194</td></tr><tr><td>2016-02-19T09:40</td><td>Grün et al. (2016)</td><td>2.40</td><td>192</td></tr><tr><td>2016-07-03T07:30</td><td>Agarwal et al. (2017)</td><td>3.32</td><td>186</td></tr></tbody></table>

<footnote><sup>a</sup>The given value is calculated considering the outburst inside the FOV. Using an additional fraction of 30 per cent of the outburst outside the FOV gives a dust mass-loss between 150 and 500 tons (Rinaldi et al. 2018). <sup>b</sup>1.6 multiplied by $2\pi$ to normalize to a hemispherical outgassing distribution. <sup>c</sup>Large spread due to the use of extreme values including uncertainty for both the dust production rate and the outburst life-time.</footnote>
<table><caption>**Table A6.** Exponential e-folding time $\tau$ values from equation (7) for all events considered in this study.</caption><thead><tr><td>Event date (UTC)</td><td>Reference</td><td>τ<sub>min</sub> (s)</td><td>τ<sub>max</sub> (s)</td></tr></thead><tbody><tr><td>2015-05-23T12:42</td><td>Feldman et al. (2016)</td><td>1050</td><td>1550</td></tr><tr><td>Summer fireworks</td><td>Vincent et al. (2016)</td><td>40</td><td>260</td></tr><tr><td>2015-09-13T13:36</td><td>Rinaldi et al. (2018)</td><td>200</td><td>520</td></tr><tr><td>2015-09-14T18:47</td><td>Rinaldi et al. (2018)</td><td>100</td><td>390</td></tr><tr><td>2015-11-07T17:32</td><td>Noonan et al. (2021)</td><td>130</td><td>520</td></tr><tr><td>2016-01-06T12:00</td><td>Agarwal et al. (2017)</td><td>120</td><td>590</td></tr><tr><td>2016-02-19T09:40</td><td>Grün et al. (2016)</td><td>800</td><td>1300</td></tr><tr><td>2016-07-03T07:30</td><td>Agarwal et al. (2017)</td><td>120</td><td>590</td></tr></tbody></table>

This paper has been typeset from a *TEX/LATEX* file prepared by the author.