High-fidelity subsurface thermal model as part of a Martian atmospheric column model

As the Martian atmosphere is observed in ever greater detail, more realistic computer models are required to interpret these measurements. Physical exchange processes between the atmosphere’s lower boundary and the surface are often simplified. This is because the atmospheric calculations can describe the behaviour of the atmosphere accurately in many cases and simplifying the boundaries saves computing resources. However, the vertical heterogeneity of the subsurface (such as the presence of dust and ice layers) will interact via heat and mass transfer with the atmosphere. Here a new realistic numerical thermal conductivity scheme is introduced for use with a 1-D atmospheric column model useful for investigating the subsurface for layered material and to provide more accurate modelling of the Martian atmosphere. The model with the updated scheme produces results that are identical to the previous versions of the model in identical (non-layered) conditions. The updated model fits well to Viking 1 temperature data from the atmosphere using realistic thermal parameters. Introducing layered material, with different thermal properties, produces noticeable changes in the maximum and diurnal temperatures when changing the thickness of its top layer. The time of maximum surface temperature is only significantly changed when the thickness of the top layer is a moderate fraction of the top layer’s skin depth.


Introduction
In the early days of Martian research, before the days of robotic space exploration, the planet received a lot of attention due to its dynamic surface features as observed through telescopes as reviewed by Sheenan (1996).We now know, from spacecraft observations, that these changing features are due to clouds and dust in the atmosphere (Ruff and Christensen, 2002).These are driven by the thermal and mechanical coupling between the atmosphere and surface (Sagan and Pollack, 1967) and not, for example, by vegetation (Sinton, 1958).It is important to determine the processes that drive the coupling between the surface and atmosphere to understand the water cycle, habitability, evolution of the climate and hazardous conditions for in situ exploration (International Medusa Team et al., 2011).The surface-atmosphere interactions are driven by the solar insolation together with the atmosphere whose properties, such as pressure and temperature, vary on daily and seasonal time scales much more than on Earth (e.g.Harri et al., 1998).Investigations of these Martian surface-atmosphere interactions call for tight coupling of in situ observations and modelling efforts (e.g.Paton, 2012;Harri, 2005).
Layered subsurface thermal schemes have been incorporated in atmospheric models to study the coupling of the Martian atmosphere and it subsurface.These have been used for a variety of purposes such as to understand the global distribution of subsurface water ice, the role of the regolith in controlling the abundance of water vapour in the atmosphere and the constraining of the horizontal/vertical heterogeneity of the Martian surface (Putzig and Mellon, 2007).The models have been used to investigate large scale properties and trends of the subsurface and atmospheric properties using mostly spacecraft observations from orbit around Mars.With an increasing number of in situ surface and near surface temperature measurements by landers (Harri et al., 1998) it will be possible to characterise the thermal response Published by Copernicus Publications on behalf of the European Geosciences Union.
of a variety of surface and subsurface types.This will be useful for interpreting local heterogeneities on near-surface atmospheric measurements, for improving the accuracy of forecasts and as ground truths for interpreting measurements made from orbit.
To more effectively investigate the coupling of the atmosphere and subsurface, we have updated the University of Helsinki (UH) 1-D column atmospheric model.This has been successfully used for characterisation of local atmospheric behaviour from in situ measurements by landers on Mars (Savijärvi, 1999;Savijärvi and Kauhanen, 2008).The model has been used extensively for simulating the Martian atmosphere.It includes predictive equations for the wind components, temperature, specific humidity and ice mixing ratio (Savijärvi, 1999).It includes a long-wave and a shortwave radiation scheme and can take in the effect of CO 2 -H 2 O and dust in the atmosphere.A turbulence scheme is included based on Monin-Obukhov and mixing-length approaches for the lower layers and higher layers, respectively.The thermal diffusion within the subsurface utilises originally a five-layer Crank-Nicholson method and the energy balance at the surface to predict the surface temperature.Sublimation and condensation is modelled from the surface only using a constant soil moisture fraction.
Here a high fidelity numerical thermal conductivity scheme is added that allows thin layers to be modelled each with individual thermal properties.The scheme is designed in such a way that temperature dependent thermal properties can also be simulated.The numerical properties of the thermal scheme and the validation of the column model are investigated for the Viking 1 landing site.This is because the Viking 1 site is well characterised and a long time series of temperature measurements were made at this location.

