<figure><img src="image_1.png" /></figure>

HAL
open science

# Permafrost modelling with OpenFOAM®: New advancements of the permaFoam solver

L. Orgogozo, T. Xavier, H. Oulbani, C. Grenier

► To cite this version:

L. Orgogozo, T. Xavier, H. Oulbani, C. Grenier. Permafrost modelling with OpenFOAM®: New advancements of the permaFoam solver. Computer Physics Communications, 2023, 282, pp.108541. 10.1016/j.cpc.2022.108541. hal-03787933

HAL Id: hal-03787933

https://hal.science/hal-03787933

Submitted on 27 Sep 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Permafrost modelling with OpenFOAM®: New advancements of the permaFoam solver ⭐, ⭐☆

<figure><img src="image_3.png" /></figure>

L. Orgogozo a, *, T. Xavier a, H. Oulbani b, C. Grenier b

<figure><img src="image_4.png" /></figure>

## ARTICLE INFO

**Article history:**

Received 7 September 2021
Received in revised form 25 August 2022
Accepted 6 September 2022
Available online 14 September 2022

**Keywords:**

Permafrost
OpenFOAM
Cryohydrogeological modelling
Freeze/thaw
Evapotranspiration
High performance computing

## ABSTRACT

Permafrost, i.e., soil that is year-round frozen in depth, is covering a quarter of the continents of the northern hemisphere. It currently experiences fast changes due to climate change at global scale and technogenic perturbations at local scale, and the assessment and anticipation of these changes are of primary importance for many environmental and engineering applications in cold regions. To these ends, permafrost modelling is required, while the strong couplings and non-linearities involved in the physics at stake make it highly challenging, especially from a computational point of view. In this work we present a new solver for permafrost hydrology developed in the framework of OpenFOAM®, allowing to benefit from its advanced high-performance computing capabilities. The solver is tested for realistic, field-based cases, and its parallel performances are characterized up to ~16 000 cores on IRENE supercomputer (TGCC, CEA).

## Program summary

**Program Title:** permaFoam
**CPC Library link to program files:** https://doi.org/10.17632/swp88cvpwb.1
**Developer's repository link:** https://develop.openfoam.com/Community/hydrology/
**Licensing provisions:** GPLv3
**Programming language:** C++

**Nature of problem:** This software solves the coupled equations that govern water flow and heat transfer in variably saturated and variably frozen porous media, for transient problems in Three-dimensional, heterogeneous domains. The equation for water flow is Richards equation, which is a very popular model for water transfer in variably saturated porous media (e.g.: soils), and the equation for heat transfer is a Fourier equation including advection and the freeze/thaw of the pore water. The solver is designed to take advantage of the massively parallel computing performance of OpenFOAM®. The goal is to be able to model natural hydrosystems of cold regions on large temporal and spatial scales.

**Solution method:** For each equation a mixed implicit (FVM for Finite Volume Method in the object oriented OpenFOAM framework) and explicit (FVC for Finite Volume Calculus in the object oriented OpenFOAM® framework) discretization with backward time scheme is embedded in an iterative linearization procedure (Picard algorithm). The coupling between the two equations is performed through an operator splitting approach. The implementation has been carried out with a concern for robustness and parallel efficiency.

**Additional comments including restrictions and unusual features:** This version of permaFoam has been tested with OpenFOAM_v1912, v2106, v2112 and v2206, thus everything might not work with other (especially older) versions of OpenFOAM®. When using permaFoam, one should be careful to use fine enough spatial and temporal discretisations where and when steep fronts (freeze/thaw fronts, imbibition/drainage fronts) occur, otherwise numerical stability problems might arise.

© 2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

<footnote>a GET (Géosciences Environnement Toulouse), UMR 5563 CNRS/UR 234 IRD/UPS, Observatoire Midi Pyrénées, Université de Toulouse, Toulouse, France</footnote>

<footnote>b Laboratoire des Sciences du Climat et de l'Environnement, Université Paris Saclay, IPSL/LSCE, UMR 8212 CNRS CEA UVSQ, Orme des Merisiers, Gif sur Yvette Cedex, France</footnote>

<footnote>☆ The review of this paper was arranged by Prof. Hazel Andrew.</footnote>

