EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 638 (June 2020) A&A, 638 (2020) A106 Full HTML
Open Access
Issue
A&A
Volume 638, June 2020
Article Number A106
Number of page(s) 16
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201936655
Published online 23 June 2020

© M. Hoang et al. 2020

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1 Introduction

The Rosetta mission studied the coma of comet 67P/Churyumov-Gerasimenko (67P) from August 2014 to September 2016. The in situ measurements revealed heterogeneities in terms of the three most abundant volatile species, H2O, CO2, and CO (Hässig et al. 2015; Hoang et al. 2017). A comparison between the illumination conditions during the mission and the change in relative abundance of CO2/H2O and CO/H2O ratios indicated a dichotomy between the northern hemisphere and the southern hemisphere. Furthermore, a dust mantle covering some parts of the surface of 67P surface was reported in Capaccioni et al. (2015), Sierks et al. (2015) and El-Maarry et al. (2015). Keller et al. (2015a, 2017) partly explained those observations by the extreme seasonal changes induced by the high obliquity (52°) of the rotational axis. The nucleus’s northern hemisphere experiences a long and soft summer while the southern hemisphere experiences a warmer but shorter summer, resulting in an intense sublimation of ices, erosion, and ejection of dust particles in the coma, which are partly redeposited in the northern hemisphere and form an insulating layer covering the surface.

Comets contain pristine material, as they are thought to have formed at the beginning of the solar system. Those that orbit in the inner solar system – such as 67P – are subject to solar heating and lose their outer layers through sublimation and erosion. Details of the comet’s nucleus remain poorly constrained, in particular regarding its composition, structure, and the depth of the pristine material (most probably present under a layer of alterated material). The measurements collected by Rosetta inside the coma provide information concerning the species inside the nucleus but the composition of the outgassing and the composition of the nucleus may be significantly different. Inferring the properties of the interior of the nucleus requires the use of numerical models or laboratory experiments. Numerical models havebeen developed to simulate the heat propagation and the gas diffusion through a porous matrix. Thermal models have improved considerably since the first 1D models, which considered a spherical body as a “fast rotator” (e.g. Espinasse et al. 1991; review in Prialnik et al. 2004 and Huebner et al. 2006). For example, Lasue et al. 2008 presented a quasi 3D model that took into account the latitudinal and longitudinal variation of illuminationbut no lateral heat transfer, and that assumed that the comet’s nucleus conducts very poorly at the scale considered. Lateral heat transfers were later introduced in a 3D thermal evolution model in Guilbert-Lepoutre et al. (2011).

The case of 67P was modelled before the Rosetta mission in De Sanctis et al. (2010), where the authors applied a quasi 3D model to the approximated shape of 67P, based on ground-based observations. They observed that both dust mantling and activity strongly depend on the shape and the obliquity of the object. Capria et al. (2017) applied the model to the real shape of 67P and illumination geometry on two locations of the big lobes, and showed that the activity remains in the first metre, which is consistent with the diurnal and orbital skin depth determined by the Microwave Instrument for the Rosetta Orbiter (MIRO; Gulkis et al. 2015). Hu et al. (2017a) analysed the influence of dust mantling and dust-to-ice ratio on the water outgassing. Rusol & Dorofeeva (2018) and Marov et al. (2019) studied the thermal evolution of 67P to investigate the depth of heated layers for twenty orbital cycles, specifically for the Ma’at region. They observed that the heat remains in a thin subsurface layer, even during a prolonged stay in the comet’s current orbit.

The Rosetta Orbiter Spectrometer for Ion and Neutral Analysis (ROSINA) experiment studied the composition and gas pressure of the 67P gas during the entire Rosetta mission, using three instruments: the Double Focusing Mass Spectrometer (DFMS), the Reflectron-type Time-Of-Flight mass spectrometer (RTOF) and the COmet Pressure Sensor (COPS; Balsiger et al. 2007). Coma models have been previously used to fit the ROSINA data (see Bieler et al. 2015b; Fougere et al. 2016), while Läuter et al. (2019) applied an inverse coma model to study the spatial distribution of gas sources from ROSINA DFMS/COPS measurements. Marschall et al. (2019) proposed a comparison of the ROSINA, MIRO, Visible InfraRed Thermal Imaging Spectrometer (VIRTIS) and Optical, Spectroscopic, and Infrared Remote Imaging System (OSIRIS) observations from May 2015 and an insolation-driven model constrained by MIRO’s measurements. Most of these models focused on the dominant species, H2O, and showed that a purely insolation driven model could not reproduce all coma measurements, leading to potentially heterogeneous activity depending on the regions. These models did not include the nucleus processes at work that induce the surface outgassing fluxes. Moreover, no model was compared with the ROSINA/RTOF dataset, whose characteristics allow for a study of the coma composition at a high time resolution.

In this work, we make a direct comparison between ROSINA’s RTOF coma measurements and the calculated surface fluxes from the thermophysical model introduced in Lasue et al. (2008) and De Sanctis et al. (2010), applied to the case of 67P. We study a large temporal period from 1 September 2014 until 15 April 2015. This approach allows us to link the nucleus processes to the coma measurements, but cannot include the gas expansion coma processes. This is why we also compare these results with the inverse coma model used by Läuter et al. (2019), who initially compared their model with ROSINA DFMS/COPS data. This enables a derivation of surface emissions rates from RTOF measurements without considering the complex nucleus processes.

2 Observational dataset

The RTOF was a time-of-flight mass spectrometer with a one-metre drift tube and an integrated reflectron. It had a mass resolution of m/δm = 500 at 50% level, a temporal resolution of 200 or 400 s depending on the operational mode, and a mass range of 1–300 u (Balsiger et al. 2007).

We used the densities of the main volatiles (H2O, CO2 and CO) recorded by RTOF during two pre-perihelion periods: Period A from 1 September to the 30 October 2014 and Period B from the 1 November 2014 until the 15 February 2015. During these periods, the comet approached the Sun from 3.45 to 3.05 au (during Period A) and from 3.05 until 2.4 au (during Period B). The illumination conditions slightly changed, with an averaged sub-solar latitude of 41° for Period A and of 33° for Period B (see detailed temporal and spatial evolution as well as more details about RTOF in Hoang et al. 2017, 2019). A comparison between this RTOF dataset and the thermophysical model will be presented in Sect. 4.3. In addition, we used the RTOF dataset to constrain the gas production on the surface of a cometary nucleus using an inverse coma model. This is presented in Sect. 4.4.

Moreover, we used some of the results obtained by other experiments (MIRO, VIRTIS, OSIRIS) to define the physical parameters of the comet and of cometary material needed for the nucleus modelling, as described in Sect. 3.2.1 and in Table 1. The MIRO gas velocity results (Gulkis et al. 2015) were also used for a comparison with the velocities derived from the surface temperature modelled (Sect. 4.3.1).

3 Thermophysical modelling

3.1 Description of the model

We simulated the thermophysical evolution of the cometary nucleus surface by using numerical simulations similar to the ones developed and described in Espinasse et al. (1991), Prialnik et al. (2004) and Huebner et al. (2006). The model is a 1D thermophysical model for simulating the evolution and the differentiation of an icy body, developed at the Istituto di Astrofiscica e Planetologia Spaziali (INAF-IAPS) in Rome and described in detail in Lasue et al. (2008) and De Sanctis et al. (2010). This mode assumed very low thermal inertia (no lateral heat transfer).

The model calculates the surface equilibrium temperature and mass balance considering the solar radiation reaching the surface, infrared energy emission, the sublimation of ices, and the transport of gases through the porous cometary material. The heat wave propagation in the interior leads to different processes: sublimation of ices, advection by gases, amorphous to crystalline ice transition, recondensation, and heat and mass transfer. Over time, as the more volatile species are sublimated, the initially homogeneous nucleus presents a layering delineated by the sublimation front of the different species and the deposit of large dust particles that are too slow to escape the gravity field on its surface. This evolution may happen until the comet loses all its volatiles or stops being active due to its insulating dust cover.

The nucleus is initially described as a homogeneous mixture of dust particles and ices, where water ice is initially in an amorphous form. The dust particles are embedded in a porous matrix of water ice and other minor species, like carbon dioxide and carbon monoxide. Species are possibly trapped in the amorphous ice and released during the crystallisation (transition from amorphous to crystalline ice). The 1D radial grid uses layers of increasing thickness with depth to ensure a good description of steep thermal gradients close to the surface. For a typical time step of 20–40 s, the code calculates the amount of ice sublimated and recondensed in each layer, by solving the coupled equations of heat transfer (Eq. (1)) and gas diffusion (described by the mass conservation Eq. (2)) using the finite difference method (see details in Capria et al. 1996, 2001; De Sanctis et al. 1999, and Huebner et al. 2006). If the amount of material in a layer is lower than a given value, the layer is added to the layer immediately below: a new layer with a new temperature, thickness, and composition is generated. Except for the first layer of the grid, which has a limit resolution of about 10−6 metre, the limit resolution of the layers underneath is one centimetre.

