Free Access
Issue
A&A
Volume 655, November 2021
Article Number A20
Number of page(s) 14
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/202039824
Published online 03 November 2021

© ESO 2021

1 Introduction

One way to investigate the history of the Solar System is through the study of its main components: planets, moons, and asteroids. However, comets, among other small bodies in the solar system, are considered to be special because their structure and composition may not have changed significantly since their formation, which is widely assumed to have occurred about 4.6 billion years ago (Bouvier & Wadhwa 2010). Therefore, the study of comets and their evolution may be key to understanding the earliest stages of our Solar System’s evolution.

The Rosetta mission accompanied comet 67P/Churyumov-Gerasimenko (hereafter 67P/CG) from August 2014 to September 2016, observing gas and dust activity in its coma. Thanks to the extensive payload complement, Rosetta collected significant data on the abundances of the major gas species near the nucleus as well as highly resolved spectra of the coma in the infrared and microwave regimes. The interpretation of these data requires a comprehensive study of the gas emission process and the gas dynamics of the inner coma. This allows us to constrain the physical processes governing the generation of the comae of comets and to link comae measurements with the near-surface structure and composition of the nucleus.

We have a basic understanding of the way gases are released from the nucleus in order to form the gas and dust comae as the comet approaches the Sun. We know that the production of these gases is driven by the incident solar radiation on the nucleus, and this leads to the sublimation of cometary ices. This, in turn, creates a gas flow that can drag dust particles away and into interplanetary space. The composition of the coma depends on the composition and thermal properties of the nucleus. In general, comets are often assumed to be a mixture of ice and dust (Whipple 1950; Mendis & Brin 1977; Fanale & Salvail 1984), however, there is an overall lack of knowledge on the actual internal structure of comets both at macroscopic and microscopic levels. For comet 67P/CG, the Rosetta Orbiter Spectrometer for Ion and Neutral Analysis (ROSINA) determined H2O, CO2, CO and O2 to be the most abundant gases in the coma (Hässig et al. 2015; Le Roy et al. 2015; Bieler et al. 2015a), which gives a clear indication of the composition of dominant ices within the nucleus. How these different ices are physically related to each other and distributed within the nucleus is a more complex issue for which we have no clear answer at this time.

Previous studies have used numerical data inversion methods or kinetic and hydro-dynamical approaches to derive the H2O activity distribution and evolution of densities observed by ROSINA in the coma of comet 67P/CG during different stages of the Rosetta mission (Bieler et al. 2015b; Biver et al. 2015; Marschall et al. 2016, 2019; Fougere et al. 2016a,b; Hoang et al. 2017; Combi et al. 2020). These studies have identified the heterogeneous distribution of H2O outgassing from the cometary surface, with strong sources localised in the Hapi region (Marschall et al. 2016), but it has been confirmed that at scales greater than ~400 m (Marschall et al. 2020), the outgassing of H2O can be assumed to be driven to a large extent by insolation with no significant thermal lag (Bieler et al. 2015b). However, CO2 does not seem to follow the same principle, because there are significant changes in the CO2/H2O mixing ratio in the coma (Bockelée-Morvan et al. 2015; Hässig et al. 2015; Fink et al. 2016; Fougere et al. 2016a; Migliorini et al. 2016). Bockelée-Morvan et al. (2015) suggested that observations of the infrared spectrometer for the Rosetta mission (VIRTIS-H) were consistent with the nightside outgassing of CO2, for example,which suggests that heat is still available to produce sublimation in the absence of insolation, on rotational timescales.

The Optical, Spectroscopic and Infrared Remote Imaging System (OSIRIS) on-board Rosetta has observed a very dark surface that, together with COSIMA dust measurements, suggested the presence of carbon-rich material mixed with an almost equal amount of silicates (Bardyn et al. 2017). The OSIRIS observations were used to estimate a geometric albedo of 6.5 ± 0.2% for its surface (Fornasier et al. 2015), which is in the same range as the albedo estimate of ~4% for other comets (Huebner et al. 2006). The Comet Nucleus Sounding Experiment by Radio wave Transmission (CONSERT) provided measurements that suggest a porosity between 75% to 85% (Kofman et al. 2015; Brouet et al. 2016). This is important because the conduction of heat into the interior decreases in highly porous materials (Seiferlin et al. 1996; Shoshany et al. 2002) and it can be used to constrain the properties of the dust mantle in thermophysical models.

Thermal conductivity is a key physical property for the heat exchange process within any material. In the case of comets, the thermal conductivity has been inferred to be extremely low (Huebner et al. 2006; Fernández et al. 2013). The degree to which a body stores or releases heat from its surface is quantified using the thermal inertia (Γ) and it is defined as the square root of the thermal conductivity (κ), the specific heat (c), and the bulk mass density (ρ): (1)

The units of thermal inertia are J m−2 K−1 s−1∕2, which we abbreviate here to thermal inertia units (TIU). Laboratory experiments for cometary analogues estimated thermal inertia values around 73 TIU for a dry and highly porous matrix and 258 TIU for a porous matrix with water vapour (Seiferlin et al. 1995, 1996). For comets 103P/Hartley 2 and 9P/Tempel 1, spacecraft observations suggest a thermal inertia lower than 250 TIU and 45 TIU, respectively (Groussin et al. 2013). However, other studies estimate a thermal inertia for comet 9P/Tempel 1 between 50–150 TIU for scarped/pitted terrain and up to 200 TIU if the terrain is relatively flat (Davidsson et al. 2013). In the case of comet 67P/CG, the Multipurpose Sensors for Surface and Sub-Surface Science (MUPUS) on the Philae lander was used to estimate a thermal inertia of 85 ± 35 TIU on the Abydos area of the nucleus (Spohn et al. 2015). Further analysis with the Microwave Instrument for the Rosetta Orbiter (MIRO) data has produced TI values from 10 to 50 TIU and suggests a highly insulating role for the dust at the surface (Schloerb et al. 2015). More recent studies using data from MIRO’s mm-channel imply this value is ≤80 TIU, while the Visible and InfraRed Thermal Imaging Spectrometer (VIRTIS) data suggests a broader range between 40 and 160 TIU (Marshall et al. 2018). Such low thermal inertia values would be consistent with the high porosity observed by CONSERT.

Although ices have been detected at the surfaces of cometary nuclei through infrared remote-sensing (Sunshine et al. 2006; Oklay et al. 2016; Barucci et al. 2016; Filacchione et al. 2016) and inferred from visual broad-band imaging (Pommerol et al. 2015; Fornasier et al. 2016), the areal extent of the ice exposed on the surface is totally insufficient to explain the total volatile mass loss (Sunshine et al. 2006; Fougere et al. 2016a,b; Marschall et al. 2016). Furthermore, measured surface temperatures under strong illumination far exceed the free sublimation temperature of water ice (~200K) (Emerich et al. 1988; Li et al. 2007; Groussin et al. 2013; Tosi et al. 2019). In addition, the mass loss per unit surface area is typically well below the free sublimation rate from a pure ice surface leading to the introduction of concepts such as the “effective active fraction” (EAF) to reduce the total sublimation rates in models to match those observed (A’Hearn et al. 1995; Lamy et al. 2001). In this way, the EAF can be interpreted either as the percentage of the surface area that is active or as the magnitude of the influence of sub-surface sublimation through a porous layer of relatively low thermal inertia (Skorov & Rickman 1995; Skorov et al. 1999, 2001, 2002; Thomas 2009).

We can construct a simple model to replicate cometary outgassing and perform a parameter study of the most important physical properties that influence the energy exchange in the sub-surface. The starting point of this work is a thermal model that neglects thermal inertia. This has been used before to replicate fairly accurately the H2O outgassing from 67P/CG (Marschall et al. 2016, 2017, 2019). However, this is probably not applicable to CO2 emissions. In order to include CO2 activity into the problem, we need to take into account the effect of thermal conductivity and the transport of heat to deeper layers of the nucleus. Thus, we have included this in our thermal model for different values of thermal inertia and we have simulated sublimation fronts at different depths for the two molecules.

Here, we study the gas emission distribution and combine this with gas dynamics simulations to investigate the effects on the coma distributions of H2O and CO2. It is important to clarify that the current work does not intend to actually fit Rosetta data, but only to understand the processes triggering nightside activity. These findings will be used in future studies to investigate extensively mixed gas emissions from the real comet’s shape and compare them with observations. In the following sections, we detail the applied method (Sect. 2) and we provide a description of the selected cases (Sect. 3) to model cometary activity for a simplified nucleus shape. We analyse the main findings of this work (Sect. 4) and we discuss their relevance for the comparison of numerical modeling methods with multi-instrument data. Finally, we present our conclusions (Sect. 6).

2 Method

