Gas flow in near surface comet like porous structures: Application to 67P/Churyumov-Gerasimenko

We performed an investigation of a comet like porous surface to study how sub-surface sublimation with subsequent ﬂ ow through the porous medium can lead to higher gas temperatures at the surface. A higher gas temperature of the emitted gas at the surface layer, compared to the sublimation temperature, will lead to higher gas speeds as the gas expands into the vacuum thus altering the ﬂ ow properties on larger scales (kilometres away from the surface). Unlike previous models that have used modelled arti ﬁ cial structures, we used Earth rock samples with a porosity in the range 24 – 92% obtained from X-ray micro computed tomography (micro-CT) scans with resolution of some μ m. Micro-CT scanning technology provides 3D images of the pore samples. The direct simulation Monte Carlo (DSMC) method for the rare ﬁ ed gas dynamics is directly applied on the digital rock samples in an unstructured mesh to determine the gas densities, temperatures and speeds within the porous medium and a few centimetres above the surface. The thicknesses of the rock samples were comparable to the diurnal thermal skin depth (5cm). H 2 O was assumed to be the outgassing species. We correlated the coma temperatures and other properties of the ﬂ ow with the rock porosities. The results are discussed as an input to analysis of data from the Microwave Instrument on Rosetta Orbiter (MIRO) on the 67P/Churyumov- Gerasimenko.


