# Subsurface Thermal Modeling of Oxia Planum, Landing Site of ExoMars 2022

Authors: M. Formisano; M. C. De Sanctis; C. Federico; G. Magni; F. Altieri; E. Ammannito; S. De Angelis; M. Ferrari; A. Frigeri
Journal: Advances in Astronomy
 (2021)

---

## Abstract

Numerical simulations are required to thermophysically characterize Oxia Planum, the landing site of the mission ExoMars 2022. A drilling system is installed on the ExoMars rover, and it will be able to analyze down to 2 meters in the subsurface of Mars. The spectrometer Ma_MISS (Mars Multispectral Imager for Subsurface, Coradini and Da Pieve, 2001) will investigate the lateral wall of the borehole generated by the drill, providing hyperspectral images. It is not fully clear if water ice can be found in the subsurface at Oxia Planum. However, Ma_MISS has the capability to characterize and map the presence of possible ices, in particular water ice. We performed simulations of the subsurface temperatures by varying the thermal inertia, and we quantified the effects of self-heating. Moreover, we quantified the heat released by the drilling operations, by exploring different frictional coefficients and angular drill velocities, in order to evaluate the lifetime of possible water ice.


---

## Body

## 1. Introduction

The target of the ExoMars 2022 mission is to investigate an ancient location with high potential for past habitability or for chemical biosignatures—indicators of extant life—as well as abiotic/prebiotic organics [1–3]. This location has been identified in Oxia Planum, a 200 km-wide low-relief terrain [4], at a latitude of about 18∘. Oxia Planum shows evidence of subaqueous episodes being characterized by hydrous clay-bearing bedrock units and layered outcrops [4–6]. The entire area of Oxia is undergoing erosion, attested by morphology (inverted features) and crater statistics [4]. This site is also important for the complex environment reconstruction (see, for example, [7]). Moreover, the surface rocks of Oxia are less exposed to oxidation and radiation compared with other sites since they have been exposed by the cosmic bombardment only recently [3, 4].The drilling system onboard the ExoMars rover will be able to penetrate down to 2 meters in the subsurface of Oxia Planum [3, 8]. The main aim of the drill is to acquire subsurface samples with high astrobiological potential to be delivered to the analytical lab for a detailed investigation, but before the sample extraction, the sample context is analyzed by Ma_MISS (Mars Multispectral Imager for Subsurface) that will observe the lateral wall of the borehole generated by the drill, with the production of hyperspectral images [9, 10]. Ma_MISS will provide information about the mineralogy, oxidation state, and hydration state of the sample before the extraction and crushing of the sample, an important aspect since deep cold samples may be altered by the drilling activity and also after their extraction. Moreover, Ma_MISS has been validated by laboratory tests [9, 11], and it can infer information on grain size using the analysis of specific spectral parameters. Also, another ExoMars instrument, such as Close-Up Imager (CLUPI) [12] can give information of the fines that drill will extract.Buried ice in the subsurface is not expected at the Oxia Planum landing site; however, recent results from the Fine-Resolution Epithermal Neutron Detector (FREND) instrument on board Trace Gas Orbiter (TGO) have pointed out the presence of ice permafrost oases near the Martian equator [13]. It is important to understand the thermal state of the subsurface at Oxia Planum and the influence of the drilling activity. The characterization and the mapping of the organic and volatile content in the subsurface have a great scientific impact [14]: therefore, numerical simulations are required to produce, for example, maps of temperature of the surface/subsurface of Oxia Planum. Also, evaluating the heat contribution released by the drilling system, we will establish if hypothetical water ice in the subsurface could be preserved or not [15, 16], having in mind that the low atmospheric pressure on Mars leads to very fast sublimation of any ice-rich deposit in direct contact with the drill tip [3].The estimation of the thermal inertia of Oxia Planum is a matter of debate. In [17], a high-resolution map of thermal inertia derived from the observations of planetary brightness temperature by the Mars Global Surveyor (MGS) Thermal Emission Spectrometer (TES) is presented; in that work, Oxia Planum exhibits thermal inertia between 300 and 400 thermal inertia units, hereafter TIU. Jones et al. [18] suggested a thermal inertia of about 260 TIU, interpreting the surface as consisting of dark fines, with some coarse sand and duricrust but very little dust, with grain size below 3 mm. However, in [4], for the clay-bearing unit, thermal inertia in the range of 550–650 TIU was reported, and this probably represents the likely landing site; river delta, instead, is characterized by a low thermal inertia (<100 TIU). Finally, in [19], the clay unit is characterized by a thermal inertia of 370–380 TIU while the river delta by a value of 300 TIU.In this work, from a thermophysical point of view, we characterized the landing site Oxia Planum. In particular, we investigated firstly the thermal response of the subsurface to different thermal inertia values, taking into account the considerations discussed above and after we evaluated the heat released by the drilling operations. We also estimated the lifetime of a hypothetical icy subsurface deposit. The paper is structured as follows: in Section2, we report the numerical model adopted in this work, in Section 3, the results obtained by the numerical simulations are given, and finally, in Section 4, the discussion and the conclusion are given.