The heat transfer equation is ρc∂T∂t=∇(κ.∇T)+QH2O+QCO2+QCO+Qtr,\begin{equation*}\rho c \frac{\partial T}{\partial t} = \nabla(\kappa . \nabla T) + Q_{\textrm{H}_{2}\textrm{O}} + Q_{\textrm{CO}_{2}} + Q_{\textrm{CO}} + Q_{\textrm{tr}} ,\end{equation*}(1)

where ρ is the density, c the specific heat, T the temperature, t the time, κ the thermal conductivity, Qi the energy transfer due to sublimation and re-condensation of species i, and Qtr is the energy released at the amorphous to crystalline ice transition.

The mass conservation equation is given by ∂ρ∂t=−∇ϕ+Qgas,\begin{equation*}\frac{\partial \rho}{\partial t} = -\nabla \phi + Q_{\textrm{gas}} ,\end{equation*}(2)

where ϕ is the gas mass flux and Qgas is the gas sublimation term depending on the latent heat of sublimation for each species (see e.g. Huebner et al. 2006, and references therein). Numerical values for the physical parameters are discussed and given in Sect. 3.2.1, Table 1.

The solar input energy arriving at the surface is used in different processes: reflection, re-emission in the infrared, eventual sublimation of ices if present on the surface, and propagation in the interior. From Eq. (1), the penetration of the heat wave depends on the properties of the material (described by ρ, c, and κ). The surface of the nucleus experiences the highest variations of temperature, which decrease exponentially with the depth. The depth where the amplitude of the heat wave is attenuated by a factor 1/e is the thermal skin depth, and is defined by Eq. (3): D=κProtπρc,\begin{equation*}D = \sqrt{\frac{\kappa P_{\textrm{rot}}}{\pi \rho c}} ,\end{equation*}(3)

where Prot is the rotation period of the comet (to obtain the diurnal skin depth) or the revolution period of the comet (to obtain the orbital skin depth). MIRO estimated a diurnal skin depth of 67P of 1–2 cm and its orbital skin depth of about 1 m (Gulkis et al. 2015).

The temperature on the surface is calculated from the balance between the solar input and the four processes. The boundary conditions are determined by the equilibrium between solar illumination, conditions of temperature in the deep interior, and the gas pressure at the surface. The gas pressure condition at the surface is given by either being zero if a thick dust layer is covering the surface or being equal to the saturation pressure Psat for a surface of ice in equilibrium with vapour as described in Skorov et al. (2011). At the bottom of the grid, corresponding to a very large depth inside the nucleus, where the heat wave cannot reach, the variation in temperature is assumed to be zero.

For the surface layer, we evaluated, in addition, the quantity of dust particles that remain on the surface. Sublimation leads to gas and dust flows, which cause erosion of the surface and the formation of a dust mantle when large dust particles accumulate at the surface, with a thickness that evolves freely as a function of time. The accumulation of dust grains on the surface leads to the formation of a crust, which reduces the local activity and may extinguish it. The model allows different size distributions for dust particles typically chosen to be silicates. A critical radius, calculated from the equilibrium between gravity, drag, and centrifugal forces, determines if a given dust particle is ejected or redeposited on the surface of the comet. The thermal conductivity of the deposited dust crust is taken from the model of heat transport through porous material of Gundlach & Blum (2012). We applied this thermophysical model to the case of 67P to study the evolution of the surface fluxes of the main volatiles (H2O, CO2, and CO) along its orbit.

Table 1

Physical parameters constrained by the Rosetta’s instruments.

3.2 Application of the model

3.2.1 Physical parameters of the cometary material

The evolution and differentiation of the nucleus depend on a series of parameters that describe the initial conditions of the cometary material. Some of them have been measured by the Rosetta mission and are listed in Table 1.

In all computations discussed below and for each chosen location, the initial parameters are set as follows: the dust/ice ratio is four, the initial temperature is 30 K (compatible with amorphous ice and observations of very low temperature species detected by ROSINA, such as N2 (Rubin et al. 2015)and O2 (Bieler et al. 2015a)), the dust distribution size range is 10−6 –10−2 m with an average dust density of 103 kg m−3, and the initial pore radius is 0.0001 m, with a maximum size fixed at 0.05 m, based on GIADA and OSIRIS measurements (Rotundi et al. 2015; Fulle et al. 2015, 2016).

The purpose of this paper is to provide some constraints on the comet structure and properties for physical parameters that were poorly constrained by other instruments. The research is based on a comparison between a nucleus model and ROSINA/RTOF coma measurements. The parameters listed in Table 1 were constrained by various Rosetta instruments, and were considered as constants for our study. The main parameter influence that would be interesting to evaluate amongst them would be the dust to ice ratio, which was debated and may significantly influence the coma measurements. However, such a study is beyond the scope of this paper and will be included in future studies.

In this paper, we focus on the role of three distinct parameters that the Rosetta mission did not directly determine but did help to constrain: the ratios of the volatiles (CO2 /H2O and CO/H2), the thickness of the dust mantle, and the presence of trapped CO in the ice. We describe them below and report in Table 2 a list of the tested parameters used in different computations.

The first studied parameter is the ratios of the volatiles. The composition of the cometary gas certainly differs from the bulk composition, so we tested different values of CO/H2O ratio (0.0005, 0.002, and 0.02) to observe how they influence the CO2 /H2O and CO/H2O ratio of the outgassing.

The second parameter is the thickness of the dust mantle covering the northern hemisphere. The OSIRIS and VIRTIS results described a dark surface largely covered by organic-rich material, with only a small percentage of the surface covered by ice patches in the northern hemisphere (Capaccioni et al. 2015; Sierks et al. 2015; De Sanctis et al. 2015; El-Maarry et al. 2015). Biele et al. 2015 estimated a granular layer of about 20 cm at the first touchdown site, Agilkia. A dust deposit of about 1 metre resulted from the significant erosion of the southern hemisphere around perihelion (Hu et al. 2017b). Keller et al. (2017) described a seasonal mass transfer, which originates from the significant erosion of the southern hemisphere and a redeposition as a granular material layer in the northern hemisphere. We studied the influence of imposing a dust mantle from the start (it then evolves freely) with different thicknesses: 0 m (dust 0/1), 0.07 m (dust 2), 0.1 m (dust 3b/c), and 0.2 m (dust 3a) covering the northern hemisphere. To simulate the OSIRIS observations of a relatively thick layer of dust deposited on the northern hemisphere, the locations with positive latitudes allow the growth of a dust layer on the surface of the nucleus, while the locations in the southern hemisphere and at the equator do not allow the apparition of a dust mantle (cases ice 0/1/2).

The last studied parameter is the possibility that CO may be present as a trapped gas in the amorphous water ice. In this case, the trapped CO molecules would be released at the same time that the water ice sublimates, unlike frozen CO that has a sublimation temperature around 50–70 K (Huebner et al. 2006). We chose to include water as completely in amorphous ice at the beginning of the computation. Volatiles may be trapped in the amorphous ice and released during the transition from amorphous ice to crystalline ice. The trapping mechanism in cometary ices is still debated, as minor species could be trapped either in amorphous ice (Bar-Nun et al. 1988), or in clathrate hydrates (Delsemme & Swings 1952). Evidence for the presence of amorphous ice is the extremely low thermal conductivity (Kouchi et al. 1992), the enrichment of CO over N2 (Notesco & Bar-Nun 1996; Bar-Nun et al. 2007), and the observed outbursts probably caused by energy release during the phase transition, which is an exothermic reaction (Prialnik & Bar-Nun 1992; Jewitt 2009; Agarwal et al. 2017). Crystallisation may be endothermic, as in the case of impure amorphous ice (Kouchi & Sirono 2001). Luspay-Kuti et al. (2016) proposed, however, that the nucleus of 67P contains crystalline ice, clathrate, and other ices. They argued that the volatile gases trapped in amorphous ice are released simultaneously, which does not seem to be the case for the minor species observed by ROSINA/DFMS.

We selected 24 locations over the surface of the nucleus, to represent the different parts of the nucleus (northern hemisphere, equator, and southern hemisphere) and different geomorphological regions (identified and named in El-Maarry et al. 2016). In Fig. 1 we indicate the positions of the locations on a 2D map showing the regions, and we list the coordinates, the regions of the locations, the corresponding hemisphere, and the computed cases in Table 3. As reported, in our computations – listed in Table 2 –, we imposed a variable thickness of the initial dust layer on the surface of the northern hemisphere (cases dust 0/1/2 for all the northern locations and dust 0/1/2/3a/3b/3c for location 1), and maintained the equator and the southern hemisphere without any dust mantle (ice 0/1/2).

Table 2

Set of initial parameters of the computations.

thumbnail Fig. 1

Two-dimensional latitude-longitude map of 67P with the regions (El-Maarry et al. 2016) and the position of the 24 points studied, as detailed in Table 3. The latitude zero is the equatorial plane of the comet, separating the northern and the southern hemispheres. The centre of the map, i.e., point of latitude 0 and longitude 0, is the extremity of the small lobe. The big lobe area is located at the borders of the maps.

