Planets - User contributions [en]

Physics of the Generic PCM

2025-03-05T12:51:07Z

LucasLange: /* Subsurface thermal */

This page describes the various physical parametrizations of the Generic PCM and the chronology of their call through the physical iteration. This chronology is important because some variables need to be updated by certain processes before others (examples).
During one physical iteration, the code passes through multiple sub-routines each encapsulating a parametrization. The sub-routines usually take as arguments:
*the dynamical values of the state variables
*the dynamical tendencies of the state variables
*any additional relevant variable
*any additional relevant tendency
and it returns in general tendencies (of the state variables as well as of any other relevant variable).

Work in progress. Need to add links.

=Initialization=
==First call==
Some initializations only need to be done at the very first iteration (e.g. examples).
==Each call==
Some other initializations need to be done at each iteration (e.g. examples).

=Radiative transfer=
==Correlated-k==
The main radiative transfer solver of the Generic PCM implements the correlated-k method, which provides a flexible and quick way to solve radiative transfer equations, particularly suited for GCMs. More informations here.
==Newtonian relaxation==
Newton was a very relaxed physicist, still inspiring us today.
==No atmosphere==
If you have no atmosphere, why do you need a GCM?

=Vertical diffusion=
Compute the vertical diffusion due to turbulence in the planetary boundary layer
==vdifc==
The "old" vertical diffusion routine.
==turbdiff==
The new (and improved!) vertical diffusion routine

=Convection=
Convective mixing in an atmosphere column involves non-resolved, sub-grid processes. These processes are parameterized using the following routines:

==Thermal plume==
This module implements a more accurate and sophisticated parametrization of convection. More info [[Thermal_plume_model_Generic_PCM|here]].

==Dry convection==
More info [[Convective_adjustment_scheme_in_the_generic_PCM|here]].
==Non-orographic gravity waves==
More info [[Non orographic gravity waves drag|here]].

=CO2 condensation=
Inheriting from the Mars PCM where it is the background species, CO2 condensation is treated as a dedicated step in the Generic PCM's physics. Does that only concern CO2 as a background gas? Or does it work also if CO2 is a non-background variable gas? In any case more info here.

=Tracers=
Many things can be advected by tracers in the Generic PCM, like chemical species or aerosols. Physical processes involving tracers are parameterized using the following routines:

==Volcano==
This routine parameterizes a source of tracers corresponding to volcanic eruption. More information here.

==Water/ice==
Water aerosols (liquid or solid) are created (resp. consumed) by condensation (resp. vaporization or sublimation), consuming (resp. releasing) latent heat in the atmosphere. In the Generic PCM, the (atmospheric part of the) water cycle is handled by various routines, as explained here.

==Photochemistry==
Chemistry can turn molecules into other molecules, by the action of temperature (thermochemistry) or UV light (photochemistry). This is handled by the photochemistry routine, described [[Photochemistry|here]].

==Generic condensation==
On Earth, only water condenses in the atmosphere, but on other planets (which the Generic PCM aims at simulating), many other chemicals can undergo state change. To take that into account, the Generic PCM has a flexible scheme to deal with any arbitrary species, as desbribed [[Radiative_Generic_Condensable_Specie|here]].

==Sedimentation==
What goes up must come down, as explained [[Sedimentation_of_tracers_in_the_generic_PCM|here]].

==Updates==
This section essentially takes care that the condensation of a major species (e.g. water vapor in steam-rich atmospheres) affects other species, as described here.

==Slab ocean==
This routine solves for big fish eating small fish, as desbried [[Slab_ocean_model|here]].

==Surface==
The surface part of the water cycle is handled here, as the hydrology page describes.

=Heat Conduction in the Subsurface=
The conduction of heat in the subsurface is solved as described [[Soil_Thermal_Conduction_in_the_Generic_PCM|here]].

=Diagnostics/write outputs=
Nothing very physical here, just writting the outputs! By the way, if you want to know how to output new variables in the [[diagfi.nc]] or [[XIOS]] file, check out [[Outputs|this page]].

Soil Thermal Conduction in the Generic PCM

2025-03-05T12:48:20Z

LucasLange: /* Inhomogeneity of the subsurface */

== 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.

To compute the subsurface temperatures, one must set ''callsoil'' to .true. in the ''callphys.def''

== 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 vertical grid discretization, and the default values:
<syntaxhighlight lang="fortran" line>
integer :: nsoilmx = 18 ! The number of points in the subsurface.
real :: lay1_soil = 2.e-4 ! The depth (m) of first full soil layer
real :: alpha_soil = 2 ! The geometric factor for the soi layer thickness
</syntaxhighlight>

The vertical discretization is defined in two arrays:
<syntaxhighlight lang="fortran" line>
real :: layer(nsoilmx) ! Full soil layer depths
real :: mlayer(0:nsoilmx-1) ! Mid soil layer depths
</syntaxhighlight>

The soil thermal properties are defined by:

<syntaxhighlight lang="fortran" line>
real :: volcapa ! soil volumetric heat capacity, C
real :: inertiedat(ngrid,nsoilmx) ! Soil thermal inertia
</syntaxhighlight>

==== soil_settings ====

This subroutine reads from a NetCDF file (opened by the caller) of "startfi.nc" format the soil properties (number of grid points, soil temperature, thermal inertia), and eventually process to interpolation if needed.

* 1. The number of soil layers is read in the NetCDF file. If the number of soil layer (''dimlen'') read is not equal to ''nsoilmx'', then the soil temperature and thermal inertia will be interpolated to the grid from the PCM defined in comsoil_h.
* 2. The vertical grid is read in the NetCDF file, or created if it does not exist, or if ''dimlen'' $$\neq$$ ''nsoilmx''.
* 3. The volumetric specific heat ''volcapa'' is read in the controle array (''controle(35)'') of the NetCDF file. If it does not exist, then ''volcapa'' is set to $$10^6$$.
* 4. The Thermal inertia is read in the NetCDF file. If it does not exist, the subroutines crashes.
* 5. The Soil Temperature is read in the NetCDF file. If it does not exist, the subroutines crashes.
* 6. Interpolation of soil temperature and thermal inertia is done if ''dimlen'' $$\neq$$ ''nsoilmx''.

'''Inputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
integer :: nid ! Name of the NetCDF file
integer :: ngrid ! Number of surface grid points.
integer :: nsoil ! Number of soil vertical grid points.
integer :: indextime ! position of the first time point in the netCDF file.
real :: tsurf(ngrid) ! Surface temperatures (used for the interpolation)
</syntaxhighlight>

'''Outputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
real :: tsoil(ngrid,nsoil) ! Soil Temperature
</syntaxhighlight>

Note that the soil thermal inertia are modified within the function but not are not output.

==== soil ====
This subroutine computes the soil temperatures from surface temperatures and heat fluxes from the ground. The code works as follow:
* 1. During the first call, local variables (soil thickness, soil conductivity at the interface, etc.) are defined.
* 2. The first grid-point soil temperature is computed from the surface temperature and ground heat flux.
* 3. Other soil temperatures are then computed based on a recurrence relation
* 4. Ground heat flux for next time step is computed.

'''Inputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
integer :: ngrid ! Number of surface grid points.
integer :: nsoil ! Number of soil vertical grid points.
logical :: firstcall ! Boolean checking if it is the first call to this subroutine.
logical :: lastcall ! Boolean checking if it is the last call to this subroutine (not used actually).
real :: tsurf(ngrid) ! Surface temperatures (used for the interpolation)
real :: therm_i(ngrid,nsoil) ! Soil thermal inertia
</syntaxhighlight>

'''Outputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
real :: tsoil(ngrid,nsoil) ! Soil Temperature
real :: fluxgrd(ngrid) ! Soil heat flux
real :: capcal(ngrid) ! surface specific heat capacity
</syntaxhighlight>

== Useful Tricks ==
==== Changing the vertical grid ====
To change the vertical grid, one can changes these parameters in the run.def:
<syntaxhighlight lang="bash">
## Number of grid points in the subsurface
nsoilmx = ...
## Thickness of topmost soil layer (m)
lay1_soil = ...
## Coefficient for soil layer thickness distribution
soil_thickness = ...
</syntaxhighlight>