## 2. Numerical Model

We performed our simulations using a 3D finite element code (e.g., [24–28]), which solves the classical heat equation in a parallelepiped (see Figure 1) compatible with the dimension of a borehole, generated by the drilling system, installed on the ExoMars rover, in Oxia Planum.Figure 1
Domain of integration for the first part of this work: zoom of the first 5 cm from the Martian surface. A radiation boundary condition is imposed at the top (x-y plane), while zero flux is applied on the other sides.The top (x-y plane) of this domain is modeled with a Gaussian random surface in order to simulate the roughness of the surface. The dimensions of the domain are 1 cm × 1 cm × 50 cm. A depth of 50 cm has been chosen since the diurnal skin depth is reasonably of the order of tens of centimeters and also compatible with the length of a rod of the drill onboard ExoMars [29].We solve the following equation:(1)ρcp∂T∂t=K∇⟶⋅∇⟶T,where T is the temperature, t the time, ρ the density, cp the specific heat, and K the thermal conductivity. Heat transfer occurs only by conduction since convection is negligible due to the small temperature gradients involved as well as the characteristic size of the sample. At the top, we imposed a radiation boundary condition, while on the other sides, zero heat flux is imposed.In particular at the top, following, e.g., [30, 31], we solve for each facet(2)Sc1−acosZ+QSH=−Kn⟶⋅∇⟶T+εσT4,where Sc is the solar constant (scaled for Mars’ heliocentric distance) in Wm−2, a is the albedo, cosZ is the cosine of the solar incidence, ε is the emissivity, and σ is the Stefan–Boltzmann constant. The illumination conditions (i.e., cosZ) are calculated according to [24, 26, 28]. As in [26], we consider a diffuse surface, which absorbs the solar irradiation and emitted IR radiation as a grey body with an emissivity of 0.97, an approach similar to [32]. The term QSH is the term related to the so-called “self-heating.”The initial temperature is set at 200 K, which is compatible with the surface equilibrium temperature.We developed three different models, characterized by different thermal inertia: model A, characterized byI = 160 TIU; model B characterized by I = 255 TIU; and model C characterized by I = 650 TIU. Models A and B are characterized by a thermal conductivity of 0.018 Wm−1K−1 and 0.048 Wm−1K−1, respectively. These values are compatible with those provided by the InSight mission [22]; the density, in these cases, is low and set at 1700 kgm−3, which corresponds to a high-porous sedimentary rock [21]. Model C is characterized by a high thermal conductivity of 0.2 Wm−1K−1 and a density typical of the clay, i.e., 2700 kgm−3 [23]. In all the models, the specific heat is set at 800 J kg−1K−1, which is a value compatible with many materials (i.e., regolith, fine sand, and coarse sand, e.g., [20]). The albedo is 0.13 [18].The heliocentric distance explored is the aphelion since the landing is currently planned for June 10, 2023, while Mars is very close to its aphelion.All the thermophysical parameters used in this work are listed in Table1.Table 1
Thermophysical parameters adopted in the first part of this work. Models A and B correspond to a thermal inertia of 160 TIU and 255 TIU, respectively. They are characterized by a composition similar to high-porous sedimentary rocks. Model C, instead, is characterized by a high thermal inertia (650 TIU), compatible with a composition clay dominant (e.g., vermiculite).
ParameterValueReferenceDistance1.67 AUThis studyInitial temperature200 KThis studySpecific heat800 Jkg−1K−1[20]Albedo0.13[18]Latitude18°[4]Model A (I = 160 TIU)Density1700kgm−3[21]Thermal conductivity0.018Wm−1K−1[22]Model B (I = 225 TIU)Density1700kgm−3[21]Thermal conductivity0.048Wm−1K−1[22]Model C (I = 650 TIU)Density2700kgm−3[23]Thermal conductivity0.2Wm−1K−1This study