3.2.2 Illumination and orbit

To keep the code stable, our computations model a progressive approach of the comet to the Sun by using several orbits. Since we wanted to assess theinfluence of simple compositional parameters, we decided to use a small set of orbits corresponding to an orbit in the Kuiper belt (a = 50 au, e = 0.5) and then a transfer orbit (a = 25 au, e = 0.5) before inserting the comet onto its current orbit. The illumination conditions are derived from SPICE (Spacecraft Planet Instrument C-matrix shape Events) kernels applied to the shape model of 67P provided by the OSIRIS team (SHAPV). The current orbit of 67P is takenfrom Maquet (2015) and has a period of 6.44 yr, a semi-major axis of 3.5 au, and an eccentricity of 0.64 (Kelley et al. 2008; Davidsson & Gutiérrez 2005). The closest distance to the Sun (at perihelion) is about 1.3 au and the aphelion is at 5.7 au. Comet 67P has a rotation period of about 12.4 h, which slightly changes with the nucleus’ activity, and an obliquity of 52° (Lamy et al. 2007; Mottola et al. 2014 and Keller et al. 2015b).

Table 3

Coordinates of the 24 chosen locations indicated in Fig. 1, related morphological regions as defined in El-Maarry et al. (2015), the corresponding hemisphere and the used cases.

3.3 Analysis of the surface and interior temperature

Figures 2 and 3 show the variations of temperature of the surface and the near-surface interior for six locationsover the surface: two in the north (points 19 and 23 in Table 3), two at the equator (points 3 and 7), and two in the south (points 17 and 18). The temperatures are given for the surface, as well as for a depth of 0.01, 0.05, 0.1, and 0.2 m. In order to observe the profiles at two different positions on the orbit, Fig. 2 presents the variation over two rotations at the beginning of September 2014, at a heliocentric distance of 3.4 au, and Fig. 3 follows the same variations at the beginning of April 2015, at a heliocentric distance of 1.9 au.

At 3.4 au (inbound), the solar illumination heated the northern hemisphere and the equator, thus we observe a large difference of temperature over the nucleus’ surface. The northern and equator locations experienced the highest temperature, with maxima of surface temperature in the north between 230 and 250 K.

For a latitude of −40 degrees and below, while the points are not illuminated and show no variation in temperature, the interior of the comet is warmer than the surface, due to the low thermal inertia of the comet and solar heating during the previous perihelion passage. The temperatures increase by a few degrees with a very strong temperature gradient.

A clear diurnal variation of the temperature can be observed in the northern locations, equator locations, and in a location at a latitude of − 30 degrees. Solar illumination heats up the surface during the day and we observe clear temperature increases, which then decrease during the night. In particular, at the surface and at 1 cm under, we observe the most pronounced thermal variations. We observe the propagation of the heat wave into deeper layers. The maximum temperature at 1 cm below the surface is slightly delayed compared to the temperature curve of the surface. This is expected because of the thermal inertia and the time delay associated with heat propagation. The small variations at 5 cm show a faint maximum almost in phase opposition with the surface increase of temperature, that is, during the night. Underneath, the temperature is practically not affected by the variation of the illumination conditions on the surface (as expected considering the diurnal skin depth of a few centimetres). The temperatures under 10 cm are almost constant in time, but are not the same for the different locations studied. In the analysis of the interior temperature, illumination is not the only factor. When not illuminated, the temperature at the surface and at 1 cm even drops below the temperature of the interior. The variations of the temperature at 1 cm under the surface are smaller than the variations of the temperature at the surface, nevertheless, they show an amplitude of about 50 K. On the surface of 67P, frost on the surface has been found by the VIRTIS instrument (De Sanctis et al. 2015).

However, all the locations present a different pattern, even the points that have the same latitude. They endure different conditions of illumination and shadowing due to the non-spherical shape of the comet. Location 19 (longitude 30°, latitude 60°) stands as an example of bimodal peaks. As the nucleus is not spherical and presents cavities and cliffs, the temperature sometimes decreases during the day time due to shadowing by the structures on the comet itself.

Comparing Figs. 2 and 3, we observe that the temperature at different depths evolves slowly as the nucleus approaches the Sun. This is linked to the orbital skin depth, which is larger than the diurnal one.

In the northern hemisphere (locations 23 and 19), the surface and the depth of 1 cm had the same temperature profile; their temperatureincreased up to 280 K during the day, and decreased down to 110 K during the night. At larger depths, the temperatureprogressively decreased. Clear bimodal peaks are seen on the two northern points. In the southern hemisphere, the heat wave reached greater depths, which are warmer for the two points. Activity started for the point at latitude −40° as solar illumination reached those areas for the first time during the orbit, slightly before the equinox.

thumbnail Fig. 2

Evolution of the temperature at different depths for six chosen locations, over two rotations at a heliocentric distance of 3.4 au. The first line shows locations in the northern hemisphere, the second line shows the locations of the equator, and the third line the locations in the southern hemisphere. Layer 0 represents the surface, layer 1 is at a depth of 1 cm, layer 5 is at 5 cm, layer 10 is at 10 cm, and layer 20 is at 20 cm.

3.4 Evolution of the stratification

At the beginning of the computation, the nucleus is homogeneous and all the water is in an amorphous state. The homogeneous nucleus evolves into a stratified structure as the heat wave propagates inwards and leads to the sublimation of the volatiles. We observe here the evolution of the stratification on three locations over the nucleus: one point in the northern hemisphere (longitude 0°, latitude 30°), one point at the equator (0°, 0°), and one point in the southern hemisphere (0°, − 60°).

Figure 4 shows the evolution of the position of the surface, the first depth with H2O ice, the first depth with CO2 ice, the first depth with CO, and the transition front between amorphous and crystalline ice from a distance of 3.6 to 1.7 au.

On the left, the point in the north is computed with the dust 1 parameters (i.e. no initial dust layer, CO2/H2O = 0.02, and CO/H2O = 0.002). We observe the apparition of a dust mantle around 3.4 au, which reaches a thickness of 6 cm around 1.7 au. The erosion of the surface is very small, with an upper value of 0.04 cm. The front of CO remains constant at a depth of 52 m under the surface. The front of CO2 sinks from adepth of 32 cm under the new surface at 3.6 au to a depth of 48 cm at 1.7 au, as well as the front with amorphous ice, which sinks from 67 cm under the new surface at 3.6 au to 1.23 m at 1.7 au.

The point at the equator shows a very similar evolution of the interior structure. The important difference is that the computation (ice 1, i.e. no initial dust layer, CO2/H2O = 0.02, and CO/H2O = 0.002) does not allow the apparition of a dust mantle. Thus, the water front is at the surface, which has been eroded by 6.5 cm at the end of the run, an order of magnitude larger than the point at the north.

The point in the southern hemisphere, computed with the ice 1 parameters, is shown on the right of Fig. 4. As for the other locations, the depth of CO is found deep inside the nucleus, while the other fronts are very close to the surface. The southern hemisphere stayed in the shadow during the first months of the computation. We observe the beginning of activity at this location around March 2015, at 2.2 au from the Sun. As a consequence, we observe the front of amorphous ice dropping quickly from a depth of 0.9 cm under the new surface to 16.6 cm at the end of the computation. The depth with CO2 starts to sink at 1.7 au; at this distance, we observe a small erosion of 0.36 cm.

The evolution of the stratification of the three locations shown in Fig. 4 is representative of the evolution of the stratification for all the studied locations. Locations 14, 3, and 2 represent, respectively, the northern hemisphere, with the apparition of a dust mantle and very little erosion, the equator with a more significant erosion of the surface, and the southern hemisphere where most of the locations remain in shadow until March 2015, and experience the sinking of the crystallisation front and CO2 front later.

thumbnail Fig. 3

Evolution of the temperature at different depths for six chosen locations, over two rotations at a heliocentric distance of 1.9 au. Layer 0 represents the surface, layer 1 is at a depth of 1 cm, layer 5 is at 5 cm, layer 10 is at 10 cm, and layer 20is at 20 cm.
thumbnail Fig. 4

Evolution of the stratification for three locations: in the northern hemisphere (left), at the equator (middle), and in the southern hemisphere (right). The second and the third lines are zooms of the first line.

3.5 Surface fluxes for different initial dust mantle thicknesses

