# Assessing simplified approaches in modeling rainfall-induced landslides using Richards' equation with Biot poroelasticity

Rashad Abbasov, Marwan Fahs, Vincent Fontaine, Husam Musa Baalousha, Anis Younes, Renaud Toussaint

► To cite this version:

Rashad Abbasov, Marwan Fahs, Vincent Fontaine, Husam Musa Baalousha, Anis Younes, et al. "Assessing simplified approaches in modeling rainfall-induced landslides using Richards' equation with Biot poroelasticity." Engineering Geology, 2026, 360, pp.108481. 10.1016/j.enggeo.2025.108481. hal-05401172

HAL Id: hal-05401172

https://hal.science/hal-05401172v1

Submitted on 5 Dec 2025

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Assessing simplified approaches in modeling rainfall-induced landslides using Richards' equation with Biot poroelasticity

Rashad Abbasov<sup>a,b</sup>, Marwan Fahs<sup>a,*</sup>, Vincent Fontaine<sup>c</sup>, Husam Musa Baalousha<sup>d</sup>, Anis Younes<sup>a</sup>, Renaud Toussaint<sup>a,e</sup>

## ARTICLE INFO

**Keywords:**

Rainfall induced landslide

Richards' equation

Biot Poroelasticity

Coupled hydromechanical models

Local factor of Safety

## ABSTRACT

Coupled Hydro-Mechanical (HM) models, which involve the interaction between fluid flow and mechanical processes in porous media, are essential for simulating rainfall-induced landslides (RILS). A common strategy to investigate RILS is through the Biot poroelasticity approach combined with the Local Factor of Safety concept. While the two-phase flow and mechanical equations are well-documented for modeling HM processes in variably saturated conditions, the use of Richards' Eq. (RE) in this context remains underexplored. This paper aims to address this gap by providing a comprehensive formulation for coupling RE with mechanical equations and evaluating the impact of various model simplifications on the reliability of RILS predictions. We present the theoretical development of the original mathematical model that leads to a strong coupling between flow and mechanical equations. We also review the different simplification approaches proposed in the literature to reduce the complexity of the original model. A straightforward implementation of these approaches in COMSOL is provided, with validation against established benchmarks. Through a comparative analysis of different approaches, using both hypothetical benchmark and real case study in Azerbaijan, this study provides new insights into the balance between model complexity and prediction accuracy. While simplified approaches are commonly used in literature, this study advocates for the fully coupled model, which delivers the most reliable results. The employment of advanced numerical techniques enhances the computational performance of the fully coupled model, making it competitive with simplified approaches.

## 1. Introduction

The interactions between fluid flow and mechanical processes in a porous medium are usually referred to as hydro-mechanical (HM) processes (Feder et al., 2022; Xue et al., 2020). These processes are coupled because the hydraulic and mechanical phenomena interact, with changes in one can affect the other (Zhao, 2023). Coupled HM processes in variably saturated porous media are relevant to several engineering applications such as in underground constructions (Gong et al., 2022), enhanced geothermal energy (Jeanne et al., 2014), nuclear waste repositories (Yang and Fall, 2021), hydrogen storage (Bai and Tahmasebi, 2022) and oil and gas reservoirs management (Yan et al., 2021). The main interest of this work is the soil stability on slopes, where losses of

mechanical balance can cause landslides (Wei et al., 2025; Paronuzzi and Bolla, 2023). Intensive rainfall infiltration through the soil body is one of the key factors contributing to slope instabilities (Wu and Zhou, 2023; Wu et al., 2020; Fang et al., 2025). This phenomenon, known as rainfall-induced landslides (RILS) causes severe losses of life and property damage (Abbasov et al., 2024; Chen et al., 2021). HM processes provide deep insights into the analysis of RILS (Kafle et al., 2022). Therefore, RILS is one of the most fundamental applications of the HM processes in a variably saturated porous medium.

Physics-based models are efficient tools in applications involving HM processes in a variably saturated porous medium (Yan et al., 2022). These models are commonly used for various purposes, such as understanding and predicting processes, validating experimental results, and

<footnote><sup>a</sup> Institut Terre et Environnement de Strasbourg, Université de Strasbourg, CNRS, ENGEES, UMR 7063, Strasbourg, France</footnote>

<footnote><sup>b</sup> Azerbaijan State Oil and Industry University, Baku, Azerbaijan</footnote>

<footnote><sup>c</sup> ESIROI, University of La Réunion, South Campus, France</footnote>

<footnote><sup>d</sup> Department of Geosciences, College of Petroleum Engineering and Geosciences, King Fahd University of Petroleum and Minerals (KFUPM), Saudi Arabia</footnote>

<footnote><sup>e</sup> PoreLab, The Njord Centre, Department of Physics, University of Oslo, Oslo, Norway</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: fahs@unistra.fr (M. Fahs).</footnote>
designing systems. In physics-based models, the physical processes are described using partial differential equations that are solved numerically using computer codes. The application of physics-based HM models to RILS is receiving significant attention because they can help in predicting the risk of RILS and build reliable warning information systems (Yang and Yeh, 2023; Zhan et al., 2022; Chiu et al., 2019; Yang et al., 2017). HM processes in a variably saturated porous medium are typically modeled by coupling the two-phase flow and the mechanical equations through the Biot poroelasticity approach. The two-phase flow equations are based on the continuity equation and Darcy's law for the wetting and non-wetting phases (Asadi and Ataie-Ashtiani, 2015; Toussaint et al., 2012; Lewis and Schrefler, 1998). The mechanical equations are based on the local equilibrium law (Moradi et al., 2018; Lu et al., 2012). Solving them requires a behavior law that expresses the stress tensor as a function of the strain tensor. The linear elastic assumption is applied in several applications involving HM processes in variably saturated soils (Abbasov et al., 2024; Shao et al., 2015; Shao et al., 2014). The porous medium can be compressible, as both the matrix and solid grains are deformable (Wang et al., 2019; Eriksen et al., 2018). The pore space between particles can be variably saturated. Deformation of the porous media can occur due to porous matrix stiffness and/or fluid flow processes. The latter can also be affected by the deformation processes. The simultaneous occurrence of fluid flow and matrix deformation is represented by the theory of Biot poroelasticity (Biot, 1941).

In several applications, especially in the vadose zone, the isobaric approach gains interest over the two-phase model (Bandara et al., 2016). This approach relies on Richards' approximation, which assumes that gas mobility is infinitely large. According to this approximation, Richards' Eq. (RE) combines Darcy's law and mass conservation to describe the wetting phase processes, while neglecting the mass balance equation for the non-wetting phase due to the static gas pressure assumption. RE is plausible if the permeability of the gas phase is large and drains freely (Grunwald et al., 2022; Wang et al., 2011). Even though RE is widely used for modeling the unsaturated flow in the vadose zone (Luo et al., 2025; Jiang et al., 2025; Lan et al., 2025; Rangarajan et al., 2024), its coupling with mechanical processes is relatively less investigated in the literature. This paper aims to provide a better understanding of the coupling between the RE, mechanical equations, and Biot poroelasticity approach for modeling HM processes in variably saturated porous media and on its applications to RILS. In particular, the paper addresses the following research gaps:

i. Mathematical formulations of the HM processes using the two-phase approach are widely discussed in the literature (Moradi et al., 2024; Asadi and Ataie-Ashtiani, 2015; Schrefler and Scotta, 2001). However, the application of the RE is still not fully understood. There is currently no comprehensive study on the derivation of a mathematical model that allows for coupling poroelasticity with RE. It is also unclear how the RE can be adapted to consider the mechanical processes, and how the mechanical equations can be modified for unsaturated flow processes with static gas pressure distribution.

ii. The approach that couples the RE with mechanical equation has been used in the literature for modeling RILS under varying levels of simplification. For instance, several studies reduce the nonlinearity of the equations by neglecting pore pressure in the mechanical equations or the effect of poroelasticity on the flow processes (Moradi et al., 2018; Shao et al., 2015). Other studies have used fully coupled equations, as demonstrated by Kafle et al. (2022). Some research works simplify the governing equations by neglecting the compressibility of the solid grains (Bandara et al., 2016). However, there are currently no comparative studies that investigate the effect of these simplifications on RILS predictions.

iii. Solving the coupled unsaturated flow and mechanical equations under the poroelastic assumption is a challenging task because of

the nonlinearity of the governing equations. While several approaches have been developed in the literature, most existing studies focus on saturated porous medium (Clément et al., 2018; Ghani et al., 2013; Niebling et al., 2012; Johnsen et al., 2008). The main issues discussed in the literature include spatial discretization, time integration, and coupling between flow and mechanical processes. COMSOL is a powerful tool for Multi-physics simulations. It is a finite element package that provides advanced capabilities in simulating multi-physical problems in engineering or scientific research applications. The software's user-friendly interface and high accuracy in the calculations facilitate an efficient progression through the various modeling steps. Despite these advantageous features, COMSOL lacks direct functionality to simulate coupled HM models for variably saturated soils with RE (Rafraf and Guellouz, 2022).

The paper aims to address the gaps discussed above and has three primary objectives. First, it seeks to provide a comprehensive insight into the mathematical formulation for coupling the RE with mechanical equations to model HM processes in variably saturated porous media, and its application to RILS. The second objective is to develop a straightforward implementation of COMSOL for simulating HM processes in variably saturated porous media using the RE. Finally, the paper aims to investigate the relationship between the complexity of the different mathematical HM models and their reliability in the context of RILS.

Thus, we present a comprehensive theoretical development of the mathematical equations and discuss different simplifications adopted in the literature. Next, we introduce an efficient technique in COMSOL for modeling HM processes in variably saturated soils with the RE. We validate the results from COMSOL against benchmarks established in the literature. The concept of the local factor of safety is used to investigate the slope stability. Later, we take advantage of the flexibility of COMSOL to perform a comparative study of the different simplifications of the HM model for RILS. This is done by simulating a hypothetical benchmark problem as well as a real case study in the municipality area of Badamdar in Azerbaijan.

## 2. Governing equations of RILS: Fully coupled approach and simplifications

### 2.1. Fully coupled mathematical HM approach (FC-HM)

The continuity equation combines mass conservation equations for both fluid and solid phases. Under Richards' assumption, the mass conservation equation for the air phase can be disregarded by assuming static air pressure distribution (Asadi and Ataie-Ashtiani, 2015; Schrefler and Scotta, 2001; Lewis and Schrefler, 1998; Both et al., 2019; Bui et al., 2020). The continuity equations for the fluid and solid phases are given as follows:

$$ \frac{\partial(n\rho_w S_w)}{\partial t} + \nabla \cdot (n\rho_w S_w \mathbf{v}_w) = 0 \quad (1) $$

$$ \frac{\partial[(1-n)\rho_s]}{\partial t} + \nabla \cdot [(1-n)\rho_s \mathbf{v}_s] = 0 \quad (2) $$

where $\rho_w$ [kg.m$^{-3}$] is the density of water, $\rho_s$ [kg.m$^{-3}$] is the density of the solid grains, $S_w$ [–] is the degree of water saturation, $n$ [–] is the porosity, $\mathbf{v}_w$ [m.s$^{-1}$] and $\mathbf{v}_s$ [m.s$^{-1}$] are the absolute velocities of the fluid phase and solid phase, respectively.

The mass conservation of the water phase is reformulated with Darcy's velocity which can be expressed in terms of the relative velocity of the fluid to the solid grains, as follows:

$$ \mathbf{q}_w = nS_w \mathbf{v}_{ws} = nS_w (\mathbf{v}_w - \mathbf{v}_s) \quad (3) $$
where $q_w$ [$m.s^{-1}$] is Darcy's velocity and $v_w$ [$m.s^{-1}$] is the relative velocity of the fluid with respect to the solid grains which is calculated as the difference between the absolute velocities of fluid and solid phases.

The expanded form of the fluid continuity equation is derived from Eqs. (1)–(3), as follows:

$$nS_w \frac{\partial \rho_w}{\partial t} + \rho_w S_w \frac{\partial n}{\partial t} + n \rho_w \frac{\partial S_w}{\partial t} + \rho_w \nabla \cdot \mathbf{q}_w + \mathbf{q}_w \nabla \rho_w + n \rho_w S_w \nabla \cdot \mathbf{v}_s + \mathbf{v}_s \cdot \nabla (n \rho_w S_w) = 0 \quad (4)$$

The term $q_w \nabla \rho_w$ is neglected in comparison to the term $\rho_w \nabla \cdot q_w$ due to the difference in their magnitudes, which are determined based on the characteristic length (details are provided in Appendix A). Based on the same analysis, the term $v_s \cdot \nabla (n \rho_w S_w)$ can be neglected in comparison to the term $n \rho_w S_w \nabla \cdot v_s$.

$$nS_w \frac{\partial \rho_w}{\partial t} + \rho_w S_w \frac{\partial n}{\partial t} + n \rho_w \frac{\partial S_w}{\partial t} + \rho_w \nabla \cdot \mathbf{q}_w + n \rho_w S_w \nabla \cdot \mathbf{v}_s = 0 \quad (5)$$

By assuming slightly compressible fluid and by using the definition of fluid compressibility (see Eq. B1 in Appendix B), Eq. (5) can be rewritten by dividing all terms by $\rho_w S_w$:

