Possibility of Concentration of Non-volatile Species near the Surface of Comet 67P/Churyumov-Gerasimenko

TL;DR

This study models how non-volatile compounds may concentrate near the surface of Comet 67P/Churyumov-Gerasimenko over time, indicating potential delivery of prebiotic materials to early Earth.

Contribution

We developed the Simplified Cometary Concentration Model (SCCM) to simulate the evaporation and concentration of non-volatile species on comet surfaces, using 67P as a benchmark.

Findings

Non-volatile species concentrate within 0-100cm of the surface in 10 comet years.

The model explains observed production rates of water and amino compounds.

Concentrated non-volatile materials could be delivered to early Earth.

Abstract

The cometary materials are thought to be the reservoir of primitive materials in the Solar System. The recent detection of glycine and CH$_3$NH$_2$ by the ROSINA mass spectrometer in the coma of 67P/Churyumov-Gerasimenko suggests that amino acids and their precursors may have been formed in an early evolutionary phase of the Solar System. We investigate the evolution of cometary interior considering the evaporation process of water followed by the concentration of non-volatile species. We develop a Simplified Cometary Concentration Model (SCCM) to simulate the evaporation and concentration processes on the cometary surface.We use 67P/Churyumov-Gerasimenko as the benchmark of SCCM. We investigate the depth of the layer where non-volatile species concentrate after the numerous passages of perihelion after the formation of the Solar System. As a result, the SCCM explains the observed production rates of water and CH$_3$NH$_2$ at 100 comet years. SCCM results suggest that the non-volatile species would concentrate at depths between 0 and 100cm of comet surface within 10 comet years. Our results also suggest that the non-volatile species would concentrate several meters beneath the surface before it hit the early Earth. This specific mass of non-volatile species may provide unique chemical condition to the volcanic hot spring pools.

Tables

Table 1. Table 1 : The Parameters for SCCM Modeling

Model	Initial ratios of non-volatile species	Layer Size (cm)	Initial Porosity
Set 1	0.5	10	0
Set 2	0.5	10	0.2
Set 3	0.5	5	0
Set 4	0.5	20	0
Set 5	0.1	10	0
Set 6	0.9	10	0

Equations

$$T[n,t]=T[n=1,t]\times f[n,t],$$

$$t[n,j,t{]}_{\mathrm{evaporation}}={\nu }_0^{-1}\exp (E[j]/k_{\mathrm{B}}T[n,t]),\mathrm{where},$$

$${\nu }_0=(2n_{\mathrm{s}}E[j]/{\pi }^{2}m{)}^{1/2}$$

$$P[n,t]=\frac{V_{\mathrm{Max}}-V_{\mathrm{Real}}[n,t]}{V_{\mathrm{Max}}},$$

$$\frac{\Delta V_{\mathrm{evap}}[1,j,t]}{\Delta t}=k[1,j,t]\times V_{\mathrm{Real}}[1,j,t].$$

$$\frac{\Delta V_{\mathrm{evap}}[n+1,j,t]}{\Delta t}=P[n,t]\times k[n+1,j,t]\times V_{\mathrm{Real}}[n+1,j,t](n\geqq 1).$$

$$\Delta V_{\mathrm{direct}}[n+1,j,t]=(\prod _{k=1}^{k=n}P[k,t])\times \Delta V_{\mathrm{evap}}[n+1,j,t].$$

$$\begin{array}{r}\Delta V_{\mathrm{Real}}[n,j,t]=-\Delta t\times (\Delta V_{\mathrm{evap}}[n,j,t]/\Delta t\\ -\Delta M_{\mathrm{evap}}[n+1,j,t]/\Delta t+\Delta V_{\mathrm{direct}}[n,j,t]/\Delta t).\end{array}$$

$$V_{\mathrm{Real}}[n,j,t+\Delta t]=V_{\mathrm{Real}}[n,j,t]+\Delta V_{\mathrm{Real}}[n,j,t].$$

$$P[n,t+\Delta t]=\frac{V_{\mathrm{Max}}-{\Sigma }_j{V}_{\mathrm{Real}}[n,j,t+\Delta t]}{V_{\mathrm{Max}}}.$$