Thermal properties of the Martian subsurface
The temperature of the Martian surface is controlled by the energy balance at the surface.Heat is transferred upwards and outwards into the environment above the surface, i.e. the "sky", that includes the atmosphere and perhaps obscuring objects that may be in the vicinity (e.g.rocks or parts of a spacecraft structure).The temperature will also be dependent on heat transferred downwards into the subsurface that may be composed of materials that vary in their thermal properties with depth gradually, as with a mixture of materials, or change abruptly, as with layered material.On Mars the dominant heat transfer mechanisms into the "sky" are the emission and absorption of radiation because convection is a relatively inefficient mechanism for the transfer of heat away from the surface due to the low density of the Martian atmosphere.At the Pathfinder site it was found the wind-dependent sensible heat flux was found to be only a few percent of the net outgoing heat flux from the surface (Savijärvi and Määttänen, 2010).At the Viking 1 landing site convection was found to be three times more vigorous than on Earth but still only around 15 % of the net outgoing heat flux from the surface (Sutten et al., 1978).These figures are for relatively calm conditions.However, under certain conditions, such as during dust storm formation, localised convective activity may be more vigorous and large wind speeds, in places, would increase the convective transport of heat from the surface (Fernandez, 1998).
The dominant heat transfer mechanism into the subsurface will be conduction although other mechanisms (like convection and latent heat) will also contribute, depending on the subsurface porosity and the amount of volatiles present.The surface energy balance equation we use is then as follows: where S is the incoming solar flux at the current Mars-Sun distance, a is the albedo, i is the solar incidence angle, ε is the surface emissivity, σ is the Stefan-Boltzmann constant, T 0 is the surface temperature, k T ,z is the temperature and depth dependent thermal conductivity, dT /dz is the vertical temperature gradient at the surface, L is the latent heat of CO 2 frost, dm/dt is the CO 2 frost deposition/sublimation rate, ρ is the density of the atmosphere at the surface, c is the heat capacity of the atmosphere, v is the velocity of the wind and T a is the temperature of the atmosphere.The efficiency of radiative heat transfer is dependent on the optical properties of the surface, such as albedo and emissivity in Eq. ( 1), which are in turn dependent on the regolith physical properties such as composition and microstructure (e.g.Pitman, 2005).The efficiency of heat conducted into the subsurface will depend on the thermal properties, such as conductivity and volumetric heat capacity, which in turn depends on regolith composition and microstructure (e.g.Paton et al., 2012).The structure of the surface material can be consolidated (rock) or granular (sand), which may contain a large amount of voids.In granular material heat transfer may be inhibited or enhanced by material filling in the voids such as gases and liquids (Seiferlin et al., 2008).It is also possible cementing agents could bond the grains together increasing the heat transfer by thermal conductivity (Seiferlin et al., 2003;Piqeux and Christensen, 2009).
The thermal inertia of the regolith is a key parameter that drives the surface-atmosphere exchange processes.It is defined as follows, I = (ρck) 0.5 , where ρ is the density, c is the heat capacity and k is the thermal conductivity.The thermal inertia has been found to be strongly correlated to terrain type (Mellon et al., 2000).In Fig. 1 the surface temperatures are plotted versus the corresponding thermal inertia on the Martian surface calculated using UH 1-D atmospheric model at the location of the Viking 1 site.The daily and seasonal variations of insolation will cause a planetary surface temperature to vary.On Earth the diurnal temperature variation may 20 Figure 1.Maximum, minimum and mean annual surface temperatures at the Viking 1 site plotted against a range of thermal inertia in SI units of J m -2 K-1 s -1/2 .The thermal inertia ranges, represented by A, B and C, correspond approximately to terrain types (r).The main terrain types on Mars are fine dust with few rocks (A), coarse loose particles with scattered rocks (B), coarse sand with strongly crusted fines, abundant rocks and scattered bedrock (C).
Above a thermal inertia of about 1500 J m -2 K -1 s -1/2 the terrain type would most be probably be continuous bedrock or ice.vary a few degrees in regions, surrounded by ocean, and up to 50 K in desert regions (Price, 1977).On Mars, where much of the planet is covered by a desert, the daily temperature variations are larger than on Earth because of the thinner atmosphere.Daily variations in temperature of around 100 K have been measured at an altitude of 1.5 m above the surface at the Mars Pathfinder site (Savijärvi, 1999).The ground, consisting of dust, sand, rock and ice, will experience a similar range of temperatures, if not larger.The thermal inertia is the square root of the product of conductivity and volumetric heat capacity and is a thermal prop- erty that can be estimated remotely from orbiting spacecraft.It is a measure of how quickly and by how much the temperature of a material can be changed for a given amount of energy input.A high thermal inertia indicates a material whose temperature is difficult to change, i.e. it takes a long time to reach a maximum temperature which is relatively low.Conversely, a low thermal inertia indicates a material whose temperature is easy to change, i.e. it takes a short time to reach a maximum temperature which is relatively high.A single temperature measurement can be used to derive the thermal inertia by fitting the points to reference curves from climate models (Fergason et al., 2006).The thermal inertia of the Martian surface shows a clear temperature dependence (e.g.Piqueux and Christensen, 2011).
The temperature dependence of the conductivity of particulate basalt in a Lunar environment is significant, however on Mars the gas phase tends to reduce this effect (Fountain and West, 1970).The variation of heat capacity at low temperatures for basalt is significant in a Lunar environment (Robbie et al., 1970).The gas phase will have a negligible effect in the Martian environment so will contribute more to the temperature dependence of the thermal inertia in a Martian environment.This may cause errors in interpretation of observations of surface temperatures on Mars (Piqueux and Christensen, 2011).
The Martian surface is covered by dust and sand that is composed largely of a silicate framework that is combined with carbonate, sulphate, pyroxene and olivine.Water and carbon dioxide ices dominate the surface in the polar regions.The thermal inertia of the surface of Mars ranges from 30 to 3000 J m −2 K −1 s 0.5 (Jakosky et al., 2000) which represents surface types, in order of increasing thermal inertia, of dust, sand, solid rock and ices.Over seasonal time scales the heat will be conducted into the surface down to maybe 50 cm for dust and several metres for ices.It is known that the near sub-surface of Mars is not vertically homogeneous as evident from excavations by the Viking landers, the MER rovers and Phoenix.For example the Viking lander dug trenches with  depths up to 23 cm and encountered cross bedding, layers of crusts and blocky slabs during excavations (Moore et al., 1982).The MER rover dug trenches with depths of between 6 and 11 cm with its wheels exposing sulphate deposits (Wang et al., 2006).The Phoenix lander dug trenches to depths of a few centimetres in a polygon and dug a trench down to 18 cm in a trough between polygons.Trenches in the polygon uncovered water-ice bearing soils under crusty to cloddy soils but no icy soils were found in the trough between polygons (Ardvidson et al., 2009).Both Phoenix and Viking 2 (48 • N) observed thin layers of frost forming on top of the regolith (Svitek and Murray, 1990;Smith et al., 2009).It may be possible that the topmost surface layer, such as dust on rock, masks the underlying material when using remote thermal measurements to determine the bulk surface composition.Water ice is believed, from interpretation of orbital observations, to exist under the surface beyond the polar regions and to within 25 • of the equator (Vincendon et al., 2010).

