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

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. &e 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.


Introduction
e 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][2][3]. is 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 claybearing bedrock units and layered outcrops [4][5][6]. e entire area of Oxia is undergoing erosion, attested by morphology (inverted features) and crater statistics [4]. is 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].
e drilling system onboard the ExoMars rover will be able to penetrate down to 2 meters in the subsurface of Oxia Planum [3,8].
e 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. e 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]. e 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) ermal 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. e paper is structured as follows: in Section 2, 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.

Numerical Model
We performed our simulations using a 3D finite element code (e.g., [24][25][26][27][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. e top (x-y plane) of this domain is modeled with a Gaussian random surface in order to simulate the roughness of the surface. e 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: where T is the temperature, t the time, ρ the density, c p 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 where S c is the solar constant (scaled for Mars' heliocentric distance) in Wm − 2 , a is the albedo, cos(Z) is the cosine of the solar incidence, ε is the emissivity, and σ is the Stefan-Boltzmann constant. e illumination conditions (i.e., cos(Z)) 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]. e term Q SH is the term related to the so-called "self-heating." e 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 by I � 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 − 1 K − 1 and 0.048 Wm − 1 K − 1 , respectively. ese 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 − 1 K − 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 − 1 K − 1 , which is a value compatible with many materials (i.e., regolith, fine sand, and coarse sand, e.g., [20]). e albedo is 0.13 [18]. e 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 Table 1.

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 selfheating. We recall that the thermal inertia is defined as follows: 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.

Influence of the ermal Inertia.
In the top panels of Figure 2, we report the temperature plot vs. time at different depths from the surface for models A, B, and C (respectively, from left to right). e 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.
In the bottom panels of Figure 2, 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. ese 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 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.

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 Figure 4(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 . e 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.

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 Figure 4(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.

Heat Released by the Drilling.
e second part of this work concerns the estimation of the heat released during the drilling operations. is contribution is important to evaluate the lifetime of possible volatile species in the subsurface of Mars. e geometry adopted in this second part of the work is shown in Figure 6.
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 Figure 6), of our domain of integration.
is heat flux is defined by the following equation (e.g., [33,34]):   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. e time reported on the x-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.

Advances in Astronomy
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). e 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. e 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]. e 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 velocity For the domain of integration (i.e., the Martian subsurface), we used the values of model C, i.e., K � 0.2 Wm − 1 K − 1 , ρ � 2700 kgm − 3 , and c p � 800 Jkg − 1 K − 1 ; however, also using the parameters of models A and B, the results are very close.
In Table 2, we report the thermophysical parameters adopted in the second part of this work.
In Figure 8, 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. e 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. e 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.

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]: where Γ is expressed in kg m − 2 s − 1 . R is the universal gas constant, while μ is the water molecular weight. e water ice saturation pressure is defined as in [37]: 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 (r p /h), as in [38]. We set r p � 3 mm as the expected grain  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. e depth is assumed to be equal to 50 cm, compatible with the length of a rod [29]. On mode Off mode Figure 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." e values on the y-axis represent a sort of efficiency. e heat flux described in equation (5) is multiplied for this function: there is a gradual transition between on mode and off mode. 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 − 2 s − 1 , while in the "coldest case" (η � 0.3 and rpm � 30), the sublimation rate is about 0.2 kg m − 2 s − 1 .
As in [39], we can calculate the lifetime of an ice deposit of mass m 0 . For simplicity, we use a spherical form with initial radius equal to the radius of the drill tip. e variation of the radius as a function of the time is related to the sublimation rate by the following formula: where for ρ ice , we used 950 kg m − 3 . After integrating equation (8) and using the Taylor series around r 0 , we obtain the relative loss of ice mass at time t after the deposit is raised to temperature T: where t 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.

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 with K � 0.018 Wm − 1 K − 1 ; model B characterized by a thermal inertia of 255 TIU, compatible with K � 0.048 Wm − 1 K − 1 ; and model (C) characterized by a thermal inertia of 650 TIU, compatible with K � 0.2 Wm − 1 K − 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). e 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. e 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. e 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. e difference between the case with self-heating and without self-heating is about 10 K: we used model A as a test case. e 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. e 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. ese 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. e 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.

Data Availability
e data used to support this study are included at https:// github.com/MiFormisano/Oxia_Planum.  Advances in Astronomy