In order to reproduce the collisional behaviour of the coma, we used a C++ parallel Direct Simulation Monte Carlo (DSMC) code, called ultraSPARTS (developed by Wu & Lian 2003; Wu et al. 2004; Wu & Tseng 2005). The DSMC code is a powerful method proposed by Graeme Bird (Bird 1994) that allows us to use a molecular model to find a probabilistic solution to the Boltzmann equation for every flow regime from continuum to free molecular flow. This method is very computationally demanding, but has been used for cometary studies in the past (Combi & Smyth 1988; Crifo et al. 2002, 2003, 2005; Zakharov et al. 2009; Finklenburg & Thomas 2014a; Finklenburg et al. 2014b; Liao et al. 2016; Marschall et al. 2016; Liao et al. 2018). The DSMC code solves the average properties of the gas flow around any type of 3D shape for given boundary conditions. In our study, this is particularly important as the gas field coming from the sublimation of ices at the surface of the nucleus encompasses different flow regimes. In the following, we discuss the assumptions we have made in order to simulate the coma around the nucleus, as well as the additional tools needed for the task.

2.1 Model assumptions

In order to model gas activity in the inner coma, we have chosen a case in which we assume the nucleus is a sphere with a 2 km radius. Although the considered spherical shape is highly simplified compared to the complex shape of for example 67P/CG’s nucleus, this approach reveals the physical interaction of the different volatile species mixed within the flow, while avoiding theadditional effects of irregular geometry. We have selected the outer limit of our simulation domain to be at 10 km from the centre of the nucleus, as beyond this distance the gas flow is supersonic (Mach number 2 < M0 < 5 in our simulations) and has been estimated to have a radiality1 between 96 and 100%.

We also neglect the effects of gravity in our study. This is based on the fact that the gravity field of comet 67P/CG is very weak (~10−4 m s−2 at the surface) and the speed at which gas particles leave the surface is much higher than the average escape velocity of 0.81 m s−1 (Thomas et al. 2015). Under the given illumination conditions most of the gas molecules are ejected from the surface at speeds above 100 m s−1. We estimatethat with constant speed, gas particles reach the edge of the simulation domain in less than 2 min. This is only 0.27% of the rotation period of comet 67P/CG (Keller et al. 2015), which means we can, in principle, neglect the influence of the nucleus rotation in the coma and use a set of steady DSMC solutions for the flow that correspond to boundary conditions for specific points in time. However, the nucleus rotation is taken into account in the thermal model used to compute the surface boundary conditions as explained later in Sect. 2.4.

We assume that the activity is mainly driven by H2O and, at a smaller percentage, by CO2, as these molecules have been determined as two of the main volatile species of 67P/CG coma (Le Roy et al. 2015; Combi et al. 2020). We use mixing ratios CO2/H2O in terms of mass loss of about 7% and 14%. We selected these values based on the mean estimates of the gas density above the nightside needed to fit ROSINA data (Bieler et al. 2015b). The day-to-nightside brightness ratio of the dust coma also indicates nightside outgassing (Gerig et al. 2020).

2.2 Collision model

For the production rates in the present work, the mean free path (λ) of the flow is larger than 10 m at 8 km above the sub-solar point, while above the anti-solar point, it can reach values greater than 10 km at the same distance. The local Knudsen number varies between 0.01 < Kn < 1000 from day to nightside, such that it covers flow regimes from the transitional to the free-molecular (Bird 1994), in which case the kinetics of molecular collisions need to be considered.

The DSMC method uses simulation particles, each of which represent a large number of molecules of the real gas. Our code tracks the random interaction of collision pairs (H2O–H2O, CO2 –CO2 and H2O–CO2) and the interaction of molecules with the boundaries in small time steps Δt, which have to be smaller than the mean collision time. We use the viscosity indexes (0.923 and = 0.93) and the molecular diameters (4.82 × 1010 m at 300 K and 5.54 × 1010 m at 273 K) of the Variable Soft Sphere (VSS) collision model (Smith & Callendar 1924; Bird 1994; Alexeenko et al. 2009). For H2O–CO2 collisions, the code uses the averaged parameters between the two particle species. Gas particles that interact with the surface boundary are diffusely reflected. The code stores the position, velocity, and rotational energy of each simulation particle to calculate the averaged macroscopic properties of the flow field. The energy redistribution applied between the translational and rotational degrees of freedom is the Larsen–Borgnakke model, for which the number of collisions necessary to reach equilibrium was selected as equal to 1 (Crifo et al. 2002). The rotational degrees of freedom are 3 and 2 for H2O and CO2, respectively,while the translational degrees of freedom are 3 for both molecules. This process is repeated multiple times until the variation in the number of particles, and in the average gas properties of the system, is very small (Wu & Tseng 2005; Scanlon et al. 2010).

2.3 Mesh generation

Although homogeneous properties have been set for the nucleus composition, the insolation conditions at the surface vary. This produces an anisotropic distribution of the flow that, for H2O, is locally denser at regions close to the sub-solar point. In order to track the interactions between simulated gas particles in the flow, the simulation domain has been split into volume elements that have to have a length dimension (Δs) of the same order as λ, which mathematically can be expressed through the equation (2)

Our DSMC simulations consequently require an unstructured mesh that takes into account the density variations within the simulated gas-flow. It is generated using PointwiseTM software and has 323 508 surface facets with an averaged Δs and angular size of 20 m and 1.96 arcsec, respectively. In total, our mesh has ~8 million tetrahedral volume elements (voxels), in which the code follows the simulation particles and calculates the average temperatures, velocities, and number densities of the gas.

2.4 Thermal model

In the simplest case of sublimation from the surface, we use a thermal model to calculate the energy balance at each surface facet, i, of the nucleus given the chemical composition of a single species gas. This is done using the surface energy balance equation at depth z = 0: (3)

The term to the left side of Eq. (3) is the input solar energy, where S is the solar irradiance at 1AU, A is the directional-hemispherical albedo, rh is the heliocentric distance and θi is the angle between the surface normal of facet i and the vector pointing to the sun. Whenever 90° < θi < 270°, it is assumed that the facet is not illuminated and the solar input is zero. The first term on the right side of Eq. (3) is the thermal radiated emission from the surface, where ϵ is the thermal emissivity, σ is the Stefan-Boltzmann constant, and T0,i is the surface temperature. The second term is the loss of energy through sublimation of ice per unit time, where L is the latent heat of the ice sublimation and dmi∕dt is the mass loss rate per unit time per unit area evaluated at z. The last term is the effect of thermal conduction and influences the amount of heat transported to depth, where κ is the thermal conductivity coefficient and dTi∕dz is the temperature gradient. The value of all parameters used in this model are listed in Table 1. For a superficial sublimation of ices, the effect of the heat conduction of the dust mantle in Eq. (3) is set to zero, because there is no desiccated dust mantle on top of the sublimation front.

For cases with sub-surface sublimation, the heat balance at a given sublimation front is given by (4a) (4b)

where z and z+ indicate the side before and after the correspondent sublimating interface for each gas species, respectively, for which the temperature gradient is calculated assuming a lower boundary condition for the temperature of TN = 50 K (Bonev et al. 2007, 2013). In the decoupled mode, Eqs. (4a) and (4b) are run separately (independent energy balance for each species). This corresponds to cases when there is an heterogeneous distribution of ices within the nucleus. In the coupled mode, the two species influence each other and both equations are used simultaneously to calculate the temperature gradient. In this case, the distribution of ices is more similar to an intimate mixture. At the same time, the balance equation at the surface layer is given by: (5)

Figure 1 shows vertical temperature profiles of two cases with: (1) sublimation fronts at the same depth for H2O and CO2 and (2) sublimation fronts at different depth for each species using the decoupled and coupled modes. Each of the panels in Fig. 1 represent different illumination conditions every 90° in longitude. The calculation using the decoupled mode retrieves two different temperature profiles for the same surface element depending on the type of molecule used at the sublimation front. Given that the internal structure of the cometary nucleus is not well understood and suspected to be very complex (Brandt 2007; Levasseur-Regourd et al. 2009), we can treat every surface element as a combination of different thermal structures. There might be region richer in H2O than in CO2 ice sources; in such cases, using the decoupled mode could still give a good estimate of the global distribution of activity. When the sublimation fronts are set at different depths, the temperature profile in coupled mode (solid pink line) is strongly influenced by the presence of CO2-ice at , decreasing the energy budget close to the surface compared to the case with pure CO2 sources at shallower depths (dashed-dotted pink line). This also depends on the local illumination conditions, such that the minimum temperatures close to the surface are produced at the sunrise terminator, while the maximum temperatures are located close to midday.

The mass loss rate dmi∕dt for a molecule of mass m is calculatedusing the Hertz–Knudsen equation: (6)

where Ps is the sublimation pressure calculated using empirical fits to experimental data (Huebner et al. 2006), kB is the Boltzmann constant, and ng and vg are the number density and the velocity of the gas, respectively.