$$\begin{array}{r}r[n,j,t+\Delta t]=V_{\mathrm{Real}}[n,j,t+\Delta t]/V_{\mathrm{Real}}[n,t+\Delta t]\\ =V_{\mathrm{Real}}[n,j,t]+\Delta V_{\mathrm{Real}}[n,j,t]/(1-P[n,t+\Delta t])\times V_{\mathrm{Max}}.\end{array}$$

$$\begin{array}{r}{P}_{\mathrm{new}}[n,t]=\frac{V_{\mathrm{Max}}-(V_{\mathrm{Real}(\mathrm{before})}[n,t]+V_{\mathrm{Real}(\mathrm{before})}[n+1,t])}{V_{\mathrm{Max}}}\\ =(P_{\mathrm{before}}[n,t]+P_{\mathrm{before}}[n+1,t])-1.\end{array}$$

$$\begin{array}{r}{r}_{\mathrm{new}}[n,j,t]\\ =\frac{(1-P_{\mathrm{before}}[n,t])r[n,j,t]+(1-P_{\mathrm{before}}[n+1,t])r[n+1,j,t]}{1-P_{\mathrm{new}}[n,t]}.\end{array}$$

$$\begin{array}{r}\mathrm{Production}[j,t]=\rho \times S\times (\frac{\Delta V_{\mathrm{evap}}[n=1,j,t]}{\Delta t}\\ +{\Sigma }_{n\geqq 1}\frac{\Delta V_{\mathrm{direct}}[n\geqq 1,j,t]}{\Delta t})\times d\text{leavevmode}\mathrm{kg}\text{leavevmode}{\mathrm{s}}^{-1},\end{array}$$

$$\mathrm{Erosion}[j]={\Sigma }_t\frac{\mathrm{Production}[j,t]\Delta t}{\rho }\text{leavevmode}\mathrm{m}$$

Full text

*Aims. *We investigate the evolution of cometary interior considering the evaporation process of water, followed by the concentration of non-volatile species.

*Methods. *We develop a Simplified Cometary Concentration Model (SCCM) to simulate the evaporation and concentration processes on the cometary surface. We use 67P/Churyumov-Gerasimenko as the benchmark of SCCM. We investigate the depth of the layer where non-volatile species concentrate after the numerous passages of perihelion after the formation of the Solar System.

*Results. *We find that the SCCM explains the observed production rates of water and CH3NH2 at 100 comet years. SCCM results suggest that the non-volatile species would concentrate at depths between 0 and 100 cm of comet surface within 10 comet years. Our results also suggest that the non-volatile species would concentrate several meters beneath the surface before it hit the early Earth. This specific mass of non-volatile species may provide unique chemical condition to the volcanic hot spring pools.

1 Introduction

Since the pioneering suggestion by Oró (1961), many researchers discussed the possibility of external delivery of organic materials by comet or meteorites as the supplier of organic materials to the early Earth. Ehrenfreund et al. (2002) strengthened the importance of comet for the sources of organic materials compared to other possible sources on the Earth, such as UV photolysis or the electron discharging in the ancient atmosphere. If this is the case, an attempt to understand the origins of life must necessarily begin with detailed studies of the formation and evolution of complex organic molecules. These organic molecules are the products of complex chemistry that most likely starts in molecular clouds and continues within the protoplanetary disk, which is the birthplace of comets and asteroids. The detections of various kinds of organic materials in the Murchison meteorites, such as sugars and amino acids (Cooper et al., 2001; Engel & Nagy, 1982), support this theory.

If this is the case, the chemical evolution of material should start from the parent cloud of the Solar System. Recent observations improved our understanding of the chemical compositions in the protoplanetary disk and star-forming regions. The high sensitivity observation of ALMA achieved the detections of complex organic molecules CH3CN (Öberg et al., 2015; Loomis et al., 2018; Bergner et al., 2018); CH3OH (Walsh et al., 2016); and HCOOH (Favre et al., 2018). Besides, observations toward low-mass star-forming regions are much easier than protoplanetary disks due to high beam-averaged column densities of molecules. More complex molecules such as glycolaldehyde (HCOCH2OH) (Jørgensen et al., 2012), methyl isocyanate (CH3NCO)(Martin-Domenech et al., 2017), and CH3Cl (Fayolle et al., 2017) were reported towards low-mass star-forming regions. Glycine (NH2CH2COOH), the simplest and only non-chiral member out of 20 standard amino acids, has gathered the attention of astronomers; however, none of the observations were successful toward any star-forming regions (e.g., Ceccarelli et al., 2000). The astrochemical modeling of interstellar glycine predicted that CH3NH2 and CH2NH would be potential precursors of glycine (Garrod, 2013; Suzuki et al., 2018b). The detections of CH3NH2 and CH2NH were reported towards various high-mass star-forming regions (Suzuki et al., 2018a; Ohishi et al., 2019). For glycine itself, Kuan et al. (2003) claimed the first detections of glycine towards high-mass star-forming regions, but their detections were refuted by Snyder et al. (2005).