A numerical scheme for thermal modelling of the subsurface
Here the numerical thermal conductivity scheme for the subsurface, in the UH 1-D atmospheric model, is updated to enable the modelling of composite materials.The ability to take into account the temperature dependence of their thermal properties is also included.The new scheme is based on the control volume method, described in Patankar (1980).A schematic of the control volumes is shown in Fig. 2. Each control volume can have unique thermal properties of density (ρ), heat capacity (c) and conductivity (k) assigned to them.The thermal conductivity through the boundaries of the control volumes is calculated using a floating mean.For example the thermal conductivity of the upper wall of the control volume centred on z P in Fig. 2 is defined as where k P is the conductivity in the centre of the control volume centred on z P , k N is the conductivity at the centre of the upper control volume and k S is the conductivity at the centre of the lower control volume.The Comparison of surface temperature calculated with and without the atmosphere.The e shows the surface temperature when there is an atmosphere present.The solid line surface temperature for when there is no atmosphere present.1-D heat conduction equation solved by our scheme is the following: Equation ( 2) with S = 0 is integrated over the control volumes to obtain the following expression: The fully implicit version (f = 1) is chosen over the explicit (f = 0) scheme and Crank-Nicholson (f = 0.5) scheme as it allows for larger time steps.It is also easier to implement composite materials and temperature dependent thermal properties.With f = 1 Eq.( 3) becomes the following expression: where However, the explicit method is included as an option within the program.This option was found to be useful as a check for the implicit scheme when using large time steps.The explicit scheme is more sensitive to instabilities and produces easily recognisable unphysical results when run at time steps that are too large.It is easy to become overconfident when using the implicit scheme as the instabilities produced at large time steps by it are not so large and could be misinterpreted as physical properties or processes.
The surface boundary condition is calculated from the net energy flux at the surface, E sf , given by the atmospheric model, and the temperature gradient in the top layer (i.e.across the uppermost control volume boundary): where T sf is the surface temperature and T sub is the temperature of the first control volume below the surface.The deep subsurface boundary is set at a constant value that is calculated from the average yearly temperature which is assumed to be representative for the deep subsurface (e.g.Bense and Kooi, 2004).
An example of the results from calculations, using the new thermal scheme, is shown in Fig. 3.The figure shows a comparison between a homogeneous subsurface and a subsurface with a dust layer on top of a rocky bedrock.While the temperature variation at the surface may be similar for both examples, the thermal behaviour is different in the subsurface.In Sect. 4 we will examine in detail the effect of subsurface layers on the thermal signature of the surface and nearsurface atmospheric temperature.