Introduction
While the European Space Agency (ESA) Rosetta mission to comet 67P/Churyumov-Gerasimenko (hereafter 67P) has formally been completed, the data analysis and interpretation continues. Rosetta has provided an in-depth investigation, for the first time, orbiting a comet for an extended period of time. For example, the high resolution scientific imaging system on board, OSIRIS, revealed interesting and unexpected surface phenomena (Thomas et al., 2015b). One of the main goals of Rosetta was to study how gas release from comets and the inner coma dynamics are affected by the surface properties. Outgassing from a cometary surface is a complex phenomenon as it takes place in an almost zero gravity regime with rapid expansion form the coma around the nucleus (Gombosi et al., 1985;Kitamura et al., 1985). The details of the coma properties near the nucleus are still under investigation Marschall et al., 2016Marschall et al., , 2017. Production from 67P is dominated by H 2 O while some other major species include CO 2 and CO (Balsiger et al., 1986). This has been confirmed by the mass spectrometry instruments of ROSINA (H€ assig et al., 2015). Outgassing from 67P takes place mostly from the illuminated part of the comet (Filacchione et al., 2016). Activity increases as the comet approaches the perihelion with a peak water production occurring 18-22 days after perihelion (Snodgrass et al., 2016). Ice sublimation from the nucleus requires as input a considerable amount of energy because the latent heat of sublimation of H 2 O is relatively large (Huebner, 2000). One outstanding question is the position of the sublimation front relative to the visible surface. It has been established that the surface area covered by pure (and therefore bright) water ice on comets is much smaller than required to explain the observed total gas production rate (Sunshine et al., 2006). There are two obvious possibilities. Either the water ice is intimately mixed with non-volatile dark material thereby masking its presence in the infrared (Yoldi et al., 2015) or the water ice is mostly below the surface and sublimed gas passes through a porous layer before reaching the visible surface and expanding into the inner coma.
Prior to Rosetta, the evidence was that cometary nuclei had densities lower than that of water ice (Davidsson and Guti errez, 2005) implying significant porosity. Radio science measurements on Rosetta essentially confirmed the low density (P€ atzold et al., 2016) and the porosity is now assumed to be 70-75% on average. Beside the instruments on-board Rosetta, the Philae lander was designed to provide surface data with several sensors. However, because the harpoon anchors failed to fire the collected data are limited. The Multipurpose Sensors for Surface and Sub-Surface Science (MUPUS) aimed to extract important data for the near surface layer (Spohn et al., 2015). Spohn et al. concluded that Philae could not penetrate into the subsurface because of a near-surface layer with a uniaxial compressive strength of the material about 2-4 MPa. This high strength was supported by an observation of (at least locally) a fine-grained porous medium material with porosity of 30-65% (Spohn et al., 2015). These observations have motivated us to study the influence of porosity on the outgassing process. Investigation of the near-surface porous medium layer has been investigated previously (Skorov and Blum, 2012;Skorov et al., 2011Skorov et al., , 2016. In these papers, the porous media layer and targeted porosity were simulated by computer using modelled non-volatile random porous media aggregates. Gundlach et al. performed a detail experimental investigation for the temperature dependent sublimation properties of a hexagonal water ice and the gas diffusion through a dry dust layer covering the ice surface (Gundlach et al., 2011). It was found that the relative permeability of the dust layer is inversely proportional to its thickness. By performing laboratory experiments and an investigation on the heat transport in porous dust layers in vacuum it was found that the head conductivity of porous material in vacuum is determined by the head conduction through solid material and heat transport due to radiation within pores of the material (Gundlach and Blum, 2012).
One of the key reasons for this type of investigation is to study the influence of the porous layer on the production rate and gas temperature at the surface. The thermal conductivity of the porous layer and the resulting higher temperature of the visible surface throttles the gas production. However, a second effect on the coma from the porous layer is the possibility that energy can be transferred from the porous layer to the gas as it passes through the medium. Heating the gas increases the energy available for expansion. The H 2 O outgassing speed has been modelled in the past (Gombosi et al., 1985) and measurements are now available from the MIRO instrument on Rosetta. Speeds vary between 400 and 700 m/s ( In the present paper, we consider H 2 O as the dominant gas for comet outgassing through a porous surface layer. Different porous medium structures with porosity ranging from 24 to 92% are used. However, unlike previous studies, the porous medium structure is directly taken from samples obtained from different Earth basins and measured using Xray micro computed tomography (micro-CT) technology. We use these samples within a gas dynamics code to provide more realistic studies of the influence of porosity on the gas outflow parameters. It should also be noted that, unlike the numerically modelled surface layers, our samples have a structural strength that can be compared to those measured by MUPUS and are thus possibly a better representation of the cometary subsurface.
In the following sections, we shall describe our method of defining porous structures by using micro-computed tomography. We shall then describe the numerical approach for the gas flow including showing results. We conclude with a short set of conclusions.

X-ray micro computed tomography (micro-CT) images
Micro-CT technology offers a three-dimensional representation of a rock with high resolution on the scale of μm. It also offers a highly accurate method to measure rock porosity. The use of X-Ray microcomputed tomography imaging has grown immensely in the past years. The reconstruction of three-dimensional (3D) objects is of great interest in a variety of engineering and medical applications. Initially, the first report of these X-Ray images was made in 1989 by Feldkamp for the representation of bone (Feldkamp et al., 1989). The resolution at the time was considerably low. Micro-CT technology attracted interest in the last decade from the petroleum industry for representation of natural gas reservoirs porous media (Dong and Blunt, 2009;Okabe and Blunt, 2004). It has now become the state-of-the-art method to represent porous medium structure with a resolution of a few microns. Micro-CT scanners use different viewing angles to represent and reconstruct images of the rock (Fig. 1). For the analysis and representation of rock samples micro-CT images are binarised by setting a grey-level threshold, above which voxels are considered as voids and below which voids are considered as porous rock (Campbell and Sophocleous, 2014).
A typical method to study gas transport in near-surface porous layers of a comet is to generate porous media with monodisperse spheres using either random ballistic deposition (RBD) or random sequential packing (RSP) ). In the current study, we use instead micro-CT technology which can offers far better representation of actual porous medium structures with relatively high resolution and with precise results. All images that are used in the current study are high resolution images at a scale of μm. While it is accepted that 67P near surface porous layer consists of high porosity medium structures we extend our investigation to lower porosity levels. Studies of lower porosity allow us to compare results with some of the unexpected observations reported on the 67P (Jia et al., 2017). These observations, suggest that lateral coma winds can move decimetre sized particles above the surface to form dune-like features. Table 1 shows the rock properties for all rock samples used here. Rock porosity has been calculated by the micro-CT scanner. The error on the porosity is below 0.1%.

Gas flow through comet analogue porous medium
Observations of lower than expected porosity close to the surface was supported by Lethuillier et al. who suggested a near-surface porosity of at most 55 per cent and concluded that a porosity increasing with depth could explain the observations (Lethuillier et al., 2016). The fact that Philae failed to penetrate the comet surface may also be attributed to a near-surface layer with relatively low porosity (Spohn et al., 2015). On the other hand, P€ atzold et al. concluded that at least the interior of the "head" of the nucleus has a porosity of 75-85% with the lower value (75%) to be more likely (P€ atzold et al., 2016).
Cometary nucleus consists mostly of dust and water ice (Whipple, 1950). The maximum water production rate as measured by the various Rosetta instruments is 0.5 Â 10 29 molecules m À2 s À1 . Thus, the uniform gas production rate on the sunlight side can reach Z ¼ 3 Â 10 22 molecules m À2 s À1 . For our simulations, we study the comet outgassing near perihelion. The comet surface temperature in illuminated regions at perihelion ranges from~120 to 330 K (Gulkis et al., 2015). The temperature of the rock may be higher than the initial gas temperature if the gas emission is from the sub-surface (Lien, 1990). Surface temperatures in the range (272-336K) were previously reported for 9P on the sunlit hemisphere (Groussin et al., 2007).
For our simulations, we assume our initial sublimation gas temperature as, T gas ¼ 200 K (consistent with the free sublimation temperature of water ice) while the rock takes values T rock ¼ 300 and 330 K. The size of the micro-CT images that are used in the current study have a height, H ¼ 5 cm. These dimensions are comparable to the near surface boundary layer (thermal skin depth). Fig. 2 illustrates the typical set-up of the simulations. The sublimation is taking place 5 cm beneath the surface. Sanctis et al. who studied subsurface temperatures reported that the temperature decreases with depth and they show that nucleus temperature is roughly constant at a distance of 4 cm beneath surface (De Sanctis et al., 2015). The effect of dust mantle thickness was studied previously and it was concluded that the water production of 67P can be attributed to sublimation water ice beneath a dust mantle several millimetres thickness (Hu et al., 2017). By increasing dust mantle thickness, sublimation flux decreases due to the reduce of mantle permeability and thereby, inhibits outgassing (Gundlach et al., 2011).
We evaluate gas properties at the surface and for the next~25 cm where rarefied conditions are found. The degree of rarefaction can be described by the dimensionless Knudsen number (Knudsen, 1909). The Knudsen number is defined as the mean free path (λ) divided with the characteristic length (D), Kn ¼ λ/D. According to the Knudsen classification, the fluid flow can be divided into four regions: continuum flow regime (Kn < 0.001), slip flow regime (0.001 < Kn < 0.1), transition flow regime (0.1 < Kn < 10) and free molecular regime (Kn > 10). For the present study the characteristic length may be considered as the average diameter of the pore of each structure. The corresponding mean free path (λ) for Z ¼ 3 Â 10 22 molecules m À2 s À1 is equal to 2 Â 10 À2 m and the average Knudsen number is equal to 5 which reduces to 0.4 if based on the rock thickness. For a production rate of Z ¼ 3 Â 10 19 molecules m À2 s À1 , the mean free path (λ) is 20 m and the average Knudsen number is 6000 which reduces to 400 if based on the rock thickness. All simulations are therefore in the transition and free molecular regimes.
The direct simulation Monte Carlo method (DSMC) which mimics the Boltzmann equation is the dominant method for simulating rarefied gas flows (Bird, 1994). DSMC is a valid approach for all flow regimes while other methods like Navier-Stokes equations cannot be applied beyond the slip flow regime (0.001 ! Kn 0.1) . DSMC can provide us with gas macroscopic properties such as temperature, velocity and pressure. For our simulations, we use an unstructured mesh with over one million computational cells and over twenty million DSMC particles for each case. The DSMC is applied directly on the computational mesh from the micro-CT images. The direct application of DSMC onto the micro-CT images has previously been reported to be able to describe the flow properties in a wide range of Knudsen number (Christou and Dadzie, 2015). The particles are inserted below the porous medium at a position where the sublimation is assumed to take place. The intermolecular collision model that is used is the Variable Soft Sphere, VSS, which has advantages over the traditional Variable Hard Sphere, VHS, in polyatomic gases (Shen, 2006). The gas-surface interaction model is a full diffuse wall reflection. It has been previously reported that the fully diffuse reflection wall model is better than a specular reflection model (Christou and Dadzie, 2015;Shahabi et al., 2017). The DSMC code implemented in the OpenFOAM framework is used in the current study (Scanlon et al., 2014). The non-destructive micro-CT technology captures the pore space in detail. Modelling flow and transport through the reconstructed 3D porous media image is the next step. Other techniques to model a flow through porous media includes a network modelling approach which is achieved by semi-analytical methods (Dong and Blunt, 2009).
Here, the transport properties are investigated directly using a Direct Simulation Monte Carlo (DSMC) code on the micro-CT images with an unstructured mesh.

Results
Gas properties are calculated as an average over the whole considered volume. Fig. 3 depicts temperature profiles in the normal outgassing direction (X-Axis) for all rock samples while in Table 2 we report the corresponding surface gas overall temperature above the comet surface. The gas surface temperature increases as the rock porosity decreases. At higher porosities, the outgassing surface temperature remains nearly constant over a longer distance above the surface than at low porosities, φ < 45%. Lower porosity rocks are associated with higher surface gas temperature as the gas interacts more strongly with the warmer surface. This may be explained by the fact that, as low porosity implies low level of void, gas passes through rock with low porosity is heated more effectively by the higher temperature rock and expands quickly as it exits. We note that porosity in the range 88-90% records surface gas temperature below 200K for a T rock ¼ 300K and just above 200K in the case for a T rock ¼ 330K. Fig. 3 illustrates in both production rates the rapid decrease of the temperature for the two lowest porosity cases (24% and 45%) at a distance of 15 cm above surface. This is due to the fact that for porosity <45% few molecules are emitted after the medium. As a result, few collisions take place and the energy transferred is less than in the cases with higher porosities. Thus, for low porosity medium we observe higher temperature at the surface and rapid decrease of it at 15 cm above surface due to the small number of particles. This result indicates the importance of the position (above surface) that gas properties are being measured from spacecraft instruments. In Fig. 3 we notice the presence of oscillations which are more noticeable for lower porosities. Such oscillations are produced due to the irregular shape of the rock geometry and the significantly reduced number of intermolecular collisions (especially for lower porosity cases).
Surface gas temperature ranges from 172 to 225K and from 189 to 271K for a medium temperature of 300K and 330K respectively ( Table 2). Shi et al. modelled surface gas temperature and reported large variations depending on rotational period (Shi et al., 2016). For the shadowed it was computed that subsurface temperature is higher than the surface. Surface temperature depicts values between 130K in the shadowed area and 290K in the illuminated area.
Comparing the outgassing temperature at the surface based on the two different rock temperatures there is a clear trend of varying gas temperature as function of porosity. The corresponding fitting profiles are depicted in Fig. 4. We also observe that the highest surface temperature occurs for the maximum rock temperature (330K) at the lowest porosity.
Using the DSMC method and various porosities for different rock samples we examine the average gas pressure at the surface over the whole volume for two different rock temperatures (Fig. 5). For both cases we observe a decrease in gas pressure as rock porosity decreases. The porous medium temperature slightly affects the outgassing pressure as we observe minor increase in outgassing pressure as porous medium temperature increases, from 300K to 330K. Outgassing pressure profiles (Fig. 5) show rapid decrease in pressure for porosities less than 45%. Porous media with low porosity prevent gas molecules to escape to the surface and as a result we observe very low pressures for such cases. Similarly, Skorov et al. showed that the porous media has a diffusion resistance that rapidly grows with increasing thickness (Skorov et al., 2017). It also concluded that only 10% of gas sublimating molecules Table 2 Surface gas temperature as function of rock porosity.

Porosity (%)
Overall T (K) at 1 cm above surface with T rock ¼ 300K Overall T (K) at 1 mm above surface with T rock ¼ 300K Overall T (K) at 1 cm above surface with T rock ¼ 330K Overall T (K) at 1 mm above surface with T rock ¼ 330K   92  172  186  177  189  90  175  193  183  204  88  185  195  187  210  80  195  215  212  227  45  222  242  237  267  24  225  244  enter the coma from a layer (porous medium) with dimensionless thickness 20 R A (R A ¼ Agglomerate Radius). For current configuration and for a constant layer thickness (5 cm) we conclude similar behaviour and that for porosities < 45% very few molecules enter coma. Fig. 6 shows the surface pressure as function of the porosity. Higher pressure is associated with higher porosity. Porosity plays a critical role in the surface pressure and number density. From Table 3, at a rock temperature of 330K the ratio between surface pressure at the 92% porosity rock and surface pressure at 45% porosity rock is about 50. Minor differences exist between the cases with rock temperatures of 300K and 330K. For the low porosities samples, Table 3 and Fig. 5 indicate that they can act as seal rock (or caprock) by preventing the migration of gas molecules to the surface. Seal rock is also found in many geological systems in Earth basins (Watts, 1987). Seal rocks exhibit low porosity (<5%) and low permeability (%1 Â 10 À20 m 2 ) . For the present configuration a height of 5 cm, a porous media with porosity <50% acts as seal rock.
The number density profiles in the flow direction are shown in Fig. 7. The results illustrate that for high porosity rocks the outgassing is denser than for in low porosity rocks. A porosity difference of 12% (92% and 80%) reduces the gas density by a factor of 3. We also observe that the outgassing density drops significantly for porosity <80%. In principle, a medium with a porosity <45% acts as a seal and the gas that passes through it has low density. For the two lowest porosity cases, the density is close to zero and not any variations are observed. This result indicates the importance of measuring both temperature and density from spacecraft instruments in order to be able to calculate the porosity of comets.
Velocity profiles along the outgassing direction are presented in Fig. 8. We observe that for both rock temperatures (T rock ¼ 300K & T rock ¼ 330K) the highest velocity occurs for the highest porosity. For the two highest porosity rocks, we calculate the maximum outgassing gas velocity as~715 m/s. Similar velocity range have been reported previously by other studies (Fougere et al., 2016). In addition, outgassing velocity between the two lowest porosity rock samples, 24% and 45% shows similarities in terms of the decreasing trend after few centimetres away from the surface. The difference in the profiles is smaller within the group 88%, 90%, 92%, than for the group 24% and 45%. Comparing the velocity profiles with Fig. 5 it is noted that the lowest outgassing pressure and velocity occur for the two lowest porosities. Bieler et al. (2015) reported similar velocity with BATS-R-US and is higher than the one from Gulkis et al. (2015) which is equal to 680 m/s. Like in the temperature profiles, the velocity also shows a rapid decrease at a distance of 15 cm above surface for lower porosities.

Discussion
In the first section we demonstrated the capabilities of the micro-CT technology. This technology allows us to represent real Earth core samples with high resolution and is a relatively new tool for the simulation of gas flows through porous media. De Sanctis et al. reported that water vapour temperature near the illuminated area of 67P is 200K (De Sanctis et al., 2015). This value served as the initial temperature of the gas that enters through our micro-CT images. To evaluate the porous media heat transfer process two rock temperatures were used which are similar to what were used in previous experiments (Gundlach and Blum, 2012). Our discussion of the results is presented for the distance up to 25 cm above the surface. There is clearly a heat transfer process from the medium to the passing gas (water vapour) and it is proportional to the rock temperature. As this heat transfer phenomenon occurs, for lower porosity in the medium the higher gas temperature is observed at the surface. Gas passing through high porosity media receives less heating, and goes through an expansion and therefore cooling before exiting the rock.  Similar results were reported by Skorov and Blum (2012). Such observation is noticeable in both cases ( Fig. 3a and b). In general, the temperature profile for the first 25 cm above the surface remains flat and with small fluctuations when the porous medium porosity is above 45%. The fluctuations observed in the temperature profiles for the low porosity cases are due to the low number of intermolecular collisions. Micro-CT samples with porosity <45% prevent molecules from properly migrating to the surface and then results in high oscillations. Based on this, the temperature profile is independent of orbital periods. Outgassing pressure increases with the rock porosity in both cases. Jia et al. reported that in the near surface boundary layer a pressure gradient drives a tangential flow from warm, high-pressure region, to cold, low-pressure regions (Jia et al., 2017). From Table 3, we may conclude that an existence of tangential porosity gradient may also explain strong tangential flows in the boundary layer. Such observation supports those reports and adds a new perspective to the explanation of coma inhomogeneities. In all surface gas property profiles, we observe a sudden drop for the low porosity cases. The sudden drop can be described by the fact that low porosity materials may be treated as seal rocks and prevent molecules to properly migrate to the nucleus surface. Interestingly, we observe that the velocity differs from the temperature profile and the highest velocity is observed for the highest porosity. Comparing Figs. 8 and 3 a negative correlation occurs among all porous samples. In other words, for the highest temperature (low porosity cases) we observe the lower velocity and vice versa. This is the result of an isentropic expansion where temperature and expansion velocity are inversely proportional.

Conclusions
We presented simulations for water vapour flow through porous medium in rarefaction conditions by varying their porosity. Simulations conditions are similar to 67P over the period near perihelion. Porous media are obtained through micro-CT technology. We present details of outgassing properties for a production rate that corresponds to a heliocentric distance of 1.24 au.
We have shown how subliming water molecules passing through a porous medium are heated when the medium is warmer than the sublimation temperature. Temperature rises at a given distance from the source after water vapour travels through the medium and is related to the porosity. The influence of the medium is demonstrated when the temperature of the medium is changed. If the porous medium has a higher temperature, the gas temperature directly above the surface also increases. With a medium with porosity less than 45% the outgassing density nearly vanishes as few molecules are being emitted indicating seal rock features. Effects of the porosity is non-linear. A reduction of 2% in the porosity produces an average surface temperature 6% higher while a reduction by 12% in the porosity can give surface gas temperature 18% higher on average. A reduction by 47% results in an increase of surface gas temperature by 36%. The peak surface gas temperature occurs at a very low porosity (24%) with a peak temperature of 271K. The lowest surface temperature is found for a porosity of 92%, as T ¼ 186K. The current range fells well into those that have previously been modelled for 67P (Shi et al., 2016).
Most importantly, coupled phenomenon exists in terms of temperature and density for porous medium with low porosity. We observe that for the porosity <45% surface temperature is on the highest levels and it rapidly decreases in a distance of 10 cm above surface. However, density profile for φ < 45% is near zero up to a distance of 25 cm above surface. The phenomenon in low porous medium where low density profile and high temperature profile occurs at the same time is due to the fact that limited number of molecules are emitted from such rock samples. On the other hand, density for high porous medium reduce gradually as we move away from the surface. High porosity samples result lower surface temperature, which however remains constant for at least 25 cm above surface. Results from current work suggest that future missions should be equipped with instruments in order to measure density and temperature of the gas at various distances on and above surface. Such measurements may allow scientists to estimate the porosity of the medium based on current results.