It is thought that the materials in star-forming regions and protoplanetary disks could be incorporated into the comet. Walsh et al. (2014) performed the chemical modeling study for the protoplanetary disk and showed that grain-surface fractional abundances of simple and complex organic molecules (relative to water ice) for the outer disk are consistent with abundances derived for comets. Biver et al. (2015) also claimed good agreement of the chemical compositions between the cometary coma and warm cores.

The amino acids are also a major target for the in-situ or sample return mission missions of comets. Elsila et al. (2009) claimed the detection of glycine in the Stardust sample with its carbon isotopic ratio $\delta^{13}C=+29\pm 6$ ‰, suggesting the extraterrestrial origin. The latest paradigm shift of our knowledge was brought about by the detailed analysis of comet 67P/Churyumov-Gerasimenko (hereafter 67P) by Rosetta mission, enabling us the direct measurement of its chemical composition. This mission confirmed all the known chemical species at that time (Le Roy et al., 2015), and new detections of various complex species as well (e.g., CH3NCO, CH3COCH3, C2H5CHO, and CH3CONH2) through the direct analysis of COmetary Sampling And Composition experiment COSAC (Goesmann et al., 2015), though the detection of CH3NCO was disputed by the proceeding study (Altwegg et al., 2017). Recent measurement suggested the existence of ammonium salt in this comet, which could be the substantial reservoir of nitrogen (Altwegg et al., 2020; Poch, 2020). In addition, the volatile glycine and its precursor, CH3NH2, were detected in the coma of comet 67P by the ROSINA (Rosetta Orbiter Spectrometer for Ion and Neutral Analysis) mass spectrometer (Altwegg et al., 2016). Hadraoui et al. (2019) argued that glycine would originate from sublimating water ice from dust particles that are ejected from the nucleus.

Though this scenario is fascinating, there would be room for discussing the cometary surface’s chemistry and physics. As Capaccioni et al. (2015) reported the very low albedo on the surface of 67P, water rapidly sublimates from the surface, and refractory materials are left behind. Capaccioni et al. (2015) also suggested that refractory organic materials cover the comet surface. The Deep Impact mission (Mumma et al., 2005) reported the significant increase of a volatile species, C2H6, after impact to comet 9P/Tempel 1. The sudden increase of C2H6 implies that the volatile species were lost from the comet surface due to thermal processing, while it is still frozen below the surface of the comet. This observation also emphasizes the importance of thermal processing to discuss the chemical compositions close to the comet surface. During the thermal processing by the Sun, volatile species would be easily lost, while non-volatile species would concentrate at the surface, making the non-homogeneous distributions of molecules.

Comets would have hit the early Earth providing the materials. Pierazzo & Chyba (1999) performed hydrocode simulations of the comet impact and studied the thermal decomposition rates of amino acids based on the time evolution of the temperature after the impact. As a result, they showed that the opposite side to the collision point is a good place where amino acids could survive the high-temperature condition during the impact. If this is the case, it will provide the high-density organic mass without diluting to the volcanic hot spring pools, which is suggested as the candidate place for the origin of life (Damer & Deamer, 2020). As the concentration of non-volatile species gives the unique chemical composition to the volcanic hot spring pools, it is interesting to investigate the chemical distribution in the cometary interiors. Besides, a previous modeling study by Keller et al. (2015) overestimated the water production rate of 67P. The concentration of the non-volatile species may reduce the evaporation of water and other species from the comet’s interior, which may improve the agreement between the modeled and the observed water sublimation rate.

In this paper, we develop a new physical model for the evolution of the cometary interior, including the evaporation processes of water, CH3NH2, and glycine. In Section 2, we describe the details of the Simplified Cometary Concentration Model (SCCM). Then, in Section 3, we benchmark SCCM with the production rates of water and CH3NH2 and the layer’s thickness composed of the non-volatile species. We summarize our work in Section 4.