Numerical limits of the model
Figure 4 shows the sensitivity of the temperature calculated by the model on different parameters such as the subsurface boundary depth and layer thickness.A control model was used that calculated the temperature of the surface over a period of one Martian year.Each parameter was varied in turn and the results of the altered model were subtracted from the control model.The most important model parameters, in order to enable the most efficient use of computational resources, is the number of levels and slab thickness.
The depth of the boundary condition for our high fidelity thermal scheme was investigated by varying the depth of the lower boundary condition.The depth of the boundary level was fixed at the seasonal thermal skin depth.The skin depth is the depth of penetration of a sinusoidal temperature change at the surface.At depths on the order of the thermal skin depth, the temperature will remain constant, essentially representing the average annual temperature of the surface.At much larger depths, the temperature will depend on internal heat sources and subsurface evolution.27 ace temperature differences for temperature dependent and temperature rmal properties for dust, sandy and rocky surfaces.For a model using constant es there is a single value of thermal inertia in the legend.This corresponds to m.For a model using temperature dependent thermal inertia there are two gend.The minimum value corresponds to the minimum temperature and the corresponds to the maximum temperature.The units of thermal inertia for the end are J m -2 K -1 s -0.5 .
Fig. 8. Surface temperature differences for temperature dependent and temperature independent thermal properties for dust, sandy and rocky surfaces.For a model using constant thermal properties there is a single value of thermal inertia in the legend.This corresponds to the value at 2 a.m.For a model using temperature dependent thermal inertia there are two values in the legend.The minimum value corresponds to the minimum temperature and the maximum value corresponds to the maximum temperature.The units of thermal inertia for the values in the legend are J m −2 K −1 s −0.5 .calculated using the following formula: where α is the thermal diffusivity and P is the period of Mars rotation.The thermal diffusivity is defined as follows, α = k/(ρc).The thermal diffusivity of the Viking 1 lander site is most likely somewhere around 2 × 10 −7 m 2 s −1 (Savijärvi, 1995) which is representative of rocky sand.For a Martian year (668 sols), the thermal skin depth would thus be approximately 200 cm for the annual period (and 8 cm for the diurnal period).
The depth of the deep subsurface boundary temperature was increased in the model until the surface temperatures were within 0.1 K of the surface temperature from a control model whose deep level boundary temperature was set at a depth of the seasonal skin depth (200 cm).This condition was found to be met at a depth of 100 cm.The temperature of this boundary condition was set to 220 K which is the average annual surface temperature at the latitude of the Viking 1 landing site (22 • N) and is representative of the lower boundary temperature.
Subsurface thermal characterisation requires thin slabs (i.e.fine grid) for accurate physical modelling of the large temperature gradient near surface and also deeper down in Table 1.Thermal properties used in the simulations of dust layers on rock.The microstructural parameter (Hertz factor) value of 0.09 was used to calculate the dust thermal conductivity.A porosity of 0.5 was used to calculate the dust density.

