© ESO, 2016

1. Introduction

Although a little consensus exists on what constitutes a cometary outburst, it is phenomenologically defined as a sudden and strong emission of gas and dust from the nucleus. Observationally, the outbursts are usually recorded in the visible wavelengths and recognized with an increase in brightness within a relatively short time with subsequent fading to the pre-outburst background level (within days/hours depending on the observing context). The different outbursts of different comets are recorded to be phenomenologically distinct, from the large-scale dust jets (e.g. Comets Wild 2; Sekanina et al. 2004 and 9P/Tempel 1; Farnham et al. 2007), the periodic outburst events of P/Schwasmannn-Wachmann 1 (Hughes 1990; Trigo-Rodríguez et al. 2008), to the massive outburst of 17P/Holmes (Montalto et al. 2008). Similarly, explaining the outbursts in general terms proves to be difficult; however, the immediate cause(s) may be roughly separated into external sources (impacts), physicochemical exothermic processes inside the nucleus (phase changes, etc.), or the potential energy stored in the internal structure of the material (Gronkowski & Sacharczuk 2010). A discussion contrasting some of the mechanisms belonging to these general categories is provided in Sect. 4.

The Jupiter Family comet 67P/Churyumov-Gerasimenko (hereafter 67P) is the target of the ESA Rosetta comet rendezvous mission whose goal was to capture the onset and evolution of the comet activity during its orbit in unprecedented spatial and temporal resolution. Currently, 67P is characterized by an eccentricity e = 0.64102, an inclination i = 7.0405°, a perihelion distance q = 1.2432 AU, and a semi-major axis of 3.4630 AU resulting in a period P = 6.44 yr. The nucleus rotates with a slightly decreasing period of about 12 h around an axis that is inclined by 52° relative to the comet orbit plane (obliquity) (Keller et al. 2015). As a consequence, mainly the northern hemisphere was illuminated when Rosetta arrived at the comet and the nucleus could be resolved by the science cameras OSIRIS (Keller et al. 2007) in July 2014. The nucleus is characterized by a bi-lobate shape where the bigger and the smaller lobes are connected by a relatively small contact area forming a neck. The overall dimensions along its principal axes are 4.3 × 2.6 × 2.1 km³ resulting in a volume of 18.8 km³ (Preusker et al. 2015). The porosity of this comet is estimated to be approximately 70 to 75% depending upon the exact mean density value and dust-to-ice ratio. The distribution between micro and macro density remains open.

The neck forms a cavity oriented towards the north pole located on the big lobe. Above this region, named Hapi (see Thomas et al. 2015 and El-Maarry et al. 2015b for the designation of the regions), the first signs of cometary dust activity were observed while the comet’s heliocentric distance was still beyond 3.5 AU (Vincent et al. 2015). The morphology of the northern hemisphere is rather diverse. The terrains often look smooth, sometimes aeolian-like landforms are visible suggestive of “air fall” (Thomas et al. 2015). Loose refractory material (dust) covers most near horizontal areas (relative to local gravity) forming terraces and alcoves. Slopes steeper than the angle of repose, often extended cliffs, show consolidated material quite often scarred by fractures extending up to hundreds of meters. Terraces and pits are numerous, as are boulders (Pajola et al. 2015) often concentrated in talus field. Many of the pits show dust activity originating from their walls or bottoms and rims (Vincent et al. 2015).

The first signs of dust activity reported by the Rosetta mission were associated with an outburst-like event on April 30, 2014, occurring in the region of the north pole (Tubiana et al. 2015). The gas production for this event could not be determined at the time, although the mm/submm wavelength instrument (MIRO) detected water in the coma near illuminated surfaces (Gulkis et al. 2015; Lee et al. 2015) at heliocentric distance of 3.9 AU. Water outgassing was only detected once at a larger heliocentric distance (de Val-Borro et al. 2014). Later H₂O, CO₂, and CO were detected by the ROSINA mass spectrometer (Altwegg et al. 2015; Hässig et al. 2015). The most recent measurements also indicate a significant spatial/temporal variability of the relative abundances of volatiles in the coma (Luspay-Kuti et al. 2015). Nevertheless, the relative abundance of CO/CO₂, CO₂/H₂O, and CO/H₂O are not fully consistent among the different instruments (Morse et al. 2015).

Although outbursts were detected during the arrival stage (Tubiana et al. 2015) and thought to originate from collapsing sink holes forming the numerous surface pits (Vincent et al. 2015), no outburst-like event was observed until about 30 days prior to the perihelion passage based on the NAVCAM images and published literature. On the other hand, the ubiquitous dust activity in the form of periodic occurrence of jet-like shapes were studied extensively on 67P. First, by Earth observers (Lara et al. 2015) and by Rosetta (Vincent et al. 2016; Lin et al. 2016). These studies concluded that although influenced by the local topography, these dust jets are primarily driven by the propagation of the thermal heat wave into subsurface ice layer. This mechanism is also reported to account for such dust jets originating from the shadow regions right after a local sunset (Shi et al. 2016).

In this paper we describe a new mechanism for short-lived cometary outbursts (<1 h) observed by the cameras on Rosetta toward the perihelion approach. The predicted duration of these outbursts is on the order of ≈30 min according to our model. The NAVCAM images are sampled from 1 to 5 h in the time near the perihelion approach and cannot directly confirm this prediction. Nevertheless, the OSIRIS results do observationally support our prediction, at least in the July−August 2015 time frame (A’Hearn et al. 2016). The strength of the outburst from the NAVCAM images is qualitatively estimated by a) the released dust cone being clearly visible on a linear brightness scale for which the ubiquitous dust-jet activity is hardly perceivable; and b) the lowest part of the outburst jet is brighter than the reflected light from the nucleus, indicating large optical thickness. The currently suggested mechanisms of cometary outburst cannot properly explain such short-lived outbursts, as is discussed below.