## 3. Results

We start the discussion about the results of part I, i.e., the influence of the thermal inertia on the subsurface temperatures of Oxia Planum, evaluating also the effects of self-heating. We recall that the thermal inertia is defined as follows:(3)I=Kρcp.Part II, instead, is related to the estimation of the heat released in the subsurface of Oxia Planum by the drilling operations. We finally provide an estimation of the lifetime of the hypothetical ice-rich sphere in contact with the tip of the drill.


### 3.1. Influence of the Thermal Inertia

In the top panels of Figure2, we report the temperature plot vs. time at different depths from the surface for models A, B, and C (respectively, from left to right). The temperatures are calculated along a line in the z-direction perpendicular to the surface passing through the center. In model A, the surface temperature oscillates between 175 K and 265 K, while in model B, it oscillates between 185 and 255 K, a result very similar to model A since thermal conductivities of these models are of the same order of magnitude. Conversely, in model C, the oscillations of the surface temperature are reduced because of the higher thermal conductivity and consequently the higher thermal inertia: in this case, the surface temperature is in the range 200–245 K.Figure 2
In panels (a–c), we show the temperature profiles versus time for the three models developed in this work. Different depths from the surface are explored. Note that, for model C, more depths are reported than models A and B since the skin depth is greater. The time reported on thex-axis is that required to reach the steady state, so we report only the last planetary rotations. In panels (d–f), we show the temperature profile vs. depth at the particular moment of the Martian day: p1 represents the minimum point in the “day-night” temperature plot, p3 the maximum point, and p2 and p4 the inflexion points before and after the maximum, respectively. Also, see Figure 3.
(a)(b)(c)(d)(e)(f)In the bottom panels of Figure2, we show the temperature vs. depth profile at different moments of the Martian day: we select the point of minimum temperature (p1) in the “day-night temperature profile,” the inflexion point (p2) before the maximum (p3), and the inflexion point (p4) after the maximum. See Figure 3. These plots provide information about the skin depth, i.e., the depth at which the amplitude of the thermal wave is attenuated by a factor 1/e. We recall that the diurnal skin depth is defined as(4)δ=IρcpPπ,where P is the rotational period. Since the temperature profiles converge at about 30 cm in models A and B, we can deduce that this value represents the skin depth; in model C, instead, the skin depth is about 40 cm.Figure 3
Points to which the profiles of temperature vs. depth in Figure2 refer to. The points correspond to the minimum (p1), the inflexion point (p2) before the maximum (p3), and the inflexion point (p4) after the maximum. This day-night profile is merely indicative.


### 3.2. Influence of the Density

In this section, we would like to demonstrate that, fixed the thermal conductivity (and the specific heat), if we use different values for the density, the surface temperature does not change in a significant way as shown in Figure4(a). As a test case, we adopt model C: we compare the surface temperature in case of density of 1900 kgm−3, 2700 kgm−3, and 3000 kgm−3. The value 1900 kgm−3 could correspond, for example, to hydrated vermiculite while 3000 kgm−3 to dehydrated nontronite or chlorite [23]. We observe no differences between the case with 2700 kgm−3 and 3000 kgm−3 and a slight difference with the case with 1900 kgm−3 in the maximum and minimum values reached: for instance, the maximum values differ of less than 5 K. We conclude that the density affects the surface temperatures only in a marginal way, unlike the thermal conductivity as previously shown since these parameters can change in a very significant way.Figure 4
(a) Temperature dependence on density: we used model C as the test case. (b) Evaluation of self-heating contribution in case of model A used as the test case: red plot includes self-heating, while the blue plot does not.
(a)(b)


### 3.3. Influence of Self-Heating