Thermal property Dust Rock
Conductivity (W m −1 K −1 ) 0.18 2 Density (kg m −3 ) 1500 3000 Heat capacity (J kg −1 K −1 ) 533 533 Volumetric heat capacity (J m −3 K −1 ) 8 × 10 5 16 × 10 5 Thermal Inertia (W m −2 s 1/2 K −1 ) 380 1788 our case for modelling possible fine structures in the subsurface structure.There will be a lower limit to ensure numerical stability of these thin slabs.With the thermal and insolation characteristics at the Viking 1 landing site, a slab thickness of 2 mm was found to be the lowest value for numerical stability and 9 mm is the upper value for a physically realistic model.The stability criterion for the explicit scheme is t<ρc( x) 2 /2k where ρc is the volumetric heat capacity, x is the level thickness and 2 k is the conductivity multiplied by 2. For example, for a volumetric heat capacity of 0.8 × 10 5 J m −3 K −1 and a thermal conductivity of 0.18 W m −1 K −1 the time step for a thickness of 2 mm is 9 s.This result was used to initialise the simulations featured in Fig. 4, except the lower boundary which was set to a depth of 2 m.This then allowed us to explore the effect of increasing the layer thickness and the depth of the lower boundary.Consequently all the simulations of the Martian atmosphere and subsurface in this paper were made with a layer thickness of 2 mm, allowing us to explore the effect of dust layers on rock, and to maintain stability when simulating the atmosphere.
Figure 5 shows the diurnal surface temperatures calculated with the atmosphere present (dotted line) and without the atmosphere present (solid line).The atmospheric heat transfer is calculated by the 1-D column model.The relative difference between the results is of the order of a few degrees (< 5 K).The calculated temperatures results in a reduced amplitude compared to the results under the assumption that no atmosphere is not present.This is to be expected, because the atmosphere absorbs some of the solar radiation and also some of the emitted infrared radiation.