As specified in Table 2, the effects of dust mantling, with a thickness of 0.07 m, and CO ratios and trapping were simulated at all the surface points. To investigate more precisely the influence of the dust mantle thickness on surface fluxes red of the northern hemisphere of the comet, we selected a representative location (location 1 (45°, 30°)) on a flat surface and calculated the average fluxes for H2O, CO2 and CO for the four different initial dust layers (see values in Table 2). This location was chosen in the northern hemisphere on a flat surface as representative of the general activity of that part of the comet, so that we can expect the results obtained on this region to be generally applicable to the northern hemisphere. As shown in Fig. 5, the impact is very clear on the H2O outgassing (upper panel). The increase of the mantle thickness reduces considerably the average H2O fluxes. The evolution of the fluxes with no initial dust mantle is not continuously increasing (shift in the blue curve), due to the discrete description of the material instead of a continuous description. At each step, the model adapts the structure of the grid to take into account the loss of mass, therefore some values may change abruptly. The initial dust mantle of 0.2 m reduces the fluxes of the case with no dust layer by one order of magnitude. The effect is very important at 3.4 au with a reduction of 5 × 102 between the two cases. The difference is smaller at 1.9 au. Comparing the mean fluxes over the computations, it appears that the mean outgassing decreases exponentially with the thickness of the layer, at least for the four points of reference. However, the increase in the dust mantle thickness does not induce large changes in the outgassing of CO2 and CO, which will have important consequences on the interpretation of the RTOF maps (see Sect. 4.3).

3.6 Surface fluxes with or without trapped CO

The presence of trapped species in the amorphous ice has an impact on the outgassing composition. If most volatile species’ fronts of sublimation drop in the interior and thus disappear in the near-surface layer (see CO front in Fig. 4) when the nucleus is heated, minor species may be trapped in the amorphous ice and be released at the sublimation temperature of the amorphous ice.

To investigate the influence of CO trapped in the amorphous ice, we compared two runs for point 1 (45°, 30°) with the same set of parameters, except for the absence (dust 3b) or the presence (dust 3c) of trapped COgas. The effect on the outgassing is shown in Fig. 6. We note that CO2 may be found as a trapped species (Kouchi & Yamamoto 1995) but this is not investigated in this study.

For location 1, we observe no large change in the fluxes of H2O and CO2, while the CO outgassing is five orders of magnitude larger when CO is trapped in the amorphous ice. The fluxes of CO are constant due to the depth of the crystallisation front, while the fluxes of H2O and CO2 show variations because they are located closer to the surface and are influenced by night-day and seasonal variations in illumination. We note that some variations are not physical, but due to the computed grid that evolves to adjust the mass loss.

The presence of trapped CO in the amorphous ice induces a flow of CO coming from a depth much closer to the surface around the water ice crystallisation front. Indeed, the sublimation front of CO is located deeper (52.45 m) than the crystallisation front (1.27 m), due to the high volatility of CO (Prialnik 2004).

thumbnail Fig. 5

Evolution of fluxes from point 1 (longitude 45°, latitude 30°), for H2O (top), CO2 (middle), and CO (bottom), for case dust 1 (initial dust layer of 0 m of thickness), case dust 2 (0.07 m), case dust 3b (0.1 m), and case dust 3a (0.2 m), from about 3.45 au to 2.3 au.

4 Study of the spatial distribution of fluxes

4.1 Visualisation and interpolation

We computed the thermophysical model surface results on the 24 locations defined in Table 3 and interpolated the results to the entire surface of the nucleus with an interpolating method using Green’s functions for splines in tension according toWessel & Becker (2008). The adding of tension to the spline reduces unwanted oscillations between the points that appear with the use of cubic splines. An example of interpolation is shown on Fig. 7. The upper map shows the temperatureof 24 locations on the surface of 67P shape (European Space Agency) for one example of illumination. The middle map is the result obtained after applying the interpolating method with a tension fixed at five to the temperature of the 24 locations. The two bottom pictures present the resulted interpolation on the 3D shape of the nucleus, for two orientations of the nucleus.

thumbnail Fig. 6

Evolution of fluxes from point 1 (45°, 30°), for H2O (top), CO2 (middle), and CO (bottom), in case dust 3b (no trapped CO) and case dust 3c (presence of trapped CO), from about 3.45 to 2.3 au. A 10 cm dust layer is imposed in the northern hemisphere for both cases.
thumbnail Fig. 7

Example maps of interpolated surface temperature [K] using the Green’s functions for splines in tension. We show a map of the 24 studied locations (top), a 2D visualisation (middle), where the white facets are due to the non-convex shape of the nucleus, and a visualisation in which the locations are projected onto a 3D shape (bottom figures) after interpolation over all the surface regions. Topography lines are overplotted on the maps in black based on the 3D shape model provided by European Space Agency/Rosetta/MPS for the OSIRIS team. In this example, colours represent the surface temperature of the nucleus at around 3.3 au.

4.2 Model fluxes map

We calculated the average H2O, CO2, and CO surface fluxes for each studied location over Period A and B. The fluxes computed by the model are normalised for a surface of 1 m2 and given in molecules per m2 per second. We projected averaged fluxes on 2D maps and on the 3D shape of 67P for the two periods. The computed H2O, CO2, and CO maps are visible, respectively in the first, second, and third line of Figs. 8 and 9. Figure 8 presents the case dust 0–ice 0, computed with the dust 0 set of parameters for the northern points and the ice 0 set of parameters for the equator and southern points, while Fig. 9 shows the case dust 2–ice 2 (using the dust 2 set of parameters for the northern points and the ice 2 set of parameters for the equator and southern points). For each volatile, we present the 3D and 2D visualisation for Period A (left panels) and Period B (right panels). At the sametime, we list in Table 4 the median value of each volatile flux, for each presented map.

In Fig. 8, maps reveal a very active northern hemisphere, as it is the illuminated hemisphere during that part of the orbit. For all volatiles we observe the consequences of the different illumination conditions between Period A and Period B, with higher outgassing fluxes from the southern hemisphere, especially from the big lobe during the second period. We observe the effect of self-shadowing due to the topography of the nucleus, in particular in regions Khonsu, Wosret, and Imhotep (also seen in the temperature map in Fig. 7). These regions are located at the extremity of the lobes, around latitude −30°, and are shadowed by the lobes themselves.

The H2O appears clearly dominant with an averaged flux of about 1018 m−2 s−1. The highest fluxes of CO2 and CO have similar origins on the surface, with larger values of CO, due to the presence of trapped CO in the amorphous ice and its release at the transition to crystalline ice.

In Fig. 9, the case dust 2–ice 2 shows that a 0.7 m thick layer of dust reduces the outgassing of H2O in the northern hemisphere and around the equator. The origin of the fluxes seems more widespread over the surface and the dichotomy between the hemispheres seen in the case dust 0–ice 0 is reduced as well. The effects of self-shadowing are noticed in regions Khonsu and Anhur. We observe less shadowing than during Period A, as the sub-solar latitude changed from 41° to 33° and the Sun illuminated slightly more the shadowed regions around latitude −30°.

The CO2 fluxes appear clearly dependent on the initial CO/H2O ratio. For the same initial CO2/H2O ratio as in the case dust 0–ice 0 (set at 0.02), but a CO/H2O ratio divided by four, the CO2 averaged fluxes obtained with dust 2–ice 2 gained 1.88 × 1017 m−2 s−1. For the same heating, we observe a higher sublimation of CO2 as the available energy is less used for the sublimation of CO.

The third line shows a very weak activity for CO, of about 1012 –1013 m−2 s−1. This is a direct consequence of the absence of trapped CO in the amorphous ice.

4.3 Global comparison with RTOF

4.3.1 Deriving surface production rates from RTOF

In order to convert the RTOF densities into fluxes, we need information concerning the ejection velocities. From the surface temperature T given by the model, the mean speed v of gas molecules of mass m may be described in a first approximation by the speed distribution of a Maxwell–Boltzmann distribution (Huebner et al. 2006): v=8kTπm.\begin{equation*}v=\sqrt{\frac{8 k T}{\pi m}} .\end{equation*}(4)

We analysed the surface temperature computed from the thermophysical model and compared the calculated velocities with the velocities measured by MIRO and COPS onboard Rosetta. As an example, from the temperature computed for August 2014, we observe that velocities in the northern hemisphere varied between 420 and 500 m s−1, and in the southern hemisphere, between 325 and 370 m s−1. From MIROdata, Gulkis et al. (2015) analysed the temperature of the subsurface at 1.6 mm and 0.5 mm and derived an ejection velocity between 550 and 750 m s−1 for water from the data collected from 7 until 9 August 2014. The corresponding analysis with Eq. (4) gives a velocity from 385 to 425 m s−1, which underestimates the MIRO measurements. This difference is due to the approximation used to derive the velocities and, in particular, the equation does not account for the gas expansion in the coma. The COPS measurements and derived gas velocities were studied in Tzou (2016): she compared the COPS derived velocities with a direct simulation Monte Carlo (DSMC) model from the 1 to 5 October 2014. The DSMC model takes into account the dynamics of particles but Bieler et al. (2015b) showed that the cometary volatiles are barely accelerated. They computed a velocity range from 280 to 750 m s−1 for COPS compared to 300–500 m s−1 for the DSMC. Our model yields a velocity range from 390 to 435 m s−1 (for the days from 1 to 5 October 2014) in good agreement with ranges from MIRO, COPS, and DSMC (cf. Fig. 4 in Hansen et al. 2016; we note the large range values of COPS, indicating large uncertainties in this estimation).