So far, we have assumed that 100% of the icy-surfaces are active. However, the sublimation rate of each gas species can be modulated by an effective active fraction (EAF) such that the global production rate for each type of gas is (7)

where ai is the area per surface facet and Ns is the total number of facets.

This approach is a modified version of the “standard thermal model” for slow-rotators (Lebofsky & Spencer 1989; Festou et al. 2004; Huebner et al. 2006). It calculates the insolation condition at a certain heliocentric distance for one point on the surface after Nrot nucleus rotations on its axis. Figure 2 shows the calculation of the temperature at the surface, the temperature at the sublimation front, and the sublimation rate of (2a) H2O and (2b) CO2 as a function of Nrot. Both variables clearly increase with the number of rotations reaching an asymptotic state at about 20–40 nucleus rotations. We chose Nrot = 20 (equivalent to 248.12 h) for our model calculations, so that we can ignore changes in the insolation conditions caused by the obliquity of 52° (Lhotka et al. 2015) of the comet. In this way, we estimated the diurnal change in temperature and mass loss rate from a sub-surface layer. Different illumination conditions are shown with different symbols, from which we can notice that on the sunrise terminator the thermal calculation takes longer to reach a stable value for both variables, given that at this region the comet is poorly illuminated and it had enough time to cool down from the previous rotation. We also notice in Figs. 1 and 2 that in poorly illumination areas, the temperatureat the sublimation front is larger than the temperature at the surface, which is caused by the rapid lost of heat at the dusty surface combined with the heat transport into the interior. In this work, we assume κ, ρ, and c to be constant across any internal boundaries. This is a major simplification over the 1D analysis seen in Huebner et al. (2006), where these variables are temperature-dependent. Furthermore, heat transport by gas and the effect of increased subsurface gas pressure are ignored. This allows us to compute results for all surface facets (3D) of the comet at the cost of simplification of a complex problem.

Figure 3 illustrates a possible distribution of H2O, CO2 and the non-volatile material within the nucleus, with the presence of a desiccated dust mantle at the surface. Sublimation fronts are shown at different depths ( and ) under the assumption that CO2 is depleted close to the surface. Shadowed areas at the sublimation fronts are an interpretation of the EAF as a pattern of activity at the surface that is linked with the composition and activity of deeper layers of the nucleus. EAF is not implemented exactly as shown in Fig. 3, however, this illustration is used to display its meaning when scaling the sublimation rate in Eq. (7). Hence, dmi∕dt is throttled by the presence of a desiccated surface layer, which reduces the energy input to the sub-surface sublimation front. The sub-surface layer cannot be so deep that EAF exceeds 1, setting a natural limit on the depth of the sublimation front. The values of EAF used for H2O and CO2 are shown in Table 2. Because of flux conservation, the production rate at the sublimation front and at the surface has to be the same. Therefore, we assume the desiccated layer does not provide resistance to the gas flow.

Other authors (Skorov & Rickman 1995; Skorov et al. 1999, 2001, 2002; A’Hearn et al. 1995) have produced mathematical models of the surface layer to account for the effect of finite porosity on the gas flow from a sub-surface sublimation front. Thesemodels are, however, 1D, whereas our aim is to provide, for the first time, 3D solutions for the gas flow and, ultimately, for complex shapes as well. This additional complexity stands in the way of more sophisticated numerical solutions. We have ignored the influence of heating of gas through collision with a hotter dust mantle, such as Christou et al. (2018), in the present model, considering that these additional effects probably only result in relatively small changes in the position of the sublimation front due to the exponential change in sublimation rate with temperature. However, we have set a test with larger surface temperatures as an approximation for such physical processes.

Table 1

Physical parameters used in the model.

thumbnail Fig. 1

Vertical profiles of the temperature with depth for different insolation conditions on the equatorial plane of a comet with a spherical nucleus. Coupled cases with sublimation fronts at the same depth ( cm) are labelled as “Coupled 1,” while coupled cases with sublimation fronts at different depths (1.4 cm and 5.6 cm) are labeled as “Coupled 2.” Decoupled mode calculations are shown with dashed and dashed-dotted lines for H2O and CO2, respectively.

thumbnail Fig. 2

Decoupled calculation of the temperature at the surface (in black), at the sublimation front (in red) and the sublimation rate (dm/dt) of H2O and CO2 (with an EAF = 100%) as a function of the number of rotations of the nucleus for a thermal inertia value of 40 TIU for both species. Each symbol indicates a different observation geometry during one day. Circles indicate facets located at sunrise, squares are facets at midday, triangles are facets at sunset, and crosses are facets at midnight.

thumbnail Fig. 3

Sketch of the energy exchanges (heat sources and sinks) at the surface and sub-surface of a cometary nucleus: is the vector normal to the surface and θi is the solar incidence angle; T0 is the temperature at the surface of the nucleus; and are the temperatures at the sublimation fronts for the respective type of ice; and TN is the lower boundary temperature of the nucleus. The coupling mode is represented as a switch in between the H2O and CO2 sublimation fronts. The shadowed layers at and represent the selected effective active fractions (EAF) used at the sublimation front of each molecule.

Table 2

List of tested parameters for model calculations.

3 Boundary conditions

In this section, we describe the selection of models with different boundary conditions (listed in Table 2) based on previous estimates of the thermal inertia of comet 67P/CG. The listed free parameters have been varied to: first, show the effect of thermal inertia in the model and second, show the influence of including a small amount of CO2 in the gas flow field. All cases listed in Table 2, except case B2(C), are set using the decoupled mode of the thermal model in which the thermal properties of each front are independent from the other sublimation front. This is justified by the fact that the EAFs are generally small compared to the total area so that the H2O sublimation front can be independentof the CO2 front. Case B2(C) uses a coupled mode to account for the effectof including CO2 on H2O sublimation and vise versa. In this case, the sublimation of both species influences the temperature gradients in the system and, thus, the production rates of each species. We find, predictably, that the EAF for both species has to be slightly larger in order to compensate the different energy distribution between the two sublimation fronts compared to the decoupled mode. The H2O production rate is set to 50 kg s−1, which is a mean estimate of the global production rate measured by ROSINA for comet 67P/CG during its spring equinox at 1.67 AU.

We started out by neglecting the thermal conductivity term and studied a case in which the thermal inertia is set to zero (case A0). It is used as a standard case to compare with results of cases with different thermal inertia values. In the case where the thermal inertia is zero, the model requires the absence of a dust mantle on the top of the sublimation front of H2O, which leads to a very low EAF compared to the other cases where a dust mantle is introduced. The EAF to match the required production rate depends not only on the depth of the sublimation front, but also on the thermal conductivity of the nucleus as indicated by cases B2 and C1 in Table 2. The thermal skin depth (δ) is related to the thermal diffusivity (α) through the equation: (8)

where α measures the rate of transfer of heat from the surface to the interior of the nucleus and P is the rotation period of the nucleus, which value is listed in Table 1.

Sublimation fronts at different depth have been selected for H2O and CO2 molecules with the assumption that CO2 has to be more depleted close to the surface because of its lower free sublimation temperature. The chosen values are of course not unique, but are intended to investigate the range of possibilities.

We have used the temperature at the sublimation front of each molecule ( and ) as the boundary condition for our DSMC simulations, so that these are the temperatures at which the sublimating flow escapes from the (sub)surface into a vacuum. However, the dust mantle on top of the fronts is not transparent to gas emissions and there should be some heat exchange between the dust-matrix and the gas particles which would be responsible for higher temperatures of the gas leaving the nucleus (Christou et al. 2018) and ultimately results in higher terminal velocities for the flow. Consequently, we also used the surface temperature as a boundary condition for a variation of case B2 and studied how it changes the average characteristics of the gas flow.

4 Results

In the first part of this section, we discuss the outcome of the thermal model to study the daily evolution of the temperatures at the sublimation front and at the surface for all the cases listed in Table 2. In the second part, we use the slices through the DSMC domain along the equatorial plane to study the role of thermal inertia in modifying surface source distributions and the effect of CO2 activity in the flow dynamics.

4.1 Diurnal evolution of temperatures at the surface

Figures 4, 5, and 6 illustrate the diurnal variation of the temperatures at the surface, the H2O and CO2 sublimation fronts, respectively, for all cases listed in Table 2. Black lines indicate the planetographic longitude from 0 to 350° in steps of 10° and latitude from 0 to 90° in steps of 20° when looking above the northern hemisphere of the comet. The sub-solar point in all cases is in the equatorial plane at 0° longitude, while the sunset and sunrise terminators are at 90° and 270° in longitude, respectively. We notice that the only case without thermal inertia (A0) has a very abrupt change between the dayside-nightside temperatures, while in cases that include thermal inertia, the variation of temperatures is more gradual. In terms of surface temperature (Fig. 4), cases with dust mantles and low thermal inertia values (40 TIU) have higher temperatures on the dayside than cases with larger thermal inertia values (80 TIU). On the nightside, however, the effect seems to be the opposite. This means that for low thermal inertia values, the surface layers heat very quickly when illuminated, but also cool faster when solar radiation decreases. The surface temperature decreases with longitude, so that between 180° and 270°, they have cooled down to approximately the same temperature at the sunrise terminator.