<footnote>★★ This paper and its associated computer program are available via the Computer Physics Communications homepage on ScienceDirect (http://www.sciencedirect.com/science/journal/00104655).</footnote>

<footnote>* Corresponding author.
E-mail address: laurent.orgogozo@get.omp.eu (L. Orgogozo).</footnote>
# 1. Introduction

The permafrost areas of the northern hemisphere are currently warming twice time faster than the global average [17]. Currently, fast and important changes of both hydrological [58] and thermal [34] states of the northern continental surfaces are observed in response to permafrost thaw. In that context arctic and sub-arctic areas are prone to major biogeochemical and ecological transformations due to permafrost thaw, with strong associated feedbacks on greenhouse gas cycling (e.g.: [25,59]). Fast-increasing economic activity in cold regions is also impacted by climate change [14,51,46,2,23,50,24], and may generate itself important impacts on permafrost at local scale [57,47,33]. Overall the evolution of permafrost in changing environmental conditions is of primary importance [38]. Thus, the development of advanced numerical modelling tools of permafrost dynamics is an active research field (e.g.: [37,9,26,5,20,32,60,63,67]), in order to help predicting the impact of permafrost thaw on arctic thermo-hydrologic functioning and providing mechanistic understanding of Arctic environmental changes under climate warming. This knowledge is necessary to further understand carbon cycling and contaminant/nutrient transport, to assess risks and opportunities for water resources, sustainable urbanization, agriculture and general sustainable development of the (sub-)Arctic.

Since permafrost is year-round frozen soil, complemented by an active layer close to the surface freezing and thawing seasonally, dealing with its dynamics implies considering both water and heat transfers within soils. However, these two phenomena are non-linear and strongly coupled, which make their numerical resolution highly challenging. For instance, the soil apparent heat transfer properties, such as its apparent thermal conductivity or heat capacity, are impacted by its water content. On the other hand, the freezing of the pore water also induces a dramatic decrease in the apparent water conductivity. These couplings and non-linearities prevents the complex set of equations describing the problem from being solved analytically without strong hypothesis [31]. For the same reasons, the numerical approach represents a huge challenge, since it requires both a good stability of the solvers and a large amount of resources to perform the computation [58]. That is the reason why the need of high-performance computing in the field of cryohydrogeological simulation is widely acknowledged [43], and several permafrost modelling tools with High Performance Computing capabilities have been recently developed (e.g.: [26,45]).

Our contribution to the development of permafrost modelling is an OpenFOAM® cryohydrogeological solver, named "permaFoam". The permaFoam solver presented here is based on a previous version published as electronic supporting information in Orgogozo et al. [41] and on the RichardsFoam3 solver [39,40,42]. By solving the coupled equations of heat transfer and water flow in variably saturated porous media with freeze/thaw of the pore water and evaporation, it enables to model permafrost subsurface dynamics. Its development in the framework of OpenFOAM® allows to use efficiently modern supercomputing infrastructures (e.g.: [48,11]), and to benefit from the momentum of a large community including both academic and industrial users.

# 2. Physical equations

The physical problem at stake is the coupled transfer of water and heat in a variably saturated porous medium, taking into account the freezing of pore water and the evapotranspiration water uptake. As a cryohydrogeological solver, permaFoam only deals

with sub-surface processes, modelling permafrost itself. Nonetheless permafrost dynamics are strongly coupled with surface processes such as snowpack dynamics (e.g.: [21]) or evapotranspiration (e.g.: [41]), which control surface energy and mass balance and thus heat and water fluxes through the top boundary of the permafrost. Thus in many cases, in order to build proper top boundary conditions when using permaFoam, field measurements, assumptions on surface processes or even coupling with adequate surface processes models will be needed. Since permaFoam, as a cryohydrogeological solver, simulates the heat and water transfers within permafrost, this paper focuses on the sub-surface processes. As said above the coupling with surface processes models such as for instance Surface Energy Balance models (e.g.: [64,52,66]) may be necessary in some cases for building top boundary conditions for permaFoam. One can refer to Endrizzi et al. [9] or to Painter et al. [45] for examples of permafrost modelling tool with integrated coupling of surface and subsurface processes. Thus permaFoam focuses on modelling permafrost at Darcy scale, dealing with transfers in porous media that may include four different phases: the solid phase, the liquid water phase, the frozen water (ice) phase, and the air phase. The soil is considered as an undecformable porous media. Moreover, cryosuction and ice volumetric expansion are neglected. One can refer to Sturop et al. [53], for a modelling work about cryosuction, and to McKenzie et al. [37] and Westermann et al. [61] for examples of permafrost modelling works in which ice volumetric expansion is neglected. Due to these different approximations, cryo-mechanical processes are not considered. One can refer to Huang and Rudolph [18] and to Huang et al. [19] for modelling works dealing with coupled thermo-hydro-mechanical phenomena in variably saturated porous media with freeze/thaw. The mathematical description adopted in the present work has been widely used in the literature (e.g.: [30,13]), and recently discussed in Orgogozo et al. [41]. It involves basic coupling between water flow and heat flow, in the sense of Yu et al. [65]. The two primary variables that describe the physical state of water into the soil are the generalized water pressure $h$ [m] (positive in saturated conditions, negative in unsaturated conditions – see [49]) and the soil temperature $T$ [K] (defined by assuming local thermal equilibrium – see [15] for a discussion of this assumption). The two equations governing these primary variables are a modified Richards equation with a source term accounting for actual evapotranspiration (Eq. (1)) for the generalized water pressure $h$ and a macro-scale heat transfer equation for the porous medium with a term of latent heat exchange (Eq. (2)) for the soil temperature $T$.

$$C_H(h) \frac{\partial h}{\partial t} = \nabla \cdot (K_H(h, T) \nabla(h+z)) + Q_{AET}(h, T) \quad (1)$$

$$\frac{\partial \left( (C_{T,eq}(h, T) + L \frac{\partial \theta_{ice}(h, T)}{\partial T}) T \right)}{\partial t} + \nabla \cdot (V(h, T) C_{T,liquid} T) = \nabla \cdot (K_{T,eq}(h, T) \nabla T) \quad (2)$$

Then the parameterizations and constitutive laws that are needed to estimate the transfer properties embedded in the three main equations above may be found below (Eq. (3) to (14)). The different variables appearing in the equations are listed and briefly described in Table 1.

- Generalized Darcy's law:

$$V(h, T) = -K_H(h, T) \nabla(h+z) \quad (3)$$
<table><caption>**Table 1**
Variables and parameters of the system of equations solved by permaFoam (in the source code the used units are the SI base units).</caption><tr><td>C<sub>H</sub></td><td>Capillary capacity [m<sup>-1</sup>]</td></tr><tr><td>C<sub>T,air</sub></td><td>Heat capacity of air [J m<sup>-3</sup> K<sup>-1</sup>] (SI base units: [kg m<sup>-1</sup> s<sup>-2</sup> K<sup>-1</sup>])</td></tr><tr><td>C<sub>T,eq</sub></td><td>Equivalent heat capacity [J m<sup>-3</sup> K<sup>-1</sup>] (SI base units: [kg m<sup>-1</sup> s<sup>-2</sup> K<sup>-1</sup>])</td></tr><tr><td>C<sub>T,ice</sub></td><td>Heat capacity of liquid water [J m<sup>-3</sup> K<sup>-1</sup>] (SI base units: [kg m<sup>-1</sup> s<sup>-2</sup> K<sup>-1</sup>])</td></tr><tr><td>C<sub>T,liquid</sub></td><td>Heat capacity of liquid water [J m<sup>-3</sup> K<sup>-1</sup>] (SI base units: [kg m<sup>-1</sup> s<sup>-2</sup> K<sup>-1</sup>])</td></tr><tr><td>C<sub>T,soil</sub></td><td>Heat capacity of the soil [J m<sup>-3</sup> K<sup>-1</sup>] (SI base units: [kg m<sup>-1</sup> s<sup>-2</sup> K<sup>-1</sup>])</td></tr><tr><td>h</td><td>Generalized water pressure head (primary variable of flow equation (1)) [m]</td></tr><tr><td>K<sub>freezing</sub></td><td>Relative hydraulic conductivity with respect to water freezing [-]</td></tr><tr><td>K<sub>freezingmin</sub></td><td>Min. value of relative hydraulic conductivity with respect to water freezing [-]</td></tr><tr><td>K<sub>H</sub></td><td>Apparent hydraulic conductivity [m s<sup>-1</sup>]</td></tr><tr><td>K<sub>rel</sub></td><td>Relative hydraulic conductivity with respect to water saturation [-]</td></tr><tr><td>K<sub>s</sub></td><td>Saturated hydraulic conductivity [m s<sup>-1</sup>]</td></tr><tr><td>K<sub>T,air</sub></td><td>Thermal conductivity of air [W m<sup>-1</sup> K<sup>-1</sup>] (SI base units: [kg m s<sup>-3</sup> K<sup>-1</sup>])</td></tr><tr><td>K<sub>T,eq</sub></td><td>Equivalent thermal conductivity [W m<sup>-1</sup> K<sup>-1</sup>] (SI base units: [kg m s<sup>-3</sup> K<sup>-1</sup>])</td></tr><tr><td>K<sub>T,ice</sub></td><td>Thermal conductivity of ice [W m<sup>-1</sup> K<sup>-1</sup>] (SI base units: [kg m s<sup>-3</sup> K<sup>-1</sup>])</td></tr><tr><td>K<sub>T,liquid</sub></td><td>Thermal conductivity of liquid water [W m<sup>-1</sup> K<sup>-1</sup>] (SI base units: [kg m s<sup>-3</sup> K<sup>-1</sup>])</td></tr><tr><td>K<sub>T,soil</sub></td><td>Thermal conductivity of the soil [W m<sup>-1</sup> K<sup>-1</sup>] (SI base units: [kg m s<sup>-3</sup> K<sup>-1</sup>])</td></tr><tr><td>L</td><td>Volumetric latent heat of fusion of ice [10<sup>3</sup> KJ kg<sup>-1</sup>] (SI base units: [kg m<sup>-1</sup> s<sup>-2</sup>])</td></tr><tr><td>n</td><td>Second Mualem-van Genuchten fitting parameter [-]</td></tr><tr><td>PET</td><td>Potential evapotranspiration [8,64,10<sup>7</sup> mm d<sup>-1</sup>] (SI base units: [m s<sup>-1</sup>])</td></tr><tr><td>Q<sub>AET</sub></td><td>Actual evapotranspiration volumetric rate [s<sup>-1</sup>]</td></tr><tr><td>Q<sub>PET</sub></td><td>Potential evapotranspiration volumetric rate [s<sup>-1</sup>]</td></tr><tr><td>S</td><td>Storage coefficient [m<sup>-1</sup>]</td></tr><tr><td>S<sub>soil</sub></td><td>Total surface of the soil [m<sup>2</sup>]</td></tr><tr><td>t</td><td>Temporal coordinate [s]</td></tr><tr><td>T</td><td>Temperature (primary variable of heat transfer equation (2)) [K]</td></tr><tr><td>T<sub>freeze</sub></td><td>Temperature of starting of the freezing of the water phase [K]</td></tr><tr><td>V</td><td>Darcy velocity [m s<sup>-1</sup>]</td></tr><tr><td>V<sub>evapotranspiration</sub></td><td>Total volume of the evapotranspiration layer [m<sup>3</sup>]</td></tr><tr><td>ω</td><td>Fitting parameter for the freezing function [K]</td></tr><tr><td>z</td><td>Vertical coordinate, oriented upward [m]</td></tr><tr><td>α</td><td>First Mualem-van Genuchten fitting parameter [m-1]</td></tr><tr><td>Δt</td><td>Length of the current time step of numerical solution [s]</td></tr><tr><td>θ</td><td>Soil volumetric water content [-]</td></tr><tr><td>θ<sub>ice</sub></td><td>Volumetric frozen water content [-]</td></tr><tr><td>θ<sub>liquid</sub></td><td>Volumetric liquid water content [-]</td></tr><tr><td>θ<sub>s</sub></td><td>Volumetric water content at saturation [-]</td></tr><tr><td>θ<sub>r</sub></td><td>Volumetric residual water content [-]</td></tr><tr><td>θ<sub>WP</sub></td><td>Soil volumetric water content at the wilting point [-]</td></tr><tr><td>X<sub>evapotranspiration</sub></td><td>Indicatrix of the evapotranspiration layer (= 1 inside; = 0 outside) [-]</td></tr><tr><td>Ω</td><td>Empirical thermal impedance factor [-]</td></tr></table>

- Apparent hydraulic conductivity [37]:

$$K_H(h, T) = K_s K_{rel}(h) K_{freezing}(T) \quad (4)$$

- Mualem-van Genuchten parameterization of soil retention curve for the pressure-based formulation of Richards equation [54,28,4]:

$$\theta(h) = \begin{cases} \theta_s & \text{if } h \ge 0 \\ \theta_r + (\theta_s - \theta_r) (1 + (-\alpha h)^n)^{-(1-\frac{1}{n})} & \text{if } h < 0 \end{cases} \quad (5)$$

$$C_H(h) = S \frac{\theta(h)}{\theta_s} + \frac{\partial \theta(h)}{\partial h} \quad (6)$$

$$K_{rel}(h) = \begin{cases} K_s & \text{if } h \ge 0 \\ K_s \left( (1 + (-\alpha h)^n)^{-(1-\frac{1}{n})} \right)^{\frac{1}{n}} \\ \times \left( 1 - \left( 1 - \left( (1 + (-\alpha h)^n)^{-(1-\frac{1}{n})} \right)^{\frac{n}{n-1}} \right)^{1-\frac{1}{n}} \right)^2 & \text{if } h < 0 \end{cases} \quad (7)$$

- Relative hydraulic conductivity with respect to water freezing [13]:

$$K_{freezing}(T) = \max(10^{-\Omega\theta_{ice}}; K_{freezingmin}) \quad (8)$$

- Soil Freezing Curve function [55]:

$$\theta_{liquid}(h, T) = \begin{cases} \theta(h) & \text{if } T > T_{freeze} \\ \theta_r + (\theta(h) - \theta_r) \exp\left(-\left(\frac{T-T_{freeze}}{\omega}\right)^2\right) & \text{if } T \le T_{freeze} \end{cases} \quad (9)$$

One can refer to Dall'Amico et al. [8] and Painter and Karra [44], for an alternative approach of handling of the freezing behaviour of soils.

- Ice volume fraction calculated from the total water volume fraction:

$$\theta_{ice}(h, T) = \theta(h) - \theta_{liquid}(h, T) \quad (10)$$

One should note here that since we consider a non-deformable porous medium, we neglect the volumetric expansion of ice [30].

- Estimator of actual evapotranspiration [40,41]:

$$Q_{AET}(h, T) = \begin{cases} Q_{PET} & \text{if } \left(\frac{\theta_{liquid}(h, T) - \theta_{WP}}{\Delta t} - Q_{PET}\right) \ge 0 \\ \frac{\theta_{liquid}(h, T) - \theta_{WP}}{\Delta t} & \text{if } -Q_{PET} < \left(\frac{\theta_{liquid}(h, T) - \theta_{WP}}{\Delta t} - Q_{PET}\right) < 0 \\ 0 & \text{if } \left(\frac{\theta_{liquid}(h, T) - \theta_{WP}}{\Delta t} - Q_{PET}\right) \le -Q_{PET} \end{cases} \quad (11)$$
# PermaFoam numerical scheme

<figure><img src="image_5.png" /></figure>

Implemented within the 3D finite volumes HPC framework OpenFOAM

<caption>Fig. 1. Numerical strategy used for the resolution of the two non-linear, coupled equations of permafrost dynamics considered in permaFoam.</caption>

$$Q_{PET} = PET \left( \frac{S_{soil}}{V_{evapotranspiration}} \right) X_{evapotranspiration} \quad (12)$$

- Equivalent thermal conductivity of the medium - geometric mean [10]:

$$K_{T,eq} = (K_{T,liquid})^{\theta_{liquid}(h,T)} (K_{T,ice})^{\theta_{ice}(h,T)} (K_{T,soil})^{1-\theta_s} \times (K_{T,air})^{\theta_s - \theta_{liquid}(h,T) - \theta_{ice}(h,T)} \quad (13)$$

- Equivalent heat capacity of the medium - arithmetic mean [27]:

$$C_{T,eq} = C_{T,liquid}\theta_{liquid}(h, T) + C_{T,ice}\theta_{ice}(h, T) + C_{T,soil}(1 - \theta_s) + C_{T,air}(\theta_s - \theta_{liquid}(h, T) - \theta_{ice}(h, T)) \quad (14)$$

Seven non-linearities and six couplings occur in these two coupled equations, making them especially difficult to solve numerically. This is the main reason why high-performance computing is required here, and thus why we choose to develop permaFoam in the framework of OpenFOAM®. Please note that the numbering of these equations is also used in the comments of the main source file of permaFoam permaFoam.C, for the sake of clarity. Please note also that the units used in permaFoam sources are the SI based units presented in Table 1.

## 3. Numerical methods

### 3.1. Resolution strategy

The numerical methods used to solve these two coupled, non-linear equations are classical: DIC (Diagonal Incomplete Cholesky) pre-conditioning and PCG (Preconditioned Conjugate Gradient) linear solver for the equation of water transfers, DILU (Diagonal Incomplete Lower Upper) pre-conditioning and BiPCG (Preconditioned Bi-Conjugate Gradient) linear solver for the heat transfer equation, two Picard loops (fixed point method), one for the linearization of each equation, sequential operator splitting for the handling of the couplings between the two equations, explicit resolution of the phase change (freeze/thaw of pore water) and implicit/backward Euler method for temporal integration with adaptive time stepping (see [41] for more numerical details). Fig. 1

presents a simple schematization of the adopted numerical strategy.

The use of an operator-splitting approach for dealing with coupled equations is widely adopted in porous media numerics [7,22,35], although it leads to splitting errors that come in addition to the truncature errors related to the finite volume discretisation. Since splitting errors are generally increasing while coarsening the used spatio-temporal discretisations, a peculiar attention should be paid to convergence studies for numerical resolutions of coupled problems based on an operator splitting approaches. A possible alternative to operator splitting would be to deal with a single block coupled discretisation of the equations (e.g.: [3,6]). The adaptive time stepping approach is a so-called Empirically Based Adaptive Strategy approach [62]; it is analogous to the one presented in Orgogozo et al. [39].

Given the strong non linearities and couplings encountered, few references analytical solutions exist for the considered set of equations, and only for very simplified problems. Thus to perform code-to-code intercomparison exercises is a commonly accepted method for validating numerical simulators for cryohydrogeology, and so has been done for the validation of permaFoam implementation for the simulation of coupled heat transfer and water flow with freeze/thaw in saturated porous media [13]. Considering water flow in variably saturated porous media, permaFoam inherits of the validation that has been done for RichardsFoam [39] since it has the same numerics for Richards equation resolution.

### 3.2. New features embedded in this release of permaFoam

This permaFoam package encompasses the new features that have been introduced in RichardsFoam3 [42] along with other specific developments:

(i) better handling of memory use during computation;

(ii) possibility of on-the-fly modification of the control parameters of the linearization procedure;

(iii) new boundary conditions to be compiled along with the solver itself, and dedicated to the simulation of rain infiltration (rainFallFlux) or no infiltration (noRainFlux) conditions while allowing exfiltration with OpenFOAM® stand-alone (i.e. without mandatory use of swak4foam as it was the case in RichardsFoam2 – Orgogozo [40] – see section 3.3);
<figure><img src="image_6.png" /><figcaption>Fig. 2. Seasonally variable thermal boundary condition used for modelling the North Aspected Slopes of the Kulingdakan watershed (see also section 4.1).</figcaption></figure>

(iv) computation of the field of piezometric head along simulation. NB: in the demonstration cases, a way to implement a constant piezometric head boundary conditions with OpenFOAM® is proposed.

(v) new functions that allow to represent seasonal variability of boundary conditions in an easy way with OpenFOAM® stand-alone, using the 'seasonal' (see Section 3.4) function object, and including the 'thawing' option (see section 3.5).

(vi) use of postprocessing procedures with OpenFOAM® stand-alone (i.e. without mandatory use of swak4foam as it was the case in Richards Foam2 - Orgogozo [40]);

(vii) update of the capillary capacity parametrization according to the pressure-based formulation of Richards equation proposed in Kavetski et al. [28] (see equation (6) and associated references).

(viii) several other minor rewrites/cleanings of the code.

### 3.3. The rainFallFlux and the noRainFlux boundary conditions

These two peculiar boundary conditions deal with the hydrologic relations between the sub-surface and the surface in a very simple, sub-surface oriented way, since permaFoam is a cryohydrogeological solver. The goal is to enable to impose an infiltration flux (>0 with rainFallFlux, =0 with noRainFlux) through a boundary (typically the top surface of the soil) while prescribing that in case of saturation of the boundary cells the imposed flux condition is turned into an imposed pressure head condition, with a pressure head equal to the atmospheric pressure. This is equivalent to the top pressure boundary conditions constructed with swak4foam in Orgogozo et al. [41] – one can see the discussion of the application of this kind of boundary conditions to cryohydrogeological modelling in this previous paper. In such a way, the computed infiltration does take account of the fact that no rain water may come in the soil if the soil surface is already saturated, and exfiltration fluxes (water coming out of the soil surface due to subsurface flow) may occur, and be quantified.

### 3.4. The 'seasonal' function object

The seasonal function object has been developed specifically for enabling easy implementation of simple seasonal variation of boundary conditions that are very useful for modelling permafrost, meanwhile it could be applied in other contexts when time varying boundary conditions are needed. It allows to mix two time evolving modes, and each mode can be a sine, a cosine or a square function. Then the two modes are concatenated according to a user specified periodic time frame, using again a basic time dependent function (dualMode true; modeSelector square1;). In the example below (see Fig. 2), extracted from the

initial temperature volScalarField 0/T of the demoCaseSinusoidalClimateForcings (see the tutorials in the permaFoam package), two sinusoidal modes are concatenated using a simple mode selector making mode 1 used when modeSelector=0 and mode2 used when modeSelector=1.

The implementation in the dictionary of the patch top (representing the soil surface) within the 0/T file of the demoCaseSinusoidalClimateForcings tutorial is the following:

```json
'top
{
    type uniformFixedValue;
    uniformValue seasonal;
    uniformValueCoeffs
    {
        dualMode true;
        modeSelector square1;
        modeSelectorCoeffs
        {
            mark 2.35872e+07;
            space 7.8624e+06;
            period constant 3.1536e+07;
            scale constant 1;
            level constant 0;
        }
        mode1 sine1;
        mode1Coeffs
        {
            t0 -2.1681e+07;
            clip (258.65 1000);
            period constant 3.1536e+07;
            scale constant 18.5;
            level constant 267;
        }
        mode2 sine1;
        mode2Coeffs
        {
            t0 -1.5768e+07;
            clip (258.65 1000);
            period constant 3.1536e+07;
            scale constant 13;
            level constant 259;
        }
    }
    value uniform 258.65;
}
```

This syntax encoded the following mathematical function:
Top temperature: $T = \max(\text{clip}; (\text{level} + \text{scale}*\sin((t - t_0)/\text{period})))$
with, for $t \le 23673600$ s,

clip = 258.65 K; level = 267 K; scale = 18.5 K;
t0 = -2.1681*107 s; period = 3.1536*107 s

and for $t > 23673600$ s,

clip = 258.65 K; level = 259 K; scale = 13 K;
t0 = -1.5768*107 s; period = 3.1536*107 s

Note the use of the clipping, for imposing a controlled lower boundary to the time dependent evolution.

There is a strong interest of using continuous time dependent evolutions for boundary conditions, because in case of discontinuities steep fronts occur resulting in a substantial increase of numerical difficulties and hence computing time. In the meantime, using only one purely sinusoidal evolution for catching seasonality may be a problematic over-simplification for some analysis. The seasonal function object has been developed to propose a trade-off between these two concerns.

### 3.5. The 'thawing' option

Additionally, there is a peculiar feature that has been added for establishing soil surface boundary conditions in cold environment: the thawing option for the seasonal function object. When thawing or freezing occurs at the top surface, one may want to start or to stop a given process (for instance, infiltration modelled with a rainFallFlux boundary condition) depending on the thermal status of the faces of the top surface patch. Generally, the permaFoam solver deals natively with such threshold effects by considering a temperature of freezing prescribed in the scalar field Tmelt. Nevertheless, since the apparent hydraulic conductivity of a frozen soil is very low, imposing a water flux with for instance rainFallFlux on a frozen soil may lead to very high pressure heads, and so may generate numerical difficulties. This may occur for instance when a snowpack melts above a frozen soil. Since in such case the amount of infiltrated water is very low, it may be of interest to directly turn the infiltration flux to zero as long as the top soil is frozen. In some cases one may even want to slightly shift the threshold of temperature for easing the numerical resolution, for instance in case of thick cells at the top surface or fast-changing conditions at the patch. This is allowed by the thawing option, which is an indicatrix that has the value 1 when $T >= Tmelt + Tshift$, and the value 0 when $T < Tmelt + Tshift$. One can find an example of use of the thawing options in the 'permaFoam_demoCase_sinusoidalClimateForcings' test case provided in the permaFoam package, in the top patch describing the boundary field of the initial pressure file "0/psi". In this patch, the temperature of start of infiltration is fixed at $Tmelt + Tshift = 273.15 + 0.5$ K:

```R
'top
{
  type rainFallFlux;
  [...]
  thawingActive true;
  rate seasonal;
  rateCoeffs
  {
    [...]
  }
  thawing
  {
    type thawingField;
    Tface T_F;
    Tmelt Tmelt;
    Tshift 0.5;
  }
}
```

```R
}
[...]
}'
```

Such a feature may help to ease significantly the numerical resolution of cases in particularly steep conditions at the period of initiation of thawing or freezing. The approximation related with the use of the 'thawing' option may generate additional errors in the numerical resolutions, so problem-wise trade off must be made between the accuracy of the simulation and the decrease of the computation time. For instance turning off the infiltration when the top soil temperature is close to 0°C may have significant impacts on the modelling of hourly to daily time scale variations of the active layer hydrological and thermal state during snowmelt [16]. Meanwhile, if the interest is in seasonal or multi-annual trends, such errors on high-frequency variabilities might be acceptable, depending of the application at stake.

## 4. Examples of applications

### 4.1. 3D computation in Kulingdakan watershed, Central Siberia

The Kulingdakan watershed has been studied numerically with a special focus on evapotranspiration fluxes in Orgogozo et al. [41]. In this previous study 2D computations of thermo-hydrological fluxes in the slopes of the catchments were carried out. More recently 3D simulations of the permafrost of this watershed have been started, as preliminary works for simulations under climate change to come. One can see in Fig. 3 a brief presentation of the Kulingdakan watershed along with views of the simplified representation used for the spatial discretization of the 3D simulations.

In order to simulate the seasonal permafrost dynamics under current climate conditions, the permafrost is simulated for multi-annual means of climate forcings (top soil temperature, liquid rain rate) with a cycling of simulations along several years to reach a kind of dynamical equilibrium under current conditions (spin-up). One can see on Fig. 4 the computed seasonal evolution of the temperature profiles in the middle of the South aspected slopes and in the middle of the North aspected slopes.

A good agreement is obtained with respect to the field observations. For instance, the computed maximum thicknesses of the active (thawed in summer) layer are 0.8 m in NAS and 1.1 m in SAS, while the observed maximum active layer thicknesses are 0.6 m and 1.2 m in SAS [12]. A slightly better agreement was found on that metric with the 2D computations of Orgogozo et al. [41], likely due to the fact that the spin-up (cycling of several years of modelling for estimating the initial conditions) was better completed in this previous study. In 3D these simulations are costly in terms of CPU hours, requiring several hundreds of hours on OCCIGEN supercomputer using parallel computation on 4000 cores simultaneously, thus we did not push the cycling as far as in this previous study. The used mesh was composed of $57.10^6$ cells, while the average time step was 303 s (minimum time step: 0.76 s; maximum time step: 600 s). This illustrates the need of High Performance Computing for permafrost modelling.

### 4.2. 2D computation in Syrdakh watershed, Eastern Siberia

The Syrdakh study site is located roughly 100 km NE of Yakutsk in Central Yakutia (Eastern Siberia, Russia). A cross section to a river has been instrumented for thermal and hydrological measurements from 2012 on. This landscape unit is a base element of a water catchment (see Fig. 5). Monitoring of temperature is implemented in the river and below (Fig. 5a locations as red crosses). Ground temperatures are monitored at various depths and distances from the river (Fig. 5a - F8, F4, F3 provide monitoring
<figure><img src="image_7.png" /><figcaption>**Fig. 3.** Localization and main features of the Kulingdakan watershed (adapted from Viers et al. [56]), along with the simplified 3D mesh (57.10⁶ cells) considered for its numerical simulation.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 4.** Temperature profiles in North and South Aspected Slopes – good agreement with field observations in terms of active layer depths.</figcaption></figure>

depths at 1, 2, 3, 4 meters, F13 and F14 monitor depths of 30 cm, 1.2 m, 2 and 3 m). Additionally, soil temperature close to the ground across the valley are monitored to provide with surface boundary conditions for the numerical simulation of coupled TH processes (3 zones are discriminated, shadowed and mixed zone, river bed, see Fig. 5a). The hydrological variables monitored are river levels, volumetric water contents at various depths two distances from the river (VWC1 at High and Low locations, displaying VWC at 10, 40 and 70 cm depths). Field studies, generally carried out in September, provide access to topography, soil properties, vegetation cover, supra-permafrost water levels and active layer depths across the valley. Thermal and hydrological parameters are measured from pits. More in Grenier et al. in prep.

The numerical simulation presented here is a preliminary run aiming at providing a basic representation of coupled TH processes and identify data or process gaps. After a spin-up, a converged mean year representing conditions encountered from 2012 to 2019 is obtained. The full set of coupled TH equations including vegetation transpiration and evaporation is solved (Eq. (1) and (2)).

The modelled domain represents the left bank of the valley containing the denser monitoring network: 25 m wide for 5 meters of ground (Fig. 5a). The surface thermal boundary conditions are imposed temperature signals, issued from near surface temperature monitoring (3 zones - shady part, mixed zone and river bed). Surface hydrological conditions are precipitation over land and imposed head for river bed. All other boundary conditions are zero-flux (temperature and water). Evapotranspiration source term is introduced as a term varying exponentially with depth and according to seasons. Average year conditions are repeated from initial conditions until convergence of temperature and liquid volumetric water content fields (spin-up). Temperature and volumetric liquid water content fields are represented in Fig. 5b displaying conditions at mid-September: the deepest active layer zone of the year (unfrozen region above permafrost) is obtained and the influence of the river recharging the supra-permafrost aquifer is visible leading to high liquid water contents as the vegetation has stopped pumping at this time of the year.

Results from this average year are then compared with averaged monitoring series over the 2012-2019 period for all monitoring lo-
<figure><img src="image_9.png" /><figcaption>**Fig. 5.** Overview of the valley, model geometry and simulation maps: (a) view of the modelled valley with simulation domain boundaries, boundary conditions zonation (Shade, Mixed, River) and monitoring locations (F8, F4, F3, F13&14 monitoring temperature - VWC1 at high and low locations monitoring volumetric water content). (b) and (c) simulation maps for respectively temperature and volumetric liquid water content fields at mid-September. (For interpretation of the colours in the figure(s), the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 6.** Inter-comparison of monitoring (solid lines) and simulation results (dotted lines) for temperature evolutions at F4 (1, 2, 3 and 4 m depth) and volumetric liquid water content evolutions at VWC1_H (10, 40, 70 cm depth), both locations are visible in Fig. 5a.</figcaption></figure>

cations (Fig. 5a). Results in Fig. 6 show that seasonal evolution is well captured when it comes to temporalities and trends, for all monitoring depths. The propagation of freezing and thawing are qualitatively well represented while vegetation and evaporation uptake are active during summer, competing with vertical infiltration from precipitations and lateral river inflows. However, the temperature and volumetric water content time series are not quantitatively well-matched requiring future parameter adjustments (esp. freezing curve and vegetation uptake level and history). This preliminary simulation demonstrates the ability to represent realistically the seasonal evolution of a river-valley system with the ad hoc set of complex processes.

## 5. Parallel performances

In this section, we evaluate the numerical performances of permaFoam, and especially its scalability up to a large number of cores. The physical problem used for this study is analogous to the one described by the "democase_sinusoidal_forcing" test case provided together with the solver. Only the spatial scales, and

especially the horizontal extensions of the considered domains, change compared to the basic test case. Three geometries have been used to evaluate performances: they are described by an inclined square with either 102 m, 307 m or 700 m long sides and a thickness of 1 meter, forming a 3D-volume very similar in shape to one side of the geometry studied in the previous section. Each lateral extension corresponds to a different mesh size, from 20.10⁶ to 10⁹ cells, so that the refinement of the discretization is kept unchanged between cases. The boundary conditions include Rain-FallFlux and noRainFlux conditions, while the initial condition is chosen so that thawing occurs in the domain (initial time during June). In this way, this configuration is considered representative of real 3D cases and takes into account all the physical processes that permaFoam models. Performances are evaluated by simulating 1 min of physical time by repeating 60 time steps of 1 second.

The first mesh we study is a 512x512x85 grid, referred as "20M" mesh. We show in Fig. 7 that the topological distribution of decomposed domains associated with each MPI-process on the grid can dramatically affect the performances. The domain is split following one, two or the three directions of space, as each curve
<figure><img src="image_11.png" /><figcaption>**Fig. 7.** Execution times obtained using different MPI decomposition strategies. Each curve corresponds to specific splitting direction(s) indicated on the curve label. Calculations have been performed on OLYMPE supercomputer (Intel® SKYLAKE processors).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 8.** Strong scaling of permaFoam on three different HPC centers. Ideal scalability is displayed with dashed line. Three meshes are tested, with 20M, 200M and 1G cells. Each scalability curve has its own reference core numbers.</figcaption></figure>

label on Fig. 7 indicates; for example, in the (xz) case, the subdomains are created by cutting the geometry in the x and z directions. For each case, we perform the test on a range of cores from 36 to 1800 cores. The dimensionless execution times are calculated by dividing the results by the maximum execution time. A large difference between the different strategies is obtained, since one order of magnitude is separating the fastest decomposition strategy from the slowest one. The best performances are obtained using the (xy) decomposition and so affecting to each MPI-process a subdomain covering the whole thickness of the domain. This result is expected as the main computational charge is located on the thaw front, i.e. close to the surface of the soil. Thus, any decomposition in the (z) direction deteriorates the simulation load-balancing.

Using this decomposition strategy, we evaluate the parallel performances of permaFoam with a strong scaling study on three different HPC centers: CALMIP (Olympe supercomputer), CINES (Occigen supercomputer) and TGCC (IRENE supercomputer, ROME partition), and for three different problems of increasing sizes: one with a 20.10⁶ cells mesh, noted 20M, one with a 200.10⁶ cells mesh, noted 200M, and one with 10⁹ cells mesh, noted 1B. This comparative study gives insights into the parallel performances by including both different hardware environments (processors, architectures, ...), cores number (up to 16384 cores) and computational loads. Results are shown in Fig. 8. PermaFoam demonstrates excellent parallel performances on every tested supercomputers. For instance it's showing results beyond ideal scalability up to 6000 cores for the 20M case on every used supercomputers. This superscalable behaviour, already noticed and discussed in Orgogozo et al. [39], may come from memory cache effects, as well as from a dependence of linear solvers operation conditions to the size of the subdomains, although no effect on results have been observed yet.

When considering a relatively small mesh of 20.10⁶ cells, scalability tends to decrease with higher number of cores, even if a 50% parallel efficiency is still observed at the highest number of cores (16384 cores on IRENE-ROME – Fig. 8 (c))). With large amounts of cores, the communication costs relative to the mesh size are being expected to affect the performances. Indeed, at 4000 cores, where the scalability is starting to decrease, 70% of the cells in memory of a given MPI-process is a ghost-cell, used for communication or boundary condition.

To evaluate this effect, a 200M cells mesh is created by increasing the cell number in both the x and y directions (1536 cells in each direction), as well as the size of the surface plan (307 m long side). In this way, the number of cells is increased without any change in the physical solved or in the size of the mesh cell. With this mesh, scalability of permaFoam is maintained ideal, even with the use of 7000 cores (Fig. 8 (b)). Since the amount of cells per core increased in the 200M case compared to the 20M case, we observe as expected better parallel performances at high number
of cores in the former than in the later. For the 200M test case, linear to super-linear behaviour is observed on all considered supercomputers for the whole range of used MPI processes, up to 16 384 on IRENE-ROME (Fig. 8 (c)). One should note that using a large amount of cores with a large mesh represents itself a particular challenge in terms of numerical resources required, especially for the mesh partition step. For instance 3 days of run are required to decompose the domain of the 200M mesh on 7000 MPI processes using the decomposePar OpenFOAM utility on a fat node with 3TB RAM (Intel® Skylake Xéon Platinum 8176@2.1GHz) available on Occigen supercomputer.

We also demonstrate the capacity of permaFoam to address huge sized problems, since we performed a scalability study up to 1800 cores on a 1billion cells grid (Fig. 8 (a)), with excellent stability and good scalability, with the only need to recompile OpenFOAM to allow the use of large integer numbers (label encoding with 64 bits).

Overall permaFoam, thanks to its development within OpenFOAM, allows to benefit from good parallel performances on different computing centers with different supercomputers. When considering big meshes, the scalability remains excellent even when considering the largest number of cores that are usable in practice for our applications (~10 000) given the usual Quality Of Services (wall time) of the fat nodes available on the used supercomputers for preprocessing operations (especially mesh partitioning).

## 6. Conclusion and perspectives

In this paper we presented permaFoam, an OpenFOAM® solver for permafrost hydrology. Thanks to the good parallel capabilities of OpenFOAM®, and to the inclusion of the main processes involved, permaFoam allows to deal with permafrost modelling for large systems and on multi-annual time scales. It will be soon applied to the study of climate change impacts on the permafrost of several environmentally monitored boreal watersheds in the framework of the HiPerBorea project (hiperborea.omp.eu). One should note that such a tool may also be applied in other contexts where freeze/thaw in porous media occurs, such as for artificial ground freezing [1], geothermal energy [29] or geohazards in high altitude environments [36].

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data availability

The source files of permaFoam are submitted for publication alongside the manuscript. It is open source under the GPL_v3 licence.

## Acknowledgements

The authors acknowledge the help of ESI Group - OpenCFD Ltd, and especially of Mark Olesen, for the set up of this version of permaFoam. This work has been funded by the French Research Agency ANR (grant n° ANR-19 CE46-0003-01), and benefited from access to the supercomputers of CALMIP (project p12166) and GENCI (project A0060110794, CINES and TGCC).

## References

[1] M.A. Alzoubi, M. Xu, F.P. Hassani, S. Poncet, A.P. Sasmito, Tunn. Undergr. Space Technol. 104 (2020) 103534, https://doi.org/10.1016/j.tust.2020.103534.

[2] A. Bartsch, G. Pointner, I. Nitze, A. Efimova, D. Jakober, S. Ley, E. Högström, G. Grosse, P. Schweitzer, Environ. Res. Lett. 16 (2021) 115013, https://doi.org/10.1088/1748-9326/ac3176.

[3] P. Beckstein, V. Galindo, V. Vukčević, J. Comput. Phys. 344 (2017) 623-646, https://doi.org/10.1016/j.jcp.2017.05.005.

[4] N. Brandhorst, D. Erdal, I. Neuweiler, Hydrol. Earth Syst. Sci. 25 (2021) 4041-4059, https://doi.org/10.5194/hess-25-4041-2021.

[5] M.T. Bui, J. Lu, L. Nie, A review of hydrological models applied in the permafrost-dominated Arctic region, Geosciences 10 (2020) 401, https://doi.org/10.3390/geosciences10100401.

[6] C. Busse, I. Tsivilskiy, J. Hildebrand, J.P. Bergmann, Comput. Fluids 216 (2021) 104807, https://doi.org/10.1016/j.compfluid.2020.104807.

[7] J. Carrayrou, R. Mosé, P. Behra, J. Contam. Hydrol. 68 (2004) 239-268, https://doi.org/10.1016/S0169-7722(03)00141-4.

[8] M. Dall'Amico, S. Endrizzi, S. Gruber, R. Rigon, Cryosphere 5 (2) (2011) 469-484.

[9] S. Endrizzi, S. Gruber, M. Dall'Amico, R. Rigon, Geosci. Model Dev. 7 (6) (2014) 2831-2857.

[10] J.M. Frederick, B.A. Buffet, J. Geophys. Res., Earth Surf. 119 (2014) 106-122, https://doi.org/10.1002/2013JF002987.

[11] M. Garcia-Gasulla, F. Banchelli, K. Peiro, G. Ramirez-Gargallo, G. Houzeaux, I. Ben Hassan Saidi, C. Tenaud, I. Spisso, F. Mantovani, Int. J. Comput. Fluid Dyn. 34 (7-8) (2020) 508-528, https://doi.org/10.1080/10618562.2020.1778168.

[12] N. Gentsch, Permafrost Soils in Central Siberia: Landscape Controls on Soil Organic Carbon Storage in a Light Taiga Biome, Akademische Verlagsgemeinschaft München, 2011.

[13] C. Grenier, H. Anbergen, V. Bense, et al., Adv. Water Resour. 114 (2018) 196-218, https://doi.org/10.1016/j.advwatres.2018.02.001.

[14] P. Groisman, H. Shugart, D. Kicklighter, et al., Prog. Earth Planet. Sci. 4 (2017) 41, https://doi.org/10.1186/s40645-017-0154-5.

[15] T. Heinze, Water Resour. Res. 57 (2021) e2021WR030067, https://doi.org/10.1029/2021WR030067.

[16] K.M. Hinkel, S.I. Outcalt, A.E. Taylor, Can. J. Earth Sci. 34 (1997) 667-678.

[17] G. Hu, L. Zhao, T. Wu, X. Xu, H. Park, A. Fedorov, Y. Wei, R. Li, X. Zhu, Z. Sun, J. Ni, D. Zou, Sci. Total Environ. 791 (2021) 148358, https://doi.org/10.1016/j.scitotenv.2021.148358.

[18] X. Huang, D.L. Rudolph, Adv. Water Resour. 154 (2021) 103945, https://doi.org/10.1016/j.advwatres.2021.103945.

[19] X. Huang, D.L. Rudolph, B. Glass, Water Resour. Res. 58 (2022) e2021WR030933, https://doi.org/10.1029/2021WR030933.

[20] A. Jan, E.T. Coon, S.L. Painter, Geosci. Model Dev. 13 (5) (2020) 2259-2276.

[21] A. Jan, S.L. Painter, Environ. Res. Lett. 15 (2020) 084026, https://doi.org/10.1088/1748-9326/ab8ec4.

[22] D. Jacques, J. Šimůnek, D. Mallants, M.Th. van Genuchten, J. Contam. Hydrol. 88 (2006) 197-218, https://doi.org/10.1016/j.conhyd.2006.06.008.

[23] H.-J. Jin, Q.-B. Wu, V.E. Romanovsky, Adv. Clim. Change Res. (February 2021), https://doi.org/10.1016/j.accre.2021.01.007.

[24] L. Jungsberg, L. Byskov Herslund, K. Nilsson, S. Wang, S. Tomaškovičová, K. Madsene, J. Scheerd, T. Ingeman-Nielsen, Polar Geogr. 45 (1) (2022) 58-76, https://doi.org/10.1080/1088937X.2021.1995067.

[25] J. Karlsson, S. Serikova, S.N. Vorobyev, G. Rocher-Ros, B. Denfeld, O.S. Pokrovsky, Nat. Commun. 12 (2021) 825, https://doi.org/10.1038/s41467-021-21054-1.

[26] S. Karra, S.L. Painter, P.C. Lichtner, Cryosphere 8 (5) (2014) 1935-1950.

[27] M. Kaviani, Principles of Heat Transfer in Porous Media, second edition, Springer-Verlag New York, Inc., 1995.

[28] D. Kavetski, P. Binning, S.W. Sloan, Adv. Water Resour. 24 (2001) 595-605.

[29] J.-H. Kupfernagel, J.C. Hesse, M. Schedel, B. Welsch, H. Anbergen, L. Müller, I. Sass, Geotherm. Energy 9 (2021) 24, https://doi.org/10.1186/s40517-021-00206-y.

[30] B.L. Kurylyk, K. Watanabe, Adv. Water Resour. 60 (2013) 160-177, https://doi.org/10.1016/j.advwatres.2013.07.016.

[31] B.L. Kurylyk, K.T.B. MacQuarrie, J.M. McKenzie, Earth-Sci. Rev. 138 (2014) 313-334.

[32] P. Lamontagne-Hallé, J.M. McKenzie, B.L. Kurylyk, J. Molson, L.N. Lyon, WIREs Water 2020 (2020) e1467, https://doi.org/10.1002/wat2.1467.

[33] H. Li, Y. Lai, L. Wang, X. Yang, N. Jiang, L. Li, C. Wang, B. Yang, Cold Reg. Sci. Technol. 157 (2019) 171-186, https://doi.org/10.1016/j.coldregions.2018.10.014.

[34] M.M. Loranty, B.W. Abbot, D. Blok, et al., Biogeosciences 15 (2018) 5287-5313.

[35] R. Lu, T. Nagel, J. Poonoosamy, D. Naumov, T. Fischer, V. Montoya, O. Kolditz, H. Shao, Comput. Geosci. 163 (2022) 105106, https://doi.org/10.1016/j.cageo.2022.105106.

[36] F. Magnin, J.-Y. Josnin, J. Geophys. Res., Earth Surf. 126 (2021) e2021JF006394, https://doi.org/10.1029/2021JF006394.

[37] J.M. McKenzie, C.I. Voss, D.I. Siegel, Adv. Water Resour. 30 (2007) 966-983, https://doi.org/10.1016/j.advwatres.2006.08.008.

[38] J.M. McKenzie, B.L. Kurylyk, M.A. Walvoord, V.F. Bense, D. Fortier, C. Spence, C. Grenier, Cryosphere 15 (2021) 479-484.

[39] L. Orgogozo, N. Renon, C. Soulaine, et al., Comput. Phys. Commun. 185 (2014) 3358-3371, https://doi.org/10.1016/j.cpc.2014.08.004.
[40] L. Orgogozo, Comput. Phys. Commun. 196 (2015) 619-620, https://doi.org/10.1016/j.cpc.2015.07.009.

[41] L. Orgogozo, A.S. Prokushkin, O.S. Pokrovsky, C. Grenier, M. Quintard, J. Viers, S. Audry, Permafr. Periglac. Process. 30 (2019) 75-89, https://doi.org/10.1002/ppp.1995.

[42] L. Orgogozo, Comput. Phys. Commun. 270 (2021) 108182, https://doi.org/10.1016/j.cpc.2021.108182.

[43] S.L. Painter, J.D. Moulton, C.J. Wilson, Hydrogeol. J. 21 (2013) 221-224, https://doi.org/10.1007/s10040-012-0917-4.

[44] S.L. Painter, S. Karra, Vadose Zone J. 13 (4) (2014).

[45] S.L. Painter, E.T. Coon, A.L. Atchley, M. Berndt, R. Garimella, J.D. Moulton, et al., Water Resour. Res. 52 (8) (2016) 6062-6077.

[46] E. Parfenova, N. Tchebakova, A. Soja, Environ. Res. Lett. 14 (2019) 065004, https://doi.org/10.1088/1748-9326/ab10a8.

[47] T.V. Ponomareva, K.Y. Litvintsev, K.A. Finnikov, N.D. Yakimov, A.V. Sentyabov, E.I. Ponomarev, Forests 12 (2021) 994, https://doi.org/10.3390/f12080994.

[48] D. Rettenmaier, D. Deising, Y. Ouedraogo, E. Gjonaj, H. De Gersem, D. Bothe, C. Tropea, H. Marschall, SoftwareX 10 (2019) 100317, https://doi.org/10.1016/j.softx.2019.100317.

[49] L.A. Richards, Physics 1 (1931) 318.

[50] T. Schneider von Deimling, H. Lee, T. Ingeman-Nielsen, S. Westermann, V. Romanovsky, S. Lamoureux, D.A. Walker, S. Chadburn, E. Trochim, L. Cai, J. Nitzbon, S. Jacobi, M. Langer, Cryosphere 15 (2021) 2451-2471, https://doi.org/10.5194/tc-15-2451-2021.

[51] D.A. Streletskiy, L.J. Suter, N.I. Shiklomanov, B.N. Porfiriev, D.O. Eliseev, Environ. Res. Lett. 14 (2019) 025003, https://doi.org/10.1088/1748-9326/aaf5e6.

[52] S.M. Stuenzi, J. Boike, W. Cable, U. Herzschuh, S. Kruse, L.A. Pestryakova, T. Schneider von Deimling, S. Westermann, E.S. Zakharov, M. Langer, Biogeo- sciences 18 (2021) 343-365, https://doi.org/10.5194/bg-18-343-2021.

[53] J.C. Stuurup, S.E.A.T.M. van der Zee, C.I. Voss, H.K. French, Adv. Water Resour. 149 (2021) 103846, https://doi.org/10.1016/j.advwatres.2021.103846.

[54] M.T. van Genuchten, Soil Sci. Soc. Am. J. 44 (5) (1980) 892.

[55] P. Vidstrand, S. Follin, J-O. Selroos, J-O. Näslund, Hydrogeol. J. 22 (2014) 1251-1267, https://doi.org/10.1007/s10040-014-1164-7.

[56] J. Viers, A.S. Prokushkin, O.S. Pokrovsky, A.V. Kirdyanov, C. Zouiten, J. Chmeleff, M. Meheut, F. Chabaux, P. Oliva, B. Dupré, Geochem. Trans. 16 (2015) 3, https://doi.org/10.1186/s12932-015-0018-0.

[57] I. Vladislav, K. Pavel, S. Dmitrii, in: A. Petriaev, A. Konon (Eds.), Transportation Soil Engineering in Cold Regions, Volume 1, in: Lecture Notes in Civil Engineering, vol. 49, Springer, Singapore, 2020.

[58] M.A. Walvoord, B.L. Kurylyk, Vadose Zone J. 15 (6) (2016), https://doi.org/10.2136/vzj2016.01.0010.

[59] M.A. Walvoord, R.G. Striegl, Front. Clim. 3 (2021) 730402, https://doi.org/10.3389/fclim.2021.730402.

[60] J. Wang, N. Kumar Shrestha, M. Aghajani Delavar, T. Worku Meshesha, S.N. Bhanja, Water 13 (2021) 518, https://doi.org/10.3390/w13040518.

[61] S. Westermann, M. Langer, J. Boike, M. Heikenfeld, M. Peter, B. Etzelmüller, G. Krinner, Geosci. Model Dev. 9 (2016) 523-546, https://doi.org/10.5194/gmd-9-523-2016.

[62] G.A. Williams, C.T. Miller, Adv. Water Resour. 22 (8) (1999) 831.

[63] X. Yang, J. Hu, R. Ma, Z. Sun, Front. Earth Sci. 9 (2021) 721009, https://doi.org/10.3389/feart.2021.721009.

[64] J. Yao, L. Gu, C. Yang, Y. Ding, R. Li, L. Zhao, J. Shi, C. Chen, Int. J. Climatol. 40 (2020) 5783-5800, https://doi.org/10.1002/joc.6551.

[65] L. Yu, Y. Zeng, Z. Su, Hydrol. Earth Syst. Sci. 24 (2020) 4813-4830, https://doi.org/10.5194/hess-24-4813-2020.

[66] J.M. Wani, R.J. Thayyen, C. Shekhar Prasad Ojha, S. Gruber, Cryosphere 15 (2021) 2273-2293, https://doi.org/10.5194/tc-15-2273-2021.

[67] C. Zheng, J. Šimůnek, Y. Zhao, Y. Lu, X. Liu, C. Shi, H. Li, L. Yu, Y. Zeng, Z. Su, J. Hydrol. 598 (2021) 126250, https://doi.org/10.1016/j.jhydrol.2021.126250.