From estimated ejection velocities, we converted RTOF densities into approximated surface production rates (allowing for a comparison with the model derived production rates) following Eq. (5), using the simple assumption of a spherical body of a surface S (4.97 × 107 m2 for 67P) and an isotropic expansion in r2 (where r is the radial distance from the nucleus centre), for a gas of density n (reported to the surface, i.e. at r = 1.99 km) and an outgassing velocity v: Q=v.n.S.\begin{equation*}Q=v.n.S .\end{equation*}(5)

One limit of this approach is the spherical assumption that allows us only to observe the large scale pattern (e.g. from which hemisphere the flux originates), as well as compare the order of magnitude of the surface production rates estimated by RTOF measurements and by the model. Moreover, the nucleus model being only 1D, considering a spherical shape should have a reduced impact. A full description of non-spherical shape effects would require 3D nucleus modelling combined with a 3D coma model, such as a DSMC model including collisions.

thumbnail Fig. 8

Representations in 3D and 2D of the H2O (first line), CO2 (second line), and CO (third line) averaged surface fluxes (in molec m−2 s−1) for Period A (six left panels) and Period B (six right panels) periods. The dust 0 set of parameters is used to compute the northern points and the ice 0 set of parameters is used for the equator and southern points.
thumbnail Fig. 9

Representations in 3D and 2D of the H2O (first line), CO2 (second line), and CO (third line) averaged surface fluxes (in molec m−2 s−1) for Period A (six left panels) and Period B (six right panels) periods. The dust 2 set of parameters is used to compute the northern points and the ice 2 set of parameters is used for the equator and southern points.
Table 4

Median values of H2O, CO2, and CO surface fluxes for periods A and B, computed with the cases dust 0–ice 0 and dust 2–ice 2.

4.3.2 Comparison between RTOF measurements and modelled fluxes

As an indication of the agreement between the computed and the RTOF derived production rates, we report the median values for each volatile measured by RTOF and computed by the model in Table 5. For RTOF production rate, we calculated the median surface production rates of RTOF derived from the densities of H2O, CO2, and CO measured during Periods A and B. For rate computed by the model, we calculated the median of the production rate averaged over the whole interpolated surface during the same periods.

For Period A, an initial dust layer of 0.07 m (dust 2–ice 2) reduces the mean water fluxes and provides the best prediction of the H2O fluxes measured, while no initial dust layer (dust 0–ice 0 and dust 1–ice 1) overestimates the observations by a factor of ten. The median fluxes of CO2 are in best agreement with the observations when the initial ratios CO2/H2O and CO/H2O are fixed at 0.02 (dust 0–ice 0). It appears that the trapped CO (dust 0–ice 0) is probably needed in the model to approach the RTOF observations for CO. It enables us to obtain a comparable outgassing of CO2 and CO.

Similarly, for Period B, the initial dust layer of 0.07 m (dust 2–ice 2) gives the best estimations of H2O observations, the initial ratios CO2/H2O and CO/H2O fixed at 0.02 m (dust 0–ice 0) are the closest to the observations of CO2, and the observations of CO are approached only with the presence of trapped CO. Nevertheless, the model reproduces closely the surface production rate for Period A, which is not the case for Period B. The evolution of the species density from Period A and Period B is apparently not explained if we only take into account different illumination conditions.

Figure 10 shows the flux maps derived from the RTOF densities for the Periods A (left) and B (right), with modelled ejection velocities of 420 m s−1 for H2O, 270 m s−1 for CO2, and 330 m s−1 for CO. Maps are produced assuming the origin of the fluxes to be given by the nadir off-pointing latitude and longitude at the moment of the detection. Also, we scaled the data to the surface of the nucleus (r = 1.98798 km), assuming a 1/r2 dependence. The resolution of the maps is 5° × 5°. Those measurements from RTOF present a behaviour that is not explained by the illumination changes: the decrease of CO2 and CO in the northern hemisphere from Period A to Period B.

The estimation of the velocities and the approximation on the radial gas expansion do not allow for a very precise comparison with the modelled maps. Nevertheless, such a comparison brings first order constraints on the amount of ejected ices and the depth of their sublimation fronts, the dust layer thickness, and the trapped CO required to explain globally the temporal and spatial variations measured by RTOF.

The first computation (dust 0–ice 0) had no initial dust layer but a thin dust layer naturally appeared in the end, with initial CO/H2O and CO2/H2O ratios of 0.02. For Period A, the comparison with the RTOF flux maps shows that this model overestimates the water outgassing in the northern hemisphere and underestimates the outgassing in the southern hemisphere. This difference may be due to the approximation of the nadir origin of the volatile species measured by RTOF: during most of the mission, the whole comet was in the field of view of RTOF. We note that RTOF data at locations situated at latitudes larger than 60 degrees latitude are not available. As seen by RTOF, the fluxes of CO are larger than those of CO2 and they originate from the same locations. They are both underestimated by the model globally over the nucleus. For Period B, H2O fluxes are quantitatively overestimated by the computation, with, however, a similarly increasing outgassing. At the equator, the model predicts H2O outgassing on the head and the big lobe, where RTOF recorded the highest densities. CO2 and CO maps from the model are completely opposite to RTOF maps. The modelled fluxes are largely overestimated in the northern hemisphere and underestimated in the southern hemisphere.

The dust 2–ice 2 computation has an initial CO/H2O ratio of 0.005 and CO2/H2O of 0.02, with an initial dust layer of 0.07 m. The addition of an initial dust mantle strongly reduces the H2O outgassing and the water flux provides a good estimation, but the more active areas are shifted compared to the ones seen by RTOF. This can be explained partly by the nadir approximation but also by the discrete number of studied locations and the surface interpolation. In this run, the detection of CO in the southern hemisphere was underestimated due to the low CO ratio. The decrease of the initial CO/H2O ratio to 0.005 strongly decreases the CO outgassing, with an averaged production rate of about 1020 s−1.

Moreover, the dichotomy seen on RTOF CO2 and CO maps is not observed in the modelled maps (using the dust 0–ice 0 or the dust 2–ice 2 parameters). In the RTOF maps, the spatial distribution evolves from Period A to Period B but does not follow the illumination variations. Except for the H2O maps, which show strong outgassing coming from the illuminated northern hemisphere, the CO2 and CO maps appear depleted in the northern hemisphere during Period B. The southern hemisphere is even more active in terms of the two volatiles, while being the less illuminated part at the time. With the set of input parameters used in our computation, this dichotomy is not observed. The evolution of the modelled maps from Period A to Period B is influenced mainly by the changes of sub-solar latitude (on average 41° for Period A and 33° for Period B).

We observe the influence of different initial conditions of the thermophysical nucleus model on the surface fluxes of the main volatiles of H2O, CO2, and CO, and on the change of spatial distribution over two periods of time. As the illumination condition slightly changed from Period A (1 September to 30 October 2014) to Period B (1 November 2014 to 15 February 2015), the sub-solar point moved from around 41° latitude during Period A to around 33° during Period B. This change in the illumination conditions had a direct incidence on the surface temperature of the nucleus. Local changes in temperature may modify the outgassing for a species, especially if the surface or internal temperature is close to the sublimation temperature. A map of the difference in surface temperature computed with the model between Period A and Period B is shown on Fig. 11. This representation indicates the areas that experienced a change of temperature, positive or negative. The southern hemisphere and the big lobe gained a few degrees, with an increase of more than 20 K around the Khonsu region, while the northern hemisphere remained constant except at high northern latitudes on the head, which lost a few degrees. From this map, we do not observe the hemispherical dichotomy measured by RTOF. The changes in illumination (and thus in surface temperature) between the two periods are not sufficient to explain the dichotomy. Direct illumination is obviously not the only factor that influenced the sublimation. Self-shadowing, self-heating, and local heterogeneities also influence locally the activity (as discussed earlier), but these local processes can hardly explain the global north-south dichotomy recorded by RTOF.

Table 5

Surface production rates [s−1] from RTOF (see Sect. 4.3.1) from the model set-ups.

thumbnail Fig. 10

Reflectron-type time of flight mass spectrometer surface mean fluxes (for a sphere of r = 1.98798 km) for Period A (left) and Period B (right), for H2O (top), CO2 (middle), and CO (bottom).
thumbnail Fig. 11

Map of the difference of average surface temperature [K] between Period A and Period B.

4.4 Comparison with a coma model