A more quantitative overview of the variation of temperatures over one comet rotation at the equator is shown in Fig. 7. A sub-solar point of 0° is equivalent to midday at ~3.1 h, while 90° and 270° longitudes correspond to the sunset terminator (~6.2 h) and the sunrise terminator (0 h), respectively. Figure 7a illustrates the daily evolution of the surface temperature, the modulation of which is similar on the dayside for all cases, with its maximum shortly after midday. Case A0 has no dust mantle at the top of the sublimation front, therefore its maximum temperature is below 200 K (at the free sublimation temperature of H2O) and its maximum is at midday. The amplitude of the diurnal temperature curves is larger for cases with low thermal inertia values. On the other hand, cases with larger thermal inertia values (C1 and C2) have lower surface temperatures on the dayside.

The evolution of the temperature at the sublimation front of H2O is shown in Fig. 7b. The maximum temperature is always on the dayside. After the sunset terminator the temperature decreases slowly with time, except for the coupled case B2(C), in which the temperature at the H2O-front decreases much faster due to the presence of the CO2 outgassing at deeper layers of the nucleus. The case with a larger ρc value (B4) has the smallest range of H2O temperatures during day because the thermal diffusivity α is smaller in this case in comparison to case B2, which has the same thermal inertia value and depth for both sublimation fronts. Therefore, the lower the value of α, the slower the rate of transfer of heat within the material. Case C2 also has a small H2O temperaturedifference between dayside and nightside, but in this case it is caused by the combination of a larger thermal inertia value and sublimation fronts two times deeper. We see in these models expected trends. Temperatures at an H2O sublimation front close to the surface (within 1 thermal skin depth δ) suggest dayside sublimation of H2O with only a small phase shift towards the terminator. In the higher thermal inertia case with larger front depth, the phase shift is larger.

Finally, the temperature evolution of CO2 is shown in Fig. 7c. It shows that CO2 activity is coming exclusively from the dayside if the sublimation front is very close to the surface (B1) for a low thermal inertia value (40 TIU). Once we increase the depth of the sublimation front from 1 cm to 5 cm and double the thermal inertia value (C1), the CO2-ice also has its maximum activity on the dayside, but it is shifted about 40 min towards sunset compared to case B1. For all other cases, the heat wave is maximum around sunset, triggering higher temperatures and fairly large production rates on the nightside. The evolution of these cases (except B4), follow the same trend, however, a nucleus with higher thermal inertia increases the thermal energy at the front, which translates into larger temperatures for case C2. We can therefore say that the temperature at the sublimation front of CO2, well below the nucleus surface, maximises close to the sunset terminator or even beyond for thermal inertia values expected at 67P/C-G.

thumbnail Fig. 4

Polar view of the map of the surface temperature for all cases listed in Table 2. Diagram at the bottom-right of the figure indicates the planetographic coordinates.

thumbnail Fig. 5

Polar view of the map of the temperature at the sublimation front of H2O for all caseslisted in Table 2. Diagram at the bottom-right of the figure indicates the planetographic coordinates.

thumbnail Fig. 6

Polar view of the map of the temperature at the sublimation front of CO2 for all cases listed in Table 2. Diagram at the bottom-right of the figure indicates the planetographic coordinates. Cases A0 and B0 have no CO2 emissions, therefore they are shown in white.

thumbnail Fig. 7

Diurnal change of the temperature at the equator of the surface (a), at the sublimation front of H2O-ice (b) and at the sublimation front of CO2-ice (c). The parameters of the different models are presented in Table 2. The dotted vertical lines indicate midday, the sunset terminator and midnight.

4.2 Diurnal evolution of flux at the surface

Maps of the production rates per unit area at the sublimation fronts are shown in Figs. 8 and 9 for H2O and CO2 gas sources, respectively. We can see there is a strong correlation between the flux and the temperatures at the sublimation fronts (Figs. 5 and 6) and that a relatively small change in temperature can trigger substantial activity for H2O and CO2 towards the nightside of the comet compared to a case that does not include the effect of thermal conductivity. We also notice the strong decrease in the flux with latitude, such that at the poles, the H2O and CO2 activity decreases by several orders of magnitude.

The maximum activity is produced on the dayside with relatively small differences for all H2O cases, while for CO2 activity is only strong at midday when ice sources are very close to the surface (B1) or when the thermal inertia of the terrain is 80 TIU (C1). In both cases, the energy budget at the sublimation front of CO2 is larger compared to other cases and, therefore, they have slightly larger fluxes. In all other cases, CO2 activity is maximum at the sunset terminator at 90° longitude or beyond. A larger CO2 flux in the B3 case has been set on purpose in order to study the effect of larger CO2 to H2O mixing ratios, although we do see from Fig. 6 that the temperature was kept the same as in the B2 case. A decrease in the thermal diffusivity α of the nucleus shifts the maximum CO2 activity close to midnight, making B4 the only case with strong CO2 emission exclusively on the nightside of the comet.

thumbnail Fig. 8

Polar view of the map of the production rate per unit area at the sublimation front of H2O for all caseslisted in Table 2. Diagram at the bottom-right of the figure indicates the planetographic coordinates.

thumbnail Fig. 9

Polar view of the map of the production rate per unit area at the sublimation front of CO2 for all caseslisted in Table 2. Diagram at the bottom-right of the figure indicates the planetographic coordinates.

4.3 Calculation with the direct simulation Monte Carlo method

4.3.1 Gas flow fields

The DSMC results presented in Figs. 10 and 11 are slices through the 3D simulation domain on the xy-plane (i.e. equatorial plane) for the number densities for all tested cases for H2O and CO2, respectively.The view is from above the north pole of the nucleus. The comet’s rotation is anti-clockwise and the sub-solar point is located at 0° longitude, with all other longitudes being indicated in the bottom right side of both figures. Both figures show the changes in the distribution of the flow field around the nucleus as a consequence of inter-molecular interactions within the mixture of gases. We can also see the changes in the gas distribution away from the sub-solar point and therefore the effect of the thermal inertia on the sublimation of H2O-ice and CO2-ice. In order to quantify such changes, we have also extracted the number density along the circle at the edge of the simulation domain in steps of 10° on the xy-plane of the DSMC results to obtain its diurnal variation in the equatorial plane as shown in Fig. 12. We notice, for example, that compared to the reference case (A0), the single-species case (B0) with a thermal inertia of 40 TIU shows a very similar outgassing pattern, but has a shift by about 20° in longitude with respect to the sub-solar point. This is also true for all multi-species cases with non-zero thermal lag.

In terms of CO2 number densities, however, the distribution of activity is more heterogeneous and can influence the number of H2O particles transported towards the nightside of the comet. Cases B1 and C1, for example, show that if CO2 outgassing is predominantly from the dayside, then H2O can be easily transported towards the nightside and the abundances measured above the anti-solar point of the nucleus could be equal to or higher for H2O than for CO2 compared to the other multi-species cases. However, if the sublimation front of CO2 is set at deeper layers for the same thermal inertia value (cases B2–4), the maximum activity of CO2 is shifted beyond the sunset terminator and the H2O number densities can decrease by at least one order of magnitude in the region above the anti-solar point. If we increase the production rate mixing ratio by a factor of 2 (case B3), we get a larger region on the nightside in which H2O is effectively “blocked” compared to case B2. However, in terminator orbits similar to those frequently used by the Rosetta spacecraft at 67P/CG, these differences would not have been detectable.

In Figs. 10 and 11, we can compare the effect of doubling the depth of the sublimation fronts for the high thermal inertia cases. For H2O the effect is minimal on the dayside, however, on the nightside, the lower number density values are linked to the maximum abundance of CO2 in this region. In other words, we obtain a stronger shift in the maximum CO2 outgassing for fronts at deeper layers, which influence the distribution abundance of H2O on the nightside of the comet. This is intuitively correct based on our understanding of how thermal inertia and the depth of the sublimation front influence the surface production.

thumbnail Fig. 10

Slice of the 3D simulation domain on the xy-plane, with information of number density within the flow for all cases listed in Table 2. The dash-dotted lines indicate the position of the field of view along the 0°, 90°, 180°, and 270° in longitude.

thumbnail Fig. 11

Slice of the 3D simulation domain on the xy-plane, with information of number density of CO2 for all cases listed in Table 2. The dash-dotted lines indicate the position of the field of view along the 0°, 90°, 180°, and 270° in longitude.

thumbnail Fig. 12