$$nC_w \frac{\partial p_w}{\partial t} + \frac{\partial n}{\partial t} + \frac{n}{S_w} \frac{\partial S_w}{\partial t} + \frac{1}{S_w} \nabla \cdot \mathbf{q}_w + n \nabla \cdot \mathbf{v}_s = 0 \quad (6)$$

where $C_w$ [$kg^{-1}.m^1.s^2$] is compressibility of water and $p_w$ [$kg.m^{-1}.s^{-2}$] is the water pressure.

In the same manner, mass balance equation for the solid phase can be expanded as follows:

$$(1-n) \frac{\partial \rho_s}{\partial t} - \rho_s \frac{\partial n}{\partial t} + \mathbf{v}_s \cdot \nabla [(1-n)\rho_s] + [(1-n)\rho_s] \nabla \cdot \mathbf{v}_s = 0 \quad (7)$$

The term $\mathbf{v}_s \cdot \nabla [(1-n)\rho_s]$ is neglected in comparison to the term $[(1-n)\rho_s] \nabla \cdot \mathbf{v}_s$ due to the difference in their magnitudes, which are determined based on the characteristic length (details are provided in Appendix A). Eq. (7) is divided by $\rho_s$ to rewrite new form of mass balance equation for solid phase:

$$\frac{(1-n)}{\rho_s} \frac{\partial \rho_s}{\partial t} - \frac{\partial n}{\partial t} + \nabla \cdot \mathbf{v}_s - n \nabla \cdot \mathbf{v}_s = 0 \quad (8)$$

The conservation equation for total mass is developed by adding mass balance equations for fluid and solid phases (Eqs. (6) and (8)) together and by multiplying the resulting equation by $S_w$ (Grunwald et al., 2022; Varela, 2018; Bandara et al., 2016; Asadi and Ataie-Ashtiani, 2015; Schrefler and Scotta, 2001; Lewis and Schrefler, 1998).

$$nS_w C_w \frac{\partial p_w}{\partial t} + n \frac{\partial S_w}{\partial t} + \nabla \cdot \mathbf{q}_w + S_w \frac{(1-n)}{\rho_s} \frac{\partial \rho_s}{\partial t} + S_w \nabla \cdot \mathbf{v}_s = 0 \quad (9)$$

The porosity is time-dependent, however combining continuity equations for both the fluid and solid phases eliminates the time derivative of the porosity, as shown in Eq. (9).

The time derivative of the solid density can be reformulated as follows (see Eqs. B2-B6 in Appendix B):

$$\frac{\partial \rho_s}{\partial t} = \frac{\rho_s C_s}{1-n} \left( -\frac{\partial \sigma}{\partial t} - n S_w \frac{\partial p_w}{\partial t} - n p_w \frac{\partial S_w}{\partial t} \right) \quad (10)$$

where $C_s$ [$kg^{-1}.m^1.s^2$] is compressibility of solid grains, $\sigma$ [$kg.m^{-1}.s^{-2}$] is isotropic component of the total stress.

The isotropic component of the total stress ($\sigma$) can be expressed as a function of the isotropic component of the effective stress (see Eq. B7 in Appendix B). Thus, the time derivative of the total stress can be expressed as follows:

$$\frac{\partial \sigma}{\partial t} = \frac{\partial \sigma_e}{\partial t} - \alpha p_w \frac{\partial S_w}{\partial t} - \alpha S_w \frac{\partial p_w}{\partial t} \quad (11)$$

where $\sigma_e$ [$kg.m^{-1}.s^{-2}$] is isotropic component of the effective stress and $\alpha$ $[-]$ is the Biot's coefficient, which is defined as follows:

$$\alpha = 1 - \frac{C_s}{C_m} \quad (12)$$

We should mention that Eq. (11) is derived by assuming the air pressure as zero. The pressures and stresses are measured with respect to it. This can be released in a more general expression of Eq. (11) by replacing $p_w$ by $p_w - p_a$. The time derivative of isotropic component of the effective stress can be expressed as a function of the compressibility of the medium and divergence of the velocity of the solid phase:

The time derivative of isotropic effective stress can be expressed as a function of the compressibility of the medium and divergence of the velocity of the solid phase:

$$\frac{\partial \sigma_e}{\partial t} = \frac{\partial \sigma_e}{\partial \varepsilon_v} \frac{\partial \varepsilon_v}{\partial t} = \frac{1}{C_m} \nabla \cdot \mathbf{v}_s \quad (13)$$

where $\varepsilon_v$ $[-]$ is the volumetric strain and $C_m$ [$kg^{-1}.m^1.s^2$] is the compressibility of the medium (i.e. matrix).

In turn, the divergence of the velocity of the solid phase has a relation with the displacement field:

$$\nabla \cdot \mathbf{v}_s = \frac{\partial \varepsilon_v}{\partial t} = \frac{\partial}{\partial t} (\nabla \cdot \mathbf{u}) \quad (14)$$

Substituting Eqs. (10)–(14) in Eq. (9) leads to the following equation (Both et al., 2021; Varela, 2018; Asadi and Ataie-Ashtiani, 2015):

$$[(\alpha - n)C_s S_w^2 + n C_w S_w] \frac{\partial p_w}{\partial t} + [(\alpha - n)C_s S_w p_w + n] \frac{\partial S_w}{\partial t} + \alpha S_w \frac{\partial}{\partial t} (\nabla \cdot \mathbf{u}) + \nabla \cdot \mathbf{q}_w = 0 \quad (15)$$

where $\mathbf{u}$[$m$] is the displacement field.

Besides the mass conservation equations, modeling fluid flow in porous media also requires a linear momentum conservation equation that is based on Darcy's law for the fluid phase. This equation is given as follows:

$$\mathbf{q}_w = - \mathbf{k} \frac{k_{rel}}{\mu_w} (\nabla p_w - \rho_w \mathbf{g}) \quad (16)$$

$\mathbf{k}$[$m^2$] is intrinsic permeability tensor of the soil, $k_{rel}$ $[-]$ is the relative permeability, $\mu_w$ [$kg.m^{-1}.s^{-1}$] is the dynamic viscosity of the fluid, and $\mathbf{g}$[$m.s^{-2}$] is the acceleration vector due to gravity.

The Eqs. (15) and (16) are closed with the soil-water retention model, which expresses water saturation as a function of capillary pressure (i.e. $p_c = -p_w$). Brooks and Corey model is used in this work:

$$S_e = \frac{S_w - S_r}{1 - S_r} = \begin{cases} \left| \frac{p_w}{p_{ent}} \right|^{-n_{BC}} & \text{if } p_w < p_{ent} \\ 1 & \text{if } p_w \ge p_{ent} \end{cases} \quad (17)$$

$$k_{rel} = (S_e)^{\ell_{BC} + 2 + \frac{2}{n_{BC}}} \quad (18)$$

where $S_e$ $[-]$ is the effective saturation, $S_w$ $[-]$ is the degree of water saturation, $S_r$ $[-]$ is the residual water saturation, $p_{ent}$ [$kg.m^{-1}.s^{-2}$] represents air entry pressure (i.e. $p_{ent} = \rho_w \mathbf{g} / \alpha_{BC}$) in which $\alpha_{BC}$ [$m^{-1}$] is the air entry pressure head, $n_{BC}$ $[-]$ and $\ell_{BC}$ $[-]$ characterize pore size distribution in soil and tortuosity, respectively.

The mechanical equilibrium equation describes soil deformation processes. This equation should be coupled with the fluid flow equations through water saturation, as changes in water content directly affects the gravity term of the mechanical equilibrium equation and the Bishop's parameter (Terzaghi, 1925; Ben-Zeev et al., 2023; Clément et al., 2018; Moradi et al., 2018; Shao et al., 2015):
$$ \nabla \cdot \sigma + (\rho_s(1-n) + n\rho_w S_w)g = 0 \quad (19) $$

where $\sigma$ [$kg.m^{-1}.s^{-2}$] is the total stress tensor.

In compressive soils, total stress is derived by subtracting pore pressure from effective stress. Biot's coefficient represents the effect of poroelasticity on the stress tensor. Bishop's parameter describes the impact of variable saturation on pore pressure.

$$ \sigma = \sigma_{\text{eff}} - \alpha\chi(S_w)p_w I \quad (20) $$

where $\sigma_{\text{eff}}$ [$kg.m^{-1}.s^{-2}$] is the effective stress tensor and $I$ [-] is the second order identity tensor, $\chi$[-] is Bishop's parameter.

Different models can be used for the Bishop's parameter. In this work it is defined as the degree of saturation (Lu et al., 2010; Lu and Likos, 2006). The effective stress tensor can be expressed as function of the strain tensor using Hooke's law:

$$ \sigma_{\text{eff}} = C : \varepsilon \quad (21) $$

where $C$ [$kg.m^{-1}.s^{-2}$] represents the elastic stiffness tensor and $\varepsilon$[-] is the strain tensor.

The final equation to be solved is derived using the displacement field as the primary variable. This is done using the relationship between the strain tensor and the displacement field:

$$ \varepsilon = \frac{1}{2}(\nabla u + \nabla u^r) \quad (22) $$

The final system to be solved consists of Eqs. (15)–(22). This system is referred in this paper as a Fully Coupled HM approach (FC-HM).

For the slope stability analysis, we use the concept of Local Factor of Safety (LFS). LFS[-] is defined as the ratio of stress for the potential failure state ($\tau^*$ [$kg.m^{-1}.s^{-2}$]) to stress of the current state ($\tau$ [$kg.m^{-1}.s^{-2}$]) (Li et al., 2025; Moradi et al., 2018; Lu et al., 2012).