2 Simplified Cometary Concentration Model (SCCM)

2.1 Calculation Method of SCCM

In the SCCM, we assume a comet’s layer structure. Our layer structure is depicted in Figure 1. The thickness of the layer is the parameter that determines the resolution of the vertical structure. While the smaller size of the layer improves the accuracy of the model, it requires higher computational time.

In the SCCM model, we calculate the temperature of the nth layer at the time $t$, $T[n,t]$, with the following procedure. $t$ is the “comet years”, where one comet year is 6.45 years according to JPL Small-Body Database Browser111https://ssd.jpl.nasa.gov/sbdb.cgi, corresponding to one orbital period of the comet 67P. First, we simulate the orbit of 67P with its orbital parameters to calculate the surface temperature at time $t$, $T[n=1,t]$. Then, we assume the radiative equilibrium between the irradiation from the Sun and the blackbody radiation of the comet. The obtained temperature is in good agreement with that of Guilbert-Lepoutre et al. (2016) for the surface. For simplicity, we calculate the surface temperature only from zero to 0.5 comet year and assume the periodic time variation after that. We note that this symmetric temperature variation is not strictly correct, as the actual peak activity of 67P is observed two weeks after the perihelion (Snodgrass et al., 2016). We use the temperatures of 260, 240, 190, and 120 K for the depths of 0, 10, 50, and 200 cm, respectively, at the perihelion ($t$=0), following the calculation of Guilbert-Lepoutre et al. (2016). The peak temperature of other depth, $T[n,t=0]$, were linearly interpolated. Following Guilbert-Lepoutre et al. (2016), we assume that the temperature is 120 K for all layers at 0.5 comet years, when the comet is the most distant from the Sun. Then, for the depth of less than 200 cm, we calculate the temperature of the nth layer at the time $t$ as

[TABLE]

where $f[n,t]=2\times(1-q[n])\times t+q[n]$ and $q[n]=T[n,t=0]/T[n=1,t=0]$. We neglect the time variation of the temperature for the layers deeper than 200 cm. We fix the temperature of 120 and 30 K at a depth of 200 and 1000 cm, respectively, and use linear interpolation at other depths. There is no temperature gradient inside the single layer.

We show our temperature distribution in Figures 2. The right figure shows the time dependency of the temperature distribution inside the comets, while the left figure represents the vertical temperature distribution at the different depth from the surface at the different distance from the Sun.

With the temperature, we calculate the evaporation rate from the comet surface with the equation (2) of Hasegawa et al. (1992), which is frequently used in the astrochemical models (Ruaud et al., 2016):

[TABLE]

[TABLE]

The inverse of the evaporation timescale is the evaporation rate $k[n,j,t]$ (i.e., $k[n,j,t]=t[n,j,t]_{\rm evaporation}^{-1}$). $E[j]$ is the binding energy, which is determined by the strength of interaction of species $j$ with water ice, $n_{\rm s}$ is the surface density of the site, which is 1.5$\times 10^{15}$ cm*-2* (Tielens & Allamandola, 1987), $m$ is the mass of the molecules, and $k_{\rm B}$ is the Boltzmann constant. The evaporation rate of species $j$ from $n^{th}$ layer at time $t$, “$k[n,j,t]$”, is proportional to exp($-E[j]/k_{\rm B}T[n,t]$). Hence, the binding energy and the temperature are essential to determine $k[n,j,t]$. We assume that all molecules can evaporate from the surface without being trapped inside. While the molecules in the deeper layer have to move to upper layers one by one through this usual evaporation process (equation (5)), they can escape directly from the comet from any layer via a direct evaporation process as described below (equation (6)). The non-volatile species should be mainly minerals and insoluble organic matter (Fray et al., 2017). The binding energies for H2O, CH3NH2 are 5700 and 6500 K, respectively, based on the theoretical prediction by Wakelam et al. (2017). We assume that thermal evaporation is negligible for non-volatile species. For glycine, the binding energy would be 13000 K (Tzvetkov et al., 2004; Suzuki et al., 2018b), and the thermal evaporation is not dominant in the orbit of 67P. Therefore, we treat glycine as a part of non-volatile species.

In SCCM modeling, we define the porosity $P[n,t]$ as

[TABLE]