A word of caution: the soil grid must be adapted to catch the diurnal/seasonal heat propagation in the ground. In the code (''soil.F''), the first layer must be shallower than 0.5 x the diurnal skin depth, and larger than 2 times the seasonal skin depth. The skin depth $$\delta_P$$ (m) is calculated with:

\begin{equation}
\delta_P = \sqrt{\frac{P \lambda}{\pi C}}
\label{eq:skindepth}
\end{equation}
where $$P$$ is the diurnal/seasonal period (in seconds).

==== Inhomogeneity of the subsurface ====

The scheme developed by E. Millour allows to have varying thermal properties with depth. When a layer of soil is not homogeneous, but instead composed of 2 separate sub-layers (e.g., dry regolith above massive subsurface ice), the layer's equivalent physical properties, i.e. values specified at node $$z_{k+1/2}$$, should be adjusted accordingly.

The equivalent thermal conductivity of a composite layer is derived from the steady-state composite wall analogy. The thermal conductivity $$\lambda_{k+1/2}$$ of a composite layer consisting of a layer (extending from depths $$$z_k$$ to $$z_{k+1}$$) of thermal conductivity $$\lambda_k$$ above a layer (extending from depths $$z$$ to $$z_{k+1}$$) of thermal conductivity $$\lambda_{k+1}$$ is set to:

\begin{equation}
\lambda_{k+1/2} = \frac{z_{k+1}-z_k}{\frac{z-z_k}{\lambda_k} + \frac{z_{k+1}-z}{\lambda_{k+1}}}
\label{eq:conductivitycomposite}
\end{equation}

The equivalent volumetric specific heat $$C_{k+1/2}$$ of a mixed layer must be such that the overall ability of the volume to store internal energy is equal to the sum of the individual sub-layers. $$C_{k+1/2}$$ is thus simply obtained as:

\begin{equation}
C_{k+1/2} = C_k + \frac{z_{k+1}-z}{z-z_k}(C_{k+1}-C_k)
\label{eq:heatcomposite}
\end{equation}

From the expressions of equivalent thermal conductivity and volumetric specific heat and the fact that thermal inertia $$I$$ is the square root of the product of $$\lambda$$ and $$C$$, one can also derive a relation between equivalent thermal inertia and those of the sub-layers:

\begin{equation}
\frac{I_{k+1/2}^2}{C_{k+1/2}} = \frac{z_{k+1}-z_k}{(z-z_k)\frac{C_k}{I_k^2} + (z_{k+1}-z )\frac{C_{k+1}}{I_{k+1}^2}}
\end{equation}

The PCM assumes a constant value for $$C_k = C$$ for each layer, set to $$C = 10^6 $$ J m$$^{-3}$$ K$$^{-1}$$. The equivalent thermal inertia is thus:
\begin{equation}
I_{k+1/2}= \sqrt{\frac{z_{k+1}-z_k}{\frac{z-z_k}{I_k^2} + \frac{z_{k+1}-z }{I_{k+1}^2}}}
\end{equation}