Temporal and longitudinal variation in the H2O and CO2 number densities obtained with DSMC calculations at 8 km from the surface in the equatorial plane. The dotted dark vertical lines indicate the location of the midday, the terminator with the sunset at 90° and midnight.

4.3.2 Diurnal evolution of mixing ratios

Number densities in Fig. 12 can be used to define a CO2 to H2O number density mixing ratio (CO2/H2O), as illustrated in Fig. 13 for all multi-species cases used herein. The maximum mixing ratio is close to local midnight in all cases. Case B1 only shows a weak variability with local time as one goes around the nucleus – but even here, the maximum mixing ratio is on the nightside.

Cases with the same boundary conditions, but set with different coupling modes (B2 and B2(C)) have the same evolution trend. However, the overall mixing ratio is smaller for the coupled case than for the decoupled case. Using Fourier’s law of thermal conduction in its one-dimensional form, we can calculate the heat flux density (qz) at each sublimation front: (9)

and estimate the heat distribution between the two layers for both coupling modes. Figure 14 shows the diurnal evolution of qz at the H2O and CO2 sublimation fronts for cases B2 and B2(C), from which we can notice how the presence of CO2-ice at deeper layers decreases the heat flow at the H2O sublimation front around midday, and therefore, the energy budget for H2O sublimation. On the nightside, the surface layer gets colder faster compared to the subsurface layers of the nucleus, therefore, the variation of the heat flux becomes negative. Similarly to what happens around midday, the difference in heat flux between the decoupled and the coupled mode for the H2O sublimation front increases about 1h after sunset.

By increasing the thermal inertia value from 40 TIU to 80 TIU, we see that for the same front depths (C1 compared to B2), the amplitude of the mixing ratio variation is smaller and occurs earlier, because the CO2 outgassing has a slightly higher production rate on the dayside of the comet. However, if one also doubles the depth of the sublimation front for both molecules, the maximum mixing ratio is after midnight due to the strong delay in the thermal wave reaching deeper layers.

thumbnail Fig. 13

Variation in the CO2 to H2O number density mixing ratio obtained with DSMC calculations at 8 km from the surface with longitude at the equator for all cases with multi-species gas sources. The dotted dark vertical lines indicate the location of the midday, the terminator with the sunset at 90°, and midnight.

thumbnail Fig. 14

Diurnal variation of the heat flux density at the H2O and CO2 sublimation fronts for the decoupled case B2 and the coupled case B2(C). The dotted dark vertical lines indicate the location of the midday, the terminator with the sunset at 90°, and midnight.

thumbnail Fig. 15

Slice of the 3D simulation domain on the xy-plane, with information of temperature (upper row), speed (middle row), and stream traces (bottom row) within the flow for model A0, B0, and model B2. The dash-dotted lines indicate the position of the field of view along the 0°, 90°, 180°, and 270° in longitude (see Figs. 10 and 11).

4.3.3 Diurnal evolution of gas temperatures and velocities in the coma

Figure 15 shows a polar view of the gas temperature (upper row), speed (middle row), and the stream traces (bottom row) of the flow field, with all data points (x,y,0) in the simulation domain for cases A0, B0, and B2. A simple comparison of cases A0 and B0 shows that including thermal inertia in the model also affects the temperature and velocity distributions of H2O about the nucleus. Similarly to the results in Fig. 10, in terms of temperature, H2O responds strongly to illumination, but there is a small shift in the same direction of the nucleus’ rotation when including thermal inertia. This can result in different temperatures very close to the nucleus when taking measurements in terminator orbits.

In terms of flow velocities, we see that when CO2 interacts with H2O (case B2), the velocities decrease quickly (by ~200 m s−1) close to the 130° and 260° longitude positions. Results for both molecules show a very similar behaviour towards the dayside, which is an indication of the strong effect of H2O increasing the nominal speed of CO2 in this region through inter-molecular collisions. On the nightside, however, the flow speed of H2O seems to be lower than for CO2 presumably as a result of the lateral flow and low collision frequencies. The stream traces of the flow velocity in Fig. 15 shows how the presence of CO2 affects the behaviour of the H2O gas field towards the nightside.

We can also extract information on the average temperature and speed of the whole flow at the edge of the simulation domain, as has been done before for the number densities. Figure 16 illustrates both variables at 8 km from the surface of the nucleus as a function of longitude (equivalent to local time). Cases where most of the flow is focussed towards the dayside of the comet (A0, B0, B1, and C1), have their maximum gas temperatures around midnight. This can be explained by the fact that in these cases, gas pressure is lower on the nightside and H2O gas can be transported more easily to the nightside compared to the cases with a larger CO2 presence. Therefore, the flow in cases that have a stronger dichotomy between day and night sides interact around midnight, which slows down the speed of the flow and increases the thermal energy in this region. All other cases have lower gas temperatures, which are linked to larger abundances of CO2 at midnight (see Fig. 13). Although larger temperatures on the nightside of the nucleus are not intuitive, this is a consequence of the smaller number of gas particles and the presence of CO2, which is heavier thanH2O and, therefore, is less effective transforming thermal energy into kinetic energy. Similarly to the effect we see in the mixing ratios, comparing cases B2 and B3 shows that the increase of a small amount of CO2 can also decrease the temperature of the flow by about 10 K on the nightside. This is caused by a smaller abundance of H2O gas in this region for case B3. The difference between cases in decoupled and coupled mode (B2 and B2(C)) is minimal in terms of temperature differences. We can also see these two cases have a very similar trend to case C2, showing that there are multiple combinations of model parameters than can lead to similar temperature results. Hence, it is difficult to unambiguously establish the physical properties of the nucleus from these quantities through measurements. The case with a larger ρc value (B4) has a smaller variation of temperature with time at 8 km from the surface compared to case B2, which has the same thermal inertia value and front depths. The difference is minimal on the illuminated areas of the comet, however, bigger differences arise closer to the sunrise terminator because the initial temperatures at the CO2 sublimation front have a larger delay with longitude.

In the panel on the bottom of Fig. 16, we see the variation of the speed of the mixed flow at 8 km from the surface over one comet rotation. On the dayside, all the cases have very similar speeds. In contrast, on the nightside, most multi-species cases have fluxes with speeds between 100 and 200 m s−1 slower than the single-species cases due to the higher abundance of CO2, which has a molecular mass more than two times bigger than H2O. For case B1, the speed of the flow on the nightside does not decrease as much because H2O molecules can travel more easily to this region. For case C1, the range in which the flow speed is small is much thinner than the other cases, which can be understood by the distribution of CO2 molecules around the nucleus that is shown is Fig. 11 for this case. Once again, a delay of between 20–30° is observed for case B4 when looking at the speed of the flow.

thumbnail Fig. 16

Averaged temperature and velocity of the gas mixture obtained with DSMC calculations at 8 km from the surface on the xy-plane (Fig. 15) as a function of the longitude for each of the cases (Table 2). The vertical dashed lines represent the location of the terminator with the sunset at 90° and the sunrise at 270°.

thumbnail Fig. 17

Ratio between the rotational temperature Trot and the translational temperature Ttrans of H2O for all cases listed in Table 2. The dash-dotted lines indicate the position of the field of view along the 0°, 90°, 180°, and 270° in longitude.

4.3.4 Non-thermal equilibrium factor Trot/Ttrans

We calculated the ratio between the rotational to translational temperatures (Trot/Ttrans) for all the cases in order to quantify the degree of equilibrium inside the gas flow of H2O in Fig. 17 and CO2 in Fig. 18. Here, Trot/Ttrans = 1 indicates regions with translation-rotational (T-R) equilibrium. Intermediate and high production rates are significantly better at transforming the energy, therefore, we see that this is only true on the dayside of the comet. On the nightside of the comet, the rotational energy dominates, which manifests a TR nonequilibrium of the gas flow.

Given that the production rate of H2O is higher than the production rate of CO2, it is reasonable to assume that the region of thermal equilibrium for CO2 shown in Fig. 18 is also on the dayside. However, this is smaller than for H2O, because the CO2 molecules are heavier, which decreases the energy transfer between gas particles. On the nightside of the comet, CO2 is more abundant (see Fig. 11), which slightly increases the energy exchange in the field of the gas mixture.

Either in single gas fields or in gas mixtures, we notice a dichotomy between the sunrise and sunset terminators in terms of energy exchange. The sunset terminator is in local thermal equilibrium (LTE) compared to the sunrise terminator, where the translational temperature is significantly lower than the rotational temperature. This has certain implications for the interpretation of the gas temperatures along different lines of sight, namely, that regions in non-LTE are more likely to experience rough changes in temperature in relatively small altitude ranges.

thumbnail Fig. 18

Ratio between the rotational temperature, Trot, and the translational temperature, Ttrans, of CO2 for all cases listed in Table 2. The dash-dotted lines indicate the position of the field of view along the 0°, 90°, 180°, and 270° in longitude.