In order to evaluate the contribution due to self-heating among the facets of the shape model under study, we carried out a further simulation using model A as the test case but without self-heating. In Figure4(b), we show the result, plotting only the surface temperature: we observe that the peak difference is about 10 K. Self-heating is taken into account with the algorithm of the hemicube, implemented in the COMSOL Multiphysics code used for these simulations. In Figure 5, we report the 3D Martian surface temperature map for model A in order to point out how the temperature can increase where the surface roughness is high. Maps refer to one planetary rotation.Figure 5
3D temperature maps along a Martian day for model A. Maps refer to the Martian surface. The time on the top left is time elapsed from the beginning of simulation.


### 3.4. Heat Released by the Drilling

The second part of this work concerns the estimation of the heat released during the drilling operations. This contribution is important to evaluate the lifetime of possible volatile species in the subsurface of Mars. The geometry adopted in this second part of the work is shown in Figure6.Figure 6
Geometry adopted in the second part of this work. A cylinder with a radius of 13 mm (the thickness of the drill tip [29]) simulated the hole of the drill. The depth is assumed to be equal to 50 cm, compatible with the length of a rod [29].To take into account this contribution in a simple way, we imposed a heat frictional flux on the sides, “in contact” with the drill (represented by the cylinder of Figure6), of our domain of integration. This heat flux is defined by the following equation (e.g., [33, 34]):(5)Φ=ηFnωrπrhWm−2,where η is the coefficient of friction, ω is the rotational velocity, r is the radius of the “hole,” i.e., the thickness of the tip, and h is the height of the domain of integration (50 cm). The coefficients of friction explored are 0.3 and 0.9. Coefficient of friction represents a sort of heat transfer efficiency, from the drill to the soil. The value of 0.3 is compatible with the ones used in [35]: we also test a frictional heating coefficient close to the unit, i.e., 0.9. We select two values for the rotation for minutes (hereafter rpm): 30 and 60 [29]. The total integration time is 90 min, while the time step is 10 s. We simulate a “drilling window” as in Figure 7: we start with 30 min in “on mode,” followed by 30 min in “off mode” and, finally, other 30 min in “on mode.” Some assumptions are made:(i)
Instantaneous drilling (drill excavation is not modeled)(ii)
Constant thrust(iii)
Constant rotational velocityFigure 7
Drilling window adopted in this work: first 30 minutes in “on mode,” followed by 30 minutes in “off mode” and, finally, other 30 minutes in “on mode.” The values on they-axis represent a sort of efficiency. The heat flux described in equation (5) is multiplied for this function: there is a gradual transition between on mode and off mode.For the domain of integration (i.e., the Martian subsurface), we used the values of model C, i.e.,K = 0.2 Wm−1K−1, ρ=2700kgm−3, and cp=800Jkg−1K−1; however, also using the parameters of models A and B, the results are very close.In Table2, we report the thermophysical parameters adopted in the second part of this work.Table 2
Thermophysical parameters adopted in the second part of this work.
ParameterValueReferenceThrust300 (N)[29]rpm30–60[29]Friction coefficient0.3∗∗–0.9This studyTip density3500 (kgm−3)[35]Tip thermal conductivity540 (Wm−1K−1)[35]Tip specific heat790 (Jkg−1K−1)[35]Drill activity90 (min)This study∗∗Compatible with [35].In Figure8, we report the 1D temperature profile along the y-direction at x = 0 and z = 50 cm (maximum depth; see Figure 6): in the x-axis, we report the distance from the hole. The shadowed red rectangle identifies the drill. At z = 50 cm, we want to point out that the solar input is negligible, so the initial temperature at that temperature, before the drilling, is reasonably close to the equilibrium temperature, i.e., 200 K. We selected four times (20, 43, 65, and 90 minutes). We observe that, in the case of μ=0.3 and rpm = 30, the maximum increase in temperature due to the drilling operations is about 55 K. The increase results in about 75 K if we use rpm = 60. If we adopt a frictional heating coefficient of 0.9, the increase in temperature is about 90 K for rpm = 30 and very high, more than 100 K, for rpm = 60.Figure 8
(a)μ=0.3 and rpm = 30; (b) μ=0.3 and rpm = 60; (c) μ=0.9 and rpm = 30; (d) μ=0.9 and rpm = 60. For the Martian subsurface, we used K = 0.2 Wm−1K−1, ρ=2700kgm−3, and cp=800Jkg−1K−1.
(a)(b)(c)(d)


#### 3.4.1. Lifetime of Subsurface Ices