== 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. [https://theses.hal.science/tel-04824230 Link to the manuscript]

[[Category:Generic-Model]]

Soil Thermal Conduction in the Generic PCM

2025-03-05T12:43:57Z

LucasLange: /* Useful Tricks */

== 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.

To compute the subsurface temperatures, one must set ''callsoil'' to .true. in the ''callphys.def''

== 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 vertical grid discretization, and the default values:
<syntaxhighlight lang="fortran" line>
integer :: nsoilmx = 18 ! The number of points in the subsurface.
real :: lay1_soil = 2.e-4 ! The depth (m) of first full soil layer
real :: alpha_soil = 2 ! The geometric factor for the soi layer thickness
</syntaxhighlight>

The vertical discretization is defined in two arrays:
<syntaxhighlight lang="fortran" line>
real :: layer(nsoilmx) ! Full soil layer depths
real :: mlayer(0:nsoilmx-1) ! Mid soil layer depths
</syntaxhighlight>

The soil thermal properties are defined by:

<syntaxhighlight lang="fortran" line>
real :: volcapa ! soil volumetric heat capacity, C
real :: inertiedat(ngrid,nsoilmx) ! Soil thermal inertia
</syntaxhighlight>

==== soil_settings ====

This subroutine reads from a NetCDF file (opened by the caller) of "startfi.nc" format the soil properties (number of grid points, soil temperature, thermal inertia), and eventually process to interpolation if needed.

* 1. The number of soil layers is read in the NetCDF file. If the number of soil layer (''dimlen'') read is not equal to ''nsoilmx'', then the soil temperature and thermal inertia will be interpolated to the grid from the PCM defined in comsoil_h.
* 2. The vertical grid is read in the NetCDF file, or created if it does not exist, or if ''dimlen'' $$\neq$$ ''nsoilmx''.
* 3. The volumetric specific heat ''volcapa'' is read in the controle array (''controle(35)'') of the NetCDF file. If it does not exist, then ''volcapa'' is set to $$10^6$$.
* 4. The Thermal inertia is read in the NetCDF file. If it does not exist, the subroutines crashes.
* 5. The Soil Temperature is read in the NetCDF file. If it does not exist, the subroutines crashes.
* 6. Interpolation of soil temperature and thermal inertia is done if ''dimlen'' $$\neq$$ ''nsoilmx''.

'''Inputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
integer :: nid ! Name of the NetCDF file
integer :: ngrid ! Number of surface grid points.
integer :: nsoil ! Number of soil vertical grid points.
integer :: indextime ! position of the first time point in the netCDF file.
real :: tsurf(ngrid) ! Surface temperatures (used for the interpolation)
</syntaxhighlight>

'''Outputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
real :: tsoil(ngrid,nsoil) ! Soil Temperature
</syntaxhighlight>

Note that the soil thermal inertia are modified within the function but not are not output.

==== soil ====
This subroutine computes the soil temperatures from surface temperatures and heat fluxes from the ground. The code works as follow:
* 1. During the first call, local variables (soil thickness, soil conductivity at the interface, etc.) are defined.
* 2. The first grid-point soil temperature is computed from the surface temperature and ground heat flux.
* 3. Other soil temperatures are then computed based on a recurrence relation
* 4. Ground heat flux for next time step is computed.

'''Inputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
integer :: ngrid ! Number of surface grid points.
integer :: nsoil ! Number of soil vertical grid points.
logical :: firstcall ! Boolean checking if it is the first call to this subroutine.
logical :: lastcall ! Boolean checking if it is the last call to this subroutine (not used actually).
real :: tsurf(ngrid) ! Surface temperatures (used for the interpolation)
real :: therm_i(ngrid,nsoil) ! Soil thermal inertia
</syntaxhighlight>

'''Outputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
real :: tsoil(ngrid,nsoil) ! Soil Temperature
real :: fluxgrd(ngrid) ! Soil heat flux
real :: capcal(ngrid) ! surface specific heat capacity
</syntaxhighlight>

== Useful Tricks ==
==== Changing the vertical grid ====
To change the vertical grid, one can changes these parameters in the run.def:
<syntaxhighlight lang="bash">
## Number of grid points in the subsurface
nsoilmx = ...
## Thickness of topmost soil layer (m)
lay1_soil = ...
## Coefficient for soil layer thickness distribution
soil_thickness = ...
</syntaxhighlight>

A word of caution: the soil grid must be adapted to catch the diurnal/seasonal heat propagation in the ground. In the code (''soil.F''), the first layer must be shallower than 0.5 x the diurnal skin depth, and larger than 2 times the seasonal skin depth. The skin depth $$\delta_P$$ (m) is calculated with:

\begin{equation}
\delta_P = \sqrt{\frac{P \lambda}{\pi C}}
\label{eq:skindepth}
\end{equation}
where $$P$$ is the diurnal/seasonal period (in seconds).

==== Inhomogeneity of the subsurface ====

The scheme developed by E. Millour allows to have varying thermal properties with depth. When a layer of soil is not homogeneous, but instead composed of 2 separate sub-layers (e.g., dry regolith above massive subsurface ice), the layer's equivalent physical properties, i.e. values specified at node $z_{k+1/2}$, should be adjusted accordingly.

The equivalent thermal conductivity of a composite layer is derived from the steady-state composite wall analogy. The thermal conductivity $\lambda_{k+1/2}$ of a composite layer consisting of a layer (extending from depths $z_k$ to $z_{k+1}$) of thermal conductivity $\lambda_k$ above a layer (extending from depths $z$ to $z_{k+1}$) of thermal conductivity $\lambda_{k+1}$ is set to:

\begin{equation}
\lambda_{k+1/2} = \frac{z_{k+1}-z_k}{\frac{z-z_k}{\lambda_k} + \frac{z_{k+1}-z}{\lambda_{k+1}}}
\label{eq:conductivitycomposite}
\end{equation}

The equivalent volumetric specific heat $C_{k+1/2}$ of a mixed layer must be such that the overall ability of the volume to store internal energy is equal to the sum of the individual sub-layers. $C_{k+1/2}$ is thus simply obtained as:

\begin{equation}
C_{k+1/2} = C_k + \frac{z_{k+1}-z}{z-z_k}(C_{k+1}-C_k)
\label{eq:heatcomposite}
\end{equation}

From the expressions of equivalent thermal conductivity and volumetric specific heat and the fact that thermal inertia $I$ is the square root of the product of $\lambda$ and $C$, one can also derive a relation between equivalent thermal inertia and those of the sub-layers:

\begin{equation}
\frac{I_{k+1/2}^2}{C_{k+1/2}} = \frac{z_{k+1}-z_k}{(z-z_k)\frac{C_k}{I_k^2} + (z_{k+1}-z )\frac{C_{k+1}}{I_{k+1}^2}}
\end{equation}

The PCM assumes a constant value for $C_k = C$ for each layer, set to $C = 10^6 $~J~m$^{-3}$~K$^{-1}$. The equivalent thermal inertia is thus:
\begin{equation}
I_{k+1/2}= \sqrt{\frac{z_{k+1}-z_k}{\frac{z-z_k}{I_k^2} + \frac{z_{k+1}-z }{I_{k+1}^2}}}
\end{equation}

== 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. [https://theses.hal.science/tel-04824230 Link to the manuscript]

[[Category:Generic-Model]]

Soil Thermal Conduction in the Generic PCM

2025-03-05T12:42:36Z

LucasLange: /* Changing the vertical grid */

== 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.

To compute the subsurface temperatures, one must set ''callsoil'' to .true. in the ''callphys.def''

== 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 vertical grid discretization, and the default values:
<syntaxhighlight lang="fortran" line>
integer :: nsoilmx = 18 ! The number of points in the subsurface.
real :: lay1_soil = 2.e-4 ! The depth (m) of first full soil layer
real :: alpha_soil = 2 ! The geometric factor for the soi layer thickness
</syntaxhighlight>

The vertical discretization is defined in two arrays:
<syntaxhighlight lang="fortran" line>
real :: layer(nsoilmx) ! Full soil layer depths
real :: mlayer(0:nsoilmx-1) ! Mid soil layer depths
</syntaxhighlight>

The soil thermal properties are defined by:

<syntaxhighlight lang="fortran" line>
real :: volcapa ! soil volumetric heat capacity, C
real :: inertiedat(ngrid,nsoilmx) ! Soil thermal inertia
</syntaxhighlight>

==== soil_settings ====

This subroutine reads from a NetCDF file (opened by the caller) of "startfi.nc" format the soil properties (number of grid points, soil temperature, thermal inertia), and eventually process to interpolation if needed.

* 1. The number of soil layers is read in the NetCDF file. If the number of soil layer (''dimlen'') read is not equal to ''nsoilmx'', then the soil temperature and thermal inertia will be interpolated to the grid from the PCM defined in comsoil_h.
* 2. The vertical grid is read in the NetCDF file, or created if it does not exist, or if ''dimlen'' $$\neq$$ ''nsoilmx''.
* 3. The volumetric specific heat ''volcapa'' is read in the controle array (''controle(35)'') of the NetCDF file. If it does not exist, then ''volcapa'' is set to $$10^6$$.
* 4. The Thermal inertia is read in the NetCDF file. If it does not exist, the subroutines crashes.
* 5. The Soil Temperature is read in the NetCDF file. If it does not exist, the subroutines crashes.
* 6. Interpolation of soil temperature and thermal inertia is done if ''dimlen'' $$\neq$$ ''nsoilmx''.

'''Inputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
integer :: nid ! Name of the NetCDF file
integer :: ngrid ! Number of surface grid points.
integer :: nsoil ! Number of soil vertical grid points.
integer :: indextime ! position of the first time point in the netCDF file.
real :: tsurf(ngrid) ! Surface temperatures (used for the interpolation)
</syntaxhighlight>

'''Outputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
real :: tsoil(ngrid,nsoil) ! Soil Temperature
</syntaxhighlight>

Note that the soil thermal inertia are modified within the function but not are not output.

==== soil ====
This subroutine computes the soil temperatures from surface temperatures and heat fluxes from the ground. The code works as follow:
* 1. During the first call, local variables (soil thickness, soil conductivity at the interface, etc.) are defined.
* 2. The first grid-point soil temperature is computed from the surface temperature and ground heat flux.
* 3. Other soil temperatures are then computed based on a recurrence relation
* 4. Ground heat flux for next time step is computed.

'''Inputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
integer :: ngrid ! Number of surface grid points.
integer :: nsoil ! Number of soil vertical grid points.
logical :: firstcall ! Boolean checking if it is the first call to this subroutine.
logical :: lastcall ! Boolean checking if it is the last call to this subroutine (not used actually).
real :: tsurf(ngrid) ! Surface temperatures (used for the interpolation)
real :: therm_i(ngrid,nsoil) ! Soil thermal inertia
</syntaxhighlight>

'''Outputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
real :: tsoil(ngrid,nsoil) ! Soil Temperature
real :: fluxgrd(ngrid) ! Soil heat flux
real :: capcal(ngrid) ! surface specific heat capacity
</syntaxhighlight>

== Useful Tricks ==
==== Changing the vertical grid ====
To change the vertical grid, one can changes these parameters in the run.def:
<syntaxhighlight lang="bash">
## Number of grid points in the subsurface
nsoilmx = ...
## Thickness of topmost soil layer (m)
lay1_soil = ...
## Coefficient for soil layer thickness distribution
soil_thickness = ...
</syntaxhighlight>

A word of caution: the soil grid must be adapted to catch the diurnal/seasonal heat propagation in the ground. In the code (''soil.F''), the first layer must be shallower than 0.5 x the diurnal skin depth, and larger than 2 times the seasonal skin depth. The skin depth $$\delta_P$$ (m) is calculated with:

\begin{equation}
\delta_P = \sqrt{\frac{P \lambda}{\pi C}}
\label{eq:skindepth}
\end{equation}
where $$P$$ is the diurnal/seasonal period (in seconds).

== 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. [https://theses.hal.science/tel-04824230 Link to the manuscript]

[[Category:Generic-Model]]

Soil Thermal Conduction in the Generic PCM

2025-03-05T12:36:02Z

LucasLange: /* Changing the vertical grid */

== 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.

To compute the subsurface temperatures, one must set ''callsoil'' to .true. in the ''callphys.def''

== 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 vertical grid discretization, and the default values:
<syntaxhighlight lang="fortran" line>
integer :: nsoilmx = 18 ! The number of points in the subsurface.
real :: lay1_soil = 2.e-4 ! The depth (m) of first full soil layer
real :: alpha_soil = 2 ! The geometric factor for the soi layer thickness
</syntaxhighlight>

The vertical discretization is defined in two arrays:
<syntaxhighlight lang="fortran" line>
real :: layer(nsoilmx) ! Full soil layer depths
real :: mlayer(0:nsoilmx-1) ! Mid soil layer depths
</syntaxhighlight>

The soil thermal properties are defined by:

<syntaxhighlight lang="fortran" line>
real :: volcapa ! soil volumetric heat capacity, C
real :: inertiedat(ngrid,nsoilmx) ! Soil thermal inertia
</syntaxhighlight>

==== soil_settings ====

This subroutine reads from a NetCDF file (opened by the caller) of "startfi.nc" format the soil properties (number of grid points, soil temperature, thermal inertia), and eventually process to interpolation if needed.

* 1. The number of soil layers is read in the NetCDF file. If the number of soil layer (''dimlen'') read is not equal to ''nsoilmx'', then the soil temperature and thermal inertia will be interpolated to the grid from the PCM defined in comsoil_h.
* 2. The vertical grid is read in the NetCDF file, or created if it does not exist, or if ''dimlen'' $$\neq$$ ''nsoilmx''.
* 3. The volumetric specific heat ''volcapa'' is read in the controle array (''controle(35)'') of the NetCDF file. If it does not exist, then ''volcapa'' is set to $$10^6$$.
* 4. The Thermal inertia is read in the NetCDF file. If it does not exist, the subroutines crashes.
* 5. The Soil Temperature is read in the NetCDF file. If it does not exist, the subroutines crashes.
* 6. Interpolation of soil temperature and thermal inertia is done if ''dimlen'' $$\neq$$ ''nsoilmx''.

'''Inputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
integer :: nid ! Name of the NetCDF file
integer :: ngrid ! Number of surface grid points.
integer :: nsoil ! Number of soil vertical grid points.
integer :: indextime ! position of the first time point in the netCDF file.
real :: tsurf(ngrid) ! Surface temperatures (used for the interpolation)
</syntaxhighlight>

'''Outputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
real :: tsoil(ngrid,nsoil) ! Soil Temperature
</syntaxhighlight>

Note that the soil thermal inertia are modified within the function but not are not output.

==== soil ====
This subroutine computes the soil temperatures from surface temperatures and heat fluxes from the ground. The code works as follow:
* 1. During the first call, local variables (soil thickness, soil conductivity at the interface, etc.) are defined.
* 2. The first grid-point soil temperature is computed from the surface temperature and ground heat flux.
* 3. Other soil temperatures are then computed based on a recurrence relation
* 4. Ground heat flux for next time step is computed.

'''Inputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
integer :: ngrid ! Number of surface grid points.
integer :: nsoil ! Number of soil vertical grid points.
logical :: firstcall ! Boolean checking if it is the first call to this subroutine.
logical :: lastcall ! Boolean checking if it is the last call to this subroutine (not used actually).
real :: tsurf(ngrid) ! Surface temperatures (used for the interpolation)
real :: therm_i(ngrid,nsoil) ! Soil thermal inertia
</syntaxhighlight>

'''Outputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
real :: tsoil(ngrid,nsoil) ! Soil Temperature
real :: fluxgrd(ngrid) ! Soil heat flux
real :: capcal(ngrid) ! surface specific heat capacity
</syntaxhighlight>

== Useful Tricks ==
==== Changing the vertical grid ====
To change the vertical grid, one can changes these parameters in the run.def:
<syntaxhighlight lang="bash">
## Number of grid points in the subsurface
nsoilmx = ...
## Thickness of topmost soil layer (m)
lay1_soil = ...
## Coefficient for soil layer thickness distribution
soil_thickness = ...
</syntaxhighlight>

== 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. [https://theses.hal.science/tel-04824230 Link to the manuscript]

[[Category:Generic-Model]]

Soil Thermal Conduction in the Generic PCM

2025-03-05T12:35:45Z

LucasLange:

== 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.

To compute the subsurface temperatures, one must set ''callsoil'' to .true. in the ''callphys.def''

== 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 vertical grid discretization, and the default values:
<syntaxhighlight lang="fortran" line>
integer :: nsoilmx = 18 ! The number of points in the subsurface.
real :: lay1_soil = 2.e-4 ! The depth (m) of first full soil layer
real :: alpha_soil = 2 ! The geometric factor for the soi layer thickness
</syntaxhighlight>

The vertical discretization is defined in two arrays:
<syntaxhighlight lang="fortran" line>
real :: layer(nsoilmx) ! Full soil layer depths
real :: mlayer(0:nsoilmx-1) ! Mid soil layer depths
</syntaxhighlight>

The soil thermal properties are defined by:

<syntaxhighlight lang="fortran" line>
real :: volcapa ! soil volumetric heat capacity, C
real :: inertiedat(ngrid,nsoilmx) ! Soil thermal inertia
</syntaxhighlight>

==== soil_settings ====

This subroutine reads from a NetCDF file (opened by the caller) of "startfi.nc" format the soil properties (number of grid points, soil temperature, thermal inertia), and eventually process to interpolation if needed.

* 1. The number of soil layers is read in the NetCDF file. If the number of soil layer (''dimlen'') read is not equal to ''nsoilmx'', then the soil temperature and thermal inertia will be interpolated to the grid from the PCM defined in comsoil_h.
* 2. The vertical grid is read in the NetCDF file, or created if it does not exist, or if ''dimlen'' $$\neq$$ ''nsoilmx''.
* 3. The volumetric specific heat ''volcapa'' is read in the controle array (''controle(35)'') of the NetCDF file. If it does not exist, then ''volcapa'' is set to $$10^6$$.
* 4. The Thermal inertia is read in the NetCDF file. If it does not exist, the subroutines crashes.
* 5. The Soil Temperature is read in the NetCDF file. If it does not exist, the subroutines crashes.
* 6. Interpolation of soil temperature and thermal inertia is done if ''dimlen'' $$\neq$$ ''nsoilmx''.

'''Inputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
integer :: nid ! Name of the NetCDF file
integer :: ngrid ! Number of surface grid points.
integer :: nsoil ! Number of soil vertical grid points.
integer :: indextime ! position of the first time point in the netCDF file.
real :: tsurf(ngrid) ! Surface temperatures (used for the interpolation)
</syntaxhighlight>

'''Outputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
real :: tsoil(ngrid,nsoil) ! Soil Temperature
</syntaxhighlight>

Note that the soil thermal inertia are modified within the function but not are not output.

==== soil ====
This subroutine computes the soil temperatures from surface temperatures and heat fluxes from the ground. The code works as follow:
* 1. During the first call, local variables (soil thickness, soil conductivity at the interface, etc.) are defined.
* 2. The first grid-point soil temperature is computed from the surface temperature and ground heat flux.
* 3. Other soil temperatures are then computed based on a recurrence relation
* 4. Ground heat flux for next time step is computed.

'''Inputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
integer :: ngrid ! Number of surface grid points.
integer :: nsoil ! Number of soil vertical grid points.
logical :: firstcall ! Boolean checking if it is the first call to this subroutine.
logical :: lastcall ! Boolean checking if it is the last call to this subroutine (not used actually).
real :: tsurf(ngrid) ! Surface temperatures (used for the interpolation)
real :: therm_i(ngrid,nsoil) ! Soil thermal inertia
</syntaxhighlight>

'''Outputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
real :: tsoil(ngrid,nsoil) ! Soil Temperature
real :: fluxgrd(ngrid) ! Soil heat flux
real :: capcal(ngrid) ! surface specific heat capacity
</syntaxhighlight>

== Useful Tricks ==
==== Changing the vertical grid ====
To change the vertical grid, one can changes these parameters in the run.def:
## Number of grid points in the subsurface
nsoilmx = ...
## Thickness of topmost soil layer (m)
lay1_soil = ...
## Coefficient for soil layer thickness distribution
soil_thickness = ...
</syntaxhighlight>

== 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. [https://theses.hal.science/tel-04824230 Link to the manuscript]

[[Category:Generic-Model]]

Soil Thermal Conduction in the Generic PCM

2025-03-05T12:32:58Z

LucasLange:

== 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.

To compute the subsurface temperatures, one must set ''callsoil'' to .true. in the ''callphys.def''

== 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 vertical grid discretization, and the default values:
<syntaxhighlight lang="fortran" line>
integer :: nsoilmx = 18 ! The number of points in the subsurface.
real :: lay1_soil = 2.e-4 ! The depth (m) of first full soil layer
real :: alpha_soil = 2 ! The geometric factor for the soi layer thickness
</syntaxhighlight>

The vertical discretization is defined in two arrays:
<syntaxhighlight lang="fortran" line>
real :: layer(nsoilmx) ! Full soil layer depths
real :: mlayer(0:nsoilmx-1) ! Mid soil layer depths
</syntaxhighlight>

The soil thermal properties are defined by:

<syntaxhighlight lang="fortran" line>
real :: volcapa ! soil volumetric heat capacity, C
real :: inertiedat(ngrid,nsoilmx) ! Soil thermal inertia
</syntaxhighlight>

==== soil_settings ====

This subroutine reads from a NetCDF file (opened by the caller) of "startfi.nc" format the soil properties (number of grid points, soil temperature, thermal inertia), and eventually process to interpolation if needed.

* 1. The number of soil layers is read in the NetCDF file. If the number of soil layer (''dimlen'') read is not equal to ''nsoilmx'', then the soil temperature and thermal inertia will be interpolated to the grid from the PCM defined in comsoil_h.
* 2. The vertical grid is read in the NetCDF file, or created if it does not exist, or if ''dimlen'' $$\neq$$ ''nsoilmx''.
* 3. The volumetric specific heat ''volcapa'' is read in the controle array (''controle(35)'') of the NetCDF file. If it does not exist, then ''volcapa'' is set to $$10^6$$.
* 4. The Thermal inertia is read in the NetCDF file. If it does not exist, the subroutines crashes.
* 5. The Soil Temperature is read in the NetCDF file. If it does not exist, the subroutines crashes.
* 6. Interpolation of soil temperature and thermal inertia is done if ''dimlen'' $$\neq$$ ''nsoilmx''.

'''Inputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
integer :: nid ! Name of the NetCDF file
integer :: ngrid ! Number of surface grid points.
integer :: nsoil ! Number of soil vertical grid points.
integer :: indextime ! position of the first time point in the netCDF file.
real :: tsurf(ngrid) ! Surface temperatures (used for the interpolation)
</syntaxhighlight>

'''Outputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
real :: tsoil(ngrid,nsoil) ! Soil Temperature
</syntaxhighlight>

Note that the soil thermal inertia are modified within the function but not are not output.

==== soil ====
This subroutine computes the soil temperatures from surface temperatures and heat fluxes from the ground. The code works as follow:
* 1. During the first call, local variables (soil thickness, soil conductivity at the interface, etc.) are defined.
* 2. The first grid-point soil temperature is computed from the surface temperature and ground heat flux.
* 3. Other soil temperatures are then computed based on a recurrence relation
* 4. Ground heat flux for next time step is computed.

'''Inputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
integer :: ngrid ! Number of surface grid points.
integer :: nsoil ! Number of soil vertical grid points.
logical :: firstcall ! Boolean checking if it is the first call to this subroutine.
logical :: lastcall ! Boolean checking if it is the last call to this subroutine (not used actually).
real :: tsurf(ngrid) ! Surface temperatures (used for the interpolation)
real :: therm_i(ngrid,nsoil) ! Soil thermal inertia
</syntaxhighlight>

'''Outputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
real :: tsoil(ngrid,nsoil) ! Soil Temperature
real :: fluxgrd(ngrid) ! Soil heat flux
real :: capcal(ngrid) ! surface specific heat capacity
</syntaxhighlight>

== Useful Tricks ==

== 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. [https://theses.hal.science/tel-04824230 Link to the manuscript]

[[Category:Generic-Model]]

Soil Thermal Conduction in the Generic PCM

2025-03-05T11:16:52Z

LucasLange: /* Numerical approach */

== 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 vertical grid discretization, and the default values:
<syntaxhighlight lang="fortran" line>
integer :: nsoilmx = 18 ! The number of points in the subsurface.
real :: lay1_soil = 2.e-4 ! The depth (m) of first full soil layer
real :: alpha_soil = 2 ! The geometric factor for the soi layer thickness
</syntaxhighlight>

The vertical discretization is defined in two arrays:
<syntaxhighlight lang="fortran" line>
real :: layer(nsoilmx) ! Full soil layer depths
real :: mlayer(0:nsoilmx-1) ! Mid soil layer depths
</syntaxhighlight>

The soil thermal properties are defined by:

<syntaxhighlight lang="fortran" line>
real :: volcapa ! soil volumetric heat capacity, C
real :: inertiedat(ngrid,nsoilmx) ! Soil thermal inertia
</syntaxhighlight>

==== soil_settings ====

This subroutine reads from a NetCDF file (opened by the caller) of "startfi.nc" format the soil properties (number of grid points, soil temperature, thermal inertia), and eventually process to interpolation if needed.

* 1. The number of soil layers is read in the NetCDF file. If the number of soil layer (''dimlen'') read is not equal to ''nsoilmx'', then the soil temperature and thermal inertia will be interpolated to the grid from the PCM defined in comsoil_h.
* 2. The vertical grid is read in the NetCDF file, or created if it does not exist, or if ''dimlen'' $$\neq$$ ''nsoilmx''.
* 3. The volumetric specific heat ''volcapa'' is read in the controle array (''controle(35)'') of the NetCDF file. If it does not exist, then ''volcapa'' is set to $$10^6$$.
* 4. The Thermal inertia is read in the NetCDF file. If it does not exist, the subroutines crashes.
* 5. The Soil Temperature is read in the NetCDF file. If it does not exist, the subroutines crashes.
* 6. Interpolation of soil temperature and thermal inertia is done if ''dimlen'' $$\neq$$ ''nsoilmx''.

'''Inputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
integer :: nid ! Name of the NetCDF file
integer :: ngrid ! Number of surface grid points.
integer :: nsoil ! Number of soil vertical grid points.
integer :: indextime ! position of the first time point in the netCDF file.
real :: tsurf(ngrid) ! Surface temperatures (used for the interpolation)
</syntaxhighlight>

'''Outputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
real :: tsoil(ngrid,nsoil) ! Soil Temperature
</syntaxhighlight>

Note that the soil thermal inertia are modified within the function but not are not output.

==== soil ====
This subroutine computes the soil temperatures from surface temperatures and heat fluxes from the ground. The code works as follow:
* 1. During the first call, local variables (soil thickness, soil conductivity at the interface, etc.) are defined.
* 2. The first grid-point soil temperature is computed from the surface temperature and ground heat flux.
* 3. Other soil temperatures are then computed based on a recurrence relation
* 4. Ground heat flux for next time step is computed.

'''Inputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
integer :: ngrid ! Number of surface grid points.
integer :: nsoil ! Number of soil vertical grid points.
logical :: firstcall ! Boolean checking if it is the first call to this subroutine.
logical :: lastcall ! Boolean checking if it is the last call to this subroutine (not used actually).
real :: tsurf(ngrid) ! Surface temperatures (used for the interpolation)
real :: therm_i(ngrid,nsoil) ! Soil thermal inertia
</syntaxhighlight>

'''Outputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
real :: tsoil(ngrid,nsoil) ! Soil Temperature
real :: fluxgrd(ngrid) ! Soil heat flux
real :: capcal(ngrid) ! surface specific heat capacity
</syntaxhighlight>

== Useful Tricks ==

== 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. [https://theses.hal.science/tel-04824230 Link to the manuscript]

[[Category:Generic-Model]]

Soil Thermal Conduction in the Generic PCM

2025-03-05T11:10:11Z

LucasLange:

== 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 vertical grid discretization, and the default values:
<syntaxhighlight lang="fortran" line>
integer :: nsoilmx = 18 ! The number of points in the subsurface.
real :: lay1_soil = 2.e-4 ! The depth (m) of first full soil layer
real :: alpha_soil = 2 ! The geometric factor for the soi layer thickness
</syntaxhighlight>

The vertical discretization is defined in two arrays:
<syntaxhighlight lang="fortran" line>
real :: layer(nsoilmx) ! Full soil layer depths
real :: mlayer(0:nsoilmx-1) ! Mid soil layer depths
</syntaxhighlight>

The soil thermal properties are defined by:

<syntaxhighlight lang="fortran" line>
real :: volcapa ! soil volumetric heat capacity, C
real :: inertiedat(ngrid,nsoilmx) ! Soil thermal inertia
</syntaxhighlight>

==== soil_settings ====

This subroutine reads from a NetCDF file (opened by the caller) of "startfi.nc" format the soil properties (number of grid points, soil temperature, thermal inertia), and eventually process to interpolation if needed.

* 1. The number of soil layers is read in the NetCDF file. If the number of soil layer (''dimlen'') read is not equal to ''nsoilmx'', then the soil temperature and thermal inertia will be interpolated to the grid from the PCM defined in comsoil_h.
* 2. The vertical grid is read in the NetCDF file, or created if it does not exist, or if ''dimlen'' $$\neq$$ ''nsoilmx''.
* 3. The volumetric specific heat ''volcapa'' is read in the controle array (''controle(35)'') of the NetCDF file. If it does not exist, then ''volcapa'' is set to $$10^6$$.
* 4. The Thermal inertia is read in the NetCDF file. If it does not exist, the subroutines crashes.
* 5. The Soil Temperature is read in the NetCDF file. If it does not exist, the subroutines crashes.
* 6. Interpolation of soil temperature and thermal inertia is done if ''dimlen'' $$\neq$$ ''nsoilmx''.

'''Inputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
integer :: nid ! Name of the NetCDF file
integer :: ngrid ! Number of surface grid points.
integer :: nsoil ! Number of soil vertical grid points.
integer :: indextime ! position of the first time point in the netCDF file.
real :: tsurf(ngrid) ! Surface temperatures (used for the interpolation)
</syntaxhighlight>

'''Outputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
real :: tsoil(ngrid,nsoil) ! Soil Temperature
</syntaxhighlight>

Note that the soil thermal inertia are modified within the function but not are not output.

==== soil ====
This subroutine computes the soil temperatures from surface temperatures and heat fluxes from the ground. The code works as follow:
* 1. During the first call, local variables (soil thickness, soil conductivity at the interface, etc.) are defined.
* 2. The first grid-point soil temperature is computed from the surface temperature and ground heat flux.
* 3. Other soil temperatures are then computed based on a recurrence relation
* 4. Ground heat flux for next time step is computed.

'''Inputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
integer :: ngrid ! Number of surface grid points.
integer :: nsoil ! Number of soil vertical grid points.
logical :: firstcall ! Boolean checking if it is the first call to this subroutine.
logical :: lastcall ! Boolean checking if it is the last call to this subroutine (not used actually).
real :: tsurf(ngrid) ! Surface temperatures (used for the interpolation)
real :: therm_i(ngrid,nsoil) ! Soil thermal inertia
</syntaxhighlight>

'''Outputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
real :: tsoil(ngrid,nsoil) ! Soil Temperature
real :: fluxgrd(ngrid) ! Soil heat flux
real :: capcal(ngrid) ! surface specific heat capacity
</syntaxhighlight>

== Useful Tricks ==

== 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. [https://theses.hal.science/tel-04824230 Link to the manuscript]

[[Category:Generic-Model]]

Soil Thermal Conduction in the Generic PCM

2025-03-05T11:09:31Z

LucasLange: /* soil_settings */

== 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 vertical grid discretization, and the default values:
<syntaxhighlight lang="fortran" line>
integer :: nsoilmx = 18 ! The number of points in the subsurface.
real :: lay1_soil = 2.e-4 ! The depth (m) of first full soil layer
real :: alpha_soil = 2 ! The geometric factor for the soi layer thickness
</syntaxhighlight>

The vertical discretization is defined in two arrays:
<syntaxhighlight lang="fortran" line>
real :: layer(nsoilmx) ! Full soil layer depths
real :: mlayer(0:nsoilmx-1) ! Mid soil layer depths
</syntaxhighlight>

The soil thermal properties are defined by:

<syntaxhighlight lang="fortran" line>
real :: volcapa ! soil volumetric heat capacity, C
real :: inertiedat(ngrid,nsoilmx) ! Soil thermal inertia
</syntaxhighlight>

==== soil_settings ====

This subroutine reads from a NetCDF file (opened by the caller) of "startfi.nc" format the soil properties (number of grid points, soil temperature, thermal inertia), and eventually process to interpolation if needed.

* 1. The number of soil layers is read in the NetCDF file. If the number of soil layer (''dimlen'') read is not equal to ''nsoilmx'', then the soil temperature and thermal inertia will be interpolated to the grid from the PCM defined in comsoil_h.
* 2. The vertical grid is read in the NetCDF file, or created if it does not exist, or if ''dimlen'' $$\neq$$ ''nsoilmx''.
* 3. The volumetric specific heat ''volcapa'' is read in the controle array (''controle(35)'') of the NetCDF file. If it does not exist, then ''volcapa'' is set to $$10^6$$.
* 4. The Thermal inertia is read in the NetCDF file. If it does not exist, the subroutines crashes.
* 5. The Soil Temperature is read in the NetCDF file. If it does not exist, the subroutines crashes.
* 6. Interpolation of soil temperature and thermal inertia is done if ''dimlen'' $$\neq$$ ''nsoilmx''.

'''Inputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
integer :: nid ! Name of the NetCDF file
integer :: ngrid ! Number of surface grid points.
integer :: nsoil ! Number of soil vertical grid points.
integer :: indextime ! position of the first time point in the netCDF file.
real :: tsurf(ngrid) ! Surface temperatures (used for the interpolation)
</syntaxhighlight>

'''Outputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
real :: tsoil(ngrid,nsoil) ! Soil Temperature
</syntaxhighlight>

Note that the soil thermal inertia are modified within the function but not are not output.

==== soil ====
This subroutine computes the soil temperatures from surface temperatures and heat fluxes from the ground. The code works as follow:
* 1. During the first call, local variables (soil thickness, soil conductivity at the interface, etc.) are defined.
* 2. The first grid-point soil temperature is computed from the surface temperature and ground heat flux.
* 3. Other soil temperatures are then computed based on a recurrence relation
* 4. Ground heat flux for next time step is computed.

'''Inputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
integer :: ngrid ! Number of surface grid points.
integer :: nsoil ! Number of soil vertical grid points.
logical :: firstcall ! Boolean checking if it is the first call to this subroutine.
logical :: lastcall ! Boolean checking if it is the last call to this subroutine (not used actually).
real :: tsurf(ngrid) ! Surface temperatures (used for the interpolation)
real :: therm_i(ngrid,nsoil) ! Soil thermal inertia
</syntaxhighlight>

'''Outputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
real :: tsoil(ngrid,nsoil) ! Soil Temperature
real :: fluxgrd(ngrid) ! Soil heat flux
real :: capcal(ngrid) ! surface specific heat capacity
</syntaxhighlight>

== 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. [https://theses.hal.science/tel-04824230 Link to the manuscript]

[[Category:Generic-Model]]

Soil Thermal Conduction in the Generic PCM

2025-03-05T11:09:14Z

LucasLange: /* soil */

== 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 vertical grid discretization, and the default values:
<syntaxhighlight lang="fortran" line>
integer :: nsoilmx = 18 ! The number of points in the subsurface.
real :: lay1_soil = 2.e-4 ! The depth (m) of first full soil layer
real :: alpha_soil = 2 ! The geometric factor for the soi layer thickness
</syntaxhighlight>

The vertical discretization is defined in two arrays:
<syntaxhighlight lang="fortran" line>
real :: layer(nsoilmx) ! Full soil layer depths
real :: mlayer(0:nsoilmx-1) ! Mid soil layer depths
</syntaxhighlight>

The soil thermal properties are defined by:

<syntaxhighlight lang="fortran" line>
real :: volcapa ! soil volumetric heat capacity, C
real :: inertiedat(ngrid,nsoilmx) ! Soil thermal inertia
</syntaxhighlight>

==== soil_settings ====

This subroutine reads from a NetCDF file (opened by the caller) of "startfi.nc" format the soil properties (number of grid points, soil temperature, thermal inertia), and eventually process to interpolation if needed.

* 1. The number of soil layers is read in the NetCDF file. If the number of soil layer (''dimlen'') read is not equal to ''nsoilmx'', then the soil temperature and thermal inertia will be interpolated to the grid from the PCM defined in comsoil_h.
* 2. The vertical grid is read in the NetCDF file, or created if it does not exist, or if ''dimlen'' $$\neq$$ ''nsoilmx''.
* 3. The volumetric specific heat ''volcapa'' is read in the controle array (''controle(35)'') of the NetCDF file. If it does not exist, then ''volcapa'' is set to $$10^6$$.
* 4. The Thermal inertia is read in the NetCDF file. If it does not exist, the subroutines crashes.
* 5. The Soil Temperature is read in the NetCDF file. If it does not exist, the subroutines crashes.
* 6. Interpolation of soil temperature and thermal inertia is done if ''dimlen'' $$\neq$$ ''nsoilmx''.

'''Inputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
integer :: nid ! Name of the NetCDF file
integer :: ngrid ! Number of surface grid points.
integer :: nsoil ! Number of soil vertical grid points.
integer :: indextime ! position of the first time point in the netCDF file.
real :: tsurf(ngrid) ! Surface temperatures (used for the interpolation)
</syntaxhighlight>

'''Outputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
integer :: tsoil(ngrid,nsoil) ! Soil Temperature
</syntaxhighlight>

Note that the soil thermal inertia are modified within the function but not are not output.

==== soil ====
This subroutine computes the soil temperatures from surface temperatures and heat fluxes from the ground. The code works as follow:
* 1. During the first call, local variables (soil thickness, soil conductivity at the interface, etc.) are defined.
* 2. The first grid-point soil temperature is computed from the surface temperature and ground heat flux.
* 3. Other soil temperatures are then computed based on a recurrence relation
* 4. Ground heat flux for next time step is computed.

'''Inputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
integer :: ngrid ! Number of surface grid points.
integer :: nsoil ! Number of soil vertical grid points.
logical :: firstcall ! Boolean checking if it is the first call to this subroutine.
logical :: lastcall ! Boolean checking if it is the last call to this subroutine (not used actually).
real :: tsurf(ngrid) ! Surface temperatures (used for the interpolation)
real :: therm_i(ngrid,nsoil) ! Soil thermal inertia
</syntaxhighlight>

'''Outputs of the subroutine:'''
<syntaxhighlight lang="fortran" line>
real :: tsoil(ngrid,nsoil) ! Soil Temperature
real :: fluxgrd(ngrid) ! Soil heat flux
real :: capcal(ngrid) ! surface specific heat capacity
</syntaxhighlight>

== 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. [https://theses.hal.science/tel-04824230 Link to the manuscript]

[[Category:Generic-Model]]

Soil Thermal Conduction in the Generic PCM

2025-03-05T11:03:44Z

LucasLange: /* soil */

Soil Thermal Conduction in the Generic PCM

2025-03-05T10:50:11Z

LucasLange: /* Description of the soil subroutines */

Soil Thermal Conduction in the Generic PCM

2025-03-05T10:48:48Z

LucasLange: /* soil_settings */

Soil Thermal Conduction in the Generic PCM

2025-03-05T10:48:27Z

LucasLange: /* soil_settings */

Soil Thermal Conduction in the Generic PCM

2025-03-05T10:38:04Z

LucasLange: /* soil_settings */

Soil Thermal Conduction in the Generic PCM

2025-03-05T10:32:22Z

LucasLange: /* soil_settings */

Soil Thermal Conduction in the Generic PCM

2025-03-05T10:31:41Z

LucasLange: /* soil_settings */

Soil Thermal Conduction in the Generic PCM

2025-03-05T10:26:11Z

LucasLange: /* Description of the soil subroutines */

Soil Thermal Conduction in the Generic PCM

2025-03-05T10:24:21Z

LucasLange: /* Description of the soil subroutines */

Soil Thermal Conduction in the Generic PCM

2025-03-05T10:24:04Z

LucasLange: /* comsoil_h */

Soil Thermal Conduction in the Generic PCM

2025-03-05T10:23:40Z

LucasLange: /* Description of the soil subroutines */

Soil Thermal Conduction in the Generic PCM

2025-03-05T10:20:08Z

LucasLange: /* comsoil_h */

Soil Thermal Conduction in the Generic PCM

2025-03-05T10:19:29Z

LucasLange: /* Description of the soil subroutines */

Soil Thermal Conduction in the Generic PCM

2025-03-05T10:16:38Z

LucasLange: /* comsoil_h */

Soil Thermal Conduction in the Generic PCM

2025-03-05T10:15:51Z

LucasLange: /* Purpose */

Soil Thermal Conduction in the Generic PCM

2025-03-05T10:15:34Z

LucasLange: /* comsoil_h */

Soil Thermal Conduction in the Generic PCM

2025-03-05T10:15:05Z

LucasLange:

Soil Thermal Conduction in the Generic PCM

2025-03-05T10:14:28Z

LucasLange: /* Description of the soil subroutines */

Soil Thermal Conduction in the Generic PCM

2025-03-05T10:14:11Z

LucasLange:

Soil Thermal Conduction in the Generic PCM

2025-03-05T10:13:45Z

LucasLange: /* Numerical approach */

Soil Thermal Conduction in the Generic PCM

2025-03-05T10:12:55Z

LucasLange: /* Numerical approach */

Soil Thermal Conduction in the Generic PCM

2025-03-05T10:12:32Z

LucasLange: /* Numerical approach */

Soil Thermal Conduction in the Generic PCM

2025-03-05T10:10:54Z

LucasLange: /* References */

Soil Thermal Conduction in the Generic PCM

2025-03-05T10:08:59Z

LucasLange: /* Purpose */

Soil Thermal Conduction in the Generic PCM

2025-03-05T10:07:49Z

LucasLange: /* Purpose */

Soil Thermal Conduction in the Generic PCM

2025-03-05T10:02:09Z

LucasLange: /* Numerical approach */

== 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.

== 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 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 seperatly. The essential feature of the artificial uncoupling is that, when solving the heat diffusion equation in the soil at time t(i+1), the value of surface temperature T (i+1)S is already known. In short, hte algorithm works as follows: the surface temperature at i + 1 is computed with a radiative balance, knowing F (i+1)S . Then soil temperatures and coefficients β are updated. F (i+2) S is finally computed and used as an input to compute T (i+2

== 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.

[[Category:Generic-Model]]

Soil Thermal Conduction in the Generic PCM

2025-03-05T09:59:52Z

LucasLange: /* Theoretical Considerations */

== 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.

== 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 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 seperatly. The essential feature of
the artificial uncoupling is that, when solving the heat diffusion equation in the soil at time t(i+1), the value of surface temperature T (i+1)S is already known. In short, hte algorithm works as follows: the surface temperature at i + 1 is computed with a radiative balance, knowing F (i+1)S . Then soil temperatures and coefficients β are updated. F (i+2) S is finally computed and used as an input to compute T (i+2

== 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.

[[Category:Generic-Model]]

Soil Thermal Conduction in the Generic PCM

2025-03-05T09:57:03Z

LucasLange: /* Numerical approach */

== 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.

== 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 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 seperatly. The essential feature of
the artificial uncoupling is that, when solving the heat diffusion equation in the soil at time t(i+1), the value of surface temperature T (i+1)S is already known. In short, hte algorithm works as follows: the surface temperature at i + 1 is computed with a radiative balance, knowing F (i+1)S . Then soil temperatures and coefficients β are updated. F (i+2) S is finally computed and used as an input to compute T (i+2

== 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.

[[Category:Generic-Model]]

Soil Thermal Conduction in the Generic PCM

2025-03-05T09:55:34Z

LucasLange: /* Numerical approach */

== 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.

== 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 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 seper-
atly, as will be described in the next sections. The essential feature of
the artificial uncoupling is that, when solving the heat diffusion equa-
tion in the soil at time t(i+1), the value of surface temperature T (i+1)
S
is already known.

== 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.

[[Category:Generic-Model]]

Soil Thermal Conduction in the Generic PCM

2025-03-05T09:54:49Z

LucasLange: /* References */

== 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.

== 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}

== 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.

[[Category:Generic-Model]]

Soil Thermal Conduction in the Generic PCM

2025-03-05T09:51:09Z

LucasLange: /* References */

Soil Thermal Conduction in the Generic PCM

2025-03-05T09:50:37Z

LucasLange: /* Numerical approach */

== 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.

== 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}

== 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.

== References ==

[[Category:Generic-Model]]

Soil Thermal Conduction in the Generic PCM

2025-03-05T09:37:07Z

LucasLange: /* Theoretical Considerations */

Soil Thermal Conduction in the Generic PCM

2025-03-05T08:55:36Z

LucasLange: /* Theoretical Considerations */

Soil Thermal Conduction in the Generic PCM

2025-03-05T08:53:54Z

LucasLange: Created page with "== 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. == Theoret..."

Subgrid Slopes in the Mars PCM

2023-10-19T09:04:25Z

LucasLange:

== Purpose ==
In order to study surface processes in the Mars Planetary Climate Model (e.g., gully formation, glacier evolution, ...), a representation of the sub-grid slopes has been introduced in the Mars PCM. This parametrization is detailed in Lange et al., 2023, "Modelling Slope Microclimates in the Mars Planetary Climate Model" https://arxiv.org/pdf/2306.12449.pdf. Contact: lucas.lange@lmd.ipsl.fr

For each GCM mesh, we decompose the cell as a distribution of sloped terrains (defined by characteristic slopes) and a flat terrain. On each sub-grid terrain, we let the microclimates evolve so that the slope-specific features (e.g., condensation of volatiles, formation of glaciers, migration of subsurface ice, etc.) can be simulated. The portion of the atmosphere above the ground within the cell sees an average of these surface microclimates, and all sub-grid terrains see the same ”shared atmosphere”.

== Practical Implementation ==
The practical implementation is detailed in section 3. of Lange et al., 2023. In short:

* The mesh can be divided in 1, 5, or 7 sub-grid slopes that are North - South facing. Then, using MOLA data (resolution: 64 pixel per degree), the distribution of sub-slope within the mesh is computed.

* All the surface quantities have now a dimension nslope: qsurf, watercap, tsurf, tsoil, albedo, emis, inertiesoil

* The radiative transfer is computed with grid-box albedo, emissivity, temperature. Then, the values are adapted for-each sub-grid slopes following Spiga & Forget 2008 (section 2. of ''physiq_mod'', ''param_slope'').

* The turbulent exchanges (''vdifc_mod'') are solved with the predominant sub-grid surface (i.e., the sub-grid surface which has the largest cover within the grid box).

* Water and CO2 ice condensation are computed for each sub-grid slopes (''vdifc_mod'', ''co2_condens'')

* Grid box averaged values are then computed using the specific computation detailed in the dedicated section.

== How to run with sub-grid slopes ? ==

By default, the model runs with only 1 sub-slope (i.e., a flat surface). If one wants to use this sub-grid parameterization, then one must use ''newstart'' and choose the option to add sub-slopes. FYI: The CPU time is increased by ~20% with 7 sub-slopes (64x48x26 grid).

== Post Processing ==

In order to go from the sub-grid surface quantity to the grid-box value, the following computations are done:
* For albedo and emissivity:
\begin{align}
X = \sum_{i=1}^7 X_i \delta_i
\label{eq:opticalproperties}
\end{align}
where $$X$$ refers to the grid parameter, $$X_i$$ the sub-grid parameter of slope $$i$$, $$\delta_i$$ the cover fraction (''subslope_dist'' in the code).

* For surface tracers, to ensure mass conservation:
\begin{align}
X = \sum_{i=1}^7 X_i \frac{\delta_i}{\cos(\mu_i)}
\label{eq:traceurs}
\end{align}
where $$\mu_i$$ is the slope angle of the sub-grid slope (''def_slope_mean'' in the code)

* For surface temperature, we average the surface flux:
\begin{align}
\epsilon T_{\rm{surf,grid}} = \sum_{i=1}^7 \epsilon_i T_{\rm{surf,i}}^4 \delta_i
\label{eq:Tsurf}
\end{align}
where $$\epsilon$$ (1) is the grid emissivity, $$T_{\rm{surf,grid}}$$ (K) the grid surface temperature

* This new parameterization enable, for each grid point, to determine the surface fields for the seven types of slopes modeled here. The question is now: for any point on the Martian surface, knowing the characteristics of the terrain (slope and azimuth), can we determine the values of the surface quantities using the GCM outputs that are discretized in space and slopes (i.e., build high resolution map of surface quantities)?

I (LL) propose to use the following method, inspired from the MCD.
Let us consider a point (noted $M$) of the surface of Mars, placed on a surface of slope angle $$\theta$$ and azimuth $$\psi$$. We compute $\mu$ and find the index $$i,i+1$$ of the two sub-grid slopes such that:
\begin{align}
\mu_i \leq \mu < \mu_{i+1}
\end{align}

The point $$M$$ is surrounded by four GCM meshes. We can interpolate the values obtained at the GCM points for the sub-grid slope $$\mu_i$$ and $$\mu_{i+1}$$ at point $$M$$. This is performed by using a bilinear interpolation in space. After that, the value for slope $$\mu$$ at $$M$$ is obtained with linear interpolation between the fields obtained after the bilinear interpolation for $$\mu_i$$ and $$\mu_i+1$$. If $$\mu < \mu_0$$ or $$\mu > \mu_7$$, an extrapolation is computed. Thus, we have virtually access to any surface for any kind of slopped surface on Mars using this complete climate model. WARNING: 1) we assume linear interpolation on processes that are not necessarily linear (although we interpolate neighbouring points); 2) we neglect the time-shift (generally lower than 1-2~hours) for slopes oriented East-West that are not well represented by this parameterization.

I (LL) have a python algorithm to do that, feel free to ask me!

[[Category:Mars-Model]]

Mars PCM physics: CO2 cycle

2023-10-19T09:03:41Z

LucasLange: /* Surface properties */

The goal of this page is to summarize key aspects, but theoretical and practical, of the parametrizations in the Mars physics package that enable modeling the Martian CO2 cycle

== The key processes ==
* CO2 condensation: atmospheric condensation. surface condensation
* polar sub-surface ice tables as high thermal inertia "heat reservoirs"
* surface properties (albedo of the surface)

== Parametrizations and key parameters in the model ==

=== CO2 condensation in the atmosphere ===
* CO2 ice clouds. These are implemented, but not used in standard simulations.
* CO2 condensation. See routine co2condens.

=== Polar sub-surface ice tables of high thermal inertia ===
These are set via the newstart utility. The information itself is stored in the startfi file. For now, the subsurface ice depth and thermal inertia values are set to fit the CO2 cycle thought an inversion of MONS spectrometer data (see Diez et al., 2008) by the routine ''xvik''. One can also use the Planetary Evolution Model to compute the position of the ice table.

=== Surface properties ===
In callphys.def one may specify ( by setting flag ''TESicealbedo'' to ''.true.'' ; recommended) that a climatology based on TES observations be used to drive (nudge) the evolution of CO2 polar caps.
<syntaxhighlight lang="bash">
# Impose polar cap surface albedos as observed by TES?
TESicealbedo = .true.
## Coefficient for Northern cap albedoes
TESice_Ncoef=1.6
## Coefficient for Southern cap albedoes
TESice_Scoef=1.6
</syntaxhighlight>
The sub-options ''TESice_Ncoef'' and ''TESice_Scoef'' are multiplying factors applied to the TES values, meant to account for non Lambertian effects.
In practice this is coded in ...

If a alternative setup with
<syntaxhighlight lang="bash">
# Impose polar cap surface albedos as observed by TES?
TESicealbedo = .false.
</syntaxhighlight>
Then the surface albedo is based on the presence or absence of CO2 ice on the surface. Without CO2 ice the albedo is that of bare ground and if CO2 ice is present it is set to 0.6 (variable ''albedice'' in surfdat_h.F90 which stores 2 values one for north and south). In this case, the perennial CO2 ice at the South Pole is represented throught a different albedo set in the ''callphys.def'':
<syntaxhighlight lang="bash">
## Albedo of the perenial CO2 ice at the South Pole
albedo_perenialco2 = 0.85
</syntaxhighlight>
[[Category:Mars-Model]]

Mars PCM physics: CO2 cycle

2023-10-19T08:58:39Z

LucasLange: /* Polar sub-surface ice tables of high thermal inertia */

Subgrid Slopes in the Mars PCM

2023-10-19T08:57:09Z

LucasLange: /* Purpose */