Soil Thermal Conduction in the Generic PCM

Purpose

The goal of this page is to summarize key aspects, both theoretical and practical, of the computation of subsurface temperatures in the generic PCM. Three subroutines are used for that:

• The 'module' comsoil_h.F90 which contains the parameters to define the subsurface grid.
• The subroutine soil_settings.F90 which reads/initializes the soil depths and properties.
• The routine soil.F which computes soil temperatures and the flux from the ground to the surface.

Theoretical considerations

The computation of the heat conduction in the soil was based on the 1D layer soil model originally described by Hourdin et al. (2013). Since that, E. Millour modified this scheme to include the possibility of varying the thermal properties with depth.

The PCM solves the time-dependent diffusion equation: \begin{equation} C \frac{\partial T}{\partial t} = - \overrightarrow{\nabla} \cdot \overrightarrow{F}_c \label{Eq:EqChaleur} \end{equation}

where $$C$$ is the volumetric specific heat (J m$$^{-3}$$ K$$^{-1}$$) such that $$C = \rho C_p$$ ($$\rho$$ is the material's density, in kg m$$^{-3}$$ and $$C_p$$ is its specific heat, in J kg$$^{-1}$$ K$$^{-1}$$), and $$F_c$$ is the conductive heat flux: $$\overrightarrow{F}_c = - \lambda \overrightarrow{\nabla}T$$ according to Fourrier's law ($$\lambda$$ is the solid’s heat conductivity, in J s$$^{-1}$$ m$$^{-1}$$ K$$^{-1}$$).

Thermal conduction is here considered as a one-dimensional (1D) process, lateral conduction is neglected (a reasonable assumption considering typical grid size in the Generic PCM). Temperature $$T$$ of the soil is thus a function of time $$t$$ and depth $$z$$, which must satisfy the following equation:

\begin{equation} C \frac{\partial T}{\partial t} = \frac{\partial}{\partial z}\left[ \lambda(z) \frac{\partial T}{\partial z} \right] \label{Eq:Chaleur1D} \end{equation}

The Mars PCM used to have the following boundary conditions:

• At the bottom (e.g.: $$z = H$$) of the layer of soil: No outgoing (or incoming) heat flux. This boundary condition is then simply: \begin{equation} \frac{\partial T}{\partial z} _{z= H} = 0 \end{equation}
• At the surface ($$z$$ = 0): Surface temperature $$T_S (t) = T (z = 0, t)$$ can be computed from the balance of heat fluxes through the surface and cooling thereof, which leads to:

\begin{equation} F_c + \sum F^\downarrow = \epsilon \sigma T_s^4 \end{equation} where $$F_c$$ is the flux of heat conduction given by Fourier's law, $$\sum F^\downarrow$$ is the sum of fluxes reaching the surface, $$\epsilon$$ is the emissivity of the ground and $$\sigma$$ the Stephan-Boltzman constant. In the case of condensation/sublimation, a latent heat term is added to this last equation.

Numerical approach

The computation of the heat conduction, and the updates of soil temperatures, are done in the routine soil.F.

The spatial discretization of the unsteady-heat equation is based on a finite volume approach. The vertical domain is divided into nsoilmx discrete layers; $$z_k$$ denote interlayer depths, and $$z_{k-1/2}$$ ($$k = 1, . . . ,$$ nsoilmx - 1) are the computational nodes representative of the $$k^{th}$$ layer. The computational grid is obtained by imposing a geometrically stretched distribution of layers: \begin{equation} z_{k-1/2} = lay1_{soil} \times 10^{-4} \times \alpha_{soil}^{k-1/2} \end{equation}

where $$k$$ varies from 1 to nsoilmx. By default, the variables nsoilmx, $$lay1_{soil}$$, and $$\alpha_{soil}$$ are set respectively to 18, 2 $$ \times 10^{-4}$$, and 2 in the routine comsoil_h.

The time integration is done via an implicit (first order) Euler scheme which approximates the differential equation $$du(t)/dt = F (u, t)$$ as: \begin{equation} \frac{u(t+\delta t) - u(t)}{\delta t} = F (u(t+\delta t), t+\delta t)) \end{equation}

The complete description of the heat equation discretization can be found in Lange (2024) - pages 59 to 61 -.

The lower boundary condition (Eq. 3) is included as such in the solver, but the upper boundary condition (Eq. 4) is not. Surface temperature is technically linked to atmospheric and ground processes which are coupled. Rather than solving the coupled problem, it is uncoupled in a way that allows one to solve atmospheric and ground processes separately. The essential feature of the artificial uncoupling is that, when solving the heat diffusion equation in the soil at time $$t+\delta t$$, the value of surface temperature $$T_s (t+\delta t)$$ is already known. In short, the algorithm works as follows: the surface temperature at $$t+\delta t$$ is computed with a radiative balance, knowing $$ F_C (t+\delta t)$$ . Then soil temperatures are updated. $$F_C (t+2\delta t)$$ is finally computed and used as an input to compute $$T_s (t+2\delta t)$$.

Description of the soil subroutines

comsoil_h

This module defines the parameters for the soil discretization

References

Hourdin, F., Van, P. L., Forget, F., and Talagrand, O. (1993). Meteorological variability and the annual surface pressure cycle on Mars. Journal of Atmospheric Sciences, 50(21):3625 – 3640.

Lange, L. (2024). Studying and modeling the dynamics of water and CO$$_2$$ ice on Mars. PhD dissertation. Link to the manuscript