At this point, we can calculate the rate of sublimation of hypothetical subsurface ices from one basis of the cylinder (i.e., the borehole). We applied the following formula [36]:(6)Γ=Psatμ2πRT,where Γ is expressed in kg m−2s−1. R is the universal gas constant, while μ is the water molecular weight. The water ice saturation pressure is defined as in [37]:(7)Psat=exp9.550426−5723.265T+3.53068lnT−0.00728332T,valid for T>110 K.Equation (6) holds at the surface: in order to take into account the effects of crust thickness h, we multiply Γ by rp/h, as in [38]. We set rp=3 mm as the expected grain size at Oxia Planum [18]. In Figure 9, we observe that, in the “hottest case” (η=0.9 and rpm = 60), the sublimation rate is tens of kg m−2s−1, while in the “coldest case” (η=0.3 and rpm = 30), the sublimation rate is about 0.2 kg m−2s−1.Figure 9
Sublimation rate for the explored cases of Figure8. The rate of sublimation is calculated through one basis of the cylinder representing the borehole.As in [39], we can calculate the lifetime of an ice deposit of mass m0. For simplicity, we use a spherical form with initial radius equal to the radius of the drill tip. The variation of the radius as a function of the time is related to the sublimation rate by the following formula:(8)drdt=−Γρice,where for ρice, we used 950 kg m−3. After integrating equation (8) and using the Taylor series around r0, we obtain the relative loss of ice mass at time t after the deposit is raised to temperature T:(9)mtm0=1−Γtρicerp,wheret is the time expressed in seconds. By using the maximum sublimation rate Γ, we obtain that, in case of η=0.3 and rpm = 30, the fraction of initial mass remaining after 90 minutes is about 60 %, in case of η=0.3 and rpm = 60, it is about 1%, while in the other two cases (η=0.9 and rpm = 30 or 60), we have a complete loss of the initial mass.


## 3.1. Influence of the Thermal Inertia

In the top panels of Figure2, we report the temperature plot vs. time at different depths from the surface for models A, B, and C (respectively, from left to right). The temperatures are calculated along a line in the z-direction perpendicular to the surface passing through the center. In model A, the surface temperature oscillates between 175 K and 265 K, while in model B, it oscillates between 185 and 255 K, a result very similar to model A since thermal conductivities of these models are of the same order of magnitude. Conversely, in model C, the oscillations of the surface temperature are reduced because of the higher thermal conductivity and consequently the higher thermal inertia: in this case, the surface temperature is in the range 200–245 K.Figure 2
In panels (a–c), we show the temperature profiles versus time for the three models developed in this work. Different depths from the surface are explored. Note that, for model C, more depths are reported than models A and B since the skin depth is greater. The time reported on thex-axis is that required to reach the steady state, so we report only the last planetary rotations. In panels (d–f), we show the temperature profile vs. depth at the particular moment of the Martian day: p1 represents the minimum point in the “day-night” temperature plot, p3 the maximum point, and p2 and p4 the inflexion points before and after the maximum, respectively. Also, see Figure 3.
(a)(b)(c)(d)(e)(f)In the bottom panels of Figure2, we show the temperature vs. depth profile at different moments of the Martian day: we select the point of minimum temperature (p1) in the “day-night temperature profile,” the inflexion point (p2) before the maximum (p3), and the inflexion point (p4) after the maximum. See Figure 3. These plots provide information about the skin depth, i.e., the depth at which the amplitude of the thermal wave is attenuated by a factor 1/e. We recall that the diurnal skin depth is defined as(4)δ=IρcpPπ,where P is the rotational period. Since the temperature profiles converge at about 30 cm in models A and B, we can deduce that this value represents the skin depth; in model C, instead, the skin depth is about 40 cm.Figure 3
Points to which the profiles of temperature vs. depth in Figure2 refer to. The points correspond to the minimum (p1), the inflexion point (p2) before the maximum (p3), and the inflexion point (p4) after the maximum. This day-night profile is merely indicative.


## 3.2. Influence of the Density