The paper is structured as follows. In Sect. 2 we provide a more detailed account and physical characteristics of the tension fractures/cracks found on the nucleus of 67P, which are the key to our outburst model. In Sect. 3 we present the general idea of the outburst model, provide numerical experiments to test its assumptions, and demonstrate that it naturally explains the observed outburst’s characteristics. In Sect. 4 we contrast the current model with other existing ideas of sudden outbursts that can potentially provide explanations, and in Sect. 5 we summarize the main results and outline future steps in applying and validating the current model.

2. Cometary fractures: geological view

Several observations made by Rosetta near perihelion (±3 weak from August 13, 2015) revealed the presence of strong outburst activity which was unexpectedly short lived (A’Hearn et al. 2016). One such event was recorded by the NAVCAM on August 8, 2015, and is shown in Fig. 1. We estimate the source of this outburst to lie in the transition between the Aker and Anhur regions in the heavily fractured terrains in the southern hemisphere (El-Maarry et al. 2015b,a). In general, in this work we do not concentrate on any particular outburst event, but present a general framework for strong and short-lived outburst activity.

Fig. 1 Outburst event captured by the NAVCAM on August 8, 2015 (ROS_CAM1_20150808T152148). The brightness is in linear scale. The estimated source region of this particular event lies around the transition between the Aker and Anhur regions.

The surface of 67P was found to be heavily fractured, which is apparent at nearly all spatial scales on the exposed surfaces; the fractures either form in networks or occur in isolation (El-Maarry et al. 2015a). In order to estimate the typical physical characteristics of the cracks/fractures we begin our discussion with the Anuket region and argue by analogy for a similar type of surface fractures elsewhere on the nucleus. The Anuket region features a system of lineaments (a few meters wide and up to a hundred meters long) trending generally along the neck–head boundary, which were first shortly described by Basilevsky et al. (2015) and Thomas et al. (2015) as probable tension fractures. The fractures are located in the lower part of the slope of the Anuket region towards Hapi and oriented generally perpendicularly to the local slope and essentially along the axis of the nucleus neck (Fig. 2).

Fig. 2 Fragment of a NAVCAM image taken on October 18, 2014, showing lineaments that are probably tension fractures. The solid lines indicate the obvious features, while dotted lines outline the suggested features.

As previously mentioned, the lineaments’ sinuous outlines and variable widths suggest that they are tension fractures known in many localities on Earth and on other bodies in the solar system (e.g., Gudmundsson 1987; Cohn & Ahrens 1981; Tanaka & Golombek 1989). Figure 3 shows a system of tension fractures at the southwest end of the boundary fault Almannagja in the Thingvellir Graben of the Holocene rift zone of southwest Iceland. The fractures deform the Holocene pahoehoe basalt lava flow (Gudmundsson 1987).