The above comparison between RTOF gas fluxes and the nucleus-model-derived surface fluxes assumes a radial expansion and a constant gasvelocity of the gas up to the coordinates of the RTOF instrument. However, this approach – which has the great advantage of linking the nucleus processes to the coma measurements – suffers from several limitations: the nucleus model is 1D and does not include full 3D effects at the complex surface of the comet, and the complex physics of gas expansion (with an angular distribution, and the mixing of several gas sources at the location of the spacecraft) are not included. This is the reason why we compare here the previous results with another complementary approach: we use the inverse coma model developed by Kramer et al. (2017) and Läuter et al. (2019), which includes the 3D coma processes, to derive surface emission rates from the inversion of the RTOF coma measurements, without making assumptions on the link between these emission rates and the nucleus physics. In the inverse coma model by Läuter et al. (2019), the volatiles in the coma are described by a collisionless gas model, while the complex shape of the nucleus is approximated by 3996 equally sized triangular elements. Each of the elements holds one gas emitter and the gas field at a given spacecraft position around the nucleus arises from the superposition of all surface emitters. The emission rates of the surface elements are the free model parameters and in the following are adjusted to best reproduce the measured coma densities. Multiple surface locations contribute to each spacecraft measurement. Therefore the inverse model leads to different surface-emission maps compared to the nadir flux attribution of Sect. 4.3.2 where each measurement is mapped to a single point on the surface. To retrieve the surface emission, we analysed RTOF measurements (17 September 2014 to 13 October 2014) from Period A and COPS/DFMS measurements (11 January 2015 to 5 February 2015) from Period B for the gas species H2O, CO2, and CO assuming the same molecular velocities as in Sect. 4.3.2. Once the surface emitters are determined by an optimization procedure from the data set, the gas model allows us to evaluate the gas densities in the coma on a spherical surface with fixed radius (here: 30 km) around the nucleus.

With RTOF measurements from Period A taken as input, the fixed-distance density of the inverse model (Fig. 12, left column)resembles closely the nadir projection (Fig. 10, left column) based on Eq. (5) for all species, H2O, CO2, and CO.

Two localised maxima around (−120°, 60°) and (60°, 60°) are discernible. In addition, the surface-integrated production rates for both methods match (see Table 5). The inverse model for the RTOF data in Period A (Fig. 12, right column) attributes the two maxima to a w-shaped band ranging between 0° and 60° northern latitude. The w-shaped band agrees with the diurnally averaged illumination during Period A (sub-solar latitude 41°) (see Hoang et al. 2017, Fig. 8), and supports an illumination-driven thermophysical model set-up as done in the set-up dust 0–ice 0. Due to the close agreement of the inverse model and the nadir projection with respect to the surface-integrated production, similar to Sect. 4.3.2, the H2O production is best predicted by the model set-up dust 2–ice 2. Despite the overestimation of the global production rate for H2O with the set-up dust 0–ice 0, this model agrees best with the inverse model in respect to the active surface areas.

The inversion of the COPS/DFMS measurements to a spherical surface by the end of Period B is shown in Fig. 13 (left column) and again shows good agreement to the nadir flux in Fig. 10 (right column). As for Period A, the H2O densities are attributed to the w-shaped band of enhanced surface flux (Fig. 13, upper right panel), albeit with overall increased emission levels closer to perihelion. In contrast to the H2O emission, CO2 and CO emissions are shifted to the southern hemisphere (see also Läuter et al. 2019, Fig. 7) and in Period B the CO2 flux exceeds the CO contribution. The deviation of the CO2 and CO emissions from the sub-solar latitude shows that a purely illumination driven ice model (set-up dust 0–ice 0) is not adequate to describe the volatile emission of 67P. Although not as close as for Period A, the H2O production for the model set-up dust 2–ice 2 best matches (factor 2.0) the inverse model (see Table 5).

The comparison between the nucleus model and inverse coma model surface production rates thus shows an overall good agreement in terms of qualitative conclusions and orders of magnitude values. The next step could be the development of a model that would combine the nucleus-model-derived surface production rates with a coma model that would use these rates as initial conditions for the coma expansion.

thumbnail Fig. 12

Inverse model applied to the RTOF measurements between 17 September 2014 and 13 October 2014 in Period A. Left column: density distribution on a 30 km sphere. Right column: surface flux. Rows: H2O, CO2, and CO; white areas indicate missing coverage.

5 Discussion

We provided here a first order comparison between a thermophysical model of the nucleus of 67P and RTOF observations. The model allows us to reproduce closely the estimated surface production rate based on the RTOF number density, in particular for Period A, which is from the 1 September to the 30 October 2014.

The model that best compares to the observations appears to have an initial ratio of CO/H2O smaller than that of CO2/H2O with the presence of trapped CO in the amorphous ice (initial value of about 0.02 for CO2/H2O and 0.005 for CO/H2O), as well as the presence of an initial dust layer of about 0.05–0.1 m (0.07 m in our case) in the northern hemisphere. However, we assumed an uniform composition to start the computation while the comet already passed the inner solar system several times. With this assumption, we overestimated the sublimation of CO and CO2 in the near surface layers, which probably led to an underestimation of the initial ratios.

The observed ROSINA gas fluxes could probably be reproduced by several sets of thermophysical parameters. However, a number of them were constrained more or less tightly by various Rosetta measurements (see Table 1). Our analysis shows the impact of several amongst the most important parameters that were poorly constrained by Rosetta, including the dust thickness, CO trapping mechanism, and minor species ratios. Moreover, we get an overall good agreement when comparing ROSINA measurements with two different models (the nucleus model and the inverse coma model), thus making us more confident in our conclusions.

However, two observations are not reproduced by the studied computations of the model: (1) the north-south dichotomy and (2) the decrease of CO and CO2 outgassing in the northern hemisphere between Periods A and B. The first observation will either require the use of different intrinsic ratios for the two hemispheres, which will be difficult to explain in terms of comet formation (all the more colliding bodies assumed to form the comet correspond to the lobes and not the hemispheres), or the forcing of different depths of sublimation fronts for CO2 and CO. Maquet (2015) indicated that the northern hemisphere might have been depleted during previous orbits, as the comet has witnessed eight weak and long summers in the north, and eight strong and short summers in the south in its current orbit since the 1959 encounter with Jupiter.

The second observation requires either an amount of CO and CO2 to be included in the dust back fall (but the estimated size of the dust particles is not large enough to explain this inclusion (Keller et al. 2017)), or the presence of a layer of pure water ice under the layer of dust that contains water molecules (Hoang et al. 2019). The near-surface interior in the northern hemisphere would be stratified as follows: a dust layer containing water ice, then a pure water ice layer, then CO2 and CO ices. In this configuration, heating would first sublimate water ice contained in the dust layer, but also CO2 and CO ices; as the temperature decreases rapidly with depth, it is not sufficient to induce the sublimation of the pure water ice layer. Closer to the Sun, once the temperature is sufficient to sublimate the water ice, little energy is left for sublimation of CO2 and CO andwe observe a decrease in the abundance of those volatiles in the outgassing.

However, this work assumed that the RTOF measurements originated from the nadir direction and that there are no collisions considered for particles. Nevertheless, the acceleration is small according to Fougere et al. (2016) and the acceleration in the direction of the spacecraft should not play a role as long as the flux is conserved. The nadir assumption may not reproduce local changes but global changes such as the north-south dichotomy should be revealed. Moreover, the comparison in the previous section with the inverse coma model developed by Kramer et al. (2017) and Läuter et al. (2019), which includes the complex gas transport (without collisions) in the coma, leads to surface production rates that are in overall good agreement with the surface production rates given by the thermophysical model. Additionally, the thermophysical modelling of the nucleus’ surface is complicated as the two hemispheres contain features, such as ice patches on the dust covered areas (Fornasier et al. 2016; Filacchione et al. 2016). Furthermore, the surface of the nucleus is not smooth but is covered by boulders of various sizes (Lee et al. 2017). Those details are not taken into account in this work, as we use the same set of parameters for each hemisphere. Moreover, fast activity events, such as outbursts, are not taken into account and can lead to non-continuous fluxes of volatiles and dust.

Beyond the comparison of the gas fluxes, it is also important to compare the model predictions and assumptions with other Rosetta observations. In particular, the model predicted surface temperatures of 140–245 K, comparable to the diurnal temperature of 180–250 K seen by VIRTIS at 3 au (Capaccioni et al. 2015). The computed surface temperaturegives ejection velocities in agreement with COPS measurements and DSMC simulations. The MUlti Purpose Sensors for Surface and Subsurface Science instrument (MUPUS) measured a dust layer of 10–20 cm on the Philae’s landing site Abydos, in the Wosret region, which is close to the values from 0.07 to 0.2 m used in our simulations. Such an initial dust layer simulates effectively the terrain covered by dust reported by El-Maarry et al. (2015) in the northern hemisphere. In addition, we obtained a water production rate of about 2 × 1025 s−1 at 3.45 AU and 2 × 1026 s−1 at 2.3 AU, in good agreement with the production rate presented in Hansen et al. (2016), ~4 × 1025 s−1 at 3.45 AU and ~3.5 × 1026 s−1 at 2.3 AU, who provided a multi-instruments study combining the results of ROSINA/DFMS, VIRTIS, Rosetta Plasma Consortium-Ion Composition Analyser (RPC-ICA), and MIRO.

thumbnail Fig. 13

Inverse model applied to the COPS/DFMS measurements between 11 January 2015 and 5 February 2015 in Period B. Left column: density distribution on a 30 km sphere. Right column: surface flux. Rows: H2O, CO2, and CO; white areas indicate missing coverage.