where $V_{\rm Real}[n,t]=\Sigma_{j}V_{\rm Real}[n,j,t]$ represents the sum of the volumes occupied by materials of all species $j$ on $n^{th}$ layer and at $t$, while $V_{\rm Max}$ denotes maximum volume in a certain layer. From this definition, it is apparent that the term $V_{\rm Real}[n,t]=(1-P[n,t])\times V_{\rm Max}$ denotes the real volume occupied by any material in a layer. Here we introduce the ratio of certain species “$j$” in nth layer at $t$, $r[n,j,t]$ as $r[n,j,t]=V_{Real}[n,j,t]/V_{Real}[n,t]$. The sum of the ratios of all molecules, $\sum_{j}r[n,j,t]$, should always be unity.

With the term of $V_{\rm Real}[n,j,t]=V_{\rm Real}[n,t]\times r[n,j,t]=(1-P[n,t])\times V_{\rm Max}\times r[n,j,t]$, the evaporated volume amount of species $j$ from the first layer ($n=1$) during the time step $\Delta t$, $V_{\rm evap}[1,j,t]$, is given by

[TABLE]

Since the porosity $P[1,t]$ represents the percentage of cavity, $(1-P[1,t])\times V_{\rm Max}$ represent the actual amount of material in the 1st layer. The volume amount of material transferred from $n+1^{th}$ layer to $n^{th}$ layer during the time step $\Delta t$ is calculated as

[TABLE]

Although the default time step $\Delta t$ is 1000 seconds, the time step is half that value when (1) the evaporation rate of at least one species $j$, $\Delta t\times k[n,j,t]$, is less than 1, and (2) the molecular ratio of species $j$, $r[n,j,t]$, is more than 1$\%$. The first criteria will enable us to perform accurate calculation of volatile species, especially H2O. When $\Delta t\times k[n,j,t]$ exceeds 1 and $r[n,j,t]$ is less than 1$\%$, the term of $\Delta t\times k[n,j,t]=1$, corresponding to the situation that all material is transferred to the upper layer during $\Delta t$. Similar to the equation for the evaporation process from the first layer, $k[n+1,j,t]\times V_{\rm Real}[n+1,j,t]$ denotes the total volume amount of thermal evaporation for species $j$. However, this equation represents the blocking of material transfer to the upper layer by multiplying the porosity of the upper layer, $P[n,t]$ (Figure 1).

We assume that a part of the material can evaporate directly from any layer if there is enough porosity in the upper layer (Figure 1). The volume amount of direct evaporation for species $j$ from the nth layer is presented as

[TABLE]

By multiplying the term of $\prod_{k=1}^{k=n}P[k,t]$, this equation permits to evaporate materials directly from the inner layer when the porosity of the upper layers are high enough.

The evaporation processes (5) to (7) determine the volume amount of lost material from the layer. The differentiation of porosity and molecular ratios can be calculated as the result of the above processes. The differentiation of $V_{\rm Real}$ during the time step $\Delta t$ can be calculated as

[TABLE]

Then, $V_{\rm Real}$ after the time step $\Delta t$ is presented as

[TABLE]

Then, recalling the definition of the porosity, new porosity $P[n,t+\Delta t]$ is calculated as

[TABLE]

Since $V_{\rm Real}[n,j,t]=(1-P[n,t])\times r[n,j,t]\times V_{\rm Max}$, this equation is independent of $V_{\rm Max}$. The time evolution of molecular ratios is calculated as

[TABLE]

Finally, we consider merging layers when the porosity gets larger than 0.7 after the above calculations. This value of 0.7 is determined considering the observed porosity of 67P (Pätzold et al., 2016). In addition, the sum of the actual volume amount of layers, $V_{\rm Real(before)}[n,t]+V_{\rm Real(before)}[n+1,t]$, must be smaller than $V_{\rm Max}$ before the merging process so that the volume amount of new layer does not exceed the limitation after merging. Then, using the definition of porosity, the merging process is considered when $P[n,t]+P[n+1,t]>1$ is achieved. With these restrictions, new molecular ratios and the porosity for the $n^{th}$ layer can be calculated from the definition of porosity. We fix $V_{\rm Max}$ during merging but the real volume amount of material is regarded as the sum of the original two layers.

[TABLE]

Then, new molecular ratios can be calculated as

[TABLE]