In this section, we would like to demonstrate that, fixed the thermal conductivity (and the specific heat), if we use different values for the density, the surface temperature does not change in a significant way as shown in Figure4(a). As a test case, we adopt model C: we compare the surface temperature in case of density of 1900 kgm−3, 2700 kgm−3, and 3000 kgm−3. The value 1900 kgm−3 could correspond, for example, to hydrated vermiculite while 3000 kgm−3 to dehydrated nontronite or chlorite [23]. We observe no differences between the case with 2700 kgm−3 and 3000 kgm−3 and a slight difference with the case with 1900 kgm−3 in the maximum and minimum values reached: for instance, the maximum values differ of less than 5 K. We conclude that the density affects the surface temperatures only in a marginal way, unlike the thermal conductivity as previously shown since these parameters can change in a very significant way.Figure 4
(a) Temperature dependence on density: we used model C as the test case. (b) Evaluation of self-heating contribution in case of model A used as the test case: red plot includes self-heating, while the blue plot does not.
(a)(b)


## 3.3. Influence of Self-Heating

In order to evaluate the contribution due to self-heating among the facets of the shape model under study, we carried out a further simulation using model A as the test case but without self-heating. In Figure4(b), we show the result, plotting only the surface temperature: we observe that the peak difference is about 10 K. Self-heating is taken into account with the algorithm of the hemicube, implemented in the COMSOL Multiphysics code used for these simulations. In Figure 5, we report the 3D Martian surface temperature map for model A in order to point out how the temperature can increase where the surface roughness is high. Maps refer to one planetary rotation.Figure 5
3D temperature maps along a Martian day for model A. Maps refer to the Martian surface. The time on the top left is time elapsed from the beginning of simulation.


## 3.4. Heat Released by the Drilling

The second part of this work concerns the estimation of the heat released during the drilling operations. This contribution is important to evaluate the lifetime of possible volatile species in the subsurface of Mars. The geometry adopted in this second part of the work is shown in Figure6.Figure 6
Geometry adopted in the second part of this work. A cylinder with a radius of 13 mm (the thickness of the drill tip [29]) simulated the hole of the drill. The depth is assumed to be equal to 50 cm, compatible with the length of a rod [29].To take into account this contribution in a simple way, we imposed a heat frictional flux on the sides, “in contact” with the drill (represented by the cylinder of Figure6), of our domain of integration. This heat flux is defined by the following equation (e.g., [33, 34]):(5)Φ=ηFnωrπrhWm−2,where η is the coefficient of friction, ω is the rotational velocity, r is the radius of the “hole,” i.e., the thickness of the tip, and h is the height of the domain of integration (50 cm). The coefficients of friction explored are 0.3 and 0.9. Coefficient of friction represents a sort of heat transfer efficiency, from the drill to the soil. The value of 0.3 is compatible with the ones used in [35]: we also test a frictional heating coefficient close to the unit, i.e., 0.9. We select two values for the rotation for minutes (hereafter rpm): 30 and 60 [29]. The total integration time is 90 min, while the time step is 10 s. We simulate a “drilling window” as in Figure 7: we start with 30 min in “on mode,” followed by 30 min in “off mode” and, finally, other 30 min in “on mode.” Some assumptions are made:(i)
Instantaneous drilling (drill excavation is not modeled)(ii)
Constant thrust(iii)
Constant rotational velocityFigure 7
Drilling window adopted in this work: first 30 minutes in “on mode,” followed by 30 minutes in “off mode” and, finally, other 30 minutes in “on mode.” The values on they-axis represent a sort of efficiency. The heat flux described in equation (5) is multiplied for this function: there is a gradual transition between on mode and off mode.For the domain of integration (i.e., the Martian subsurface), we used the values of model C, i.e.,K = 0.2 Wm−1K−1, ρ=2700kgm−3, and cp=800Jkg−1K−1; however, also using the parameters of models A and B, the results are very close.In Table2, we report the thermophysical parameters adopted in the second part of this work.Table 2
Thermophysical parameters adopted in the second part of this work.
ParameterValueReferenceThrust300 (N)[29]rpm30–60[29]Friction coefficient0.3∗∗–0.9This studyTip density3500 (kgm−3)[35]Tip thermal conductivity540 (Wm−1K−1)[35]Tip specific heat790 (Jkg−1K−1)[35]Drill activity90 (min)This study∗∗Compatible with [35].In Figure8, we report the 1D temperature profile along the y-direction at x = 0 and z = 50 cm (maximum depth; see Figure 6): in the x-axis, we report the distance from the hole. The shadowed red rectangle identifies the drill. At z = 50 cm, we want to point out that the solar input is negligible, so the initial temperature at that temperature, before the drilling, is reasonably close to the equilibrium temperature, i.e., 200 K. We selected four times (20, 43, 65, and 90 minutes). We observe that, in the case of μ=0.3 and rpm = 30, the maximum increase in temperature due to the drilling operations is about 55 K. The increase results in about 75 K if we use rpm = 60. If we adopt a frictional heating coefficient of 0.9, the increase in temperature is about 90 K for rpm = 30 and very high, more than 100 K, for rpm = 60.Figure 8
(a)μ=0.3 and rpm = 30; (b) μ=0.3 and rpm = 60; (c) μ=0.9 and rpm = 30; (d) μ=0.9 and rpm = 60. For the Martian subsurface, we used K = 0.2 Wm−1K−1, ρ=2700kgm−3, and cp=800Jkg−1K−1.
(a)(b)(c)(d)