6 Conclusions

In this paper, we compared the Rosetta ROSINA/RTOF coma measurements of several species (H2O, CO2, and CO) with the surface production rates derived from a thermophysical model of comet 67P. This comparison, though limited by the 1D nature of the model and the radial gas expansion assumed, allows us to link the complex subsurface and internal processes to the main coma heterogeneities observed by ROSINA/RTOF. We thus provide constraints on the structure and the composition of the layers involved in producing the observed activity.

We focused in particular on the influence of several parameters that were poorly constrained by the Rosetta instruments. We show that the surface outgassing is strongly influenced by the initial conditions such as the initial dust layer, the trapping of CO, and the initial ratios of minor species. The comparison between model and RTOF maps reveals the following conclusions:

  • A smaller CO/H2O initial ratio increases considerably the CO2/H2O fluxes. More energy is available for CO2 sublimation, as less CO ice is present.

  • Making CO and CO2 fluxes similar to RTOF observations appears impossible without trapped CO in the amorphous ice. We note that Kouchi & Yamamoto (1995) mentioned the possibility of having CO trapped in CO2 ice, and experiments reveal CH4 and CO co-desorbing with CO2 at the sublimation temperature of CO2.

  • H2O and CO2 fluxes are slightly affected by the presence of trapped CO.

  • H2O fluxes are slightly affected by the initial ratios of minor species.

  • The dustmantle only slightly influences the CO2 and CO fluxes but rapidly reduces the fluxes of H2O.

  • An initial dust layer mantle of 0.05–0.1 m underestimates slightly the water fluxes but allows a better match to the relative amount of CO2 versus H2O measured by RTOF.

  • The decrease of the RTOF fluxes from Period A to Period B is not reproduced by the computations. This phenomenon is not well understood as the nucleus is approaching the Sun and the rest of the data set globally increases when the nucleus is closer to the Sun.

  • The north-south dichotomy that is clearly visible in the CO and CO2 RTOF maps cannot be reproduced by any model scenario above. In particular, a dust layer of 5–10 cm in the north does not sufficiently reduce the CO and CO2 fluxes; a much thicker layer would be needed, but this would almost extinguish the H2O fluxes in the north.

One major observation in the RTOF measurements (Hoang et al. 2017) was the clear dichotomy between thenorthern and the southern hemispheres. From our model results, further work is required to investigate this dichotomy by imposing different sublimation front depths, or setting very different CO2/H2O and CO/H2O ratios.

A comparison between the inverse coma model providing surface emission rates and the RTOF coma measurements shows a good agreement overall with the thermophysically derived emission rates. However, differences still exist that point to the need for a full coupling of a 3D nucleus model with a model of the coma, such as the DSMC model by Fougere et al. (2016), which would include both the nucleus processes and the gas transport and collisions in the coma. More work also needs to be done in terms of direct comparison with the other Rosetta instruments (like VIRTIS and MIRO) to refine the single-instrument-based interpretations.

Acknowledgement

The authors thank the following institutions and agencies, which supported this work. Work at IRAP was supported by the French space agency CNES. Work at the University of Bern was funded by the State of Bern and the Swiss National Science Foundation (200020-182418). The results from ROSINA would not be possible without the work of the many engineers, technicians, and scientists involved in the mission, the Rosetta spacecraft, and the ROSINA instrument team, whose contributions are gratefully acknowledged. All ROSINA flight data have been or will be released to the PSA archive of ESA and to the PDS archive of NASA.

References

  1. Agarwal, J., Della Corte, V., Feldman, P., et al. 2017, MNRAS, 469, s606 [NASA ADS] [CrossRef] [Google Scholar]
  2. Balsiger, H., Altwegg, K., Bochsler, P., et al. 2007, Space Sci. Rev., 128, 745 [NASA ADS] [CrossRef] [Google Scholar]
  3. Bar-Nun, A., Kleinfeld, I., & Kochavi, E. 1988, Phys. Rev. B, 38, 7749 [NASA ADS] [CrossRef] [Google Scholar]
  4. Bar-Nun, A., Notesco, G., & Owen, T. 2007, Icarus, 190, 655 [NASA ADS] [CrossRef] [Google Scholar]
  5. Biele, J., Ulamec, S., Maibaum, M., et al. 2015, Science, 349, aaa9816 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bieler, A., Altwegg, K., Balsiger, H., et al. 2015a, Nature, 526, 678 [Google Scholar]
  7. Bieler, A., Altwegg, K., Balsiger, H., et al. 2015b, A&A, 583, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Capaccioni, F., Coradini, A., Filacchione, G., et al. 2015, Science, 347, aaa0628 [Google Scholar]
  9. Capria, M., Capaccioni, F., Coradini, A., et al. 1996, Planet. Space Sci., 44, 987 [Google Scholar]
  10. Capria, M., Coradini, A., De Sanctis, M., & Blecka, M. 2001, Planet. Space Sci., 49, 907 [NASA ADS] [CrossRef] [Google Scholar]
  11. Capria, M. T., Capaccioni, F., Filacchione, G., et al. 2017, MNRAS, 469, S685 [Google Scholar]
  12. Combi, M., Shou, Y., Fougere, N., et al. 2019, Icarus, 335 113421 [Google Scholar]
  13. Davidsson, B. J., & Gutiérrez, P. J. 2005, Icarus, 176, 453 [NASA ADS] [CrossRef] [Google Scholar]
  14. De Sanctis, M., Capaccioni, F., Capria, M., et al. 1999, Planet. Space Sci., 47, 855 [Google Scholar]
  15. De Sanctis, M. C., Lasue, J., Capria, M. T., et al. 2010, Icarus, 207, 341 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. De Sanctis, M. C., Capaccioni, F., Ciarniello, M., et al. 2015, Nature, 525, 500 [NASA ADS] [CrossRef] [Google Scholar]
  17. Delsemme, A. H., & Swings, P. 1952, Ann. Astrophys., 15, 1 [NASA ADS] [Google Scholar]
  18. El-Maarry, M.R., Thomas, N., Giacomini, L., et al. 2015, A&A, 583, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. El-Maarry, M.R., Thomas, N., Gracia-Berná, A., et al. 2016, A&A, 593, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Espinasse, S., Klinger, J., Ritz, C., & Schmitt, B. 1991, Icarus, 92, 350 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Filacchione, G., Capaccioni, F., Ciarniello, M., et al. 2016, Icarus, 274, 334 [NASA ADS] [CrossRef] [Google Scholar]
  22. Fornasier, S., Mottola, S., Keller, H. U., et al. 2016, Science, 354, 1566 [NASA ADS] [CrossRef] [Google Scholar]
  23. Fougere, N., Altwegg, K., Berthelier, J.-J., et al. 2016, A&A, 588, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Fulle, M., Della Corte, V., Rotundi, A., et al. 2015, ApJ, 802, L12 [NASA ADS] [CrossRef] [Google Scholar]
  25. Fulle, M., Marzari, F., Della Corte, V., et al. 2016, ApJ, 821, 19 [NASA ADS] [CrossRef] [Google Scholar]
  26. Guilbert-Lepoutre, A., Lasue, J., Federico, C., et al. 2011, A&A, 529, A71 [Google Scholar]
  27. Gulkis, S., Allen, M., von Allmen, P., et al. 2015, Science, 347, aaa0709 [Google Scholar]
  28. Gundlach, B., & Blum, J. 2012, Icarus, 219, 618 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Hansen, K. C., Altwegg, K., Berthelier, J.-J., et al. 2016, MNRAS, 462, S491 [Google Scholar]
  30. Hässig, M., Altwegg, K., Balsiger, H., et al. 2015, Science, 347, aaa0276 [Google Scholar]
  31. Hoang, M., Altwegg, K., Balsiger, H., et al. 2017, A&A, 600, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Hoang, M., Garnier, P., Gourlaouen, H., et al. 2019, A&A, 630, A33 [Google Scholar]
  33. Hu, X., Shi, X., Sierks, H., et al. 2017a, MNRAS, 469, S295 [CrossRef] [Google Scholar]
  34. Hu, X., Shi, X., Sierks, H., et al. 2017b, A&A, 604, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Huebner, W. F., Benkhoff, J., Capria, M.-T., et al. 2006, Heat and Gas Diffusion in Comet Nuclei (Bern, Switzerland: International Space Science Institute) [Google Scholar]
  36. Jewitt, D. 2009, AJ, 137, 4296 [NASA ADS] [CrossRef] [Google Scholar]
  37. Keller, H. U., Mottola, S., Davidsson, B., et al. 2015a, A&A, 583, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Keller, H. U., Mottola, S., Skorov, Y., & Jorda, L. 2015b, A&A, 579, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Keller, H. U., Mottola, S., Hviid, S. F., et al. 2017, MNRAS, 469, S357 [CrossRef] [Google Scholar]
  40. Kelley, M. S., Reach, W. T., & Lien, D. J. 2008, Icarus, 193, 572 [NASA ADS] [CrossRef] [Google Scholar]
  41. Kofman, W., Herique, A., Barbin, Y., et al. 2015, Science, 349, aab0639 [Google Scholar]
  42. Kouchi, A., & Sirono, S.-i. 2001, Geophys. Res. Lett., 28, 827 [NASA ADS] [CrossRef] [Google Scholar]
  43. Kouchi, A., & Yamamoto, T. 1995, Prog. Cryst. Growth Character. Mater., 30, 83 [CrossRef] [Google Scholar]
  44. Kouchi, A., Greenberg, J., Yamamoto, T., Mukai, T., et al. 1992, ApJ, 388, L73 [Google Scholar]
  45. Kramer, T., Läuter, M., Rubin, M., & Altwegg, K. 2017, MNRAS, 469, S20 [CrossRef] [Google Scholar]
  46. Lamy, P. L., Toth, I., Davidsson, B. J., et al. 2007, Space Sci. Rev., 128, 23 [NASA ADS] [CrossRef] [Google Scholar]
  47. Lasue, J., De Sanctis, M., Coradini, A., et al. 2008, Planet. Space Sci., 56, 1977 [NASA ADS] [CrossRef] [Google Scholar]
  48. Läuter, M., Kramer, T., Rubin, M., & Altwegg, K. 2019, MNRAS, 483, 852 [Google Scholar]
  49. Lee, J.-C., Massironi, M., Ip, W.-H., et al. 2017, MNRAS, 462, S573 [CrossRef] [Google Scholar]
  50. Luspay-Kuti, A., Mousis, O., Hässig, M., et al. 2016, Sci. Adv., 2, e1501781 [NASA ADS] [CrossRef] [Google Scholar]
  51. Maquet, L. 2015, A&A, 579, A78 [Google Scholar]
  52. Marov, M. Y., Rusol, A., & Dorofeeva, V. 2019, Doklady Physics, 64, 34 [Google Scholar]
  53. Marschall, R., Rezac, L., Kappel, D., et al. 2019, Icarus, 328, 104 [NASA ADS] [CrossRef] [Google Scholar]
  54. Mottola, S., Lowry, S., Snodgrass, C., et al. 2014, A&A, 569, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Notesco, G., & Bar-Nun, A. 1996, Icarus, 122, 118 [NASA ADS] [CrossRef] [Google Scholar]
  56. Pätzold, M., Andert, T., Hahn, M., et al. 2016, Nature, 530, 63 [NASA ADS] [CrossRef] [Google Scholar]
  57. Prialnik, D. 2004, Comets II (Tucson: University of Arizona Press), 1, 359 [Google Scholar]
  58. Prialnik, D., & Bar-Nun, A. 1992, A&A, 258, L9 [NASA ADS] [Google Scholar]
  59. Prialnik, D., Benkhoff, J., & Podolak, M. 2004, Modeling the Structure and Activity of Comet Nuclei (Tucson: University of Arizona), 1, 359 [Google Scholar]
  60. Rotundi, A., Sierks, H., Della Corte, V., et al. 2015, Science, 347, aaa3905 [NASA ADS] [CrossRef] [Google Scholar]
  61. Rubin, M., Altwegg, K., Balsiger, H., et al. 2015, Science, 348, 232 [NASA ADS] [CrossRef] [Google Scholar]
  62. Rusol, A. V., & Dorofeeva, V. A. 2018, Open Astron., 27, 175 [Google Scholar]
  63. Schloerb, F. P., Keihm, S., von Allmen, P., et al. 2015, A&A, 583, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  64. Sierks, H., Barbieri, C., Lamy, P. L., et al. 2015, Science, 347, aaa1044 [Google Scholar]
  65. Skorov, Y. V., Van Lieshout, R., Blum, J., & Keller, H. U. 2011, Icarus, 212, 867 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  66. Spohn, T., Knollenberg, J., Ball, A., et al. 2015, Science, 349, aab0464 [NASA ADS] [CrossRef] [Google Scholar]
  67. Tzou, C.-Y. 2016, Ph.D. thesis, Universitaet Bern, Switzerland [Google Scholar]
  68. Wessel, P., & Becker, J. 2008, Geophys. J. Int., 174, 21 [Google Scholar]