After the merging process, we reassign the layer number to compensate the merged layer. In this process, the old $n+2^{th}$ layer is regarded as new $n+1^{th}$ layer by applying $r_{\rm new}[n+1,j,t]=r_{\rm before}[n+2,j,t]$ and $P_{\rm new}[n+1,t]=P_{\rm before}[n+2,t]$. The conservation of the total volume amount of material is kept with the above equations.

The timescale of SCCM is also essential since non-volatile species would efficiently concentrate after the long timescale of the heating by the Sun. Ip et al. (2016) calculated the history of orbital evolution of 67P due to the gravitational interactions between giant planets. They suggested that the perihelion of comet 67P changed suddenly in 1959 from 2.7 to 1.3 AU. Before 1959, the perihelion was between 2 and 3 au for 300 years, but it would have the perihelion of $\sim$1.2 au between about 300 about 900 years ago. Therefore, this comet would have experienced strong heat processing by the solar flux for several hundred years. To simulate such a long thermal processing timescale, we calculated the time evolution of the porosity and the molecular ratios during 100 comet years. In the discussion, we will also discuss the case where the comet 67P did not experience the heating before 1959.

Finally, our model is a very simplified one in many points. Although the photo-dissociation process may destroy the molecules, we do not include it in the SCCM modeling. According to the experimental result, molecules would be protected from solar UV photons since they can penetrate only thousands of nanometers (Barnett et al., 2012). However, solar and cosmic particles may affect our results by destroying molecules. In addition, we neglect the sideways and inward diffusion. While the diffusion sideways will be canceled by the diffusion of opposite direction, the inward diffusion may delay the concentration process. As we will see in the next section, the concentration process mainly gets completed within 10 comet years. Therefore, the delay would not be so important if we assume the longer timescale. Our model does not consider the effect of the rotation of 67P as well.

3 SCCM Results

3.1 Parameter Dependency

We performed SCCM simulations under the six different sets of chemical ratios, sizes of layer, and initial porosity. These parameters are summarized in Table 1. We include H2O and CH3NH2, and non-volatile species. As water is the dominant component of cometary gas, it is used to simulate the concentration mechanism of non-volatile species. In addition, we include CH3NH2 as representative of less-volatile species, whose binding energy is higher than water, to test its behavior during the calculation. We also note that CH3NH2 is known as a precursor of glycine. Though the list of included species is limited, they will provide us enough information to discuss the concentration mechanism.

We set the initial porosity of zero for Set 1. Though the porosity of zero is extremely lower than the observation of 67P, we use this porosity to assume an ideal comet, which is free from the heat processing, for the initial condition. We set the volume ratio of non-volatile species to be 50$\%$. Other materials are mainly water, but CH3NH2 is also included in our model to see the behavior of molecules with binding energy higher than water. Considering the simulated abundance ratio of CH3NH2/H2O on interstellar grain, we set the ratio of CH3NH2 to be 0.1$\%$ (Suzuki et al., 2018b). All layers have the same initial molecular ratios regardless of the depth.

We use Sets 2 to 6 to investigate the parameter dependency of the SCCM. We change the initial porosity from zero to 0.2 in Set 2 while the other parameters are same as Set 1. The averaged porosity and the time evolution of the chemical composition of Set 2 (Figures 3 through 6) show a good agreement with Set 1. Therefore, the initial porosity does not affect our results.

The layer thickness is another parameter of the SCCM model. We set the layer thickness to be 5 and 20 cm, respectively, in Set 3 and 4, while the other parameters are same as Set 1. As a result, the time evolution of the porosity and the chemical ratios (Figures 3 through 6) show a small difference, especially in the later phase ($\sim$ 90 comet years) of the simulation. With these results, we fix the porosity to be zero and the layer thickness to be 10 cm in the following discussion.

3.2 Simulation Results

In this section, we discuss the evolution of the comet interior under the assumption of different dust-to-ice ratio. Since it is difficult to estimate the initial dust-to-ice ratio of the comet before the thermal processing of the Sun, the initial dust-to-ice ratio of SCCM is determined from the observation of the comet 67P. Fulle et al. (2017) suggested a high dust-to-ice ratio of 7.5, while other studies reported a lower dust-to-ice ratio of between 0.6 and 7 (Choukroun et al., 2020, and references therein). Therefore, we test three models assuming the different ratios of non-volatile species. In Sets 1, 5, and 6, respectively, the ratio of non-volatile species is assumed to be 10, 50, and 90$\%$. We set the layer thickness of 10 cm and the initial porosity of zero for all cases.