4.3.5 Artificial heat transfer from the dust mantle to the gas

As mentioned in Sect. 3, one extra case was tested (with Tsurface as the initial temperature at which gas particles leave the surface) in order to include the heat exchange between the dust mantle and the gas particles released after sublimation. This was compared to case B2. Figure 19 shows how increasing the initial gas temperature at the sublimation front increases the mean free path of the gas particles close to the surface. The absolute difference between both thermal inputs varies from 30 to 1100 m, the effect being smaller around the sub-solar point (0°), than in poorly illuminated areas (90°, 180°, and 270°), where a small difference in temperature can significantly affect the way in which gas particles are released from the nucleus.

Vertical profiles of number density, Doppler velocity and the weighted mean of the rotational and translational temperatures or “overall kinetic temperature” (Bird 1994) of the gas above the sub-solar point for these cases are shown in Fig. 20a. In terms of number density, a higher initial gas temperature does not make any difference. However, the change in velocity is approximately +60 m s−1 for H2O gas. For CO2 gas, the change is +90 m s−1 close to the nucleus and 60 m s−1 at 8 km from the surface. Assuming that the heat exchange between the dust mantle and the releasedgas molecules is 100% efficient, we can also determine from these cases that the H2O gas flow canbe about 100 K warmer at 8 km from the surface when looking exactly above the sub-solar point. For CO2, however, there is a slight difference at high altitudes.

In terminator orbits, such differences are difficult to detect for the same model parameters. In Fig. 20b, we have plotted the average number density, Doppler velocity, and the overall temperature for the mixed flow at 90° and 270° longitude. Similarly to Fig. 20a, we see no change in terms of number densities for the same viewing geometry and different thermal inputs. When looking at flow velocities, we can see slightly bigger differences far from the nucleus. However, these are less than 30 m s−1 and 50 m s−1 for the sunset terminator and sunrise terminator, respectively. Finally, our results suggest that in terminator orbits, small changes in the surface temperature can generate significant differences at 8 km from it. At the sunset terminator, a temperature difference of ~10 K, can increaseto about 70 K at large altitudes. Meanwhile, for the sunrise terminator, the same difference at the surface can produce temperature differences in the coma up to 90 K at 8 km from the surface.

thumbnail Fig. 19

Polar view of the mean free path of both thermal inputs Tfront and Tsurface at the nucleus. Diagram on the right side of the figure indicates the planetographic coordinates, where 0° indicates the sub-solar point.

thumbnail Fig. 20

Vertical profiles of number density, Doppler velocity and the overall temperature with altitude for different observation geometries. Case B2 was tested with two different boundary temperatures at the surface: the temperature at the sublimation front and the temperature at the surface, which is about 100 K greater than the other at the sub-solar point.

5 Discussion

We have found that thermal inertia produces moderate changes, which are potentially measurable, to the distribution of H2O emissions. However, CO2 is more volatile and its surface emission distribution can be substantially changed by the assumed properties of the surface layer. For cases with sublimation fronts at large depths or relatively low thermal inertia values (40 TIU), CO2 activity was mainly focussed on the nightside of the comet, which is in agreement with the inferences made by Bockelée-Morvan et al. (2015) and Gerig et al. (2020). Therefore, the choice of different model parameters can impact the interpretation of remote sensing and in situ measurements significantly.

Several authors have used different approaches to reproduce local measurements of number density within the coma of comet 67P/CG taken by ROSINA. However, when studying multiple data sets from other instruments on board the Rosetta spacecraft, Marschall et al. (2019) found a discrepancy between the model results and the measurements of temperature and speed along the MIRO line of sight. Given that many of the Rosetta orbits were on the terminator plane, our results show that thermal lag can produce differences in the distribution of the gas sources at the surface, which can be detectable at large distances from the nucleus when looking above the terminator. However, a similar analysis needs to be done for the complex shape of comet 67P/CG (Preusker et al. 2017), in which in situ ROSINA/DFMS and MIRO measurements of number densities, temperatures, and speeds along the viewing geometry are used to constrain the activity distribution at the nucleus. This can also be used to reproduce patterns of dust activity triggered by the gas as it is observed in OSIRIS images.

Given the large complexity of the nucleus structure and composition, the use of a 1D thermal model for the energy transport within the nucleus only gives us a rough estimate of the distribution of ice sources. Energy losses due to lateral expansion or re-condensation would also play a role in gas diffusion through the dust matrix. However, compared to the radial diffusion of the flow, we assume their contribution to be very small. The cases in the current paper are possibly extreme scenarios of the real conditions in cometary nuclei. However, their results help us to constrain the thermophysical properties used when modelling a single- and multi-species cometary coma.

6 Conclusions

We built a model to simulate cometary outgassing, which serves as the basis for our detailed study of the influence of CO2–H2O mixtures on the distribution of gas activity around a spherical cometary nucleus. We used a simplified thermal model to provide input to a DSMC code to study the properties of the gas flow under several conditions. We see evidence of the strong effect of thermal inertia on the distribution of the gas sources at the surface. This is especially important for contributions to nightside activity. The interaction between H2O and CO2 plays a significant role in the overall flow dynamics. We enumerate the main findings of our study in the following points:

  1. Sublimation fronts for H2O and CO2 within the first cm depth and low thermal inertia values result in most of the activity being focussed towards the illuminated regions of the comet. In order to produce substantial nightside activity, CO2-ice needs to be deeper than H2O-ice sources to produce a delay in the outgassing pattern towards sunset. CO2 fronts atdepths greater than 5cm and 10cm and thermal inertia values of 40 TIU and 80 TIU, respectively, produce their maximum activity shortly after sunset, which leads to strong emission from the non-illuminated regions of the comet;

  2. Low thermal diffusivity values can delay strongly the CO2 emissions from the sub-surface of the nucleus, locating its maximum activity close to midnight when CO2-ice sources are deeper than 5 cm;

  3. H2O and CO2 activity decrease significantly with latitude. Emissions of H2O at the poles are very weak compared to emissions at the equatorial plane. For CO2 emissions, this difference is less pronounced, given that CO2-ice requires less heat to sublimate than H2O-ice, even if it is located at a significant depth;

  4. A comparison between purely insolation-driven models with different thermal inertia values shows that by increasing these values, it is possible to shift the distribution of H2O emission towards the evening terminator by about 20° in longitude from the sub-solar point. This results in considerable differences in number densities at positions far from the sub-solar point when compared to purely insolation driven models (Bieler et al. 2015b; Marschall et al. 2019). In terms of the temperature and speed of the H2O flux, no large changes are detectable. For CO2-ice, cases with different thermal inertia values and sublimation fronts can display very different outgassing distributions. These can have their maxima beyond the terminator and, thus, make CO2-ice the main source of nightside activity (as proposed by Bockelée-Morvan et al. (2015) to explain VIRTIS observations of CO2 and by Gerig et al. (2020) to explain dust observations made using OSIRIS);

  5. Number density mixing ratios show that H2O is the most abundant species on the illuminated side of the comet as expected. There is also a strong effect of CO2 activity on the distribution of the H2O flow field on the nightside, which can decrease the amount of H2O molecules per cubic metre by at least one order of magnitude compared to a pure H2O case. This arises from CO2 limiting lateral expansion of the H2O;

  6. On average, CO2 gas decreases the flow velocities on the nightside of the comet. In the cases we studied, velocities were between 100 and 200 m s−1 slower than the cases without CO2 activity. This is connected to the larger abundance in this region of CO2. However, on the dayside, CO2 has no significant effect on the temperature and speed of the flow;

  7. There isa dichotomy between the sunrise and sunset terminators in terms of thermal equilibrium once we include the effect of thermal conductivity for the distribution of ice sources. The energy transfer is less effective in CO2 than in H2O and it depends strongly on the production rates along the line of sight;

  8. Models with different energy inputs for the gas release from the surface only show significant differences at large distances from the surface when looking above the sub-solar point. Assuming that the dust mantle on top of the sublimation front can heat the gas by 100K as it leaves the nucleus, we can measure flow speeds 60–90 ms−1 faster and temperatures of H2O gas around 100K larger within the first 10 km from the nucleus. However, in terminator orbits, differences can only be seen in terms of the flow temperature, which can be between 70 and 90 K warmer. For an ideal adiabatic expansion of the gas, these differences could be detected at larger distances from the nucleus.

In summary, the CO2 emission distribution may be markedly different from H2O and can also influence the H2O gas distribution pattern. For non-zero thermal inertia and sublimation from depth, nightside gas emission will be non-negligible and is a potential source of dust above the nightside hemisphere of 67P/CG (Gerig et al. 2020).

Acknowledgements

This work has been carried out within the framework of the National Centre of Competence in Research PlanetS supported by the Swiss National Science Foundation. The authors acknowledge the financial support of the SNSF. Raphael Marschall acknowledges the support from the Swiss National Science Foundation grant P2BEP2 184 482.