Validation of the model
Figure 6 shows a comparison between the previous version of the UH 1-D model and the updated version with the modified thermal model for the lower boundary.The models were applied to the Viking 1 lander site at a latitude of 22 • N and a solar longitude of zero, L s = 0.A thermal inertia of 380 W m −2 s 1/2 K −1 and an albedo of 0.24 (Savijärvi, 1995) was used.Figure 6  with less layers.The deepest layer compared, 50 mm below the surface, shows a slight discrepancy that is probably due to the finer grid of the new model.
To test our model we then fitted it to Viking 1 temperature data available at the NASA Planetary Data System from an altitude of 1.6 m above the surface and derived the thermal properties of the Martian surface.The Viking 1 temperature data has been investigated extensively in the literature and is fairly well understood (Smith, 2008).The data was averaged over 10 min to filter out the temperature variations due to the turbulence near the surface.The model temperature was fitted to the measured temperature using the least square fitting method and varying the thermal inertia for the L s = 191.
Figure 7 shows a fit of the output of the new model, from an altitude of 2 m, to the Viking 1 temperature data at different times during the first VL-1 Martian year (i.e.summer, autumn, winter and spring).The average thermal inertia for all the fits to the measurements was calculated to be 219 W m −2 s 1/2 K −1 using an albedo, α = 0.18, and an opacity, τ = 0.6.The significant decrease in the temperature variation in the winter is due to the second dust storm that occurred in the first Martian year of Viking 1's operation (Ryan and Henry, 1979).For this we increased the opacity to, τ = 3.0, which corresponds to derived values during the known Viking dust storm.Our model produces the correct range of temperatures but there is an offset in the average temperature which could be due to cooling due to the large scale motion of the atmosphere (i.e.winds) that has not been accounted for in the 1-D model.
To further confirm the model behaviour was stable and producing realistic results, the temperature dependent thermal properties (TDTPs) were activated in the subsurface thermal numerical scheme and tested.This allowed us to explore if TDTPs present any numerical stability issues and to compare with other published works in this area.The temperature dependent thermal conductivity (TDTC) was modelled using a linear relationship that approximates the trend of increasing thermal conductivty with temperature that can be seen in the data for particulate basalt with a particle size between 37 and 62 µm and a bulk density of 1.5 kg m −3 in a simulated Martian environment obtained by Fountain and West (1970).The linear relationship for thermal conductivity is k T = T /60 000 + 1/115 which gives a conductivity value of 0.012 W m −1 K −1 at a temperature of 200 K and a conductivity value of 0.014 W m −1 K −1 at a temperature of 320 K.This amounts to about a 10 % to 20 % change in thermal conductivity over the temperature ranges we are calculating.The temperature dependent heat capacity (TDHC) was modelled using a linear relationship fitted to data obtained from by Robbie et al. (1970) for basalt in a simulated Lunar environment.The linear relationship for TDHC used here is c T = 2T + 120 which gives a TDHC of 520 J kg −1 K −1 at 200 K and 760 at 320 K.This amounts to a 25 % to 50 % change in heat capacity over the range of temperatures we are calculating.The thermal inertia, using these equations and a bulk density of 1.5 kg m −3 , ranges from about 100 to 130 J m −2 K −1 s −0.5 which would correspond to a dusty surface.
Figure 8 compares the surface temperatures from a model with and without TDTPs for dusty, sandy and rocky surfaces.To obtain thermal inertias appropriate for sandy and rocky surface, as published relationships are unavailable, the equation for TDTC was multiplied by a factor of 15 for a sandy surface and for a factor 330 for a rocky surface.The equation TDHC remained unchanged as the effect of gas phase of the specific heat is probably negligible.The model results shown using temperature independent thermal properties (TITPs) used a model with thermal properties values calculated using the equations for TDTPs but keeping the temperature fixed so the temperatures at 2 a.m.conincided.This was done so the temperature difference could be easily seen.The time was chosen because this is commonly used time for fitting model diurnal temperature curves to observed temperature from Mars to determine the surface properties such as grain size.As can be seen in Fig. 8 there is a non-negligible difference for dust and sandy surfaces.
The maximum difference between temperatures calculated using a model with TITPs and a model with TDTP is around 4-5 K for dusty and sandy surface.This agrees well with more realistic simulations by Piqueux and Christensen (2011).The model with TITPs produces higher temperatures, during the day, than produced by the model with TDTPs because as the temperature goes up the thermal inertia increases which in turn tends to cause the material to resist further temperature changes.For the rocky surface the difference between a model with TDTP and a model with TITP is about 0.2 K.This is small because the the temperature variation through the diurnal cycle is also small.We do not discuss the problem of interpretation of Martian temperature observations of the surface here and only note that there is a significant difference between models with and without TDTP.The problem of determining Martian surface properties such as grain size, fitting models with TITP to observations of the Martian surface is discussed in detail by Piqueux and Christensen (2011).They conclude that because the models used for this task use TITPs obtained from laboratory measurements of analogue materials above room temperature, the grain sizes are underestimated when interpreting spacecraft measurements from orbit.
6 Effects of layered material on the surface temperature A vertically heterogeneous subsurface will produce distinct diurnal and seasonal thermal signatures such as differences in the surface temperature variation (Putzig and Mellon, 2007) and time of maximum temperature.This depends on the amount of each material present, its depth below the surface and the thickness of the layer.For example one might expect to find a material on Mars with a low thermal inertia (dust), perhaps of around one skin depth thickness (a few cm), located on top of a slab of high thermal inertia material (rock).As most of the temperature variation will occur in the top layer one may expect to see the typical signature of a dusty material, i.e. low thermal conductivity and low volumetric heat capacity.In the morning the heat from the sun is absorbed and stored in a thin layer near the surface, due to poor conduction into the subsurface.The stored heat then increases the temperature by a relatively large amount due to the low volumetric heat capacity.There is only small time lag before the surface temperature starts to decrease after noon (time of maximum insolation) as the heat is quickly radiated away from the surface.
If the thickness of the dust layer is reduced to some small fraction of the thermal skin depth, say a few mm, the temperature variation observed on the surface will be strongly influenced by the high thermal inertia of the rock underneath.Firstly, the temperature variation will be lower than in a thick dust layer because heat is conducted into the subsurface.The volumetric heat capacity of the near subsurface rock will be high resulting in a small temperature change of the surface for a given amount of energy.The rocky layer beneath the dust will store heat for later release in the afternoon.The surface temperature will continue to rise for some time after noon (time of maximum insolation) as long as the heat conducted upwards to the surface is higher than the heat radiated into the sky.
The updated model was used to investigate the effect of a layered subsurface on the surface and near-surface temperatures.Dust layers of varying thickness were placed on top of rock that was composed of the same rocky material as the dust.Table 1 lists the thermal properties used to calculate the thermal parameters for the dust layer simulations.The dust layer was composed of modelled slabs, each 2 mm in thickness.So for a dust layer of 1 cm thickness there would be five modelled slabs in the model and for a dust layer of 2 cm thickness there would be ten modelled slabs.The rock substrate was modelled with slabs of a thickness of 5 mm.The dust layer was varied from a thickness of 0 cm to 6 cm in steps of 2 mm.Diurnal surface temperatures were then plotted for dust layers varying from 0 cm to 6 cm in 1 cm steps.The maximum temperature was plotted in 2 mm steps to make clear the variation in the lag of the maximum temperature as the dust layer decreses in thickness.
Figure 9 shows the results from the simulations with the dust layers.The figure demonstrates how the temperatures vary for a range of dust layer thicknesses on a rocky substrate underneath.Notice that the surface temperature diurnal range is greatly affected by the thickness of the dust layer while the atmospheric temperature ranges less.The lag of the temperature maximum is not significantly affected until the dust layer is about 2 cm thick.The time of maximum temperature for both the surface and the atmospheric calculations varies over a period of about 1.5 h.This is clearer in the atmospheric temperatures because the curves are closer together.
In Fig. 10 the dust-layer model where the dust layer is set to 2 cm, is compared to a homogeneous material or "rock" which has constant thermal properties in the vertical direction.Even though the amplitude of the temperature variations is similar in all cases there is a significant lag between the layered material and the solid material.This is presumably due to the larger volumetric heat capacity of the "rock" and its ability to store the heat for later release in the afternoon.