In Sets 1, 5, and 6, the averaged porosity in 0-60, 60-120, and 120-180 cm from the surface get to $\sim$0.5, $\sim$0.6, and $\sim$0.1 within 15 comet years (see Figure 3). The averaged ratios at the different depth beneath the comet surface are shown in Figures 4 to 6. At the depths of 0-60 and 60-120 cm, water evaporates so quickly. The porosity at 0-60 and 60-120 cm is almost 0.5, the ratios of H2O and CH3NH2 are zero, while the ratio of non-volatile species is unity throughout the simulation due to the rapid evaporation of H2O. The evaporation of H2O leads to the concentration of non-volatile species (Figure 6) for all sets. Though the temperature is slightly high for the upper layer, the porosity is small for the deeper layer. As a result, H2O evaporates from the deeper layer and is frozen on the upper layer when the comet is away from the Sun, which increases the H2O ratio while it decreases the porosity and the non-volatile species ratio. On the other hand, H2O ratio decreases when the comet approaches the Sun. This behavior causes the oscillation in the porosity and the ratios. In addition, the dynamic time evolution of the porosity and ratios are seen in Set 5, where the initial ratio of H2O is 0.9. The complete evaporation of H2O will lead to porosity of more than 0.7, the merging process of layers leads to a sudden change of the porosity and the ratios. On the other hand, the merging process is not essential for Set 6, where the ratio of H2O is only 0.1. With this small ratio, the evaporation of H2O has less significance for the merging process. The porosity increases slowly once the ratio of non-volatile species gets high.

The layers beneath more than 180 cm take more time to lose H2O due to the lower temperature. The averaged porosity at 240-300 cm position (the twenty-fourth to thirty-th layer) still changes even after 100 comet years. The evaporation of H2O is still triggering the merging process in this region. The difference of CH3NH2 ratio in Sets 1 and 5 at the 240-300 cm position is due to the higher abundance of non-volatile species in Set 1.

4 Discussion

4.1 Water Production Rate

Since water is the dominant compound of comets, the production rate of water would be the best benchmark of SCCM. In our model, the production rate of jth species at the time of $t$, Production$[j,t]$, is given by summing the all volume amount of evaporation from the surface ($n=1$) and the indirect evaporation from the deeper layer.

[TABLE]

where $\rho$ is the mean density of comet 67P, corresponding to 470 kg m*-3* (Sierks et al., 2015), $d$ is the layer thickness, and $S$ is the total surface area of the comet surface. Assuming that 67P is the sphere with a radius of 2 km, $S$ is 5 $\times$ 107  m2. We calculate the water production rate of Sets 1, 5, and 6 at 100 comet years to compare SCCM results with the result of the Rosetta mission. We show the water production rate of Sets 1, 5, and 6 in Figure 7 by the function of the distance to the Sun. The water production rate in August 2014 was 1.2 kg s*-1* at 3.6 au (Gulkis et al., 2015), while the maximum value is 1.73$\times$102 kg s*-1* after perihelion (Bertaux et al., 2015).

The peak production rates of water in Sets 5, 1, and 6 are 447, 73, and 7.7 kg s*-1*, while they are 11, 1.7, and 0.2 kg s*-1* at 3.6 AU. The result of Set 6 underestimates the production rate. The production rate of Set 1 shows a good agreement with the observation.

It would be worth comparing our model with the previous theoretical modeling study of Keller et al. (2015). They simulated the detailed morphology of 67P, but the thickness of the layer of non-volatile species was assumed to be the micrometer scale. According to their model B, the water production rate at 3.6 AU was 61 kg s*-1*, which is higher than the observation by Gulkis et al. (2015) by order of magnitude. Our agreement on the water production rate emphasizes the validity of SCCM. We evaluated the erosion depth with our model by summing up the production rate along with an orbit followed by the division by the density.

[TABLE]

We find the erosion depths of 0.3, 0.03, and 0.006 m for Sets 1, 5, and 6, respectively, at 100 comet years. These depths are smaller than the result of Keller et al. (2015), where the erosion depths were between 0.67 and 2.9.

4.2 Chemical Composition

We compare the water production rates obtained by Equation (14) with the observation of the Rosetta mission. The production rate ratio, CH3NH2/H2O, is 0.007, 0.004, and 0.002 for Sets 1, 5, and 6, respectively, when the distance to the Sun is more than 2 AU. These ratios agree well with the observed ratio of 0.006 (Goesmann et al., 2015).