References

  1. A’Hearn, M., Millis, R. C., Schleichter, D. G., Osip, D. J., & Birch, P. V. 1995, Icarus, 118, 223 [CrossRef] [Google Scholar]
  2. Alexeenko, A. A., Ganguly, A., & Nail, S. L. 2009, J. Pharm. Sci., 98, 3483 [Google Scholar]
  3. Bardyn, A., Baklouti, D., Cottin, H., et al. 2017, MNRAS, 469, S712 [NASA ADS] [CrossRef] [Google Scholar]
  4. Barucci, M., Filacchione, G., Fornasier, S., et al. 2016, A&A, 595, A102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Bieler, A., Altwegg, K., Balsiger, H., et al. 2015a, Nature, 526, 678 [Google Scholar]
  6. Bieler, A., Altwegg, K., Balsiger, H., et al. 2015b, A&A, 583, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Bird, G. A. 1994, Molecular Gas Dynamics And The Direct Simulation Of Gas Flows (Oxford: Clarendon Press) [Google Scholar]
  8. Biver, N., Hofstadter, M., Gulkis, S., et al. 2015, A&A, 583, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Bockelée-Morvan, D., Debout, V., Erard, S., et al. 2015, A&A, 583, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Bonev, B. P., Mumma, M. J., Villanueva, G. L., et al. 2007, ApJ, 661, L97 [NASA ADS] [CrossRef] [Google Scholar]
  11. Bonev, B. P., Villanueva, G. L., Paganini, L., et al. 2013, Icarus, 222, 740 [Google Scholar]
  12. Bouvier, A., & Wadhwa, M. 2010, Nat. Geosci., 3, 637 [NASA ADS] [CrossRef] [Google Scholar]
  13. Brandt, J. C. 2007, CHAPTER 30, Physics and Chemistry of Comets, eds. L.-A. McFadden, P. R. Weissman, T. V. Johnson, (Academic Press), 557 [Google Scholar]
  14. Brouet, Y., Levasseur-Regourd, A. C., Sabouroux, P., et al. 2016, MNRAS, 462, S89 [NASA ADS] [CrossRef] [Google Scholar]
  15. Christou, C., Dadzie, K., Thomas, N., et al. 2018, PSS, 161, 57 [Google Scholar]
  16. Combi, M., & Smyth, W. 1988, ApJ, 327 [Google Scholar]
  17. Combi, M., Shou, Y., Fougere, N., et al. 2020, Icarus, 335, 113421 [NASA ADS] [CrossRef] [Google Scholar]
  18. Crifo, J., Lukianov, G., Rodionov, A., Khanlarov, G., & Zakharov, V. 2002, Icarus, 156, 249 [NASA ADS] [CrossRef] [Google Scholar]
  19. Crifo, J., Loukianov, G., Rodionov, A., & Zakharov, V. 2003, Icarus, 163, 479 [CrossRef] [Google Scholar]
  20. Crifo, J., Loukianov, G., Rodionov, A., & Zakharov, V. 2005, Icarus, 176, 192 [NASA ADS] [CrossRef] [Google Scholar]
  21. Davidsson, B., Gutiérrez, P., Groussin, O., et al. 2013, Icarus, 224, 154 [NASA ADS] [CrossRef] [Google Scholar]
  22. Emerich, C., Lamarre, J. M., Moroz, V. I., et al. 1988, Temperature and size of the nucleus of comet P/Halley deduced from IKS infrared Vega 1 measurements, (Berlin, Heidelberg: Springer), 839 [Google Scholar]
  23. Fanale, F., & Salvail, J. 1984, Icarus, 60, 476 [NASA ADS] [CrossRef] [Google Scholar]
  24. Fernández, Y. R., Kelley, M. S., Lamy, P. L., et al. 2013, Icarus, 226, 1138 [Google Scholar]
  25. Festou, M., Keller, H., & Weaver, H. 2004, A brief conceptual history of cometary science in Comets II, (University of Arizona Press), 3 [Google Scholar]
  26. Filacchione, G., De Sanctis, M. C., Capaccioni, F., et al. 2016, Nature, 529, 368 [NASA ADS] [CrossRef] [Google Scholar]
  27. Fink, U., Doose, L., Rinaldi, G., et al. 2016, Icarus, 277, 78 [Google Scholar]
  28. Finklenburg, S., & Thomas, N. 2014a, PSS, 93, 71 [NASA ADS] [Google Scholar]
  29. Finklenburg, S., Thomas, N., Su, C., & Wu, J. 2014b, Icarus, 236, 9 [Google Scholar]
  30. Fornasier, S., Hasselmann, P. H., & Barucci, M. 2015, A&A, 583, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Fornasier, S., Mottola, S., Keller, H., et al. 2016, Science, 354, 1566 [NASA ADS] [CrossRef] [Google Scholar]
  32. Fougere, N., Altwegg, K., Berthelier, J., et al. 2016a, A&A, 588, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Fougere, N., Altwegg, K., Berthelier, J.-J., et al. 2016b, MNRAS, 462, S156 [NASA ADS] [CrossRef] [Google Scholar]
  34. Gerig, S.-B., Pinzón-Rodríguez, O., Marschall, R., Wu, J., & Thomas, N. 2020, Icarus, 351, 113968 [Google Scholar]
  35. Groussin, O., Sunshine, J., Feaga, L., et al. 2013, Icarus, 222, 580 [Google Scholar]
  36. Hoang, M., Altwegg, K., Balsiger, H., et al. 2017, A&A, 600, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Hässig, M., Altwegg, K., Balsiger, H., et al. 2015, Science, 347, aaa0276 [Google Scholar]
  38. Huebner, W., Benkhoff, J., Capria, M.-T., et al. 2006, Heat and Gas Diffusion in Comet Nuclei (Noordwijk: ESA Publication Division) [Google Scholar]
  39. Keller, H. U., Mottola, S., Skorov, Y., & Jorda, L. 2015, A&A, 579, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Kofman, W., Herique, A., Barbin, Y., et al. 2015, Science, 349, 6247 [CrossRef] [Google Scholar]
  41. Lamy, P., Toth, I., A’Hearn, M., Weaver, H., & Weissman, P. 2001, Icarus, 154, 337 [NASA ADS] [CrossRef] [Google Scholar]
  42. Le Roy, L., Altwegg, K., Balsiger, H., et al. 2015, A&A, 583, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Lebofsky, L., & Spencer, J. 1989, in Radiometry and thermal modeling of asteroids, Asteroids II, eds. R. P. Binzel, T. Gehrels, T., & M. S. Matthews, 128 [Google Scholar]
  44. Levasseur-Regourd, A., Hadamcik, E., Desvoivres, E., & Lasue, J. 2009, PSS, 57, 221 [Google Scholar]
  45. Lhotka, C., Reimond, S., Souchay, J., & Baur, O. 2015, MNRAS, 455, 3588 [Google Scholar]
  46. Li, J.-Y., A’Hearn, M., Mcfadden, L., & Belton, M. 2007, Icarus, 188, 195 [CrossRef] [Google Scholar]
  47. Liao, Y., Su, C. C., Marschall, R., et al. 2016, Earth, Moon Planets, 117, 41 [NASA ADS] [CrossRef] [Google Scholar]
  48. Liao, Y., Marschall, R., Su, C., et al. 2018, PSS, 157, 1 [Google Scholar]
  49. Marschall, R., Su, C., Liao, Y., et al. 2016, A&A, 589, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. Marschall, R., Mottola, S., Su, C., et al. 2017, A&A, 605, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Marschall, R., Rezac, L., Kappel, D., et al. 2019, Icarus, 328, 104 [Google Scholar]
  52. Marschall, R., Liao, Y., Wu, J.-S., & Thomas, N. 2020, Icarus, 346, 113742 [Google Scholar]
  53. Marshall, D., Groussin, O., Vincent, J.-B., et al. 2018, A&A, 616, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Mendis, D., & Brin, G. 1977, The Moon, 17, 359 [NASA ADS] [CrossRef] [Google Scholar]
  55. Migliorini, A., Piccioni, G., Capaccioni, F., et al. 2016, A&A, 589, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Oklay, N., Sunshine, J., Pajola, M., et al. 2016, MNRAS, 462, S394 [NASA ADS] [CrossRef] [Google Scholar]
  57. Pommerol, A., Thomas, N., Elmaarry, M. R., et al. 2015, A&A, 583, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  58. Preusker, F., Scholten, F., Matz, K.-D., et al. 2017, A&A, 607, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  59. Scanlon, T., et al. 2010, Comput. Fluids, 39, 2078 [Google Scholar]
  60. Schloerb, F., Keihm, S., Von Allmen, P., et al. 2015, A&A, 583, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Seiferlin, K., Spohn, T., & Benkhoff, J. 1995, Adv. Space Res., 15, 35 [NASA ADS] [CrossRef] [Google Scholar]
  62. Seiferlin, K., Spohn, T., Kömle, N., & Kargl, G. 1996, Planet. Space Sci., 44, 691 [CrossRef] [Google Scholar]
  63. Shoshany, Y., Prialnik, D., & Podolak, M. 2002, Icarus, 157, 219 [NASA ADS] [CrossRef] [Google Scholar]
  64. Skorov, Y., & Rickman, H. 1995, PSS, 43, 1587 [NASA ADS] [Google Scholar]
  65. Skorov, Y., Kömle, N., Markiewicz, W., & Keller, H. 1999, Icarus, 140, 173 [NASA ADS] [CrossRef] [Google Scholar]
  66. Skorov, Y., Kömle, N., Keller, H., Kargl, G., & Markiewicz, W. 2001, Icarus, 153, 180 [NASA ADS] [CrossRef] [Google Scholar]
  67. Skorov, Y., Keller, H., Jorda, L., & Davidsson, B. 2002, Earth, Moon Planets, 90, 293 [NASA ADS] [CrossRef] [Google Scholar]
  68. Smith, C. J., & Callendar, H. L. 1924, Proc. Roy. Soc. London A, 106, 83 [NASA ADS] [CrossRef] [Google Scholar]
  69. Spohn, T., Knollenberg, J., Ball, A., et al. 2015, Science, 349, 6247 [CrossRef] [Google Scholar]
  70. Sunshine, J., A’Hearn, M., Groussin, O., et al. 2006, Science, 311, 1453 [NASA ADS] [CrossRef] [Google Scholar]
  71. Thomas, N. 2009, PSS, 57, 1106 [NASA ADS] [Google Scholar]
  72. Thomas, N., Davidsson, B., El-Maarry, M. R., et al. 2015, A&A, 583, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  73. Tosi, F., Capaccioni, F., Capria, M., et al. 2019, Nat. Astron., 3, 649 [NASA ADS] [CrossRef] [Google Scholar]
  74. Whipple, F. 1950, ApJ, 111, 375 [NASA ADS] [CrossRef] [Google Scholar]
  75. Wu, J., & Lian, Y. 2003, Comput. Fluids, 32, 1133 [Google Scholar]
  76. Wu, J., & Tseng, K. 2005, Numer. Methods Eng., 63, 37 [NASA ADS] [CrossRef] [Google Scholar]
  77. Wu, J., Tseng,K., & Wu, F. 2004, Comput. Phys. Commun., 162, 166 [NASA ADS] [CrossRef] [Google Scholar]
  78. Zakharov, V., Rodionov, A., Lukianov, G., & Crifo, J. 2009, Icarus, 201, 358 [NASA ADS] [CrossRef] [Google Scholar]