Concluding remarks
A high fidelity numerical thermal conductivity scheme was included in the UH 1-D atmospheric model to enable the inclusion of composite materials in the subsurface (e.g.dustice layers) and to more accurately model the temperature change with time on the surface.The model with the new thermal scheme reproduces the results from the model with the previously established thermal scheme in identical conditions.The numerical limits of the new thermal scheme were explored and the results were found to comply with the wellknown stability limits of the applied numerical method that have been applied.The model was fitted to spacecraft (VL-1) data, over diurnal periods, with physically sensible parameters and found to agree with previously published material.
The new model was run with a layer of dust of various thicknesses on top of rock to investigate features in the temperature variation with time and determine if dust layers on rock could be distinguished.Hourly temperature measurements of the surface or near surface can be used to determine if there is a dust cover or not by virtue of a lag in the maximum temperature.
The model could be applied to spacecraft data to detect thermal signatures from dust layers on rock or ice.The model could also presumably be used to more accurately simulate and predict the atmospheric conditions from past, present and future in situ measurements.

Fig. 1 .
Fig.1.Maximum, minimum and mean annual surface temperatures at the Viking 1 site plotted against a range of thermal inertia in SI units of J m −2 K −1 s −1/2 .The thermal inertia ranges, represented by A, B and C, correspond approximately to terrain types (r).The main terrain types on Mars are fine dust with few rocks (A), coarse loose particles with scattered rocks (B), coarse sand with strongly crusted fines, abundant rocks and scattered bedrock (C).Above a thermal inertia of about 1500 J m −2 K −1 s −1/2 the terrain type would most be probably continuous bedrock or ice.