$$ LFS = \frac{\tau^*}{\tau} = \frac{\cos(\varphi') \cdot (c' + \sigma_{\text{eff}}^I \cdot \tan(\varphi'))}{\sigma_{\text{eff}}^{II}} \quad (23) $$

where $c'$ [$kg.m^{-1}.s^{-2}$] is the effective soil cohesion parameter, $\varphi'$[-] is the effective friction angle, the effective stresses $\sigma_{\text{eff}}^I$ [$kg.m^{-1}.s^{-2}$] and $\sigma_{\text{eff}}^{II}$ [$kg.m^{-1}.s^{-2}$] are the center and the radius of the Mohr circle, respectively. Center and radius of the Mohr circle can be found from minor and major effective stress eigenvalues as follows:

$$ \sigma_{\text{eff}}^I = \frac{\sigma_{\text{eff}}^1 + \sigma_{\text{eff}}^3}{2} \quad (24) $$

$$ \sigma_{\text{eff}}^{II} = \frac{\sigma_{\text{eff}}^1 - \sigma_{\text{eff}}^3}{2} \quad (25) $$

## 2.2. Governing equations: Constant solid grains density assumption (CSGD-HM)

Several studies on RILS have suggested more simplified versions of the FC-HM approach. These simplifications reduce the complexity of the governing equations and, consequently, the computational requirements of the corresponding numerical solutions. Solving the mass conservation of the solid and water phases is computationally challenging because of the high nonlinearity and the strong coupling between the flow and the mechanical processes. Solving this system of equations numerically is highly expensive in CPU time. This complexity may limit the applicability of the FC-HM approach in real applications at large time and space scales. To address this issue, several assumptions have been made to simplify the equations. In this context, many studies have adopted the assumption of incompressible soil grains ($C_s = 0$) with a deformable matrix (soil skeleton) (Rafraf and Guellouz, 2022; Bandara et al., 2016; Guellouz et al., 1997). Based on Eq. (12), this assumption leads to a Biot's coefficient equal to one (Taherdangkoo et al., 2022;

Kafle et al., 2022). In this work, this approach is referred to Constant Solid Grains Density Assumption (CSGD-HM). The continuity equation of the CSGD-HM approach is derived by substituting Eq. (14) in Eq. (9) while neglecting the time derivative of solid grain density (i.e. $\partial\rho_s/\partial t$) and assuming the Biot's coefficient equal to one:

$$ nS_w C_w \frac{\partial p_w}{\partial t} + n \frac{\partial S_w}{\partial t} + \nabla \cdot \mathbf{q}_w + \alpha S_w \frac{\partial}{\partial t} (\nabla \cdot \mathbf{u}) = 0 \quad (26) $$

Darcy's law, Brooks and Corey, mechanical equations, and the LFS evaluation for the CSGD-HM are the same as those in the FC-HM approach. As per Eq. (12), $\alpha$ is assumed to be equal to 1 in the CSGD-HM approach.

## 2.3. Governing equations: Partial decoupling between HM processes (PDC-HM)

Several studies have suggested a more simplified mathematical HM model for RILS by reducing the degree of coupling between the flow and mechanical processes. This is achieved by neglecting the time derivative of the volumetric strain term in the flow equation, thereby ignoring the impact of mechanical deformation on flow (Ireson et al., 2023; Farthing and Ogden, 2017; Fahs et al., 2009). However, this is not a fully decoupled approach because it considers the impact of fluid pressure on deformation processes (Ireson et al., 2023). This study also neglects fluid compressibility. With these assumptions, the governing equations can be solved sequentially. Flow is first solved to obtain the pressure. The latter is then used to solve the mechanical equations and update the displacement field. This partial decoupling of flow and mechanical processes is referred to as PDC-HM. The flow equation is given as follows:

$$ n \frac{\partial S_w}{\partial t} + \nabla \cdot \mathbf{q}_w = 0 \quad (27) $$

This equation is closed with Darcy's law (Eq. (16)) and the Brooks and Corey model (Eqs. (17) and (18)). The mechanical equations and the evaluation of the LFS are the same as those in the CSGD-HM and the FC-HM approaches.

## 2.4. Governing equations: Assumption of uncoupled HM processes (UNC-HM)

Several studies neglect the impact of fluid pressure on the mechanical processes while still incorporating it in the evaluation of the LFS. This approach is widely used in the literature because it allows for significant simplification of the governing equations by reducing the degree of coupling between HM processes (Moradi et al., 2018; Shao et al., 2015; Abbasov et al., 2024). This approach is implemented in the hill-slope cube module which is widely used in the literature (Yang and Yeh, 2023; Raia et al., 2014; Lu et al., 2012). Although the flow and mechanical processes are not fully uncoupled since water content is used to evaluate the gravity term in the mechanical equation the coupling is limited to the water content. This approach is referred to as UNC-HM (uncoupled HM processes). The flow equation (i.e. continuity equation) is the same as Eq. (27). The pore pressure is neglected from the mechanical equilibrium Eq. (28) ($\tilde{\sigma} = -\sigma$).

$$ \nabla \cdot \tilde{\sigma} + (\rho_s(1-n) + n\rho_w S_w)g = 0 \quad (28) $$

The impact of pore pressure is considered in the evaluation of the LFS. The LFS is calculated using Eq. (29), with $\tilde{\sigma}_{\text{eff}}^I$ and $\tilde{\sigma}_{\text{eff}}^{II}$ ($\tilde{\sigma}_{\text{eff}} = -\sigma_{\text{eff}}$) evaluated as follows (30)–(31) (Moradi et al., 2018; Lu et al., 2012):

$$ LFS = \frac{\tau^*}{\tau} = \frac{\cos(\varphi') \cdot (c' + \tilde{\sigma}_{\text{eff}}^I \cdot \tan(\varphi'))}{\tilde{\sigma}_{\text{eff}}^{II}} \quad (29) $$
<figure><img src="image_3.png" /><figcaption>Fig. 1. Schematic of the flow equations implemented in COMSOL for the four approaches: FC-HM, CSGD-HM, PDC-HM, and UNC-HM. Flow channels and solid grain geometry inspired by the “pore scale flow” model in COMSOL Library.</figcaption></figure>

$$ \tilde{\sigma}_{eff}^{I} = \frac{\tilde{\sigma}^{1} + \tilde{\sigma}^{3}}{2} - \alpha \chi(S_{w}) p_{w} \qquad (30) $$

$$ \tilde{\sigma}_{eff}^{II} = \frac{\tilde{\sigma}^{1} - \tilde{\sigma}^{3}}{2} \qquad (31) $$

### 3. Implementation in COMSOL multiphysics

We use COMSOL Multiphysics to solve the governing equations of the FC-HM, CSGD-HM, PDC-HM, and UNC-HM approaches. This software is an advanced finite element method package that supports a range of physics simulations, both individually and in combination to perform multi-physics simulations. COMSOL also includes some inherent Multiphysics modules. In COMSOL, it is possible to simulate HM processes in saturated porous media with the Biot poroelasticity model via the “Poroelasticity” module. This module is based on the coupling between the “Darcy’s Law” and “Solid Mechanics” modules. However, the current versions of COMSOL have no functionality in processing simulating HM processes in variably saturated porous media, with the Biot poroelasticity model. Analyzing the effect of Biot poroelasticity on rainfall-induced landslides is the main target of this study. This research addresses how to use COMSOL Multiphysics for modeling coupled HM processes in variably saturated soils. This is done by modifying the equations in the COMSOL interface to reach the expected functionality. The COMSOL implementation is based on coupling the “Richards’ equation” and “Solid Mechanics” modules (existing in COMSOL) and adapting both modules from the user interface. The modification process is outlined in the following stages, with more details about this implementation documented in Appendix C.

COMSOL uses pressure as the primary variable. Therefore, as the first step, we reformulated the governing equations of the four approaches (i. e. FC-HM, CSGD-HM, PDC-HM, and UNC-HM) to use water pressure as a primary variable. This is done by applying the definition of specific moisture capacity as follows:

$$ C_{moist} = n \rho_{w} g \frac{\partial S_{w}}{\partial p_{w}} \qquad (32) $$

where $C_{moist}$ [$m^{-1}$] is the specific moisture capacity.

By substituting Eq. (32) in Eq. (15), and refereeing to Eq. (14), we obtain the following equation for the FC-HM approach:

$$ \left[ (\alpha - n) C_s S_w^2 + n C_w S_w + \frac{C_{moist}}{\rho_w g} \left( 1 + \frac{(\alpha - n) C_s S_w p_w}{n} \right) \right] \frac{\partial p_w}{\partial t} + \alpha S_w \frac{\partial e_v}{\partial t} + \nabla \cdot \mathbf{q}_w = 0 \qquad (33) $$

In the same manner, Eq. (26) for CSGD-HM and Eq. (27) for PDC-HM and UNC-HM approaches are reformulated, respectively, as follows:

$$ \left( n C_w S_w + \frac{C_{moist}}{\rho_w g} \right) \frac{\partial p_w}{\partial t} + \nabla \cdot \mathbf{q}_w + S_w \frac{\partial e_v}{\partial t} = 0 \qquad (34) $$

$$ \frac{C_{moist}}{\rho_w g} \frac{\partial p_w}{\partial t} + \nabla \cdot \mathbf{q}_w = 0 \qquad (35) $$

The flow equations for the four approaches as well as the corresponding assumptions are described in Fig. 1.

#### 3.1. Implementation of the FC-HM approach in COMSOL

We began by selecting the “Richards’ equation” and “Solid Mechanics” COMSOL modules individually. We then adapted the coefficients of the governing equations to be like the FC-HM approach.

RE in COMSOL is given as follows:

$$ \rho_w \left[ S_e S_p + \frac{C_{moist}}{\rho_w g} \right] \frac{\partial p_w}{\partial t} + \nabla \cdot (\rho_w \mathbf{q}_w) = Q_m \qquad (36) $$

where $S_p$ [$kg^{-1}.m^{-1}.s^2$] is the specific storage model and $Q_m$ [$kg.m^{-3}.s^{-1}$] is the mass source term.

The first stage involves modifying the storage coefficient in COMSOL. This is done with the help of variable settings and user defined selection in COMSOL:

$$ S_p = \frac{\left[ (\alpha - n) C_s S_w^2 + n C_w S_w + \frac{C_{moist}}{\rho_w g} \frac{(\alpha - n) C_s S_w p_w}{n} \right]}{S_e} \qquad (37) $$

The second stage involves adding the volumetric strain term, which describes the impact of deformation on flow processes. This is done through mass source functionality in COMSOL by adding the following term:

$$ Q_m = - \rho_w \alpha S_w \frac{\partial e_v}{\partial t} \qquad (38) $$
In COMSOL, the main equation of the “Solid Mechanical” module is given:

$$ \nabla \cdot \boldsymbol{\sigma} + \boldsymbol{F}_v = 0 \qquad (39) $$

$\boldsymbol{\sigma}$ [kg.m$^{-1}$.s$^{-2}$] is the total stress tensor which is expressed as a function of the strain tensor using Hooke's law, $\boldsymbol{F}_v$ [kg.m$^{-2}$.s$^{-2}$] is the source for total stress.

The third stage involves adapting Eq. (39) to get a form similar to Eq. (19) by adding the gravity term. This is achieved using the existing body load functionality in the COMSOL. Thus, the source term in COMSOL is defined as follows:

$$ \boldsymbol{F}_v = [\rho_s(1-n) + n\rho_w S_w]\boldsymbol{g} \qquad (40) $$

We then include the pore pressure (i.e. if the condition is saturated) or suction stress (i.e. if the condition is unsaturated) term as in Eq. (20). This is done via existing external stress functionality in COMSOL. When the external stress is included, the mechanical equation in COMSOL becomes:

$$ \nabla \cdot (\boldsymbol{\sigma}_{\text{eff}} + \boldsymbol{\sigma}_{\text{ext}}) + \boldsymbol{F}_v = 0 \qquad (41) $$

To obtain the same form as Eq. (20), the external stress tensor can be written as follows:

$$ \boldsymbol{\sigma}_{\text{ext}} = -\alpha\chi(S_w)p_w \boldsymbol{I} \qquad (42) $$

## 3.2. Implementation of the CSGD-HM approach in COMSOL

The implementation of the CSGD-HM approach in COMSOL follows the same procedure as for FC-HM (Wu et al., 2018; Wu et al., 2017; Bandara et al., 2016). All the steps remain the same, except the definition of the storage model. Given the assumption of incompressible solid grain (i.e. $C_s = 0$), the terms $(\alpha - n)C_sS_w^2$ and $(\alpha - n)C_sS_w p_w$ cancel out. The storage term implementation, described above, should be defined as follows:

$$ S_p = \frac{nC_wS_w}{S_e} \qquad (43) $$

The Biot coefficient should also be dropped from the governing equations of the mass source term and the external stress tensor by considering it equal to one.

## 3.3. Implementation of the PDC-HM and UNC-HM approaches in COMSOL

In the PDC-HM approach, neglecting derivative of the volumetric strain term made the implementation in COMSOL more practical. All the steps remain the same as in the CSGD-HM model, except the storage and mass source term stage. To avoid convergence issues, a small value is assigned for the storage coefficient ($S_p = 10^{-6}$). The mass source term (volumetric strain term) is not added due to the uncoupled assumption. The remaining stages are the same as in the development of CSGD-HM. The implementation of the UNC-HM approach is similar to PDC-HM. The external stress is excluded (see Eq. (41)) and the LFS is calculated as in Eqs. (23)–(25) for PDC and Eqs. (29)–(31) for UNC.

## 4. Verification of the COMSOL models

Four COMSOL models are developed for the FC-HM, CSGD-HM, PDC-HM, and UNC-HM approaches. These models are created by modifying the “Richards’ equation” and “Solid Mechanics” modules in COMSOL. Although these modules are well-established and validated, the correctness of the four developed COMSOL models should be investigated, as several modifications have been performed through the user interface. Typically, verification of numerical models can be done through comparison with analytical solutions of the governing equations.

<table><caption>**Table 1**

Mechanical and hydrodynamic parameters for the Liakopoulos test.</caption><thead><tr><th>Symbol</th><th>Parameter</th><th>Units</th><th>Value</th></tr></thead><tbody><tr><td>E</td><td>Young's modulus</td><td>MPa</td><td>1.3</td></tr><tr><td>ν</td><td>Poisson's ratio</td><td>–</td><td>0.4</td></tr><tr><td>n</td><td>Initial porosity</td><td>–</td><td>0.2975</td></tr><tr><td>ρ<sub>s</sub></td><td>Density of solid grains</td><td>kg/m<sup>3</sup></td><td>2000</td></tr><tr><td>C<sub>s</sub></td><td>Compressibility of solid grains</td><td>1/MPa</td><td>10<sup>-6</sup></td></tr><tr><td>K<sub>p</sub></td><td>Intrinsic permeability</td><td>m<sup>2</sup></td><td>4.5 × 10<sup>-13</sup></td></tr><tr><td>μ</td><td>Dynamic viscosity of water</td><td>Pa·s</td><td>10<sup>-3</sup></td></tr><tr><td>ρ<sub>w</sub></td><td>Density of water</td><td>kg/m<sup>3</sup></td><td>1000</td></tr><tr><td>C<sub>w</sub></td><td>Compressibility of water</td><td>1/MPa</td><td>5 × 10<sup>-4</sup></td></tr></tbody></table>

However, analytical solutions are not available for the FC-HM equations (Bandara et al., 2016). In this case, validation can be performed by comparing different numerical solutions or by comparing them with previous studies. Therefore, this section investigates the correctness of the four developed COMSOL models (FC-HM, CSGD-HM, PDC-HM and UNC-HM) by comparing them to the results of the open-source platform OpenGeoSys (OGS) (Virupaksha et al., 2024; Bilke et al., 2019) and with previous studies (Abbasov et al., 2024; Moradi et al., 2018).

### 4.1. The verification of the FC-HM and CSGD-HM: Liakopoulos experiment

This test case is used to investigate the correctness of the FC-HM COMSOL model by comparing its results to OGS. It is based on the experiment by Liakopoulos (1964), which is commonly used in the literature for verifying single and multiphase flow processes in deformable porous media (Grunwald et al., 2022; Bandara et al., 2016; Asadi and Ataie-Ashtiani, 2015; Schrefler and Scotta, 2001). The Liakopoulos experiment is useful for the validation of fully coupled consolidation models because it is sensitive to the coupled HM processes. In this experiment, a column made in Perspex is filled with Del Monte sand and equipped with a tensiometer to measure the pore pressure along the sand column. Water is allowed to infiltrate from the top of the column, while drainage occurs from the bottom. The procedure continues until a steady-state regime is established. This step is important to identify the initial conditions of the experiment. Once the pore water pressure gradient reaches zero, water infiltration from the top is stopped, and drainage of water from the column occurs due to gravity. During the desaturation period, pore water pressure and the water outflow rate of the system are recorded at different time outputs. Certain mechanical properties were not provided in the Liakopoulos experiment (Table 1) and were determined through researchers' back estimations (Abdullah et al., 2023; Grunwald et al., 2022; Bandara et al., 2016). Additionally, researchers obtained the remaining outcomes, such as air pressure, displacement, capillary pressure, and water saturation with the use of these estimated parameters. In the Liakopoulos experiment, the soil water retention model is given via the following constitutive relations:

$$ S_w = 1 - 1.97220 \times 10^{-11} \left( \frac{p_{\text{cap}}}{[\text{Pa}]} \right)^{2.4279} \qquad (44) $$

$$ k_{\text{rel}} = 1 - 2.207(1 - S_w)^{1.0121} \qquad (45) $$

where $p_{\text{cap}}$ [kg.m$^{-1}$.s$^{-2}$] is the capillary pressure.

Grunwald et al. (2022) documented Liakopoulos's experiment for single-phase fluid flow in deformable porous media, by using Richards' approximation. The benchmark model is available in the open-source platform OpenGeoSys. In this work, we compare the results of the FC-HM model developed with COMSOL to that of OGS. We use the same configuration as in OGS. Thus, the sand-filled column extends at a height of 1 m and a width of 0.1 m. Boundary conditions specific to the study follow fixed horizontal displacement for lateral, fixed vertical, and
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** The 2D profiles of saturation of the Liakopoulos example at 0, 5, 10, 20, 40, 80, and 120 min: Results obtained with the FC-HM approach.</figcaption></figure>

horizontal displacements for the bottom, and traction-free for the top boundaries. For fluid flow, lateral boundaries are impermeable, the top boundary is exposed to the atmosphere, and the bottom side allows free drainage of water. The pore water pressure is equal to zero for the initial condition and at the bottom surface (Grunwald et al., 2022; Bandara et al., 2016; Asadi and Ataie-Ashtiani, 2015).

While some studies use quadrilateral mesh with fewer elements for Liakopoulos’s experiment (Grunwald et al., 2022), we have discretized the domain with 608 triangular elements as in (Abdullah et al., 2022). The simulations were performed to study the system over 2 h with an adaptive time-stepping strategy. The simulation results of the

Liakopoulos experiment well characterize the desaturation process in a physical sense. Fig. 2 depicts the spatial distribution of the water saturation at different time outputs, highlighting the water drainage from the soil. The results of COMSOL and OGS are compared in Fig. 3. The COMSOL results are plotted with continuous lines while the results of OGS are plotted in dashed lines. Fig. 3a shows the vertical profiles of the water saturation at different time outputs. The figure shows that water saturation is decreasing with time at the top part of the column. Fig. 3b represents the vertical displacement of the soil. It shows that the invasion of pore space by the non-wetting phase leads to the consolidation of the soil over time. The results of the capillary pressure are given in Fig. 3c. The capillary pressure increases with time. The results of COMSOL and OGS match quite well. This confirms the correctness of the developed COMSOL model for the FC-HM approach. Fig. 3d shows the comparison between the COMSOL solutions and the experimental results. In this figure, the comparison is limited to time 20, 40, 80 and 120 where perfect agreement can be observed between numerical and experimental results. At the early observation times (5 and 10 min), discrepancies between the numerical and experimental results are noted, however, such differences are consistent with results reported by other studies (Abdullah et al., 2022).

Liakopoulos’ experiment is also used to verify the correctness of the COMSOL model developed for the CSGD-HM approach. The results are depicted in Fig. 3 which shows high similarity between the results of the CSGD-HM and FC-HM models. This is coherent with the results presented in Bandara et al. (2016).

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Comparison between the results of COMSOL (continuous lines) and OGS (dashed lines) for the FC-HM approach. The figure shows also the results of COMSOL with the CSGD-HM approach (gray dashed lines), and FC-HM-OGS: (a) the vertical profile of the water saturation, (b) vertical displacement of the soil, (c) the capillary pressure, and (d) the comparison with experimental result.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 4. The RILS benchmark: Geometry of the domain, boundary conditions, and mesh structure.</figcaption></figure>

<table><caption>**Table 2**
Hydrodynamic and mechanical parameters for the soil layers.</caption><thead><tr><th>Symbol</th><th>Parameter</th><th>Units</th><th>Upper layer</th><th>Lower layer</th></tr></thead><tbody><tr><td>θ<sub>s</sub></td><td>Saturated water content</td><td>–</td><td>0.412</td><td>0.385</td></tr><tr><td>θ<sub>r</sub></td><td>Residual water content</td><td>–</td><td>0.041</td><td>0.09</td></tr><tr><td>K<sub>s</sub></td><td>Saturated hydraulic conductivity</td><td>cm/h</td><td>2.59</td><td>0.06</td></tr><tr><td>ρ<sub>w</sub></td><td>Density of water</td><td>kg/m<sup>3</sup></td><td></td><td>1000</td></tr><tr><td>C<sub>w</sub></td><td>Compressibility of water</td><td>1/MPa</td><td></td><td>5 × 10<sup>-4</sup></td></tr><tr><td>α<sub>BC</sub></td><td>Brooks-Corey fitting parameters</td><td>1/cm</td><td>0.068</td><td>0.027</td></tr><tr><td>n<sub>BC</sub></td><td></td><td>–</td><td>0.322</td><td>0.131</td></tr><tr><td>l<sub>BC</sub></td><td></td><td>–</td><td></td><td>1</td></tr><tr><td>ρ<sub>s</sub></td><td>Density of solid grains</td><td>kg/m<sup>3</sup></td><td></td><td>2636</td></tr><tr><td>α</td><td>Biot coefficient</td><td>–</td><td></td><td>1</td></tr><tr><td>C<sub>s</sub></td><td>Compressibility of solid grains</td><td>1/MPa</td><td></td><td>10<sup>-6</sup></td></tr><tr><td>E</td><td>Young's modulus</td><td>MPa</td><td></td><td>10</td></tr><tr><td>ν</td><td>Poisson's ratio</td><td>–</td><td></td><td>0.35</td></tr><tr><td>φ'</td><td>Friction angle</td><td>°</td><td></td><td>35</td></tr><tr><td>c'</td><td>Effective cohesion</td><td>kPa</td><td>3</td><td>6</td></tr></tbody></table>

## 4.2. The verification of the PDC-HM and UNC-HM models: the benchmark of RILS

Liakopoulos' experiment cannot be used for verifying the PDC-HM and UNC-HM approaches because these approaches neglect the impact of deformation on flow processes. Thus, we use the common hypothetical benchmark of RILS to validate these models by comparing the results to previous studies. The benchmark is widely investigated in the literature (Shao et al., 2015; Moradi et al., 2018, and Abbasov et al., 2024). It deals with a domain characterized by a two-layered soil, with an upper layer exhibiting sandy loam with high permeability and a lower clay layer of low permeability (Fig. 4). The topsoil surface is initially stress-free with imposed rainfall flux (i.e. when the top surface becomes fully saturated, the corresponding atmospheric boundary acts as a seepage boundary condition). The hillslope is inclined at 23 degrees. A preliminary simulation is performed over 10 years with an average infiltration of 600 mm/year. This simulation is used to determine the initial state of the aquifer. The stability of the slope is then simulated with high-intensive rainfall events of 20 mm/h over 15 h. The infiltration over the inclined surface is assumed to be less than that over the horizontal surfaces. It is reduced by multiplying it by the cosine slope angle. The bottom and left boundaries are impermeable. Also, roller

<figure><img src="image_7.png" /><figcaption>Fig. 5. Maps of the volumetric water content (top) and LFS (bottom) at the last time of simulations (i.e. 15 h). Results obtained with the UNC-HM and PDC-HM approaches.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Maps of the water content and LFS at the last time of the simulation (i.e. 15 h). Results of the hypothetical benchmark with different HM approaches.</figcaption></figure>

boundary conditions are imposed on these sides in conjunction with the right boundary. Hydrostatic pressure at the lower right part and seepage boundary conditions at the top right side are applied. We referred to the same parameters used in Shao et al. (2015) in the simulations of this benchmark (Table 2). The water table position initially is at $z = 17$ m. The domain is discretized by applying triangular elements and locally refining above $z = 17$ m (Fig. 4).

The results of the volumetric water content and Local Factor of Safety (LFS) are given in Fig. 5. The water content maps show the infiltration of the water over the top layer of soil. The presence of a dry zone at the interface between the top and bottom layers of soil, and of a fully saturated zone at the bottom part of the inclined surface indicates a preferential path of flow in this top layer due to its high permeability. Both the PDC-HM and UNC-HM approaches give the same results regarding the water content, which is expected because the governing equations of the flow model are the same in both approaches. The results of the LFS show a high risk of landslide over a thin layer on the top surface of the soil. The high zone of risk is at the lower part of the inclined surface where the soil is highly saturated. The results of the UNC-HM COMSOL model are in full agreement with the results of Shao et al. (2015), Moradi et al. (2018), and Abbasov et al. (2024). The models used in those papers are based on the same UNC-HM approach. However, there are no available results of the PDC-HM approach for this benchmark. The correctness of the PDC-HM COMSOL model regarding LFS is validated by comparing its results to the UNC-HM approach. Fig. 5 shows a slight difference between the results of these models. This

confirms the correctness of the PDC-HM COMSOL model.

## 5. Model's complexity versus reliability

As discussed previously, in the context of RILS, several assumptions have been applied to the HM models in order to improve the computational performance of the simulations and to allow for more realistic studies at field scale. These assumptions lead to the four approaches presented above (FC-HM, CSGD-HM, PDC-HM, and UNC-HM). While these assumptions may reduce the complexity of the mathematical models, they can also impact the reliability of the models' predictions. This section aims to investigate the impact of the different simplifications on the complexity and accuracy of the HM models in the case of RILS. This is done by comparing the performance (i.e. CPU time) and the results of the four approaches FC-HM, CSGD-HM, PDC-HM, and UNC-HM in two test cases dealing respectively with the hypothetical benchmark of RILS and an application to a real case in Azerbaijan. The result of the FC-HM approach is used as a reference because it does not involve any simplification. The complexity of the models is measured using the required CPU time, time step size, and number of time steps required in the simulation. The metrics used to evaluate the accuracy are the maps of water content and LFS. More reliable quantities of interest are needed to enhance the clarity and effectiveness of the comparative analysis. Thus, we use the area of the fully saturated zone and the area of the failure zone in the top layer of the soil.

The area of the fully saturated zone is defined as follows:
<figure><img src="image_9.png" /><figcaption>Fig. 7. Plots of the area of the fully saturated zone ($A^{FSZ}$) predicted by the different HM approaches (left) and the discrepancy between the simplified approaches and the FC-HM (right) over 15 h of an intensive rainfall event.</figcaption></figure>

$$A^{FSZ} = \iint_{\Omega} S_{w}^{*} \cdot d\Omega \qquad (46)$$

where $A^{FSZ}$ is the area of the fully saturated zone, $\Omega$ is the domain corresponding to the top layer of the soil and $S_{w}^{*}$ is given as:

$$S_{w}^{*} = \begin{cases} 1 & \text{if } S_{w} = 1 \\ 0 & \text{if } S_{w} < 1 \end{cases} \qquad (47)$$

The area of the failure prone zone corresponds to the area of the domain where LFS is less than one. This metric is relevant to indicate the

landslide susceptibility (Abbasov et al., 2024). It is calculated as follows:

$$A^{LFS} = \iint_{\Omega} LFS^{*} \cdot d\Omega \qquad (48)$$

where $A^{LFS}$ is the area of failure prone zone, $\Omega$ is the domain corresponding to the top layer of the soil and $LFS^{*}$ is given as:

$$LFS^{*} = \begin{cases} 1 & \text{if } LFS < 1 \\ 0 & \text{if } LFS \ge 1 \end{cases} \qquad (49)$$

<figure><img src="image_10.png" /><figcaption>Fig. 8. Time variation of the failure prone zone (a), maximum displacement magnitude (b), maximum von Mises stress (c) and maximum volumetric strain (d). The results are obtained by the four models over 15 h of an intensive rainfall event.</figcaption></figure>
## 5.1. The RILS benchmark

In this section, we compare the different HM approaches using a hypothetical benchmark of RILS. This benchmark, described in the previous section, offers deep insights into the performance and reliability of the models due to its relative simplicity and controlled conditions. The benchmark is simulated here using the four approaches: FC-HM, CSGD-HM, PDC-HM, and UNC-HM. The parameters used in the simulations are the same as in Table 2. The boundary conditions, initial conditions, and rainfall configuration are the same as in the previous section. A mesh sensitivity analysis has been performed to avoid any numerical artifacts. The mesh-independent solutions can be obtained with a computational triangular grid of about 3 K nodes. We also investigate the sensitivity of the solutions to the time step size. The results show that the solutions with the four HM modes are highly sensitive to the time step size. This sensitivity to the time step size is associated with the Richards' equation (Fahs et al., 2009). In COMSOL, an automatic time-stepping scheme is employed. At the start of the simulation, COMSOL uses small time steps to ensure the convergence of the nonlinear solver. The time step is then automatically increased through a prediction-correction process that is based on an estimation of the truncation error. However, as the simulation progresses, this automatic adjustment can result in large time steps, which may introduce numerical diffusion because the time integration is based on an implicit scheme. To mitigate excessive numerical diffusion, an upper limit on the step size time must be imposed. In this study, the maximum time step is set to 0.1 h, based on empirical testing through several numerical experiments aimed at achieving a time-step-independent solution.

The resulting maps of water content and LFS at the end of the simulation time (i.e. end of the intensive rainfall period at $t = 15$ h) are given in Fig. 6. The maps of water content show similar results across all approaches. This indicates that the simplifications in the different approaches are not impacting the unsaturated flow processes. The fact that the three simplified approaches provide similar results as FC-HM confirms the reliability of the simplifications applied to the unsaturated flow equation. As discussed previously, these simplifications assume that solid grain density is constant (CSGD-HM) and the impact of mechanical deformation on flow processes (i.e. volumetric strain term) is negligible (PDC-HM and UNC-HM).

The LFS results reveal some discrepancies between the approaches. As shown in Fig. 6, the results of the FC-HM, CSGD-HM, and PDC-HM approaches are quite similar. This means that the assumption of constant solid grain density and neglecting the volumetric strain term do not impact the prediction of the LFS maps. However, the UNC-HM approach stands out because, unlike the other models, UNC-HM underestimates the vulnerable zone. Thus, the common approach used in the literature, that neglects the pore pressure in the mechanical equations and incorporates it directly into the LFS evaluation, leads to an underestimation of the risk of landslide. It should be also noted that the UNC-HM approach, contrary to other approaches, provides LFS values close to one in the bottom part of the domain. This result is incoherent with the risk of landslide analysis.

The previous analysis, based on the maps, provides an overview of the comparison between approaches because it is qualitative, and is limited to a last time step. For deeper analysis, we discuss here the time variation of the area of the fully saturated zone ($A^{FSZ}$) and the area of the failure prone zone ($A^{LFS}$). The results are depicted in Figs. 7 and 8.

Fig. 7 shows that all approaches provide close results regarding $A^{FSZ}$. This result is coherent with that observed for the map of water content (Fig. 6). However, slight discrepancy can be observed between the approaches. CSGD-HM seems to have better accuracy than PDC-HM and UNC-HM (Fig. 7). Fig. 7 shows two notable peaks for $A^{FSZ}$ within the first six hours of the simulation. Beyond the seventh hour, the area of the fully saturated zone increases gradually. These peaks are related to infiltration and seepage processes. The infiltration processes increase

<caption>**Table 3**</caption>

The CPU time and number of time steps used with the four HM approaches for the simulation of the 15 h of intensive rainfall period. Results are presented for two configurations of time step management: i) the time step size is freely selected by the solver (left) and ii) a maximum time step size of 0.1 h is imposed (right).

<table><thead><tr><th rowspan="2">Approach</th><th colspan="2">Free maximum time step</th><th colspan="2">Maximum time step (0.1 h)</th></tr><tr><th>CPU (min)</th><th>Nb. of time steps</th><th>CPU (min)</th><th>Nb. of time steps</th></tr></thead><tbody><tr><td>FC-HM</td><td>2.3</td><td>171</td><td>12.3</td><td>942</td></tr><tr><td>CSGD-HM</td><td>4.3</td><td>332</td><td>12.6</td><td>1073</td></tr><tr><td>PDC-HM</td><td>4.5</td><td>322</td><td>19.0</td><td>1643</td></tr><tr><td>UNC-HM</td><td>3.0</td><td>272</td><td>18.1</td><td>1686</td></tr></tbody></table>

$A^{FSZ}$, while seepage decreases it. The configuration under investigation deals with the high intensity of the rainfall which is close to the hydraulic conductivity of the top layer of soil. For this reason, infiltration and seepage are alternating in short periods at the beginning. To confirm that these peaks are attributed to the rainfall intensity and soil hydraulic conductivity, we perform additional simulations with lower rainfall intensities (i.e., 16, 12, and 8 mm/h), and no peaks were recorded in these cases. After seven hours of rainfall, the discrepancy between the FC-HM and other approaches becomes more important. The fact that CSGD-HM, PDC-HM, and UNC-HM provide similar results, and these results are slightly different from FC-HM, indicates that the discrepancy between approaches is related to the assumption of constant solid grain density.

Results of the area of the failure-prone zone ($A^{LFS}$) show perfect matching between the FC-HM, CSGD-HM, and PDC-HM approaches (Fig. 8a), before 9 h of rainfall. However, afterward, CSGD-HM and PDC-HM overestimate $A^{LFS}$, compared to the reference approach. A significant discrepancy is observed between UNC-HM and other HM approaches. This indicates that neglecting the pore pressure in the mechanical equations has a significant impact on the results. The results of CSGD-HM and PDC-HM are quite similar, indicating that neglecting the volumetric strain term in the flow equations has no significant impact on the stability predictions. Fig. 8b, c and d show the variation of the maximum displacement, Von-Mises stress and volumetric strain with time. These figure show that the models exhibit similar comparative performance to $A^{LFS}$.

As discussed previously, the performance of the different approaches is assessed using the CPU required in the corresponding numerical simulations. To understand the discrepancy between the approaches in terms of CPU time, we also investigate the time step size and the number of time steps used in COMSOL. Table 3 summarizes the CPU time and number of time steps required for simulating the intensive rainfall event with the four HM approaches. The simulations are performed on a computer system with an Intel@544CoreTM i7-6700 CPU at 3.4 GHz and with 16 GB of RAM. The results are presented for two different strategies of time step management. In the first strategy, the time step size is freely selected by COMSOL without any constrain (i.e. default setting). In this case, the time step size is calculated based on an estimation of the error introduced by the time integration scheme. The time step size is evaluated to satisfy the prescribed accuracy. However, as discussed in the previous section, this technique does not provide accurate results because it leads to large time steps for which the error cannot be controlled. With the second strategy, we impose a maximum time step of 0.1 h in order to obtain accurate solutions. Thus, the time step size is updated as with the first strategy, but its value cannot exceed 0.1 h.

The FC-HM model is expected to be less efficient than the other models. However, the results in Table 3 show that, with the first strategy (i.e. free maximum time step), the CPU times are equivalent for all approaches. There is a slight advantage of FC-HM and UNC-HM, and this is because these approaches require less time steps than other approaches. The time step size used during the simulation is represented in Fig. 9a. The general behavior of the approaches is that the time step size is
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Time step size used in the simulations for all HM approaches over the intensive rainfall period. Results are presented for two strategies of time step management: a) the time step size is freely selected by the solver and b) a maximum time step size of 0.1 h is imposed (right).</figcaption></figure>

decreasing at the beginning of the simulation. This is related to the high rainfall intensity that impacts the conditioning of the nonlinear system. The solver should take small time steps to improve the conditioning of the system and to ensure convergence. After that, when the change in the solution is slow, the time step size start increasing to reach a high value where it stays constant. At the last part of the simulation, the time step size start decreasing again. This is because the water reaches the outlet of the top layer, where a seepage boundary condition is imposed. This change in the boundary condition affects the conditioning of the system and lead to decreased time step. The FC-HM approach, due its high nonlinearity, requires smaller time steps than other approaches at the beginning of the simulation. The maximum time step used with this approach is also smaller than other approaches. However, in the latter part of the simulation, the time steps used with the FC-HM are larger than other approaches. This due to the fact that the additional terms involved in the equations of this approach improve the conditioning of the system and avoid singularity issues. Due to this behavior, FC-HM require less time steps than other approaches as shown in Table 3. UNC-HM does not require very small-time steps at the beginning because the flow and mechanical processes are decoupled. For this reason, this approach outperforms CSGD-HM and PDC-HM approaches.

In the second configuration, when the maximum time step size is set

to 0.1 h, FC-HM is the most efficient approach, followed by CSGD-HM. The time steps used with all approaches are given in Fig. 9b. This figure shows that all approaches have the same behavior during the simulation, except at the end when the seepage processes start. The decrease in the time step size starts earlier with PDC-HM and UNC-HM than FC-HM and CSGD-HM. This explains why the number of time steps is higher with PDC-HM and UNC-HM than FC-HM and CSGD-HM. All approaches show a significant improvement in both CPU time and the number of time steps when using the adaptive time step strategy, with FC-HM consistently outperforms other approaches. This analysis highlights the effectiveness of adaptive time steps in improving computational efficiency across different approaches.

The results presented above demonstrate that FC-HM outperforms other approaches in terms of accuracy and efficiency. However, these findings are specific to the given values of parameters. The performance of the four models could be impacted by the parameter's values. Therefore, for a more generic analysis, it is important to compare the approaches under different sets of parameter values. We performed preliminary runs and found that the convergence and performance of the numerical models are highly impacted by the intensity of rainfall. Hence, in this section, we compare the models under different values of rainfall intensity (i.e., 16, 12, and 8 mm/h). As in the previous section,

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Area of the failure-prone zone ($A^{\text{LFS}}[m^2]$) at the last time of simulations (i.e. 15 h) for the different rainfall intensities (i.e. 20, 16, 12, and 8 mm/h).</figcaption></figure>
<table><caption>**Table 4**
The CPU time (in min) of the four HM approaches under the different rainfall intensities (i.e. 20, 16, 12, and 8 mm/h).</caption><thead><tr><th colspan="5">Maximum time step (0.1 h)</th></tr><tr><th></th><th>FC-HM</th><th>CSGD-HM</th><th>PDC-HM</th><th>UNC-HM</th></tr></thead><tbody><tr><td>20 mm</td><td>12.3</td><td>12.6</td><td>19.0</td><td>18.1</td></tr><tr><td>16 mm</td><td>3.8</td><td>8.0</td><td>6.5</td><td>7.5</td></tr><tr><td>12 mm</td><td>2.2</td><td>2.1</td><td>2.3</td><td>1.4</td></tr><tr><td>8 mm</td><td>2.0</td><td>2.3</td><td>2.1</td><td>1.2</td></tr><tr><td colspan="5">Free maximum time step</td></tr><tr><td>20 mm</td><td>2.3</td><td>4.3</td><td>4.5</td><td>3.0</td></tr><tr><td>16 mm</td><td>1.4</td><td>1.6</td><td>1.8</td><td>1.4</td></tr><tr><td>12 mm</td><td>1.0</td><td>0.9</td><td>1.0</td><td>0.6</td></tr><tr><td>8 mm</td><td>0.7</td><td>0.7</td><td>0.7</td><td>0.4</td></tr></tbody></table>

we perform runs with both strategies of time stepping (i.e. free maximum and imposed maximum).

Fig. 10 shows bar graphs depicting $A^{LFS}$ at the last time of simulations (i.e., 15 h) for four HM approaches under different rainfall intensities (i.e., 16, 12, and 8 mm/h). In this figure, all simulations are carried out using the same mesh. As previously described in the paper, the mesh is selected based on a mesh sensitivity analysis to ensure mesh-independent solutions. All soil and fluid parameters are kept identical across simulations, and the time-stepping error is controlled by imposing a maximum time step size. The only varying parameter between simulations is the rainfall intensity. The results indicate that the CSGD-HM and PDC-HM approaches achieve the same level of accuracy as the FC-HM model. In contrast, the UNC-HM model tends to either

overestimate or underestimate the landslide risk depending on the intensity of rainfall. As shown in Fig. 10, when this intensity is high, UNC-HM underestimates the risk of sliding while conversely overestimates it for low rainfall intensity. Since all models except UNC-HM yield consistent results and as numerical errors are controlled and all parameters are the same except rainfall intensity, the observed discrepancies can be attributed to the decoupling assumption inherent in the UNC-HM model. The efficiency of the four approaches is assessed based on their CPU time. The results are summarized in Table 4. This table shows that for all models, the CPU time decreases as intensity of rainfall decreases. The FC-HM outperforms other approaches when the rainfall intensity is high.

## 5.2. The case of Badamdar Lanslide

The Badamdar landslide is located in the Badamdar municipality area, near the road of the Badamdar-20th zone within the Sabail district of Baku in Azerbaijan (Fig. 11). The sliding occurs on the eastern wing of the Bibiheibat brachianticline, an area of high risk. The first landslide initiation was recorded in December 2012 (SCLCRA (State Committee for Land and Cartography of the Republic of Azerbaijan), 2025). Until now, several sliding events occurred in the zone. The boundaries of sliding zones at different times are well documented in previous studies (Salamov et al., 2020). Intensive rainfall events and infiltration of sewage water are considered as crucial factors in the slope destabilization process. In addition, shallow surface waters triggered the formation of ravines at the nearby sliding zone. Several management techniques have been implemented to limit the risks. Mitigation techniques cover the construction of drainage wells to collect sewage water and retaining walls to make the barrier to debris flow in the territory (Fig. 12). To the

<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Map of the Badamdar landslide area. The red polygon depicts the boundary of a sliding zone (after SCLCRA (State Committee for Land and Cartography of the Republic of Azerbaijan), 2025). The blue and green lines with arrows indicate the boundary of sliding in 2019 and 2014, respectively. The black line shows a cross-section of the topography. The lithological descriptions of all layers are provided in the map (after Salamov et al., 2020, see within for landslide elements and geological cross-section of the study area). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Representation of sliding zone and applied mitigation techniques to reduce the risk of possible landslide occurrence in the Badamdar municipality area of the Sabail district of Baku in Azerbaijan (photographs by authors).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 13.** The Badamdar field case: Geometry of the domain, boundary conditions, and schematic of structure.</figcaption></figure>

best of our knowledge, no HM model has been developed specifically for investigating landslide in this specific area. In this study, we use the different HM models discussed above to investigate this real case. The main goal is to compare these models under realistic conditions.

A vertical 2D section has been selected to perform the simulations with the different models (Fig. 11). Geological and geophysical investigations in the territory have revealed that the eluvial-deluvial deposits of the Quaternary period were deposited in the shallow surface of the study area. Salamov et al. (2020) conducted comprehensive geophysical investigations using the vertical electrical sounding method. Through this approach, they achieved a more precise delineation of the lithological layers. The eluvial-deluvial sediments are deposited on the top surface, followed by the clay, clay with thin sand layers, sandy clay, and clayey sand layers toward the depth (Fig. 11).

<caption>**Table 5**</caption>

Hydrodynamic and mechanical parameters for the soil layers for the case of the Badamdar area.

<table><thead><tr><th>Symbol</th><th>Parameter</th><th>Units</th><th>1st layer</th><th>2nd layer</th><th>3rd & 4th & 5th layers</th></tr></thead><tbody><tr><td>θ<sub>s</sub></td><td>Saturated water content</td><td>–</td><td>0.4</td><td>0.385</td><td>0.39</td></tr><tr><td>θ<sub>r</sub></td><td>Residual water content</td><td>–</td><td>0.06</td><td>0.075</td><td>0.1</td></tr><tr><td>K<sub>S</sub></td><td>Saturated hydraulic conductivity</td><td>cm/h</td><td>2</td><td>0.4</td><td>0.35</td></tr><tr><td>ρ<sub>w</sub></td><td>Density of water</td><td>kg/m<sup>3</sup></td><td></td><td>1000</td><td></td></tr><tr><td>C<sub>w</sub></td><td>Compressibility of water</td><td>1/MPa</td><td></td><td>5 × 10<sup>-4</sup></td><td></td></tr><tr><td>α<sub>BC</sub></td><td>van Genuchten fitting parameters</td><td>1/cm</td><td>0.008</td><td>0.006</td><td>0.01</td></tr><tr><td>n<sub>BC</sub></td><td>Density of solid grains</td><td>kg/m<sup>3</sup></td><td>2640</td><td>2680</td><td>2670</td></tr><tr><td>C<sub>s</sub></td><td>Compressibility of solid grains</td><td>1/Pa</td><td></td><td>1.2 × 10<sup>-6</sup></td><td></td></tr><tr><td>α</td><td>Biot coefficient</td><td>–</td><td></td><td>1</td><td></td></tr><tr><td>E</td><td>Young's modulus</td><td>MPa</td><td>25</td><td>8</td><td>14</td></tr><tr><td>ν</td><td>Poisson's ratio</td><td>–</td><td></td><td>0.3</td><td></td></tr><tr><td>φ'</td><td>Friction angle</td><td>°</td><td>28</td><td>22</td><td>29</td></tr><tr><td>c'</td><td>Effective cohesion</td><td>kPa</td><td>64</td><td>56</td><td>87</td></tr></tbody></table>

### 5.3. The model results for the Badamdar case

This study deals with rainfall-induced slope destabilization. Therefore, the field case deals with rainfall on a hillslope. Salamov et al. (2020) serve as a point of reference for the lithology of the cross-section of the sliding zone. The domain consists of five layers and extends approximately over 104 m (i.e. cross-section length is 95 m) and 48 m in depth with a 24-degree slope (Fig. 11, Fig. 13).

The hydro-mechanical parameters of the layers were selected precisely in coherence with the soil texture of the layers and previous investigations made by Salamov et al. (2020) in the area. All parameters are summarized in Table 5. Salamov et al. (2020) defined the properties of the first layer with geological and geophysical exploration techniques.
<figure><img src="image_16.png" /><figcaption>Fig. 14. Maps of the water content at the last time of simulation (i.e. 24 h). Results of the Badamdar field case with all four HM approaches.</figcaption></figure>

To ensure we have all the necessary properties of the layers for simulation purposes, we referred to an already existing hydro-mechanical experimental dataset in Carsel and Parrish (1988), UNSODA (UNsaturated SOil hydraulic Database) and database of Geotechdata site (Koliji, 2008; Nemes et al., 1999). We used the UNSODA viewer for convenient access to the data (Seki et al., 2023). For the first layer, data from Salamov et al. (2020) is coherent with the data gathered by Carsel and Parrish (1988). In Carsel and Parrish (1988), the experimental database is built with the use of van Genuchten's soil water retention model. Therefore, we have performed simulations with the following retention model:

$$S_e = \frac{S_w - S_r}{1 - S_r} = \begin{cases} \left[ \frac{1}{1 + \left| \frac{p_w}{p_{ent}} \right|^{n_{VG}}} \right]^{m_{VG}} & \text{if } p_w < 0 \\ 1 & \text{if } p_w \ge 0 \end{cases} \quad (50)$$

$$k_{rel} = S_e^{\ell_{VG}} \left[ 1 - \left[ 1 - \frac{1}{S_e^{m_{VG}}} \right]^{m_{VG}} \right]^{2} \quad (51)$$

where $n_{VG}[-]$ (i.e. $m_{VG} = 1 - 1/n_{VG}$) and $\ell_{VG}[-]$ characterize pore size distribution in soil and tortuosity, respectively.

The topsoil surface is initially stress-free with an imposed rainfall flux. When the topsoil surface becomes saturated, it behaves as a seepage boundary. On the inclined surface, runoff is higher, with rain flux reduced by the cosine of the slope angle. The bottom and left-top boundaries are impermeable with roller boundary conditions (no

normal displacement allowed). The right surface also has a roller boundary condition, with hydrostatic pressure at the right-lower and a seepage boundary at the right-top part. Also, hydrostatic pressure is imposed on the left-bottom boundary to keep the water level at a given inclination (Fig. 13). For the initial conditions, we assume that the water table is an inclined line between two points (0, 63) and (170, 44).

The simulations are performed for the domain given in Fig. 13. Results are assessed for the investigated area in the center as shown in Fig. 13. The domain's size exceeds the investigated area to minimize the effect of boundary conditions on the results. The simulations are performed for two rainfall regimes. The first simulation covers 5 years with an average rainfall of 220 mm/yr (i.e. Salamov et al. (2020) referred for the average rainfall amount) to achieve equilibrium conditions, and results are used as an initial condition for the second simulation which deals with intensive rainfall of 8 mm/h over 24 h. The domain is discretized by applying an extra fine mesh of around 1.8 K nodes. Adaptive time stepping is used in the simulations, without any constraint on the maximum time step size.

Fig. 14 illustrates the results of water content and LFS maps for the four HM approaches (i.e. FC-HM, CSGD-HM, PDC-HM, and UNC-HM). The distribution of the water content obtained with FC-HM and CSGD-HM are quite similar. This confirms that the assumption of constant solid grain density does not affect the water content results. As expected, the results of PDC-HM and UNC-HM are the same because both have the same governing flow equations. However, their results are different from the reference ones (FC-HM). More water storage can be observed in the soil top layer with FC-HM than PDC-HM and UNC-HM. This is due to neglecting the volumetric strain term in PDC-HM and UNC-HM
<figure><img src="image_17.png" /><figcaption>Fig. 15. Maps of the LFS at the last time of simulation (i.e. 24 h). Results of the Badamdar field case with all four HM approaches.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 16. Plots of the area of the failure-prone zone ($A^{\text{LFS}}$) predicted by the different HM approaches (left) and the discrepancy between the simplified approaches and the FC-HM (right) over 24 h of an intensive rainfall event in the Badamdar landslide case.</figcaption></figure>

approaches.

Fig. 15 shows the LFS maps obtained using the different HM approaches. The LFS results support the previous assertion that the volumetric strain term in flow equations does not significantly affect the stability predictions. UNC-HM stands out from all other HM models in the slope stability assessment. This is coherent with the results of the RILS benchmark in the previous section. The difference is because UNC-HM neglects pore pressure in the mechanical equations and incorporates

it in the equation of LFS.

The results for the failure-prone zone ($A^{\text{LFS}}$) from the FC-HM, CSGD-HM, and PDC-HM approaches show close alignment (Fig. 16). In contrast, UNC-HM shows a significant deviation from the other approaches throughout the entire 24-h period, with the discrepancy increasing as time progresses. This discrepancy indicates that the omission of pore pressure in the mechanical equations has a considerable impact on the predicted failure-prone area, particularly under
extended rainfall conditions. Moreover, the similar behavior of CSGD-HM and PDC-HM indicates that neglecting the volumetric strain term in the flow equations does not substantially affect the stability predictions in this case. However, the significant deviation of the UNC-HM approach from the others suggests that it fails to capture critical factors influencing slope stability, emphasizing the importance of considering pore pressure in the mechanical equations when modeling landslide scenarios under prolonged rainfall.

The complexity of the different approaches is measured using the CPU time required in the simulations. Equivalent CPU time is observed for FC-HM, CSGD-HM, PDC-HM. This CPU time is measured to be 2.5 min. The CPU time for UNC-HM is about 1.5 min. This result indicates that the simplifications of the mathematical models do not lead to a substantial improvement in the CPU time.

## 6. Discussions

The four models, FC-HM, CSGD-HM, PDC-HM and UNC-HM, have been implemented in COMSOL with the same numerical schemes for spatial discretization and time integration. The COMSOL models have been validated by comparing their results to experimental observations and results from other simulators. Perfect agreement has been reached, confirming the validity of the newly developed COMSOL models. All models are then used to simulate a hypothetical benchmark in order to evaluate and compare their accuracy and efficiency under controlled conditions. The results show that for the water content and area of the fully saturated zone, all models provide quite similar results. This confirms that neglecting the solid grains compressibility and impact of volumetric strain does not impact the flow processes. Concerning the risk of landslides, all models also provide similar results, except UNC-HM model. Despite this approach is widely adapted in literature, the results show that neglecting the coupled HM processes lead to overestimation of the failure risk. A careful analysis of this result shows that, depending on rainfall intensity, the UNC-HM could also underestimate the failure risk. Concerning the performance of models, comparisons are performed with both adaptive time stepping and constant time step. The results reveal that the performance of all four models was quite similar, with the FC-HM approach even outperforming the others in some cases. In fact, the analysis indicates that employing an advanced time integration technique for solving the equations can improve the computational demands of the fully coupled approach, making its performance comparable to that of the simplified approaches. The high performance of the FC-HM model is related to the fact that this approach improves the conditioning of the system and avoids singularity issues. The four models are finally evaluated using a real case study from Badamdar, Azerbaijan. The results confirm the findings obtained from the hypothetical benchmark, providing greater confidence in the conclusions under field scale and realistic conditions.

## 7. Conclusion

Richards' Eq. (RE) is widely employed to represent the flow in unsaturated soils when modeling rainfall induced landslide (RILS) with hydro-mechanical (HM) models. This work aims at providing an in-depth understanding of the application of RE with a poroelasticity approach for the simulation of RILS. While there is no comprehensive study on the mathematical derivation of the coupling between RE and mechanical equations under the Biot poroelasticity approach, this study provides the full mathematical derivation, relying on the conservation

laws and local equilibrium. Various simplifications have been proposed in the literature to reduce the complexity of the governing equations and to improve practicality. This work reviews these simplifications, along with their mathematical derivations. To the best of our knowledge, the reliability of these simplifications for RILS has never been systematically examined by confronting them with the full coupled model. This gap is addressed here by investigating the impact of various model simplifications on models' reliability and performance, based on a hypothetical RILS benchmark and a real-world application in Azerbaijan. Minor discrepancies in flow processes were observed. For stability analysis, the comparisons support the reliability of the constant solid grain density assumption and the simplification based on ignoring impact of volumetric strain on porosity. The uncoupled model, widely used in literature, could either overestimate or underestimate failure risk, depending on rainfall intensity. The results indicate that employing an advanced time integration technique for solving the equations can improve the computational demands of the fully coupled approach, making its performance comparable to that of the simplified approaches. In summary, this study demonstrates that the fully coupled model is the most suitable approach for simulating RILS because it offers the most reliable results. Although it is computationally demanding, advancements in numerical methods and computational techniques improve the efficiency of the fully coupled model and render it the most suitable approach for simulating RILS. The findings of this work indicate that HM model simplifications are not necessary for simulating rainfall-induced landslides, as the original fully coupled model provides superior accuracy and comparable computational efficiency when equations are solved with suitable numerical techniques. The present study is limited to elastic model for the description of mechanical response of soil. This model could be limited in capturing the feature of stress re-distribution in the local failure zone. A further extension of this work could be a comparative study between the different models in case of plastic deformation.

## CRediT authorship contribution statement

**Rashad Abbasov:** Writing – original draft, Validation, Software, Methodology, Investigation, Conceptualization. **Marwan Fahs:** Writing – original draft, Supervision, Methodology, Funding acquisition. **Vincent Fontaine:** Writing – review & editing, Methodology. **Husam Musa Baalousha:** Writing – review & editing, Methodology, Conceptualization. **Anis Younes:** Writing – review & editing, Software. **Renaud Toussaint:** Writing – review & editing, Supervision, Methodology, Funding acquisition.

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

## Acknowledgments

The authors would like to acknowledge the support of the state program on increasing the international competitiveness of the higher education system in the Republic of Azerbaijan for 2019-2023. RT acknowledges the support of the Research Council of Norway through its Centers of Excellence funding scheme, project number 262644, and of IRP D-FFRACT, CNRS.

## Appendix A. Evaluation of the terms of the mass conservation equation

This appendix aims at explaining how some terms in Eq. (4) can be neglected.

The term $q_w \nabla \rho_w$ in Eq. (4) is neglected in front of the term $\rho_w \nabla q_w$. This is done through the analysis of their magnitudes, based on the
corresponding characteristics lengths.

$$ \frac{\mathbf{q}_w \nabla \rho_w}{\rho_w \nabla \cdot \mathbf{q}_w} = \frac{\mathbf{q}_w \rho_{w0} C_w \nabla (p_w - p_{w0})}{\rho_w \nabla \cdot \mathbf{q}_w} \approx \frac{|\mathbf{q}_w| \rho_{w0} C_w \frac{(p_w - p_{w0})}{l}}{\rho_w \frac{|\mathbf{q}_w|}{l}} \approx C_w (p_w - p_{w0}) \quad (A1) $$

$$ C_w(p_w - p_{w0}) \approx 5 \times 10^{-10} [1/Pa] \times 10^5 [Pa] \ll 1 \quad (A2) $$

l is the characteristic length of the problem.

The term $\mathbf{v}_s \nabla [(n\rho_w S_w)]$ is neglected in front of the term $n\rho_w S_w \nabla \cdot \mathbf{v}_s$. Let's consider $\mathbf{v}_s \nabla [(n\rho_w S_w)]$ as two terms $nS_w \mathbf{v}_s \nabla \rho_w$ and $\rho_w \mathbf{v}_s \nabla (nS_w)$. We have analyzed above neglecting $\nabla \rho_w$ term. So, let's analyze $\rho_w \mathbf{v}_s \nabla (nS_w)$ term in front of $n\rho_w S_w \nabla \cdot \mathbf{v}_s$ term.

$$ \frac{\rho_w \mathbf{v}_s \nabla (n S_w)}{n \rho_w S_w \nabla \cdot \mathbf{v}_s} \approx \frac{\rho_w |\mathbf{v}_s|^{\frac{n S_w}{l}}}{n \rho_w S_w \frac{|\mathbf{v}_s|}{l}} \approx 1 \quad (A3) $$

In the given case of characteristic length $l_v \ll l_{nS_w}$, the expression will be less than one:

$$ \frac{\rho_w \mathbf{v}_s \nabla (n S_w)}{n \rho_w S_w \nabla \cdot \mathbf{v}_s} \approx \frac{\rho_w |\mathbf{v}_s|^{\frac{n S_w}{l_{nS_w}}}}{n \rho_w S_w \frac{|\mathbf{v}_s|}{l_{nS_w}}} \ll 1 \quad (A4) $$

Where $l_v$ is the characteristic length of the solid grain velocity and $l_{nS_w}$ is the characteristic length of volumetric water content.

The term $\mathbf{v}_s \nabla [(1-n)\rho_s]$ is neglected in front of $[(1-n)\rho_s] \nabla \cdot \mathbf{v}_s$ term. Let's consider $\mathbf{v}_s \nabla [(1-n)\rho_s]$ as two terms $(1-n)\mathbf{v}_s \nabla \rho_s$ and $\rho_s \mathbf{v}_s \nabla (1-n)$. First, we analyze $(1-n)\mathbf{v}_s \nabla \rho_s$ term in front of $[(1-n)\rho_s] \nabla \cdot \mathbf{v}_s$ term.

$$ \frac{(1-n)\mathbf{v}_s \nabla \rho_s}{[(1-n)\rho_s] \nabla \cdot \mathbf{v}_s} = \frac{(1-n)\mathbf{v}_s \rho_{s0} C_s \nabla (p_s - p_{s0})}{[(1-n)\rho_s] \nabla \cdot \mathbf{v}_s} \approx \frac{(1-n)|\mathbf{v}_s| \rho_{s0} C_s \frac{(p_s - p_{s0})}{l}}{[(1-n)\rho_s] \frac{|\mathbf{v}_s|}{l}} \approx C_s (p_s - p_{s0}) \quad (A5) $$

$$ C_s(p_s - p_{s0}) \approx 10^{-12} [1/Pa] \times 10^5 [Pa] \ll 1 \quad (A6) $$

Then, $\rho_s \mathbf{v}_s \nabla (1-n)$ term in front of $[(1-n)\rho_s] \nabla \cdot \mathbf{v}_s$ term.

$$ \frac{\rho_s \mathbf{v}_s \nabla (1-n)}{[(1-n)\rho_s] \nabla \cdot \mathbf{v}_s} \approx \frac{\rho_s |\mathbf{v}_s|^{\frac{(1-n)}{l}}}{[(1-n)\rho_s]^{\frac{|\mathbf{v}_s|}{l}}} \approx 1 \quad (A7) $$

In the given case of characteristic length $l_v \ll l_{(1-n)}$, the expression will be less than one:

$$ \frac{\rho_s \mathbf{v}_s \nabla (1-n)}{[(1-n)\rho_s] \nabla \cdot \mathbf{v}_s} \approx \frac{\rho_s |\mathbf{v}_s|^{\frac{(1-n)}{l_{(1-n)}}}}{[(1-n)\rho_s]^{\frac{|\mathbf{v}_s|}{l_{(1-n)}}}} \ll 1 \quad (A8) $$

## Appendix B. Equations of hydro-mechanical processes in unsaturated porous media

Definition of the fluid compressibility:

$$ \frac{\partial \rho_w}{\partial t} = \frac{d\rho_w}{dp_w} \frac{\partial p_w}{\partial t} = \rho_w C_w \frac{\partial p_w}{\partial t} \quad (B1) $$

where $C_w [kg^{-1}.m^1.s^2]$ is compressibility of water, $\rho_w [kg.m^{-3}]$ is the density of water and $p_w [kg.m^{-1}.s^{-2}]$ is the water pressure.

The density of the solid grains depends on the average pressure and the isotropic total stress. Thus, the time derivative of the solid grain density can be expressed as follows:

$$ \frac{\partial \rho_s}{\partial t} = \frac{\partial \rho_s}{\partial \sigma} \frac{\partial \sigma}{\partial t} + \frac{\partial \rho_s}{\partial p_{av}} \frac{\partial p_{av}}{\partial t} \quad (B2) $$

where $\rho_s [kg.m^{-3}]$ is the density of the solid grains, $p_{av} [kg.m^{-1}.s^{-2}]$ is the average pressure and $\sigma [kg.m^{-1}.s^{-2}]$ is isotropic total stress.

The average pressure is expressed as an average of the water and air pressures, as follows:

$$ p_{av} = S_w p_w + (1 - S_w) p_a \quad (B3) $$

where $p_a [kg.m^{-1}.s^{-2}]$ is the air pressure, which is assumed to be zero according to Richards' approximation and $S_w [-]$ is the degree of water saturation.

The time derivative of solid density can be formulated as follows:

$$ \frac{\partial \rho_s}{\partial t} = \frac{\partial \rho_s}{\partial \sigma} \frac{\partial \sigma}{\partial t} + \frac{\partial \rho_s}{\partial p_{av}} \left( S_w \frac{\partial p_w}{\partial t} + p_w \frac{\partial S_w}{\partial t} \right) \quad (B4) $$

$\partial \rho_s / \partial \sigma$ and $\partial \rho_s / \partial p_{av}$ can be expressed using the compressibility of solid grains as follows:
$$ \frac{\partial \rho_s}{\partial \sigma} = - \frac{\rho_s C_s}{1 - n} \qquad (B5) $$

$$ \frac{\partial \rho_s}{\partial p_{av}} = - \frac{\rho_s C_s}{1 - n} n \qquad (B6) $$

where $\sigma [kg.m^{-1}.s^{-2}]$ is isotropic total stress, $C_s [kg^{-1}.m^2.s^2]$ is compressibility of solid grains, and $n [-]$ is the porosity.

The isotropic total stress can be formulated as isotropic effective stress using the effective stress principle:

$$ \sigma = \sigma_e - \alpha S_w p_w \qquad (B7) $$

where $\sigma_e [kg.m^{-1}.s^{-2}]$ is isotropic effective stress and $\alpha [-]$ is the Biot's coefficient.

## Appendix C. On the use of COMSOL Multiphysics for modeling Fully Coupled hydro-mechanical (FC-HM) processes in a variably saturated porous medium

<figure><img src="image_19.png" /><figcaption>Fig. C1. On the use of COMSOL software to achieve FC-HM model.</figcaption></figure>
# Data availability

No data was used for the research described in the article.

## References

Abbasov, R., Fahs, M., Younes, A., Nowamooz, H., Jørgen Måløy, K., Toussaint, R., 2024. Modeling rainfall-induced landslide using the concept of local factor of safety: uncertainty propagation and sensitivity analysis. Comput. Geotech. 167, 106102. https://doi.org/10.1016/j.compgeo.2024.106102.

Abdullah, S., Ma, Y., Chen, X., Khan, A., 2022. A fully coupled hydro-mechanical-gas model based on mixture coupling theory. Transp. Porous Media 143, 47–68. https://doi.org/10.1007/s11242-022-01784-6.

Abdullah, S., Ma, Y., Wang, K., Subramanyam, S., Chen, X., Khan, A., 2023. Two-phase interaction in isothermal hydro-mechanical-chemical coupling for improved carbon geological sequestration modelling. Geomech. Geoen. 1–19. https://doi.org/10.1080/17486025.2023.2292160.

Asadi, R., Ataie-Ashtiani, B., 2015. A comparison of finite volume formulations and coupling strategies for two-phase flow in deforming porous media. Comput. Geotech. 67, 17–32. https://doi.org/10.1016/j.compgeo.2015.02.004.

Bai, T., Tahmasebi, P., 2022. Coupled hydro-mechanical analysis of seasonal underground hydrogen storage in a saline aquifer. J. Energy Storage 50, 104308. https://doi.org/10.1016/j.est.2022.104308.

Bandara, S., Ferrari, A., Laloui, L., 2016. Modelling landslides in unsaturated slopes subjected to rainfall infiltration using material point method. Num. Anal. Meth. Geomechan. 40, 1358–1380. https://doi.org/10.1002/nag.2499.

Ben-Zeev, S., Goren, L., Toussaint, R., Aharonov, O., 2023. Drainage explains soil liquefaction beyond the earthquake near-field. Nat. Commun. 14, 5791. https://doi.org/10.1038/s41467-023-41405-4.

Bilke, L., Flemisch, B., Kalbacher, T., Kolditz, O., Helmhig, R., Nagel, T., 2019. Development of open-source porous media simulators: principles and experiences. Transp. Porous Media 130, 337–361. https://doi.org/10.1007/s11242-019-01310-1.

Biot, M.A., 1941. General theory of three-dimensional consolidation. J. Appl. Math. 12, 155–164. https://doi.org/10.1063/1.1712886.

Both, J.W., Kumar, K., Nordbotten, J.M., Radu, F.A., 2019. Anderson accelerated fixed-stress splitting schemes for consolidation of unsaturated porous media. Comput. Math. Appl. 77, 1479–1502. https://doi.org/10.1016/j.camwa.2018.07.033.

Both, J.W., Pop, I.S., Yotov, I., 2021. Global existence of weak solutions to unsaturated poreelasticity. ESAIM: M2AN 55, 2849–2897. https://doi.org/10.1051/m2an/2021063.

Bui, Q.M., Osei-Kuffuor, D., Castelletto, N., White, J.A., 2020. A scalable multigrid reduction framework for multiphase poromechanics of heterogeneous media. SIAM J. Sci. Comput. 42, B379–B396. https://doi.org/10.1137/19M1256117.

Carsel, R.F., Parrish, R.S., 1988. Developing joint probability distributions of soil water retention characteristics. Water Resour. Res. 24, 755–769. https://doi.org/10.1029/WR024i005p00755.

Chen, X., Zhang, Lulu, Zhang, Limin, Zhou, Y., Ye, G., Guo, N., 2021. Modelling rainfall-induced landslides from initiation of instability to post-failure. Comput. Geotech. 129, 103877. https://doi.org/10.1016/j.compgeo.2020.103877.

Chiu, Y.-Y., Chen, H.-E., Yeh, K.-C., 2019. Investigation of the influence of rainfall runoff on shallow landslides in unsaturated soil using a mathematical model. Water 11, 1178. https://doi.org/10.3390/w11061178.

Clément, C., Toussaint, R., Stojanova, M., Aharonov, E., 2018. Sinking during earthquakes: critical acceleration criteria control drained soil liquefaction. Phys. Rev. E 97, 022905. https://doi.org/10.1103/PhysRevE.97.022905.

Eriksen, F.K., Toussaint, R., Turquet, A.L., Måløy, K.J., Flekkøy, E.G., 2018. Pressure evolution and deformation of confined granular media during pneumatic fracturing. Phys. Rev. E 97, 012908. https://doi.org/10.1103/PhysRevE.97.012908.

Fahs, M., Younes, A., Lehmann, F., 2009. An easy and efficient combination of the mixed finite element method and the method of lines for the resolution of Richards' Equation. Environ. Model Softw. 24, 1122–1126. https://doi.org/10.1016/j.envsoft.2009.02.010.

Fang, K., Jia, S., Tang, H., Dong, A., Ding, B., An, P., Zhang, B., Miao, M., 2025. New insights into convergent slope instability: Physical modeling and stability analysis. Eng. Geol. 356, 108295. https://doi.org/10.1016/j.enggeo.2025.108295.

Farthing, M.W., Ogden, F.L., 2017. Numerical solution of Richards' Equation: a review of advances and challenges. Soil Sci. Soc. Amer. J. 81, 1257–1269. https://doi.org/10.2136/sssaJ2017.02.0058.

Feder, J., Flekkøy, E.G., Hansen, A., 2022. Physics of Flow in Porous Media, 1st ed. Cambridge University Press. https://doi.org/10.1017/9781009100717.

Ghani, I., Koehn, D., Toussaint, R., Passchier, C.W., 2013. Dynamic development of hydrofracture. Pure Appl. Geophys. 170, 1685–1703. https://doi.org/10.1007/s00024-012-0637-7.

Gong, C., Wang, Y., Peng, Y., Ding, W., Lei, M., Da, Z., Shi, C., 2022. Three-dimensional coupled hydromechanical analysis of localized joint leakage in segmental tunnel linings. Tunn. Undergr. Space Technol. 130, 104726. https://doi.org/10.1016/j.tust.2022.104726.

Grunwald, N., Lehmann, C., Maßmann, J., Naumov, D., Kolditz, O., Nagel, T., 2022. Non-isothermal two-phase flow in deformable porous media: systematic open-source implementation and verification procedure. Geomech. Geophys. Geos-energ. Georesour. 8, 107. https://doi.org/10.1007/s40948-022-00394-2.

Guellouz, L., Bonelli, S., Débordes, O., 1997. Modélisation du couplage fluide-squelette en mécanique des sols (Application à la simulation numérique de la construction d'un barrage en terre).

Ireson, A.M., Spiteri, R.J., Clark, M.P., Mathias, S.A., 2023. A simple, efficient, mass-conservative approach to solving Richards' equation (openRE, v1.0). Geosci. Model Rev. 16, 659–677. https://doi.org/10.5194/gmd-16-659-2023.

Jeanne, P., Rutqvist, J., Vasco, D., Garcia, J., Dobson, P.F., Walters, M., Hartline, C., Borgia, A., 2014. A 3D hydrogeological and geomechanical model of an Enhanced Geothermal System at the Geysers, California. Geothermics 51, 240–252. https://doi.org/10.1016/j.geothermics.2014.01.013.

Jiang, S.-H., Yuan, Z.-R., Wang, Z.Z., Liu, X., 2025. Improved rainfall infiltration model for probabilistic slope stability assessment considering transition layers and spatial variability. Eng. Geol. 353, 108096. https://doi.org/10.1016/j.enggeo.2025.108096.

Johnsen, Ø., Toussaint, R., Måløy, K.J., Flekkøy, E.G., Schmittbuhl, J., 2008. Coupled air/granular flow in a linear Hele-Shaw cell. Phys. Rev. E 77, 011301. https://doi.org/10.1103/PhysRevE.77.011301.

Kafle, L., Xu, W.-J., Zeng, S.-Y., Nagel, T., 2022. A numerical investigation of slope stability influenced by the combined effects of reservoir water level fluctuations and precipitation: a case study of the Bianjiazhai landslide in China. Eng. Geol. 297, 106508. https://doi.org/10.1016/j.enggeo.2021.106508.

Koliji, A., 2008. Geotechdata.info website.

Lan, P., Huang, J., Su, J., Zhu, S., Chen, J., Zhang, S., Jiang, S.-H., 2025. A rapid unsaturated infiltration prediction method for slope stability analysis considering uncertainties: deep operator networks. Eng. Geol. 347, 107886. https://doi.org/10.1016/j.enggeo.2024.107886.

Lewis, R.W., Schrefler, B.A., 1998. The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media, 2nd ed. Wiley, Chichester.

Li, X., Ishikawa, T., Kasama, K., Zhang, T., 2025. Slope stability-soil runout integrated analysis considering surface flow during heavy rainfall. Eng. Geol. 356, 108268. https://doi.org/10.1016/j.enggeo.2025.108268.

Liakopoulos, A.C., 1964. Transient Flow through Unsaturated Porous Media (PhD thesis). University of California, Berkeley.

Lu, N., Likos, W.J., 2006. Suction stress characteristic curve for unsaturated soil. J. Geotech. Engenviron. Eng. 132, 131–142. https://doi.org/10.1061/(ASCE)1090-0241(2006)132:2(131).

Lu, N., Godt, J.W., Wu, D.T., 2010. A closed-form equation for effective stress in unsaturated soil: effective stress in unsaturated soil. Water Resour. Res. 46. https://doi.org/10.1029/2009WR008646.

Lu, N., Şener-Kaya, B., Wayllace, A., Godt, J.W., 2012. Analysis of rainfall-induced slope instability using a field of local factor of safety. Water Resour. Res. 48, 2012WR011830. https://doi.org/10.1029/2012WR011830.

Luo, Y., Zhang, J., Tang, C.-S., Jiang, G., Bogaard, T., 2025. Impact of dynamic desiccation cracks on hydrological processes and stability in expansive clay slopes: a coupled dual-permeability modeling approach. Eng. Geol. 357, 108377. https://doi.org/10.1016/j.enggeo.2025.108377.

Moradi, S., Huisman, J., Class, H., Vereecken, H., 2018. The effect of bedrock topography on timing and location of landslide initiation using the local factor of safety concept. Water 10, 1290. https://doi.org/10.3390/w10101290.

Moradi, S., Huisman, J.A., Vereecken, H., Class, H., 2024. Comparing different coupling and modeling strategies in hydromechanical models for slope stability assessment. Water 16, 312. https://doi.org/10.3390/w16020312.

Nemes, A., Schaap, M.G., Leij, F.J., 1999. The UNSODA Unsaturated Soil Hydraulic Database Version 2.0. US Salinity Laboratory, Riverside, CA.

Niěbling, M.J., Toussaint, R., Flekkøy, E.G., Måløy, K.J., 2012. Dynamic aerofracture of dense granular packings. Phys. Rev. E 86, 061315. https://doi.org/10.1103/PhysRevE.86.061315.

Paronuzzi, P., Bolla, A., 2023. Rainfall infiltration and slope stability of alpine colluvial terraces subject to storms (NE Italy). Eng. Geol. 323, 107199. https://doi.org/10.1016/j.enggeo.2023.107199.

Rafraß, S., Guellouz, L., 2022. Modeling of Water Flow in Deformable Unsaturated Porous Media.

Raia, S., Alvioli, M., Rossi, M., Baum, R.L., Godt, J.W., Guzzetti, F., 2014. Improving predictive power of physically based rainfall-induced shallow landslide models: a probabilistic approach. Geosci. Model Dev. 7, 495–514. https://doi.org/10.5194/gmd-7-495-2014.

Rangarajan, S., Rahardjo, H., Satyanaga, A., Li, Y., 2024. Influence of 3D subsurface flow on slope stability for unsaturated soils. Eng. Geol. 339, 107665. https://doi.org/10.1016/j.enggeo.2024.107665.

Salamov, A.M., Gadirov, A.H., Mukhtarov, A.Sh., Zamanova, A.H., Gasimov, E.E., 2020. Development of landslide processes on the eastern wing of the Bibiheybat brachianticline using the electrical survey. ANAS Trans. https://doi.org/10.33677/ggianas20200100039.

Schrefler, B.A., Scotta, R., 2001. A fully coupled dynamic model for two-phase fluid flow in deformable porous media. Comput. Methods Appl. Mech. Eng. 190, 3223–3246. https://doi.org/10.1016/S0045-7825(00)00390-X.

SCLCRA (State Committee for Land and Cartography of the Republic of Azerbaijan), 2025. Badamdar Sliding Zone.

Seki, K., Toride, N., van Genuchten, M.Th., 2023. Evaluation of a general model for multimodal unsaturated soil hydraulic properties. J. Hydrol. Hydromechan. 71, 22–34. https://doi.org/10.2478/johh-2022-0039.

Shao, W., Bogaard, T., Bakker, M., 2014. How to use COMSOL multiphysics for coupled dual-permeability hydrological and slope stability modeling. Procedia Earth Plan. Sci. 9, 83–90. https://doi.org/10.1016/j.proeps.2014.06.018.

Shao, W., Bogaard, T.A., Bakker, M., Greco, R., 2015. Quantification of the influence of preferential flow on slope stability using a numerical modeling approach. Hydrol. Earth Syst. Sci. 19, 2197–2212. https://doi.org/10.5194/hess-19-2197-2015.
Taherdangkoo, R., Nagel, T., Tang, A.M., Pereira, J.-M., Butscher, C., 2022. Coupled hydro-mechanical modeling of swelling processes in clay-sulfate rocks. Rock Mech. Rock. Eng. 55, 7489-7501. https://doi.org/10.1007/s00603-022-03039-8.

Terzaghi, K., 1925. Erdbaumechanik auf bodenphysikalischer grundlage. F. Deuticke.

Toussaint, R., Mäløy, K.J., Méheust, Y., Løvoll, G., Jankov, M., Schäfer, G., Schmittbuhl, J., 2012. Two-phase flow: structure, upscaling, and consequences for macroscopic transport properties. Vadose Zone J. 11 (vzj2011), 0123. https://doi.org/10.2136/vzj2011.0123.

Varela, J., 2018. Implementation of an MPFA/MPSA-FV Solver for the Unsaturated Flow in Deformable Porous Media (Master Thesis). The University of Bergen.

Virupaksha, A.G., Nagel, T., Lehmann, F., Rajabi, M.M., Hoteit, H., Fahs, M., Le Ber, F., 2024. Modeling transient natural convection in heterogeneous porous media with Convolutional Neural Networks. Int. J. Heat Mass Transf. 222, 125149. https://doi.org/10.1016/j.ijheatmasstransfer.2023.125149.

Wang, W., Rutqvist, J., Görke, U.-J., Birkholzer, J.T., Kolditz, O., 2011. Non-isothermal flow in low permeable porous media: a comparison of Richards' and two-phase flow approaches. Environ. Earth Sci. 62, 1197-1207. https://doi.org/10.1007/s12665-010-0608-1.

Wang, J.-H., Desai, C.S., Zhang, L., 2019. Soft soil and related geotechnical engineering practice. Int. J. Geomech. 19, 02019001. https://doi.org/10.1061/(ASCE)GM.1943-5622.0001494.

Wei, S., Ji, F., Ding, J., Guan, M., Lu, Y., 2025. Influence of extremely rapid cyclic reservoir water level fluctuations on bank slope stability: Insights from model testing and numerical simulation of a pumped storage power station slope. Eng. Geol. 357, 108380. https://doi.org/10.1016/j.enggeo.2025.108380.

Wu, L., Zhou, J., 2023. Rainfall Infiltration in Unsaturated Soil Slope Failure, SpringerBriefs in Applied Sciences and Technology. Springer Nature Singapore, Singapore. https://doi.org/10.1007/978-981-19-9737-2.

Wu, L.Z., Xu, Q., Zhu, J.D., 2017. Incorporating hydro-mechanical coupling in an analysis of the effects of rainfall patterns on unsaturated soil slope stability. Arab. J. Geosci. 10, 386. https://doi.org/10.1007/s12517-017-3147-1.

Wu, L.Z., Xu, Q., Wang, T., 2018. Incorporating Hydromechanical Coupling of Unsaturated Soils into the Analysis of Rainwater-Induced Groundwater Ponding. Int.

J. Geomech. 18, 06018010. https://doi.org/10.1061/(ASCE)GM.1943-5622.0001149.

Wu, L., Huang, R., Li, X., 2020. Hydro-Mechanical Analysis of Rainfall-Induced Landslides. Springer, Singapore.

Xue, L., Guo, X., Chen, H., 2020. Fluid Flow in Porous Media: Fundamentals and Applications, 1st ed. World Scientific, New Jersey London Singapore Beijing Shanghai Hong Kong Taipei Chennai Tokyo.

Yan, X., Sun, H., Huang, Z., Liu, L., Wang, P., Zhang, Q., Yao, J., 2021. Hierarchical modeling of hydromechanical coupling in fractured shale gas reservoirs with multiple porosity scales. Energy Fuel 35, 5758-5776. https://doi.org/10.1021/acs.energyfuels.0c03757.

Yan, C., Xie, X., Ren, Y., Ke, W., Wang, G., 2022. A FDEM-based 2D coupled thermal-hydro-mechanical model for multiphysical simulation of rock fracturing. Int. J. Rock Mech. Min. Sci. 149, 104964. https://doi.org/10.1016/j.ijrmms.2021.104964.

Yang, J., Fall, M., 2021. Hydro-mechanical modelling of gas transport in clayey host rocks for nuclear waste repositories. Int. J. Rock Mech. Min. Sci. 148, 104987. https://doi.org/10.1016/j.ijrmms.2021.104987.

Yang, Y.-S., Yeh, H.-F., 2023. Effects of hysteretic soil water characteristic curves on the hydromechanical behaviour. J. Earth Syst. Sci. 132, 136. https://doi.org/10.1007/s12040-023-02146-6.

Yang, Z., Qiao, J., Uchimura, T., Wang, L., Huang, D., Lei, X., Shi, L., 2017. Explore on hydro-mechanical threshold for early warning of rainfall induced shallow landslides. In: Mikoš, M., Casagli, N., Yin, Y., Sassa, K. (Eds.), Advancing Culture of Living with Landslides. Springer International Publishing, Cham, pp. 193-199. https://doi.org/10.1007/978-3-319-53485-5_21.

Zhan, Q., Wang, S., Guo, F., Chen, Y., Wang, L., Zhao, D., 2022. Early warning model and model test verification of rainfall-induced shallow landslide. Bull. Eng. Geol. Environ. 81, 318. https://doi.org/10.1007/s10064-022-02827-4.

Zhao, Z., 2023. Coupled Thermo-Hydro-Mechanical-Chemical Processes in Fractured Rocks: Fundamentals, Modeling and Applications. Springer Nature Singapore, Singapore. https://doi.org/10.1007/978-981-99-6210-5.