1

Radiality , where are the components of the velocity in the x, y and z axes; d is the radial distance.

All Tables

Table 1

Physical parameters used in the model.

Table 2

List of tested parameters for model calculations.

All Figures

thumbnail Fig. 1

Vertical profiles of the temperature with depth for different insolation conditions on the equatorial plane of a comet with a spherical nucleus. Coupled cases with sublimation fronts at the same depth ( cm) are labelled as “Coupled 1,” while coupled cases with sublimation fronts at different depths (1.4 cm and 5.6 cm) are labeled as “Coupled 2.” Decoupled mode calculations are shown with dashed and dashed-dotted lines for H2O and CO2, respectively.

In the text
thumbnail Fig. 2

Decoupled calculation of the temperature at the surface (in black), at the sublimation front (in red) and the sublimation rate (dm/dt) of H2O and CO2 (with an EAF = 100%) as a function of the number of rotations of the nucleus for a thermal inertia value of 40 TIU for both species. Each symbol indicates a different observation geometry during one day. Circles indicate facets located at sunrise, squares are facets at midday, triangles are facets at sunset, and crosses are facets at midnight.

In the text
thumbnail Fig. 3

Sketch of the energy exchanges (heat sources and sinks) at the surface and sub-surface of a cometary nucleus: is the vector normal to the surface and θi is the solar incidence angle; T0 is the temperature at the surface of the nucleus; and are the temperatures at the sublimation fronts for the respective type of ice; and TN is the lower boundary temperature of the nucleus. The coupling mode is represented as a switch in between the H2O and CO2 sublimation fronts. The shadowed layers at and represent the selected effective active fractions (EAF) used at the sublimation front of each molecule.

In the text
thumbnail Fig. 4

Polar view of the map of the surface temperature for all cases listed in Table 2. Diagram at the bottom-right of the figure indicates the planetographic coordinates.

In the text
thumbnail Fig. 5

Polar view of the map of the temperature at the sublimation front of H2O for all caseslisted in Table 2. Diagram at the bottom-right of the figure indicates the planetographic coordinates.

In the text
thumbnail Fig. 6

Polar view of the map of the temperature at the sublimation front of CO2 for all cases listed in Table 2. Diagram at the bottom-right of the figure indicates the planetographic coordinates. Cases A0 and B0 have no CO2 emissions, therefore they are shown in white.

In the text
thumbnail Fig. 7

Diurnal change of the temperature at the equator of the surface (a), at the sublimation front of H2O-ice (b) and at the sublimation front of CO2-ice (c). The parameters of the different models are presented in Table 2. The dotted vertical lines indicate midday, the sunset terminator and midnight.

In the text
thumbnail Fig. 8

Polar view of the map of the production rate per unit area at the sublimation front of H2O for all caseslisted in Table 2. Diagram at the bottom-right of the figure indicates the planetographic coordinates.

In the text
thumbnail Fig. 9

Polar view of the map of the production rate per unit area at the sublimation front of CO2 for all caseslisted in Table 2. Diagram at the bottom-right of the figure indicates the planetographic coordinates.

In the text
thumbnail Fig. 10

Slice of the 3D simulation domain on the xy-plane, with information of number density within the flow for all cases listed in Table 2. The dash-dotted lines indicate the position of the field of view along the 0°, 90°, 180°, and 270° in longitude.

In the text
thumbnail Fig. 11

Slice of the 3D simulation domain on the xy-plane, with information of number density of CO2 for all cases listed in Table 2. The dash-dotted lines indicate the position of the field of view along the 0°, 90°, 180°, and 270° in longitude.

In the text
thumbnail Fig. 12

Temporal and longitudinal variation in the H2O and CO2 number densities obtained with DSMC calculations at 8 km from the surface in the equatorial plane. The dotted dark vertical lines indicate the location of the midday, the terminator with the sunset at 90° and midnight.

In the text
thumbnail Fig. 13

Variation in the CO2 to H2O number density mixing ratio obtained with DSMC calculations at 8 km from the surface with longitude at the equator for all cases with multi-species gas sources. The dotted dark vertical lines indicate the location of the midday, the terminator with the sunset at 90°, and midnight.

In the text
thumbnail Fig. 14

Diurnal variation of the heat flux density at the H2O and CO2 sublimation fronts for the decoupled case B2 and the coupled case B2(C). The dotted dark vertical lines indicate the location of the midday, the terminator with the sunset at 90°, and midnight.

In the text
thumbnail Fig. 15

Slice of the 3D simulation domain on the xy-plane, with information of temperature (upper row), speed (middle row), and stream traces (bottom row) within the flow for model A0, B0, and model B2. The dash-dotted lines indicate the position of the field of view along the 0°, 90°, 180°, and 270° in longitude (see Figs. 10 and 11).

In the text
thumbnail Fig. 16

Averaged temperature and velocity of the gas mixture obtained with DSMC calculations at 8 km from the surface on the xy-plane (Fig. 15) as a function of the longitude for each of the cases (Table 2). The vertical dashed lines represent the location of the terminator with the sunset at 90° and the sunrise at 270°.

In the text
thumbnail Fig. 17

Ratio between the rotational temperature Trot and the translational temperature Ttrans of H2O for all cases listed in Table 2. The dash-dotted lines indicate the position of the field of view along the 0°, 90°, 180°, and 270° in longitude.

In the text
thumbnail Fig. 18

Ratio between the rotational temperature, Trot, and the translational temperature, Ttrans, of CO2 for all cases listed in Table 2. The dash-dotted lines indicate the position of the field of view along the 0°, 90°, 180°, and 270° in longitude.

In the text
thumbnail Fig. 19

Polar view of the mean free path of both thermal inputs Tfront and Tsurface at the nucleus. Diagram on the right side of the figure indicates the planetographic coordinates, where 0° indicates the sub-solar point.

In the text
thumbnail Fig. 20

Vertical profiles of number density, Doppler velocity and the overall temperature with altitude for different observation geometries. Case B2 was tested with two different boundary temperatures at the surface: the temperature at the sublimation front and the temperature at the surface, which is about 100 K greater than the other at the sub-solar point.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.