Fig. 2 .
Fig. 2. Control volumes for modelling the heat transfer in the subsurface.

Figure 3 .Fig. 3 .
Figure 3. Diagram showing results for temperature of the subsurface with time obtained from 12 the new scheme.On the left the entire subsurface is made of dust.On the right the first 5 mm 13 is dust with the rest being rock.14 15 16 17 18 19 20 21 22 23 24 25 26

Figure 4 .
Figure 4.The sensitivity of the surface temperature to variations in the depth of the fixed temperature lower boundary condition (left) and the sensitivity of the model on the thickness of the layers (right).The difference in surface temperature for the left hand figure was calculated by running the model with different boundary depths and subtracting the results from a control simulation whose boundary temperature was set at a depth of 2 m.The right hand figure was calculated by subtracting the results from runs of different level thicknesses from a control simulation with a thickness set at 2 mm.

Fig. 4 .
Fig. 4. The sensitivity of the surface temperature to variations in the depth of the fixed temperature lower boundary condition (left) and the sensitivity of the model on the thickness of the layers (right).The difference in surface temperature for the left hand figure was calculated by running the model with different boundary depths and subtracting the results from a control simulation whose boundary temperature was set at a depth of 2 m.The right hand figure was calculated by subtracting the results from runs of different level thicknesses from a control simulation with a thickness set at 2 mm.

Fig. 5 .
Fig. 5. Comparison of surface temperature calculated with and without the atmosphere.The dotted line shows the surface temperature when there is an atmosphere present.The solid line shows the surface temperature for when there is no atmosphere present.

Figure 6
Figure4shows the sensitivity of the temperature calculated by the model on different parameters such as the subsurface boundary depth and layer thickness.A control model was used that calculated the temperature of the surface over a period of one Martian year.Each parameter was varied in turn and the results of the altered model were subtracted from the control model.The most important model parameters, in order to enable the most efficient use of computational resources, is the number of levels and slab thickness.The depth of the boundary condition for our high fidelity thermal scheme was investigated by varying the depth of the lower boundary condition.The depth of the boundary level was fixed at the seasonal thermal skin depth.The skin depth is the depth of penetration of a sinusoidal temperature change at the surface.At depths on the order of the thermal skin depth, the temperature will remain constant, essentially representing the average annual temperature of the surface.At much larger depths, the temperature will depend on internal heat sources and subsurface evolution.The skin depth can be

Fig. 6 .Figure 7 .
Fig. 6.Comparison between the previous version of the 1-D model and the model with the updated thermal numerical scheme.The continuous lines are the previous version of the model and the crosses (+) are from the new version of the model.In the left hand figure is the surface temperature and the temperature at 2 m altitude.On the right is the surface temperature, the temperature at a depth of 14 mm and the temperature at a depth of 50 mm.

Fig. 7 .
Fig. 7. Martian atmospheric temperature measured by the Viking 1 lander compared to the modified column model.The same thermal inertia was used for each plot (219 W m −2 s 1/2 K −1 ).The thermal inertia may be expected to change during the seasons due to exchange of volatiles between the regolith and the atmosphere as well as a varying amount of dust on the surface.The amplitude of the temperature variations from the model match those of the observations quite well.The absolute temperatures also fit quite well except for Ls 282 (bottom-left) where the deep subsurface boundary was technically set at 170 K to obtain a reasonable match whereas the others where set at 220 K.

Figure 9 .Fig. 9 .
Figure6shows a comparison between the previous version of the UH 1-D model and the updated version with the modified thermal model for the lower boundary.The models were applied to the Viking 1 lander site at a latitude of 22 • N and a solar longitude of zero, L s = 0.A thermal inertia of 380 W m −2 s 1/2 K −1 and an albedo of 0.24(Savijärvi, 1995)   was used.Figure6shows that the new model produces more or less identical results to the previous version of the model