### 3.4.1. Lifetime of Subsurface Ices

At this point, we can calculate the rate of sublimation of hypothetical subsurface ices from one basis of the cylinder (i.e., the borehole). We applied the following formula [36]:(6)Γ=Psatμ2πRT,where Γ is expressed in kg m−2s−1. R is the universal gas constant, while μ is the water molecular weight. The water ice saturation pressure is defined as in [37]:(7)Psat=exp9.550426−5723.265T+3.53068lnT−0.00728332T,valid for T>110 K.Equation (6) holds at the surface: in order to take into account the effects of crust thickness h, we multiply Γ by rp/h, as in [38]. We set rp=3 mm as the expected grain size at Oxia Planum [18]. In Figure 9, we observe that, in the “hottest case” (η=0.9 and rpm = 60), the sublimation rate is tens of kg m−2s−1, while in the “coldest case” (η=0.3 and rpm = 30), the sublimation rate is about 0.2 kg m−2s−1.Figure 9
Sublimation rate for the explored cases of Figure8. The rate of sublimation is calculated through one basis of the cylinder representing the borehole.As in [39], we can calculate the lifetime of an ice deposit of mass m0. For simplicity, we use a spherical form with initial radius equal to the radius of the drill tip. The variation of the radius as a function of the time is related to the sublimation rate by the following formula:(8)drdt=−Γρice,where for ρice, we used 950 kg m−3. After integrating equation (8) and using the Taylor series around r0, we obtain the relative loss of ice mass at time t after the deposit is raised to temperature T:(9)mtm0=1−Γtρicerp,wheret is the time expressed in seconds. By using the maximum sublimation rate Γ, we obtain that, in case of η=0.3 and rpm = 30, the fraction of initial mass remaining after 90 minutes is about 60 %, in case of η=0.3 and rpm = 60, it is about 1%, while in the other two cases (η=0.9 and rpm = 30 or 60), we have a complete loss of the initial mass.


## 3.4.1. Lifetime of Subsurface Ices

At this point, we can calculate the rate of sublimation of hypothetical subsurface ices from one basis of the cylinder (i.e., the borehole). We applied the following formula [36]:(6)Γ=Psatμ2πRT,where Γ is expressed in kg m−2s−1. R is the universal gas constant, while μ is the water molecular weight. The water ice saturation pressure is defined as in [37]:(7)Psat=exp9.550426−5723.265T+3.53068lnT−0.00728332T,valid for T>110 K.Equation (6) holds at the surface: in order to take into account the effects of crust thickness h, we multiply Γ by rp/h, as in [38]. We set rp=3 mm as the expected grain size at Oxia Planum [18]. In Figure 9, we observe that, in the “hottest case” (η=0.9 and rpm = 60), the sublimation rate is tens of kg m−2s−1, while in the “coldest case” (η=0.3 and rpm = 30), the sublimation rate is about 0.2 kg m−2s−1.Figure 9
Sublimation rate for the explored cases of Figure8. The rate of sublimation is calculated through one basis of the cylinder representing the borehole.As in [39], we can calculate the lifetime of an ice deposit of mass m0. For simplicity, we use a spherical form with initial radius equal to the radius of the drill tip. The variation of the radius as a function of the time is related to the sublimation rate by the following formula:(8)drdt=−Γρice,where for ρice, we used 950 kg m−3. After integrating equation (8) and using the Taylor series around r0, we obtain the relative loss of ice mass at time t after the deposit is raised to temperature T:(9)mtm0=1−Γtρicerp,wheret is the time expressed in seconds. By using the maximum sublimation rate Γ, we obtain that, in case of η=0.3 and rpm = 30, the fraction of initial mass remaining after 90 minutes is about 60 %, in case of η=0.3 and rpm = 60, it is about 1%, while in the other two cases (η=0.9 and rpm = 30 or 60), we have a complete loss of the initial mass.