Fig. 3 System of tension fractures at the southwest end of the boundary fault Almannagja in the Thingvellir Graben of the Holocene rift zone of southwest Iceland. Courtesy of Agust Gudmundsson (http://agustgudmundsson.co.uk/wp-content/uploads/2012/04/Fig.-16-1024x705.jpg).

Gudmundsson (1987), in reference to the works of Lachenbruch (1961) and Nur (1982), argues that in evolving tension-fracture systems, the fracture depth should be the same order of magnitude as, or less than, the fracture spacing, and that the length of the fractures should usually be equal to or greater than their depth. The Anuket lineaments are arranged in-echelon system with fracture-to-fracture spacing of 15 to 70 m, and the length of individual fractures that vary from 30 to 300 m. These values were estimated from the images themselves and validated using synthetic data with the Preusker et al. (2015) shape model. Thus, it can be concluded that the depth of Anuket tension fractures ranges from 15 to 70 m. Following Gudmundsson (2011) we can also estimate a potential maximum depth of tension fractures in the consolidated material of the 67P nucleus using an alternative method. For this purpose we use the Griffith failure criterion, which describes the maximum depth of tensile fracture as d_max = 3σ_T/ (ρg), where σ_T is tensile strength, ρ is the density of the host material, and g is surface gravity acceleration (Gudmundsson 2011). According Groussin et al. (2015)σ_T ≈3−15 Pa, the density ρ ≈ 470 kg m^-3 and g ≈ 2 × 10^-4 m s^-2, then the maximum fracture depth is d_max ≈100−500 m. So, the depth of the fractures present in the Anuket region could range from 15 to 70 m, but may be as deep as 100 to 500 m (although the range of the two estimates suggests a maximum depth of 300−500 m). Such depths are significantly greater than the penetration depth of a thermal heat waves, both diurnal and seasonal as discussed in Sect. 3.2. Therefore, these fractures provide a likely connection between the deeper interior and the nucleus surface.

Thomas et al. (2015) argue that the fractures’ position close to the neck may be “indicative of stresses at the head–body interface which might ultimately lead to nucleus splitting”. El-Maarry et al. (2015a) suggested that the formation of these fractures could be attributed to rotational- or orbital-induced stresses in that specific part of the comet. Scheeres et al. (2016) conclude that the changing rotation period introduces cracks by stress at the same location where cracks are observed.

3. Gas-dust jet formation by deepening fracture

The main features of the strong short-lived outbursts are a short duration (less than an hour) and usually rather strong collimation of the jet. It should be noted that the outbursts do not necessarily coincide with highest illumination conditions (noon). A realistic outburst model must account for the above-mentioned key facts.

In this paper a novel mechanism for the transitional outburst activity is described. We propose that the deepening of a pre-existing crack (or fracture) leads to a sudden opening of a deeper interior region containing highly volatile ices (e.g., CO, CH₄). The sudden disequilibration of this thermodynamically steady state region results in a sudden and intense release of gas and dust. Because the sublimation is an endothermic process, the temperature of the active region rapidly decreases, which at some point effectively shuts down this activity, in accordance with the observations.

In this section we focus on the analysis of the physical conditions and the processes required to implement the scenario described above. The following issues are addressed in our model:

• Can the super-volatiles exist at the depth evaluated inSect. 2?

• Can the sudden disequilibration of the region, which is rich in the super-volatiles, produce reasonable amounts of dust and gas?

• Can this activity be stopped after a short time?

3.1. Depth of super-volatiles: model description

The presence of highly volatile substances in the coma is directly proved from the observations of CO₂, CO, and CH₄ by multiple instruments, and particularly of N₂ by the ROSINA instrument (Rubin et al. 2015). The authors of this article concluded that “cometary grains formed at low-temperature conditions below ≈30 K”. Two additional facts, namely a very low average thermal inertia of nucleus material determined from the MIRO measurements (Schloerb et al. 2015; Choukroun et al. 2015) and a significant loss of material due to sublimation estimated by Keller et al. (2015), also provide strong arguments that a pristine low temperature is preserved at larger depths. Although the presence of super-volatile species is an established fact, for our model we need to estimate the depth at which such ices may persist for a long time while the comet is in the inner solar system.

Because the orbital period is much greater than the nucleus rotation period (2398 vs. 0.5 days), a natural division of the relevant processes into long- and short-term exists. The difference between time scales also naturally causes differences in spatial scale since diffusion processes are considered. An analytical solution of the linear heat transfer equation shows that the temperature variations decrease exponentially with a greater depth into the subsurface, with an e-folding scale length equal to the so-called thermal skin depth (Eq. (1)) (1)where I_th is the thermal inertia of nucleus material, Ω is the period (orbital or rotational), c_p is the heat capacity, and ρ is the effective medium density. Using the values retrieved by the Rosetta instruments (ρ ≈ 470 kg m^-3, I_th ∈ (10−30) J m^-2 K^-1 s^0.5), we obtain L_Torbital ≈ 0.6 m and L_Trotational ≈ 0.01 m, which means that the orbital heat wave penetrates only about several meters into nucleus. As with many comets, we do not know how long comet 67P has been moving in its current orbit. Its dynamical evolution is rather uncertain (Groussin et al. 2007). It cannot be excluded that the comet’s orbit has changed dramatically. The perihelion distance changed after the close approach to Jupiter in 1959 (Krolikowska 2003).

To evaluate a global temperature distribution at larger depth of the nucleus, we use a thermal model that describes the heating and cooling of the nucleus located in the modern orbit. We consider a 1D heterogeneous two-layer model of nucleus, i.e., the “ice-free dust crust” is distinguished explicitly. The second layer is taken to be a porous mixture of water ice and dust. The thickness of the crust is fixed and it is a free model parameter. The sublimation of water ice is included in the model. When an ice-free porous dust layer is present on the nucleus surface, the effective sublimation rate is determined by the layer resistance. Because the ice sublimation is added to the simulation, the possibility of the corresponding nucleus erosion should also be included in the model. The 1D non-stationary non-linear heat equation is solved by a finite element method on an adaptive grid with moving boundaries of the computational domain. The boundary conditions on the nucleus surface, at the lower boundary of the dust layer, and at the depth of the isothermal core must be determined. The first two equations below express the energy balance:

here I_v is the solar irradiation, R_H is heliocentric distance, A_v is the albedo, ϵ is the emissivity, σ is the Stefan-Boltzmann constant, T_s is the surface temperature, is the sublimation rate, and is the latent heat of sublimation. The parameters λ_dust(T) and λ_core(T) is the effective heat conductivity of the dust and core layer, respectively.

The sublimation rate is given following the Hertz-Knudsen formula, Z(T) = P(T) / (0.5πν_th), where the water saturation vapor pressure is P(T) = 3.56 × 10¹² × exp(−6141.667 /T) kg m^-1 s^-2 (Fanale & Salvail 1984) and the thermal velocity m s^-1, and μ is the molar mass of water. The porous surface layer reduces the effective gas flux from under it. The permeability of the ice-free porous dust layer, is a function of the agglomerate radius R_A and the thickness of the dust layer L. For the collisionless gas the mass flow through pore volume can be calculated via a modified Clausing formula (Skorov et al. 2011).

An effective solid conductivity is retrieved from the thermal inertia. The value of thermal inertia, porosity, and material density reported in the most recent publications related to the Rosetta mission were used (Schloerb et al. 2015; Groussin et al. 2015; Preusker et al. 2015). The values of the thermal inertia obtained from MIRO correspond to temperatures near 200 K, where the radiative conductivity becomes important only for the large dust particles r> 100 μm (see Fig. 11 in Gundlach & Blum 2012). Thus the thermal conductivity retrieved from observations might be treated as an upper limit for the solid conductivity of the surface layer. We assume that the solid conductivities of the dust crust and the ice-dust mixture are the same and that the material porosity is constant everywhere (as in Skorov & Blum 2012). In addition to the solid conductivity, the radiative conductivity as a part of the effective conductivity is considered in the model (Gundlach & Blum 2012). This introduces a temperature dependence to the effective conductivity and a relative increase of this value at small heliocentric distances. Kossacki et al. (2015) assumed that the sintering of ice due to vapor re-condensation inside pores can increase the contact area between the ice particles, which is typically associated with an increase in solid thermal conductivity of icy regions. However, the ice sintering can only play a role near the surface, where water ice can sublime effectively. Because the nucleus is most likely composed mainly of dust (Rotundi et al. 2015), the significance of this effect should be small; even so, we performed calculations for the extreme model case where the solid thermal conductivity of the inner layers is about ten times higher than the solid thermal conductivity of the surface layer.

The bottom boundary condition determines the zero thermal flux on the interior boundary. If the corresponding depth is selected properly, the temperature at the interface does not change with time and equals the initial pristine temperature of the comet nucleus. This depth should be much larger than the orbital thermal skin depth and this parameter is initially evaluated following (McKay et al. 1986). The total depth of the modeled region is a free parameter that is verified by the calculation results.

3.2. Depth of super-volatiles: model results

Several numerical experiments were set up in order to confidently estimate the depth where highly volatile ices may remain in near thermodynamic equilibrium for a long time. Model A neglects sublimation process, Model B accounts for a free ice sublimation, i.e., no porous dust layer on the surface, and Model C includes a porous dust layer (restricted permeability of surface layer) above the ice layer. In addition, for each model we distinguish fast and slow rotators, determined by the value of the cosine of local solar zenith angle (SZA) (0.25 and 1 respectively). Finally, these simulations are performed for three values of the thermal inertia, I_th: low (10 J m^-2 K^-1 s^0.5), high (30 J m^-2 K^-1 s^0.5), and with inertia three times higher than the upper MIRO evaluation, i.e., 90 J m^-2 K^-1 s^0.5 (Schloerb et al. 2015). Effectively, Models A–C test the role of the different terms in Eqs. (2) and (3) on the depth where highly volatile ices may survive for a long time.

Figures 4–7 show the temperature profiles versus depth calculated for the models described above. Calculations were made for two primordial temperatures (20 K and 30 K) running the thermal model for ten full cometary orbits. A further increase of the simulation interval leads to an insignificant increase in the heated region, and does not change the results significantly. All temperature profiles presented in this section are calculated for the comet orbital position corresponding to early August 2015 (R_h ≈1.26 AU).

Fig. 4 Temperature as a function of cometary depth for the orbital position corresponding to early August 2015. A slow rotator approach is applied for Model A (surface temperatures close to 350 K) and for Model B (surface temperatures close to 200 K). Two primordial temperature values are tested: 20 K (dashed) and 30 K (solid). See Sect. 3.2 for a detailed description.

Fig. 5 Temperature as a function of cometary depth for the orbital position corresponding to early August 2015. The Model B slow rotator at two primordial temperatures is tested: 20 K (dashed lines) and 30 K (solid lines). Results are shown for three thermal inertia values 10 (thin), 30 (normal), and 90 J m-2 K-1 s0.5 (bold) lines. Generally, the higher the thermal inertia the deeper the orbital heat wave penetrates (see Sect. 3.2 for further description).

The results of calculations for Model A and Model B for a slow rotator with I_th= 30 J m^-2 K^-1 s^0.5 are shown in Fig. 4, and isolate only the effects of incoming energy on the internal temperature profile. It is clear that the absorbed energy flux does not strongly modify temperature profiles at depths greater than about 10 m. In all cases the temperature gradient becomes small (~0.2 K/m) at a depth of about 30 m, and this scale is mainly determined by the effective material conductivity. In Fig. 5 we plot the temperature profile after ten orbital revolutions for Model B, for three different values of thermal inertia and two core temperatures. The case where I_th = 90 J m^-2 K^-1 s^0.5 can be qualitatively considered as an illustration of the ice sintering effect, when conductivity grows owing to vapor re-condensation (Kossacki et al. 2015). The calculations show that for the thermal inertia consistent with the MIRO data, the penetration depth of the heat wave reaches about 30 m. However, the model with ice sintering allows the heat wave to penetrate to a depth of about hundred meters. Finally, in Fig. 6 we show the results of Model C and its effects on the heat wave penetration into the deep interior. The porous dust layer has minor effects where only the first few meters below the surface are noticeably affected, and the temperature profiles converge at a depth around 30 m, irrespective of the primordial temperature of the interior. Based on these results we conclude that neither the variation of the absorbed energy nor the energy loss due to ice sublimation (with or without the presence of a surface porous dust layer) significantly changes the penetration depth of a heat wave into the deep interior during the orbital motions of comet. For the fixed orbital parameters this depth is mainly determined by the thermal inertia of the nucleus material, and for the inertia values obtained from the MIRO measurements the heated depth does not exceed 30−40 m.

Fig. 6 Temperature as a function of cometary depth for the orbital position corresponding to early August 2015. Model C slow and fast rotator cases are studied for two primordial core temperatures, 20 K and 30 K. The curves (solid and dashed) starting at higher surface temperature correspond to a slow rotator. The MIRO mean upper value of thermal inertia of 30 J m-2 K-1 s0.5 is used (Schloerb et al. 2015).

However, there is another process that needs to be taken into account since it may modify the heat penetration depth into the deep interior, i.e., sublimation erosion (removal of ice and dust from the surface). Erosion leads to a moving boundary and, as a result, to a shrinking of the heated region (Gortas et al. 2011). This problem is commonly referred to as a Stefan problem. The effective rate of erosion was determined from the energy balance on the surface and the structure of the surface layer (Keller et al. 2015). In Fig. 7 results are presented for the case where the surface is shifted owing to reduced ice sublimation in the case of thermal inertia of 30 J m^-2 K^-1 s^0.5. Surface shifts are plotted for several cases (orbits 5 to 10). The model with material ablation causes the orbital heat wave to penetrate only to a shallower depth, about ten meters, which would mean the high volatiles can be preserved for a long time at even smaller depths.

Fig. 7 Temperature as a function of cometary depth for the orbital position corresponding to early August 2015. A slow rotator model is applied for Model C accounting for surface erosion. The depth is measured from the surface position after ten orbital periods (showing results for orbits 5−10). The primordial temperature is fixed at 30 K and the thermal inertia is set to 30 J m-2 K-1 s0.5 (upper value from MIRO analysis).

3.3. Transient processes in a gas-dust-ice layer

Next we investigate thermophysical conditions existing in the unheated interior prior to fracture penetration and the short evolution of a gas-laden ice-dust porous media just after the fracture opening. The sublimation energy of a CO ice is approximately ten times smaller than that of water ice (see Gundlach et al. 2015, and references therein). Thus, it is expected that the medium pores at temperatures above 30 K are filled by CO vapor. For example, the CO ice has saturation pressure of about 0.01 Pa at 30 K, while at 35 K it is about two and a half times higher than the saturation pressure of water ice at 200 K. The saturation pressure of CO vapor increases very quickly with rising temperature, and already at 40 K it is about 6 Pa, i.e., the pressure rises about six hundred times for a temperature increase from 30 to 40 K (see Fig. 8). The mass of CO in a gas phase inside the porous medium in a dynamic balance with ice can be easily estimated assuming that the porosity does not vary significantly with depth. When 70% of a volume is occupied by gas (porosity), then the ideal gas law estimates about 0.5 g of saturated gas per cubic meter of cometary material at 40 K. The mass ratio of volatile to non-volatile components is believed to be 1/4 (Rotundi et al. 2015) and the effective density about 470 kg m^-3 (Sierks et al. 2015), which means that one cubic meter of interior material contains a few kilograms of CO for the (mean) 5% of relative abundance to water (Morse et al. 2015). It means that the dominant part of CO remains in a solid phase with gas in a dynamic equilibrium. Taking these simple evaluations with the results presented in the previous section confirm that CO ice can exist for a long time in the unheated depth of a nucleus formed at low-temperature conditions (30−40 K).

Fig. 8 CO saturation vapor pressure (left y-axis scale, solid line) as a function of temperature. The corresponding mass loss rate is shown on the right y-axis scale (dashed line).

Next we try to evaluate the physical state of the deeper interior medium with super-volatile ices present. A gas diffusion always exists in the porous ice-free non-volatile region, bordering a medium with sublimating ice at the lower and the vacuum at the upper boundary. This diffusion occurs because of the concentration gradient. If the molecular free path is larger than the characteristic pore size, the Knudsen diffusion determines the gas density, which goes from saturated values to the outgassing density value. The characteristic escaping time of molecules Γ_esc can be evaluated using the classical formulas of free molecular diffusion into a random porous media (4)where L is the depth of sublimating ice and D_diff is the effective gas diffusivity, expressed as (5)and D_kn is the Knudsen diffusivity of gas, φ is the porosity, and ζ is the tortuosity of medium. The Knudsen diffusion depends on the pore diameter d_pore, species molecular weight μ, and temperature T: (6)If we assume that the pore diameter is about one millimeter and the tortuosity is about 3 (see references in Skorov & Blum 2012), then the expected escaping time is a few days. At the same time, diffusive resistance of the overlying layer is so high (it is proportional to the ratio of the depth to the pore diameter) that only a fraction of 100 ppm of sublimating molecules escapes, which means that the effective mass loss is very low. Therefore, we conclude that there is a stationary diffusion of the sublimation product (e.g., CO) and the macroscopic gas speed (flow speed) is close to zero at the sublimation depth (i.e., gas is at rest at the depth where ice exists).

Based on these estimates, we qualitatively analyze the scenario when the crack reaches a layer containing super-volatile ices, dust, and gas and destroys the equilibrium, which is illustrated graphically in Fig. 9. The deepening crack opens a boundary and gas starts to diffuse immediately. Assuming that the total amount of CO gas released during the outburst is about 1 kg and the mass ratio of carbon monoxide to water is about 0.05 (Morse et al. 2015), then (using the ideal gas law) the CO vapor has to be evacuated from a region of about 2000 m³ if the region T = 40 K. We can also estimate the opening area of the crack required to release this amount of gas. The diffusive resistance of a porous layer is proportional to the ratio of a layer thickness to a pore size (in our approximation particle size can be used in place of pore size): outgassing drops very quickly when the gas filling the pores is released and the saturated front goes down. For millimeter-sized particles the effective evacuation depth is smaller than a decimeter, which means that the open area should be about 20 000 m² if the layer temperature is about 40 K. Thus in the conservative scenario (grain size below one millimeter) the diffusion of existing CO vapor is not large enough to explain the observed strong outburst.

Fig. 9 Representation of the proposed model for the transient outburst activity as a result of pre-existing fracture deepening into a region of super-volatile ices. See text for a detailed description of the mechanism.

So far the vapor release from the pores has been treated as a single source of the outgassing and activity. However, we noted that CO ice sublimates actively at temperatures of 30 K and above. Therefore, a deepening crack creates not only conditions for diffusion of vapor filling the pores, but also conditions for effective mass loss owing to super-volatile ice sublimation (Fig. 9). As discussed previously, at low temperatures (T ≲ 40 K) only a small part of CO mass is in a gas phase, whereas the ice mass is about 5 kilograms per cubic meter. Thus, if about twenty thousand square meters are required to produce a kilogram of CO via vapor release from fine pores (if the region T = 40 K), a ten thousand times smaller area is needed for CO production due to ice sublimation. In that case only two square meters are required if there are millimeter-sized grains. Obviously, the sublimation of ice is a much more powerful source of outburst. Next, we address the question of what process is present and how quickly it stops the outburst.

Sublimation is a strongly endothermic process. The latent heat of sublimation of CO ice is about 3.0 × 10⁵ J kg^-1 (Brown & Ziegler 1979). It follows that without external sources of energy (we consider depths where thermal flux due to the solar illumination goes to zero) the medium is quickly cooled by the sublimation of super-volatile ices. A 1D thermophysical modeling is performed to investigate the cooling rate. The external boundary in the model is the open side of the deepening crack. The thermal re-emission and the thermal flux from the overlying region are negligible, and the initial medium state is isothermal (Fig. 9). The numerical experiment assumed primordial temperatures in the range 32−38 K, which implies a saturation pressure difference of more than factor of forty in this temperature range. The results of these simulations are presented in Fig. 10.

Fig. 10 Left panel: time evolution of external boundary temperature during active sublimation of super-volatile ice (CO) triggered by the deepening crack. Solid curves: the effective thermal conductivity corresponds to the thermal inertia equal to the most probable value (MIRO). Dashed curves: the effective thermal conductivity is increased ten times. Right panel: same for the case of relative mass loss.

The decrease in the temperature of the external boundary (crack opening interface) as a function of time is shown in Fig. 10 panel A (left). For all the considered cases of initial temperature, the final temperature drops strongly within 20 min. As expected, if the initial temperature is higher the difference between final and initial temperatures is also higher. Panel B in Fig. 10 (right) shows the corresponding (normalized) mass production as a function of time for the different initial temperatures. The mass production goes down by a factor of ten for T₀= 32 K and by a factor of one hundred for T₀= 38 K in the range of 20 min. The duration of the activity is mainly determined by the sublimation energy loss and the effective conductivity of the medium: the first determines an energy sink; the second determines an energy influx from interior. To examine the effect related to the warming of the open boundary region from the deeper interior we compare two limiting cases with low thermal inertia (I_th = 10 J m^-2 K^-1 s^0.5), and the very high inertia (I_th = 90 J m^-2 K^-1 s^0.5). When the effective conductivity is higher, the temperature cools down more slowly owing to the increased heat flux from the warmer interior. These effects are directly reflected in the mass production (Fig. 10, panel B). The relaxation time is rather insensitive to the initial temperature and in about twenty minutes the open boundary temperature becomes so low that the sublimation is effectively stopped. Therefore, the active short-lived release of gas and dust that was observed is naturally explained by the proposed model.

4. Discussion

Above we propose and describe a new thermophysical model that successfully explains the main features of the several short-lived outbursts as observed by the Rosetta NAVCAM and OSIRIS (A’Hearn et al. 2016). Studies of cometary outbursts have a long history, and below we contrast the current model with two other models which may be able to explain the observation. Nevertheless, it will be shown that they cannot fully capture all the characteristics of this outburst.

4.1. Subsurface geyser

Sekanina (1991) was one of the first to speculate that the development of cracks in the nucleus leads to the formation of under-surface cavities. Then, the dust jet is formed by ice sublimation with subsequent dust removal. In his model, the formation of the cavity is a result of a vent opening into the cometary crust, but cavities are not considered inherent to the nucleus. Yelle et al. (2004) subsequently developed this idea further, considering the collimated dust jets observed in the inner coma of Comet 19P/Borelly. They presented a thermophysical model describing the ice sublimation from the subsurface cavity with a gas escape and dust lifting through a narrow orifice to the surface, the cometary geyser formation. The porous dust layer cover induces an increase in the ice temperature and the gas pressure under the crust (called the “cooking effect”). Owing to pressure growth a vent is opened and rapid water ice sublimation starts. We should note that in these models of activity the cavity formation and the gas/dust release are both initiated by the heat flux due to the solar insolation of the cometary surface. Therefore, the suggested mechanisms are valid only for a subsurface layer in which the heat flux due to conductivity and the medium temperature are high enough to stimulate a rapid water ice sublimation. The jet activity (its duration and intensity) is also then directly connected with the solar insolation and governed by the supplied energy.

Both of those models lack applicability to the observed short-lived outbursts. The MIRO estimates of thermal inertia clearly show that water ice sublimation occurs only in proximity to the surface (see also Figs. 5−8 where temperature profiles are presented). In addition the cooking effect occurs only under a thin layer owing to the very low conductivity temperature; as the pressure starts to decrease with increasing layer thickness, the formation of vents becomes impossible. It is also natural to suppose that gas and dust release from under a thin porous crust produces a spreading plume instead of a collimated jet (Ipatov 2012). Finally, the characteristic duration of such an outburst should be similar to the rotation period of the nucleus, i.e., it is at least an order of magnitude longer than the duration of considered outburst. In the case of jets on the comet 19P/Borrelly, Yelle et al. (2004) simulated slow outbursts with a duration of months. Therefore, we conclude that the original form of this model of developing a cavity is inapplicable and cannot be used to explain specific features of the observed short-lived outburst.

An alternative view is that the existence of subsurface cavities is an inseparable characteristic of a cometary nucleus (see, e.g., Hughes 1984). This theory argues that cometary nuclei are probably formed by the aggregation of cometesimals of different sizes, which results in the creation of voids in the form of caves among the cometesimals. This hypothesis was recently discussed in Gronkowski & Weslowski (2015) and Ipatov (2012). In this scenario the outburst (dusty jet) is caused by the abrupt opening of a large subsurface cavity containing gas under high pressure. The beginning of activity is not directly connected with the solar insolation and may be triggered by external factors (e.g., by collision and abrupt cavity excavation, Ipatov & A’Hearn 2011).

Comparing the two scenarios for the formation of dust emission, Ipatov (2012) talks about rapid (typical duration of a few minutes) and slow (up to months) outbursts. Considering the observations of comet 9P/Tempel he connects the rapid outbursts only with the excavation of pre-existing cavities of different sizes placed at different depths. Thus, the opening (uncorking) of a cavity is the preferable assumption for our case, when a very short duration of the outburst is its peculiar characteristic. However, below we set up an argument that this idea is inconsistent, hence, not physical.

First, only super-volatiles can fill the cavities lying at a depth of several meters below the surface. At a temperature of about 40 K one cubic meter of nucleus volume contains less than 1 g of saturated CO vapor. Therefore, to release one kilogram of carbon monoxide, a cavity with a volume greater than 1000 m³ is required. The size of such a cavity is close to the limiting evaluation of interior heterogeneities obtained by CONSERT (Kofman et al. 2015). Second, the cavity scenario requires that the gas should remain in the pre-existing cavity for a long time, i.e., a time at least comparable to the typical lifetime of comets in the inner solar system. This is only possible if the diffusivity of media is very low, which also requires characteristic pore size to be very low. Therefore, this scenario is in a contradiction: the size of the cavity must be large (characteristic radius 10 m), but the surrounding pores must be extremely small (<1−10 μm). The presence of many large voids leads to a very rapid vapor transport to the surface and final ice depletion. Finally, because of the low thermal inertia the water ice sublimation can occur only very close to the surface (the first meter under the surface) (see Figs. 5−7 for temperature distribution). Therefore, the cavity cannot be filled with water vapor. We conclude that the short outburst from large primary pre-existing cavities filled by quasi-saturated vapor of super-volatiles is not a physically viable scenario.

4.2. Crystallization of amorphous ice and release of trapped super-volatiles

Cometary outburst formation due to crystallization of the amorphous water ice also deserves attention in connection with the discussed short-lived outbursts. The very possibility of the existence of amorphous water ice in the nucleus has been in debate for several decades (Bar-Nun et al. 1987; Prialnik & Bar-Nun 1987, 1990; Prialnik 1997; Bar-Nun & Laufer 2003 and many others). The rebuttal of the idea about the existence of amorphous water ice in comets may be found (e.g., Sandford & Allamandola 1990; Kouchi & Sirono 2001; Ivanova et al. 2011). Our focus is only to estimate whether the crystallization of amorphous ice can be a possible driver for the very short outburst activity.

Water condensing to amorphous ice can trap a quantity of molecules of other gases like CO and CO₂. Amorphous water ice is not stable and transforms into a crystalline form whenever conditions allow. The crystallization rate is an exponential function of temperature and becomes rapid at temperature above ≈140 K (Schmitt et al. 1989). It is widely accepted that the phase transition of pure water ice is an exothermic process (Ghormley 1968). In turn, it causes the release of trapped super volatiles. Klinger (1980) and Smoluchowski (1981) were the first to note that this necessarily leads to the runaway crystallization effect. Later, the problem was studied in detail by Prialnik (1997) and González et al. (2008). In the first paper the distant activity of comet (C/1995 O1) Halle-Bopp was investigated and showed that in the first cometary passage (fresh amorphous ice) the runaway crystallization subsides on a time scale of about 100 days. In the second publication, the authors demonstrated that the crystallization evolution depends strongly on the characteristics of the nucleus interior and on the crystallization energy release. They conclude that “the evolution of the crystallization front becomes continuous and almost steady, if the net energy released is half that of the pure, exothermic case, regardless of the characteristics of the nucleus interior”. They also found that “if the dust-to-ice ratio is high (larger than 1) and/or the thermal inertia is high (larger than 100), the crystallization front evolves in a continuous and smooth manner, even for pure, exothermic crystallization”.

The physics of amorphous ice crystallization as summarized above suggests that in the case of 67P these scenarios cannot explain the sudden and short-lived outburst. The dust-to-gas ratio in this comet is estimated to be about 4 (Rotundi et al. 2015), which would imply that the crystallization occurs in a continuous and smooth regime (González et al. 2008). Second, as presented in Sect. 3.2, the subsurface temperatures above 140 K are reached within a meter for the MIRO estimated thermal inertia. We conclude that the current understanding of comet activity driven by amorphous ice crystallization implies longer time scales and less rapid increase of activity than those implied by the observations.

5. Conclusion

The comet 67P experienced numerous strong and very shot-lived outbursts around perihelion passage (August 13, 2015) as documented by the Rosetta NAVCAM and OSIRIS cameras (A’Hearn et al. 2016). In this paper we propose and model a new mechanism for these short-lived outbursts of cometary activity. The presented model explains the unique outburst characteristics, and actually predicts that these events are even shorter, ≈20−30 min, which is also now supported by observations (A’Hearn et al. 2016). We also reviewed other proposed mechanisms for cometary outbursts available in the literature; however, basic physical considerations show that they fail to explain this particular outburst.

The new model is build around the geomorphological features of fractures and cracks observed on this particular cometary nucleus, which are argued to extend into the deeper interior where super-volatile ices (e.g., CO) can be preserved for a long time in steady-state conditions. We speculate as the nucleus experiences growing gravitational and non-gravitational stresses during the perihelion approach, these cracks widen and deepen (Scheeres et al. 2016). The deepening crack leads to a sudden opening and subsequent disequilibration of the volatile-ice laden region resulting in a violent release of gas and dust. The sublimation process induces a rapid decrease in temperature, effectively shutting itself down, at which point the medium slowly returns to its equilibrium (i.e., very weak diffusive sublimation). The new features of the model which are pertinent to the explanation of the observed outburst are

• the independence of external heat flux induced by solarirradiation;

• natural self-induced quenching of the outburst activity after a short time;

• in addition, the model does not require presence of large isolated cavities in the nucleus,

• and/or crystallization of amorphous ice and trapped super-volatiles release.

In the current stage, our model can be used to evaluate the expected gas and dust mass production during the outburst induced by the deepening crack. It also provides a unique opportunity to evaluate physical properties of the deep region because the freed dust and gas are evacuated from the nucleus deeper regions. Nevertheless, there are no quantitative data yet available regarding the amount of extra gas and dust released from the other instruments on board Rosetta. We plan to continue this work as the needed data are published in the future. Other important future work is to build a physical model of the cracks and their deepening mechanism, and analyze more of these outbursts that have been observed (A’Hearn et al. 2016); at least 25 of similar nature (J. B. Vincent, priv. comm., Sept. 2015). We believe that, at least, for this particular nucleus, the outbursts induced by the geomorphological activity of nucleus fractures can be a common mechanism. Furthermore, we intend to investigate whether this mechanism can explain the outbursts of main belt comets, especially understanding why no water vapor has ever been observed even with most sensitive water vapor detection instruments (O’Rourke et al. 2013; de Val-Borro et al. 2012).

Acknowledgments

We acknowledge the entire European Space Agency (ESA) Rosetta team and thank them for their contribution, without which this work could not have been done. Rosetta is an ESA mission with contributions from its member states and NASA. The NAVCAM images (ESA/Rosetta/NAVCAM – CC BY-SA IGO 3.0) are provided by ESA archive under Creative Commons Attribution-ShareAlike 3.0 IGO (CC BY-SA 3.0 IGO) licence.

References

All Figures

	Fig. 1 Outburst event captured by the NAVCAM on August 8, 2015 (ROS_CAM1_20150808T152148). The brightness is in linear scale. The estimated source region of this particular event lies around the transition between the Aker and Anhur regions.
In the text

	Fig. 2 Fragment of a NAVCAM image taken on October 18, 2014, showing lineaments that are probably tension fractures. The solid lines indicate the obvious features, while dotted lines outline the suggested features.
In the text

	Fig. 3 System of tension fractures at the southwest end of the boundary fault Almannagja in the Thingvellir Graben of the Holocene rift zone of southwest Iceland. Courtesy of Agust Gudmundsson (http://agustgudmundsson.co.uk/wp-content/uploads/2012/04/Fig.-16-1024x705.jpg).
In the text

	Fig. 4 Temperature as a function of cometary depth for the orbital position corresponding to early August 2015. A slow rotator approach is applied for Model A (surface temperatures close to 350 K) and for Model B (surface temperatures close to 200 K). Two primordial temperature values are tested: 20 K (dashed) and 30 K (solid). See Sect. 3.2 for a detailed description.
In the text

Fig. 5 Temperature as a function of cometary depth for the orbital position corresponding to early August 2015. The Model B slow rotator at two primordial temperatures is tested: 20 K (dashed lines) and 30 K (solid lines). Results are shown for three thermal inertia values 10 (thin), 30 (normal), and 90 J m-2 K-1 s0.5 (bold) lines. Generally, the higher the thermal inertia the deeper the orbital heat wave penetrates (see Sect. 3.2 for further description).

In the text

Fig. 6 Temperature as a function of cometary depth for the orbital position corresponding to early August 2015. Model C slow and fast rotator cases are studied for two primordial core temperatures, 20 K and 30 K. The curves (solid and dashed) starting at higher surface temperature correspond to a slow rotator. The MIRO mean upper value of thermal inertia of 30 J m-2 K-1 s0.5 is used (Schloerb et al. 2015).

In the text

Fig. 7 Temperature as a function of cometary depth for the orbital position corresponding to early August 2015. A slow rotator model is applied for Model C accounting for surface erosion. The depth is measured from the surface position after ten orbital periods (showing results for orbits 5−10). The primordial temperature is fixed at 30 K and the thermal inertia is set to 30 J m-2 K-1 s0.5 (upper value from MIRO analysis).

In the text

	Fig. 8 CO saturation vapor pressure (left y-axis scale, solid line) as a function of temperature. The corresponding mass loss rate is shown on the right y-axis scale (dashed line).
In the text

	Fig. 9 Representation of the proposed model for the transient outburst activity as a result of pre-existing fracture deepening into a region of super-volatile ices. See text for a detailed description of the mechanism.
In the text

Fig. 10 Left panel: time evolution of external boundary temperature during active sublimation of super-volatile ice (CO) triggered by the deepening crack. Solid curves: the effective thermal conductivity corresponds to the thermal inertia equal to the most probable value (MIRO). Dashed curves: the effective thermal conductivity is increased ten times. Right panel: same for the case of relative mass loss.

In the text