All Tables

Table 1

Physical parameters constrained by the Rosetta’s instruments.

In the text
Table 2

Set of initial parameters of the computations.

In the text
Table 3

Coordinates of the 24 chosen locations indicated in Fig. 1, related morphological regions as defined in El-Maarry et al. (2015), the corresponding hemisphere and the used cases.

In the text
Table 4

Median values of H2O, CO2, and CO surface fluxes for periods A and B, computed with the cases dust 0–ice 0 and dust 2–ice 2.

In the text
Table 5

Surface production rates [s−1] from RTOF (see Sect. 4.3.1) from the model set-ups.

In the text

All Figures

thumbnail Fig. 1

Two-dimensional latitude-longitude map of 67P with the regions (El-Maarry et al. 2016) and the position of the 24 points studied, as detailed in Table 3. The latitude zero is the equatorial plane of the comet, separating the northern and the southern hemispheres. The centre of the map, i.e., point of latitude 0 and longitude 0, is the extremity of the small lobe. The big lobe area is located at the borders of the maps.
In the text
thumbnail Fig. 2

Evolution of the temperature at different depths for six chosen locations, over two rotations at a heliocentric distance of 3.4 au. The first line shows locations in the northern hemisphere, the second line shows the locations of the equator, and the third line the locations in the southern hemisphere. Layer 0 represents the surface, layer 1 is at a depth of 1 cm, layer 5 is at 5 cm, layer 10 is at 10 cm, and layer 20 is at 20 cm.
In the text
thumbnail Fig. 3

Evolution of the temperature at different depths for six chosen locations, over two rotations at a heliocentric distance of 1.9 au. Layer 0 represents the surface, layer 1 is at a depth of 1 cm, layer 5 is at 5 cm, layer 10 is at 10 cm, and layer 20is at 20 cm.
In the text
thumbnail Fig. 4

Evolution of the stratification for three locations: in the northern hemisphere (left), at the equator (middle), and in the southern hemisphere (right). The second and the third lines are zooms of the first line.
In the text
thumbnail Fig. 5

Evolution of fluxes from point 1 (longitude 45°, latitude 30°), for H2O (top), CO2 (middle), and CO (bottom), for case dust 1 (initial dust layer of 0 m of thickness), case dust 2 (0.07 m), case dust 3b (0.1 m), and case dust 3a (0.2 m), from about 3.45 au to 2.3 au.
In the text
thumbnail Fig. 6

Evolution of fluxes from point 1 (45°, 30°), for H2O (top), CO2 (middle), and CO (bottom), in case dust 3b (no trapped CO) and case dust 3c (presence of trapped CO), from about 3.45 to 2.3 au. A 10 cm dust layer is imposed in the northern hemisphere for both cases.
In the text
thumbnail Fig. 7

Example maps of interpolated surface temperature [K] using the Green’s functions for splines in tension. We show a map of the 24 studied locations (top), a 2D visualisation (middle), where the white facets are due to the non-convex shape of the nucleus, and a visualisation in which the locations are projected onto a 3D shape (bottom figures) after interpolation over all the surface regions. Topography lines are overplotted on the maps in black based on the 3D shape model provided by European Space Agency/Rosetta/MPS for the OSIRIS team. In this example, colours represent the surface temperature of the nucleus at around 3.3 au.
In the text
thumbnail Fig. 8

Representations in 3D and 2D of the H2O (first line), CO2 (second line), and CO (third line) averaged surface fluxes (in molec m−2 s−1) for Period A (six left panels) and Period B (six right panels) periods. The dust 0 set of parameters is used to compute the northern points and the ice 0 set of parameters is used for the equator and southern points.
In the text
thumbnail Fig. 9

Representations in 3D and 2D of the H2O (first line), CO2 (second line), and CO (third line) averaged surface fluxes (in molec m−2 s−1) for Period A (six left panels) and Period B (six right panels) periods. The dust 2 set of parameters is used to compute the northern points and the ice 2 set of parameters is used for the equator and southern points.
In the text
thumbnail Fig. 10

Reflectron-type time of flight mass spectrometer surface mean fluxes (for a sphere of r = 1.98798 km) for Period A (left) and Period B (right), for H2O (top), CO2 (middle), and CO (bottom).
In the text
thumbnail Fig. 11

Map of the difference of average surface temperature [K] between Period A and Period B.
In the text
thumbnail Fig. 12

Inverse model applied to the RTOF measurements between 17 September 2014 and 13 October 2014 in Period A. Left column: density distribution on a 30 km sphere. Right column: surface flux. Rows: H2O, CO2, and CO; white areas indicate missing coverage.
In the text
thumbnail Fig. 13

Inverse model applied to the COPS/DFMS measurements between 11 January 2015 and 5 February 2015 in Period B. Left column: density distribution on a 30 km sphere. Right column: surface flux. Rows: H2O, CO2, and CO; white areas indicate missing coverage.
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.