4.3 Timescale of the Heating

If comet 67P did not experience any heating before 1959, the heating timescale of 67P would be only eight comet years. We show the water production rates at this age with Sets 1, 5, and 6 in Figure 7. At this early phase of the thermal processing, the water production rates for Sets 1 and 5 are as high as $\sim$103 kg s*-1*, since there is abundant water near the comet surface. The water production rate for Set 6 is $\sim$10 for the perihelion, while it is $\sim$0.6 at 3 au. Therefore, the SCCM model can explain the water production rate if we assume that the ratio of the non-volatile species is high for the interior of 67P. We note that CH3NH2/H2O ratio does not dramatically change with the history of the thermal processing. Without any information of H2O ratio inside 67P, it would be difficult to constrain its history of thermal processing.

4.4 Possible Effect to the Chemistry in the Coma

The comet activity may release the concentrated materials to the coma. Rosetta observed 34 outbursts within three months surrounding its perihelion passage, suggesting this event is ubiquitous at least for 67P (Vincent et al., 2016). Since the outburst is defined by the observable change of the luminosity, smaller activity should be more ubiquitous. Several mechanisms for comet activity would be crystallization of amorphous ice, enrichment of volatile inside the comet, and receding fractured cliffs (e.g., Agarwal et al., 2017; Vincent et al., 2016).

Altwegg et al. (2016) suggested the correlation of glycine abundance with the dust density, giving the possibility that glycine is not directly originated from the nucleus of the comet but evaporated from dust particles, ejected by comet activity. The dust particles are efficiently heated by solar flux (Lien, 1990), resulting in a higher temperature to release the absorbed glycine. Hadraoui et al. (2019) also showed that glycine could originate from icy dust particles with water ejected from the nucleus. Though these studies does not consider the concentration of the non-volatile species, the release of non-volatile species from the comet surface can affect the observed molecular abundances in the coma.

4.5 Effect to the Early Earth Chemistry

It is believed that the comets would have hit the early Earth to provide the materials for the origin of life (Ehrenfreund et al., 2002). Concentration of materials, perhaps in micrometeorites of cometary origin, would provide the unique chemical condition to the volcanic hot spring pools to discuss the chemical evolution toward life. However, the surviving rate of organic materials during the cometary bombardment, especially for the large comets, is still controversial. Pierazzo & Chyba (1999) investigated the thermal decomposition rates of amino acids based on the time evolution of the temperature after the impact with their hydrocode simulations. They showed that the time evolution of the temperature after the impact strongly depends on the position. The opposite side of the impact position is less heated and provides a higher possibility of survival of amino acids.

Figure 3 shows that the timescale of the concentration process for 67P is very short. The concentration of non-volatile species at depths between 0 and 100 cm would be achieved within ten comet years, which corresponds to 65.7 years. As we can see the concentration process for the deeper layers in Figure 3, the typical depth of this concentration would be several meters beneath the surface. Since this timescale is very short for the planet formation, most comets are expected to hit the early Earth after the concentration of non-volatile species. The non-volatile species would include amino acids and other complex and large molecules, which may be essential for the chemical evolution toward the origin of life. Also, the concentration of these materials at the outer layer would show a high survival rate after the cometary impact due to the lower temperature compared to the impact position. The outer layers on the opposite side of the comet would be interesting for further studies.

5 Summary

We investigated the evolution of the composition of the outer cometary layers after the heat processing by the Sun. We developed SCCM model to see the importance of the concentration process of non-volatile species by simulating the evaporation and the merging processes along with the orbit of comet 67P. We found that the SCCM model successfully reproduced the water production rate and CH3NH2/H2O ratio for 67P. SCCM model suggests that non-volatile species may concentrate quickly near the comet surface, while the initial ratio is saved at the deep inside of comets. Comets may have a layer of non-volatile species beneath the several meters from the surface when they hit the early Earth. The outer layer of the comet would be interesting for further studies as they give the unique chemical condition to the volcanic hot spring pools.

Acknowledgements.

We would like to thank the anonymous referee for constructive comments that helped to improve the manuscript. We are grateful to Dr. Takefumi Ootsubo and Dr. Hideyo Kawakita for the discussion. This work was supported by JSPS KAKENHI Grant Number JP16J09590 and 15J10864.