## 4. Discussion and Conclusions

In this work, we explored the dependence of the thermal inertia on the subsurface temperatures of Oxia Planum, the landing site of the mission ExoMars 2022, and we also evaluated the heat released by the drilling operations. In particular, in the first part, we investigated the influence of the thermal inertia on the subsurface temperatures by testing three models: model A, characterized by a thermal inertia of 160 TIU, compatible withK = 0.018 Wm−1K−1; model B characterized by a thermal inertia of 255 TIU, compatible with K = 0.048 Wm−1K−1; and model (C) characterized by a thermal inertia of 650 TIU, compatible with K = 0.2 Wm−1K−1. Models A and B are moreover characterized by a low density, i.e., 1700 kgm−3, compatible with a composition made of high-porous sedimentary rocks; model C, instead, is characterized by a density of 2700 kgm−3, compatible with a clay dominant composition (e.g., vermiculite). The heliocentric distance considered in these simulations is the aphelion since the landing is currently planned, while Mars is very close to the aphelion.In model A, the surface temperature ranges from 175 K to 265 K, while in model B, it ranges from 175 K to 255 K. The behaviour is very similar since thermal conductivities of these models are of the same order of magnitude. In these cases, the skin depth is close to 30 cm. The situation is different for model C: here, in fact, the oscillations of the surface temperature are limited between 200 K and 245 K because of the high thermal inertia of this case (i.e., 650 TIU), and the skin depth is close to 40 cm. We also verified the effects of self-heating among the facets of the shape model by modeling the surface of our domain as a Gaussian random surface. The difference between the case with self-heating and without self-heating is about 10 K: we used model A as a test case.The second part of the paper is dedicated to test the effect of the drilling activities on the subsurface thermal environment. We used different parameters to simulate the drilling characteristics as well as the subsurface conditions, including the presence of ice patches.Heat released by the drilling operations has been simulated in different cases of friction: in the case of low frictional coefficient (η = 0.3), the maximum temperature increase at a depth of 50 cm is about 55 K (with rpm = 30) and about 75 K (with rpm = 60) after 90 minutes. In the case of high frictional heating (η = 0.9), the maximum temperature increase becomes 90 K (with rpm = 30) and 130 K (with rpm = 60). We also evaluated the behaviour of a hypothetical area with water ice in the subsurface: we calculated the sublimation rate and the lifetime of a hypothetical spherical ice-rich deposit, with a radius of 13 mm (the radius of the drill tip), using the temperature profiles obtained by the simulation of the drilling activity. After 90 minutes, the remaining ice mass is about 60% in case of η=03 and rpm = 30, while in the other simulated cases, the remaining mass is negligible. However, a more complete treatment of the diffusion of possible volatiles in the subsurface of Mars will be addressed in the next paper.The simulations reported here suggest that the thermal environment of the Martian subsurface can be strongly influenced by the drilling activity, with consequent changes in the characteristics of the subsurface areas and layers in the immediate vicinity of the borehole and in contact with the drill tip and rods. These thermal changes, induced by the heat released by the drilling, can contribute to variations in the original volatile content of the subsurface layers, including the possible loss of water ice, if initially present.The results obtained by our numerical simulations offer a first depiction of thermal state, after the drilling, of the subsurface of Oxia Planum, because some assumptions were made in this work, in particular, we considered the drilling as an instantaneous process without treating excavation. Moreover, we applied the flux defined by equation (5) in an equal way to all the walls representing the contact sides with the drill. However, in order to increase the possibility of survival of any volatiles inside the subsurface of Oxia, we would reduce frictional heating or adopting appropriate window drilling. Future improvements will be made as well as a validation with laboratory experiments.



---

*Source: 9924571-2021-09-02.xml*