# Implicit EXP-RBF techniques for modeling unsaturated flow through soils with water uptake by plant roots

Mohamed Boujoudar, Abdelaziz Beljadid, Ahmed Taik

► To cite this version:

Mohamed Boujoudar, Abdelaziz Beljadid, Ahmed Taik. Implicit EXP-RBF techniques for modeling unsaturated flow through soils with water uptake by plant roots. Applied Numerical Mathematics: an IMACS journal, 2025, 208 (PartB), pp.79-97. 10.1016/j.apnum.2024.10.003. insu-04753295

HAL Id: insu-04753295

https://insu.hal.science/insu-04753295v1

Submitted on 25 Oct 2024

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
Research Paper

# Implicit EXP-RBF techniques for modeling unsaturated flow through soils with water uptake by plant roots

Mohamed Boujoudar<sup>a,b</sup>, Abdelaziz Beljadid<sup>a,c,*</sup>, Ahmed Taik<sup>d</sup>

## ARTICLE INFO

Keywords:

Richards equation

Infiltration

Root water uptake

Soil-water-plant interactions

Meshfree methods

Localized radial basis function

EXP-RBF

BDF2

## ABSTRACT

Modeling unsaturated flow through soils with water uptake by plant root has many applications in agriculture and water resources management. In this study, our aim is to develop efficient numerical techniques for solving the Richards equation with a sink term due to plant root water uptake. The Feddes model is used for water absorption by plant roots, and the van-Genuchten model is employed for capillary pressure. We introduce a numerical approach that combines the localized exponential radial basis function (EXP-RBF) method for space and the second-order backward differentiation formula (BDF2) for temporal discretization. The localized RBF methods eliminate the need for mesh generation and avoid ill-conditioning problems. This approach yields a sparse matrix for the global system, optimizing memory usage and computational time. The proposed implicit EXP-RBF techniques have advantages in terms of accuracy and computational efficiency thanks to the use of BDF2 and the localized RBF method. Modified Picards iteration method for the mixed form of the Richards equation is employed to linearize the system. Various numerical experiments are conducted to validate the proposed numerical model of infiltration with plant root water absorption. The obtained results conclusively demonstrate the effectiveness of the proposed numerical model in accurately predicting soil moisture dynamics under water uptake by plant roots. The proposed numerical techniques can be incorporated in the numerical models where unsaturated flows and water uptake by plant roots are involved such as in hydrology, agriculture, and water management.

## 1. Introduction

Understanding the dynamics of water flows in unsaturated soils and their interactions with plant root systems are critical for agricultural efficiency, environmental sustainability, and climate predictions. These interactions play a crucial role in the hydrological cycle, particularly within the vadose zone. This zone has a significant impact on subsurface water resources and acts as a reservoir for soil moisture, which offers essential water for plant growth. The mechanism by which plants absorb water from the soil, known as root water uptake, has an important influence on the distribution of soil moisture. A comprehensive understanding of the plant root water uptake process is pivotal due its applications in agriculture, irrigation, and water management. Understanding this process can

<footnote><sup>a</sup> University Mohammed VI Polytechnic, Morocco</footnote>

<footnote><sup>b</sup> Université d'Orléans, CNRS, BRGM, ISTO, UMR 7327, France</footnote>

<footnote><sup>c</sup> University of Ottawa, Canada</footnote>

<footnote><sup>d</sup> University Hassan II, Morocco</footnote>

<footnote>* Corresponding author at: University Mohammed VI Polytechnic, Morocco.
E-mail addresses: abdelaziz.beljadid@um6p.ma, abeljadi@uottawa.ca (A. Beljadid).</footnote>
aid in developing strategies to enhance crop productivity, improve water use efficiency, and better manage limited water resources in various agricultural contexts.

The plant root water uptake is influenced by soil heterogeneity and root architecture [1,2]. Furthermore, environmental variables such as humidity and temperature affect the plant's transpiration rate, which in turn impacts water uptake [3]. These complex interactions between soil properties, root characteristics, and environmental factors make root water uptake a challenging process to study and understand comprehensively. Robust coupled numerical models can be used as efficient tools to understand these complex interactions and study the impact of root water uptake on soil water distribution. These coupled models provide accurate simulations of water dynamics in the root zone and contribute to the enhancement of water resource management and informed decision-making in agriculture applications.

The infiltration of water into soils can be described at the continuum scale using Richards' equation [4]. This equation is derived from the combination of Buckingham-Darcy's law and the principle of mass conservation [5]. We note that many extensions of the Richards equation have been developed to reproduce the fingering phenomenon in unsaturated soils (see, e.g., [6–9] for an overview). Recently, Beljadid, Cueto-Felgueroso and Juanes [10] developed a continuum model of unstable infiltration in porous media, endowed with an entropy function. While many infiltration models are available, Richards' equation is extensively used in many practical applications for describing unsaturated flow in soils. In this study, we will use the Richards equation which exhibits high non-linearity due to the functions of the hydraulic conductivity and capillary pressure which are dependent on the soil water content [11–13].

Water uptake by roots is commonly described using two main approaches: macroscopic and microscopic [14–16]. Microscopic models, which are physically based, describe water extraction at the individual root level and consider the radial flow of soil water towards specific roots [14,17,18]. These approaches require detailed information about root geometry. The second category employs an empirical macroscopic approach [1,19,20], and considers the entire root system as a unified entity to account for the combined effects of individual roots. In the macroscopic approach, the plant root water uptake is considered as a volumetric sink term in the Richards equation [21,1]. Macroscopic models have been favored in several studies [16,22,23] due to their simplicity, as they do not necessitate detailed information about the root system's geometry. For instance, Molz and Remson [21] emphasized the challenge of modeling water transport in soil when considering microscopic models. The dynamic and complex geometry of the root system makes it impractical to measure accurately. Additionally, the water permeability of roots varies with their positions along the root, which further complicates the modeling process [24,16]. To address these complexities, most extraction functions have been developed using a macroscopic approach rather than a microscopic one [16,22,23].

In the present study, we use the macroscopic model proposed by Feddes et al. [1] to represent plant root water uptake. This model accounts for plant water stress, normalized root distribution, and transpiration potential. In addition, the van Genuchten model [13] is employed to describe capillary pressure, and the van Genuchten-Mualem model [25] is used to represent the relative permeability. For the validation and effectiveness of the proposed numerical model, we employed simplified models, such as stepwise and exponential forms [26] for plant root water uptake and the Gardner model for capillary pressure. The validation is established through comparison with existing analytical solutions.

The design of suitable numerical schemes for the Richards equation present significant challenges due to its highly nonlinear nature. In addition, incorporating an implicit sink term for root water uptake introduces additional numerical complexity of the governing equation. Few numerical methods have been proposed to solve the Richards equation with a sink term due to plant root water uptake. For instance, Wilderotter [27] used an adaptive finite-element method for solving the Richards equation with plant root growth. Janz and Stonier [28] applied finite difference method to solve Richards equation with the presence of simplified model for water absorption by roots. Machado et al. [29] used finite volume method to solve the Richards' equation with evapotranspiration. Efficient and robust numerical schemes are still in demand for modeling water flow in unsaturated soils while considering the root water uptake.

In terms of temporal discretization techniques, it is crucial to address the temporal derivative appropriately to ensure accurate numerical solutions for Richards' equation. This is particularly important due to the equation's highly nonlinear and the presence of stiff unsaturated soil properties [30]. Implicit or semi-implicit schemes are commonly employed for the temporal discretization to solve the Richards equation [31]. The primary reason for utilizing these schemes is to achieve stability and enable the use of practical time steps. The Backward Euler method is frequently used for temporal discretization in solving the Richards equation [32]. Few second-order temporal schemes have been proposed in the literature for solving Richards equation, including the Crank-Nicolson method [33–35] and the backward differentiation formula [36–38].

The primary linearization methods applied to Richards' equation include the Newton method, Picard method, and the modified Picard method. The Newton method, which is quadratically convergent, has been successfully applied to Richards' equation [39]. However, a significant drawback of the Newton method is its local convergence nature and the requirement for computing derivatives, which can be computationally intensive [39]. The Picard technique, despite its wide usage in addressing Richards' equation, faces challenges with convergence as highlighted in [40,41]. An improvement is proposed in [40] to obtain the modified Picard method. While this method maintains only linear convergence, it offers greater robustness compared to the Newton method [39].

In this paper, an efficient approach that combines the localized exponential RBF method [42] with a second-order backward differ-
entiation formula for temporal discretization is developed. Localized RBF methods, which employ a set of scattered nodes distributed
throughout the computational domain and its boundaries, eliminate the need for mesh generation and simplify the computational pro-
cess. This local approach generates sparse matrices, which avoids ill-conditioning problems and reduces computational time [42–44].
The selection of EXP-RBF is due to its high convergence rate, as demonstrated in numerous studies [45–47]. Additionally, its positive
definiteness has been validated in various system solutions [46,48,47]. Furthermore, the BDF2 scheme ensures stability and a rea-
sonable choice of time steps. Furthermore, it enhances computational efficiency in dealing with stiffness in the resolution of partial differential equations [49–51]. The modified Picard method [40] is used to linearize the system. This technique has been used in various studies to overcome mass conservation problems encountered in solving the Richards equation [40]. A variety of numerical schemes, such as the finite difference method [52], the finite element method [39], and the finite volume method [53], have been utilized to solve the Richards equation using the modified Picard iteration.

In this study, the developed approach is enhanced by the use of modified Picard iterations. This approach is employed to solve the Richards equation with a sink term due to plant root water uptake. The main goal is to provide accurate numerical results of soil moisture in the root zone and to investigate the influence of plant root water uptake on the distribution of soil moisture. In comparison to our previous studies [54,43,44], the model used in this study employs more nonlinear functions describing the capillary pressure and the relative permeability. This makes the development of efficient numerical techniques for solving the model more challenging. Furthermore, in the present study, we take into account the plant root water uptake as a sink term in the Richards equation, which was not considered in previous studies [54,43,44].

The paper is structured as follows: Section 2 presents the mathematical model of infiltration and plant root water uptake. In Section 3, we introduce the LRBF meshless method and linearization techniques proposed to solve the governing system. Section 4 is dedicated to numerical experiments for modeling soil moisture distribution in the root zone. Finally, Section 5 provides the concluding remarks.

## 2. Unsaturated flow and root water uptake models

### 2.1. Unsaturated flow model

Infiltration of water in soils is commonly described at the continuum scale using Richards' equation [4]. This equation can be derived by combining Buckingham-Darcy's law and the mass conservation equation [5]. The Richards equation with a sink term due to plant root water uptake and subject to initial and boundary conditions can be written as follows:

$$ egin{cases} \frac{\partial \theta}{\partial t} - \nabla \cdot [K_s(\mathbf{x})k_r \nabla \psi] - \frac{\partial [K_s(\mathbf{x})k_r]}{\partial z} = -s(\mathbf{x}, \psi), \quad \Omega \times (0, T), \\ -K_s(\mathbf{x})k_r \nabla (\psi + z).n = Q^N, \quad \partial \Omega^N \times (0, T), \\ \psi = \psi_D, \quad \partial \Omega^D \times (0, T), \\ \psi_{t=0} = \psi^0, \quad \Omega \times \{0\}, \end{cases} \quad (2.1) $$

where $\psi$ [L] is the pressure head, $\theta$ [$L^3/L^3$] is the volumetric soil water content, $k_r$ [-] is the water relative permeability, $K_s(\mathbf{x})$ [$L/T$] is the saturated hydraulic conductivity, $\mathbf{x} = (x, y, z)^T$ is the coordinate vector, $\mathbf{x}$ [L] and $\mathbf{y}$ [L] denote the horizontal dimensions and $z$ [L] denotes the vertical dimension positive upward, $T > 0$ is a fixed time, $\Omega$ is an open set of $\mathbb{R}^d$ ($d = 1, 2, 3$), $\partial\Omega = \partial\Omega^D \cup \partial\Omega^N$ is the border of $\Omega$, $\psi_D$ represents the prescribed pressure head associated with the Dirichlet boundary $\partial\Omega^D$. Similarly, $Q^N$ denotes the infiltration flux associated to the Neumann boundary $\partial\Omega^N$, which can represent precipitation, irrigation, and evaporation rates, $\psi^0$ is the initial condition associated to the pressure head and $n$ is the outward unit normal to the boundary of the domain. The sink term $s(\mathbf{x}, \psi)$ [$L^3/L^3T$] is defined as the volume of water removed per unit of time from a unit of volume of soil as a result of plant water uptake.

Here, we use the mixed form of the Richards equation (2.1), where both the pressure head and water content are considered. The Richards equation is highly nonlinear because of the nonlinear constitutive relations of the relative permeability and capillary pressure in terms of saturation [11–13].

The hydraulic conductivity $K = K_s k_r$ and water content in Equation (2.1) are estimated using empirical models, such as Gardner [11], Brooks-Corey [12], van Genuchten [13] and van Genuchten-Mualem [25,13] models. In this study, the capillary pressure is described using the van Genuchten constitutive relationship [13] model:

$$ \theta(\psi) = \begin{cases} \theta_r + \frac{\theta_s - \theta_r}{[1 + (\alpha_{vg} |\psi|)^{n_{vg}}]^{m_{vg}}}, & \text{if } \psi < 0, \\ \theta_s, & \text{if } \psi \ge 0, \end{cases} \quad (2.2) $$

and the van Genuchten-Mualem model [25,13] is used for relative permeability:

$$ K(\psi) = \begin{cases} K_s \frac{[1 - (\alpha_{vg} |\psi|)^{n_{vg}-1} [1 + (\alpha_{vg} |\psi|)^{n_{vg}}]^{-m_{vg}}]^2}{[1 + (\alpha_{vg} |\psi|)^{n_{vg}}]^{m_{vg}/2}}, & \text{if } \psi < 0, \\ K_s, & \text{if } \psi \ge 0, \end{cases} \quad (2.3) $$

where $\theta_r$ and $\theta_s$ represent the residual and saturated water contents, respectively [$L^3/L^3$], $\alpha_{vg}$ [$L^{-1}$] is related to the inverse of the air-entry pressure, $n_{vg} > 1$ is a measure of the pore-size distribution and $m_{vg} = 1 - 1/n_{vg}$.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Schematic of the plant water stress response function [1].</figcaption></figure>

## 2.2. Root water uptake model

Various mathematical models are developed to describe plant root water uptake and they are based on microscopic or macroscopic approaches [16]. In this study, we consider the macroscopic model proposed by Feddes et al. [1], which is expressed as follows:

$$s(x, \psi) = \alpha(\psi)S_p, \qquad (2.4)$$

where $S_p$ is the potential water uptake rate [$T^{-1}$] and $\alpha(\psi)$ [$-$] is the soil water stress function ($0 \le \alpha(\psi) \le 1$). Fig. 1 shows a schematic of the stress response function proposed by Feddes et al. [1], defined by the following equation:

$$\alpha(\psi) = \begin{cases} 0, & \psi \ge \psi_1, \psi \le \psi_4, \\ \frac{\psi - \psi_1}{\psi_2 - \psi_1}, & \psi_2 \le \psi \le \psi_1, \\ \frac{1}{\psi - \psi_4}, & \psi_3 < \psi \le \psi_2, \\ \frac{\psi - \psi_4}{\psi_3 - \psi_4}, & \psi_4 < \psi \le \psi_3. \end{cases} \qquad (2.5)$$

The stress response function $\alpha(\psi)$ is characterized by four critical pressure head values: $\psi_1$ denotes the pressure head level at which roots initiate water extraction from the soil, $\psi_2$ represents the pressure head level below which roots extract water at their maximum rate, $\psi_3$ indicates the limiting pressure head level below which roots cannot extract water at the maximum rate, and finally, $\psi_4$ signifies the pressure head level below which root water uptake ceases, often corresponding to the wilting point.

The variable $S_p$ denotes the water uptake rate under conditions of no water stress when $\alpha(\psi) = 1$. The potential root water uptake $S_p$ is described according to Feddes et al. [1] as follows:

$$S_p = b_r(x)T_p, \qquad (2.6)$$

where $T_p$ denotes the potential transpiration rate [$LT^{-1}$] and the function $b_r(x)$ [$L^{-3}$] describes the spatial variation of the potential extraction term. Various approaches are developed to express the spatial variation of the potential extraction term $b_r(x)$: constant, linear with the soil depth [1], or using the Hoffman and van Genuchten function [55]. For non-uniform root distributions, the function proposed by Vrugt [56] can be employed.

# 3. Numerical method

In this study, an efficient approach is proposed to solve the governing equation. This approach combines the localized EXP-RBF method with a second-order backward differentiation formula for temporal discretization. The modified Picards iteration method is employed to linearize the system. The localized EXP-RBF approach yields a sparse matrix for the global system, which reduces memory usage and computational time. The BDF2 scheme ensures the choice of reasonable time steps and enhances computational efficiency in dealing with stiffness in the resolution of the Richards equation.

## 3.1. Time discretization and linearization procedure

Here, we employ a second-order backward differentiation formula to approximate equation (2.1). The time interval $[0, T]$ is discretized into $N_t$ equally spaced intervals using a fixed time step size of $\Delta t = T/N_t$. Each discrete time point $t_n$ is defined as $t_n = n\Delta t$, where n ranges from 0 to $N_t$. The BDF2 discretization of equation (2.1) is defined as follows:

$$\frac{3\theta^{n+1} - 4\theta^n + \theta^{n-1}}{2\Delta t} - \nabla \cdot (K^{n+1} \nabla \psi^{n+1}) - \frac{\partial K^{n+1}}{\partial z} = -s^{n+1}. \qquad (3.1)$$
The nonlinear equation (3.1) is solved by computing successive approximations of $\psi^{n+1}$. If $m$ denotes the iteration level, the Picard iteration scheme can be expressed as follows:

$$ \frac{3\theta^{n+1,m+1} - 4\theta^n + \theta^{n-1}}{2\Delta t} - \nabla \cdot (K^{n+1,m} \nabla \psi^{n+1,m+1}) - \frac{\partial(K^{n+1,m})}{\partial z} = -s^{n+1,m}. \quad (3.2) $$

In this study, we employ the mixed form of the Richards equation (2.1), which incorporates both the $\psi$ and $\theta$ variables. We adopt the approach proposed by Celia et al. [40] to represent $\theta^{n+1,m+1}$ as a truncated Taylor series expansion with respect to $\psi$. The expansion can be written as:

$$ \theta^{n+1,m+1} = \theta^{n+1,m} + \frac{d\theta}{d\psi} \Big|_{\psi}^{n+1,m} (\psi^{n+1,m+1} - \psi^{n+1,m}) + O[(\delta^m)^2], \quad (3.3) $$

where $\delta^m$ is the iteration increment given by:

$$ \delta^m = \psi^{n+1,m+1} - \psi^{n+1,m}. \quad (3.4) $$

Substituting Equation (3.3) into (3.2) yields the following expression:

$$ \left(\frac{3}{2\Delta t} C^{n+1,m}\right) \delta^m - \nabla \cdot (K^{n+1,m} \nabla \psi^{n+1,m+1}) = \frac{\partial K^{n+1,m}}{\partial z} - \frac{3\theta^{n+1,m} - 4\theta^n + \theta^{n-1}}{2\Delta t} - s^{n+1,m}, \quad (3.5) $$

where $C = d\theta/d\psi$ is the specific moisture capacity function [1/L]. Equation (3.5) can be reformulated in terms of the iteration increment as follows:

$$ \left(\frac{3}{2\Delta t} C^{n+1,m}\right) \delta^m - \nabla \cdot (K^{n+1,m} \nabla \delta^m) = \nabla \cdot (K^{n+1,m} \nabla \psi^{n+1,m}) + \frac{\partial K^{n+1,m}}{\partial z} - \frac{3\theta^{n+1,m} - 4\theta^n + \theta^{n-1}}{2\Delta t} - s^{n+1,m}, \quad (3.6) $$

Equation (3.6) represents the general mixed-form Picard iteration, called the modified Picard approximation proposed by Celia et al. [40]. This representation is used in several studies for solving the Richards equation [40,57,58]. The final discrete form of the approximation can be obtained by applying finite difference method [52], finite element method [39] and finite volume method [53]. In this study, we use the localized EXP-RBF method to solve the Richards equation (3.6) based on the modified Picard iteration.

## 3.2. Localized EXP-RBF method

The localized EXP-RBF method is employed to solve the governing equation (3.6). First, it is necessary to derive the two linearized operators: $L^m$ corresponding to equation (3.6), and $B^m$ corresponding to the boundary conditions.

To simplify the three-dimensional representation of the spatial operator from Equation (3.6), we introduce the following simplified expressions:

$$ \mathcal{L}_d^m \delta_i = \frac{\partial}{\partial x^{(d)}} \left( K_i \frac{\partial \delta_i}{\partial x^{(d)}} \right), \quad (3.7) $$

and

$$ \mathcal{L}_4^m K_i = \frac{\partial K_i}{\partial z}, \quad (3.8) $$

where $(x^1, x^2, x^3) = (x, y, z)$ and $d$ represents the spatial dimension of the computation domain ($d = 1, 2, 3$).

We propose to approximate the operators $\mathcal{L}_d^m$ and $\mathcal{L}_4^m$ using the second-order standard finite difference method [40]. This approach is based on the average ratio of $K$ at soil interfaces. The resulting approximations are presented in the following equations:

$$ \mathcal{L}_d^m \delta_i = \frac{1}{(\Delta x^{(d)})^2} \left( K_{i+1/2}^{(d)} (\delta_{iR}^{(d)} - \delta_i) - K_{i-1/2}^{(d)} (\delta_i - \delta_{iL}^{(d)}) \right), \quad (3.9) $$

and

$$ \mathcal{L}_4^m K_i = \frac{1}{\Delta z} \left( K_{i+1/2}^{(d)} - K_{i-1/2}^{(d)} \right), \quad (3.10) $$

where $K_{i+1/2}^{(d)}$ and $K_{i-1/2}^{(d)}$ are given by:

$$ \left\{ \begin{array}{l} K_{i+1/2}^{(d)} = \frac{1}{2}(K_i + K_{iR}^{(d)}), \\ K_{i-1/2}^{(d)} = \frac{1}{2}(K_i + K_{iL}^{(d)}), \end{array} \right. \quad (3.11) $$

where $K_{iR}^{(d)}$ and $K_{iL}^{(d)}$ represent the right and left neighboring points along the $x^{(d)}$-axis, respectively, which are given by:
$$ K_{iR}^{(d)} = \begin{cases} K(x_{i+1}, y_i, z_i) & \text{if } d=1, \\ K(x_i, y_{i+1}, z_i) & \text{if } d=2, \\ K(x_i, y_i, z_{i+1}) & \text{if } d=3, \end{cases} \qquad (3.12) $$

and

$$ K_{iL}^{(d)} = \begin{cases} K(x_{i-1}, y_i, z_i) & \text{if } d=1, \\ K(x_i, y_{i-1}, z_i) & \text{if } d=2, \\ K(x_i, y_i, z_{i-1}) & \text{if } d=3. \end{cases} \qquad (3.13) $$

The same equations are utilized for the variable $\delta^m$ and $\psi^{n+1,m}$. The second spatial operator $\nabla \cdot (K_s k_r \nabla \delta^m)$ in Equation (3.6) can be represented as follows:

$$ \nabla \cdot (K_s k_r \nabla \delta^m) = L^m \delta^m, \qquad (3.14) $$

where $L^m = (\mathcal{L}_1^m + \mathcal{L}_2^m + \mathcal{L}_3^m)$. The same applies for $\nabla \cdot (K_s k_r \nabla \psi^{n+1,m}) = L^m \psi^{n+1,m}$.

Therefore, the Richards equation (3.6), can be expressed in the following form:

$$ \mathcal{L}^m \delta^m = a^{n+1,m}, \qquad (3.15) $$

where:

$$ \left\{ \begin{aligned} \mathcal{L}^m \delta^m &= \left( \frac{3}{2\Delta t} C^{n+1,m} \right) \delta^m - L^m \delta^m, \\ a^{n+1,m} &= L^m \psi^{n+1,m} + \mathcal{L}_4^m K^{n+1,m} - \frac{3\theta^{n+1,m} - 4\theta^n + \theta^{n-1}}{2\Delta t} - s^{n+1,m}. \end{aligned} \right. \qquad (3.16) $$

$\mathcal{L}^m$ represents the linearized operator within the domain $\Omega$ at each time step $n+1$ and iteration $m+1$.

Next, we will determine the linearized operator associated with the boundary conditions. The Dirichlet and Neumann boundary conditions (BC) are considered using the following equations:

$$ \left\{ \begin{aligned} \psi^{n+1,m+1} &= \psi_D, & \text{Dirichlet BC,} \\ -K^{n+1,m} \nabla(\psi^{n+1,m+1} + z).n &= Q^N, & \text{Neumann BC.} \end{aligned} \right. \qquad (3.17) $$

These boundary conditions can be written in terms of the variable $\delta^m$. For Dirichlet conditions, it implies that $\delta^m = 0$. For Neumann conditions, it leads to $-K^{n+1,m} \nabla(\delta^m + z).n = Q^N - K^{n+1,m} \nabla \psi^{n+1,m}.n$. As a result, the boundary operator $B^m$ can be represented as:

$$ B^m \delta^m = b^{n+1,m}, \qquad (3.18) $$

where the boundary operator is defined as follows:

$$ B^m \delta^m = \begin{cases} \delta^m, & \text{Dirichlet BC,} \\ -K^{n+1,m} \nabla(\delta^m + z).n, & \text{Neumann BC,} \end{cases} \qquad (3.19) $$

and

$$ b^{n+1,m} = \begin{cases} 0, & \text{Dirichlet BC,} \\ Q^N - K^{n+1,m} \nabla \psi^{n+1,m}.n, & \text{Neumann BC.} \end{cases} \qquad (3.20) $$

Finally, the resulting system that needs to be solved is represented by:

$$ \left\{ \begin{aligned} \mathcal{L}^m \delta^m &= a^{n+1,m}, \\ B^m \delta^m &= b^{n+1,m}. \end{aligned} \right. \qquad (3.21) $$

At each time level $t^{n+1}$, the linearized system (3.21) is solved using the localized EXP-RBF meshless method. The solution process entails iteratively solving the system at each iteration $m+1$ of the Picard iteration until the subsequent inequality is fulfilled:

$$ |\delta^m| \le Tol, \qquad (3.22) $$

where Tol is the error tolerance.

In the following, we describe the concept of the localized RBF method [59,60]. Let $\{x_s\}_{s=1}^N$ represent a set of collocation points distributed in the domain $\Omega \cup \partial\Omega$, where $N$ denotes the total number of points. For each point $x_s \in \Omega \cup \partial\Omega$, we create a corresponding subdomain $\Omega_s = \{x_k^{[s]}\}_{k=1}^{n_s}$ called influence domain, which contains the $n_s$ nearest neighboring points to $x_s$. We use the kd-tree algorithm [61] to determine the nearest neighbor points for each point $x_s$. Fig. 2 displays an example of influence domain $\Omega_s$ with seven points within a cube. The kd-tree algorithm is accurate and computationally efficient in identifying neighboring nodes of the evaluation point [62,60,54]. To approximate the variable $\delta^m$ using the localized EXP-RBF collocation method, we restrict the
<figure><img src="image_3.png" /><figcaption>Fig. 2. Example of influence domain $\Omega_s$ with seven points within a cube.</figcaption></figure>

collocation scheme to the subdomain $\Omega_s$ instead of the entire domain $\Omega \cup \partial\Omega$. The solution $\delta^m$ can be approximated as a linear combination of RBFs within each influence domain as follows:

$$ \delta^m(\mathbf{x}_s) = \sum_{k=1}^{n_s} \lambda_k^{[s]^{n+1,m+1}} \varphi(\|\mathbf{x}_s - \mathbf{x}_k^{[s]}\|), \quad (3.23) $$

where $\{\mathbf{x}_k^{[s]^{n+1,m+1}}\}_{k=1}^{n_s}$ represent the unknown coefficients, and $\varphi$ denotes the radial basis function. In this study, we use the exponential RBF, given by:

$$ \varphi(r_k) = \exp(-\epsilon^2 r_k^2), \quad (3.24) $$

where $r_k = \|\mathbf{x}_s - \mathbf{x}_k^{[s]}\|$ denotes the distance between $\mathbf{x}_s$ and $\mathbf{x}_k^{[s]}$, and $\epsilon > 0$ is a shape parameter. The appropriate selection of the shape parameter plays a crucial role in ensuring the accuracy and stability of RBF meshless methods [63]. Various studies have proposed optimal choices for the shape parameter of specific RBFs [64,65]. The localized RBF methods are less influenced by the selection of the shape parameter compared to the global RBF methods, as shown in previous studies [42,66].

The choice of the exponential RBF is motivated by its proven positive definiteness in various studies [67,47,48], which ensures that the resulting matrix is non-singular [68,69]. Additionally, a primary benefit of using the exponential RBF is the high convergence rate that can be achieved, as shown in several studies [45,46].

Based on Equation (3.23), we can represent the solution $\delta^m$ in matrix form as follows:

$$ \delta_{[s]}^{m} = \varphi^{[s]} \lambda^{[s]^{n+1,m+1}}, \quad (3.25) $$

where $\delta_{[s]}^m = [\delta^m(\mathbf{x}_1^{[s]}), \delta^m(\mathbf{x}_2^{[s]}), \dots, \delta^m(\mathbf{x}_n^{[s]})]^T$,

$\lambda^{[s]^{n+1,m+1}} = [\lambda^{[s]^{n+1,m+1}}(\mathbf{x}_1^{[s]}), \lambda^{[s]^{n+1,m+1}}(\mathbf{x}_2^{[s]}), \dots, \lambda^{[s]^{n+1,m+1}}(\mathbf{x}_n^{[s]})]^T$ and

$\varphi^{[s]} = [\varphi(\|\mathbf{x}_i^{[s]} - \mathbf{x}_j^{[s]}\|)]_{1 \le i,j \le n_s}$ is a real symmetric matrix of size $n_s \times n_s$. The unknown coefficients $\lambda^{[s]^{n+1,m+1}}$ can be determined as follows:

$$ \lambda^{[s]^{n+1,m+1}} = (\varphi^{[s]})^{-1} \delta_{[s]}^{m}. \quad (3.26) $$

Applying the linear operator $\mathcal{L}^m$ to Equation (3.10) results in the following equations for $\mathbf{x}_s \in \Omega$:

$$ \begin{aligned} \mathcal{L}^m \delta^m(\mathbf{x}_s) &= \sum_{k=1}^{n_s} \lambda_k^{[s]^{n+1,m+1}} \mathcal{L}^m \varphi(\|\mathbf{x}_s - \mathbf{x}_k^{[s]}\|) = \sum_{k=1}^{n_s} \lambda_k^{[s]^{n+1,m+1}} \Psi^m(\|\mathbf{x}_s - \mathbf{x}_k^{[s]}\|) \\ &= \vartheta_{[s]}^m \lambda_{[s]}^{n+1,m+1} = \vartheta_{[s]}^m (\varphi^{[s]})^{-1} \delta^m = \Lambda_{[s]}^m \mu_{[s]}^{n+1,m+1}, \end{aligned} \quad (3.27) $$

where $\Psi^m = \mathcal{L}^m \varphi$, $\vartheta_{[s]}^m = [\Psi^m(\|\mathbf{x}_s - \mathbf{x}_1^{[s]}\|), \dots, \Psi^m(\|\mathbf{x}_s - \mathbf{x}_n^{[s]}\|)]$ and $\Lambda_{[s]}^m = \vartheta_{[s]}^m (\varphi^{[s]})^{-1}$. To reformulate Equation (3.27) in terms of the global variable $\delta^m$ instead of the local variable $\delta_{[s]}^m$, we extend $\Lambda^m$ as the expansion of $\Lambda_{[s]}^m$ by adding zeros to the local vector where necessary. Thus, we establish the following correspondence:

$$ \mathcal{L}^m \delta_{[s]}^m(\mathbf{x}_s) = \Lambda^m \delta^m, \quad (3.28) $$

where we redefine $\delta^m$ at the discrete level by $\delta^m = [\delta^m(\mathbf{x}_1), \delta^m(\mathbf{x}_2), \dots, \delta^m(\mathbf{x}_N)]^T$.
Similarly, we apply the linear operator $\mathcal{B}^m$ to Equation (3.23) when $\mathbf{x}_s \in \partial\Omega$:

$$ 
\begin{aligned}
\mathcal{B}^m \delta_{[s]}^m(\mathbf{x}_s) &= \sum_{k=1}^{n_s} \lambda_k^{[s]^{m+1,n+1}} \mathcal{B}^m \varphi(\|\mathbf{x}_s - \mathbf{x}_k^{[s]}\|) = (\mathcal{B}^m \varphi^{[s]}) \lambda^{[s]^{m+1,n+1}} \tag{3.29} \\
&= (\mathcal{B}^m \varphi^{[s]}) (\varphi^{[s]})^{-1} \delta_{[s]}^m = \sigma^{[s]} \delta_{[s]}^m = \sigma \delta^m.
\end{aligned}
$$

Here, $\sigma^{[s]}$ is defined as $(\mathcal{B}^m \varphi^{[s]})(\varphi^{[s]})^{-1}$, and $\sigma$ represents the expansion of $\sigma^{[s]}$ obtained by introducing zeros at the appropriate positions. By combining Equations (3.28) and (3.29) into Equation (3.21), the following system is derived:

$$ 
\begin{aligned}
\mathcal{L}^m \delta^m(\mathbf{x}_s) &= \Lambda^m(\mathbf{x}_s) \delta^m = a^{n+1,m}(\mathbf{x}_s), \\
\mathcal{B}^m \delta^m(\mathbf{x}_s) &= \sigma(\mathbf{x}_s) \delta^m = b^{n+1,m}(\mathbf{x}_s).
\end{aligned}
\tag{3.30}
$$

As a result, the localized EXP-RBF meshless method applied to the governing equation yields the final system to be solved:

$$ A^m X^m = B^m, \tag{3.31} $$

$$ \text{where } A^m = \begin{pmatrix} \Lambda^m(\mathbf{x}_1) \\ \Lambda^m(\mathbf{x}_2) \\ \vdots \\ \Lambda^m(\mathbf{x}_{N_i}) \\ \sigma(\mathbf{x}_{N_i+1}) \\ \vdots \\ \sigma(\mathbf{x}_N) \end{pmatrix}, X^m = \begin{pmatrix} \delta^m(\mathbf{x}_1) \\ \delta^m(\mathbf{x}_2) \\ \vdots \\ \delta^m(\mathbf{x}_{N_i}) \\ \delta^m(\mathbf{x}_{N_i+1}) \\ \vdots \\ \delta^m(\mathbf{x}_N) \end{pmatrix} \quad \text{and} \quad B^m = \begin{pmatrix} a^{n+1,m}(\mathbf{x}_1) \\ a^{n+1,m}(\mathbf{x}_2) \\ \vdots \\ a^{n+1,m}(\mathbf{x}_{N_i}) \\ b^{n+1,m}(\mathbf{x}_{N_i+1}) \\ \vdots \\ b^{n+1,m}(\mathbf{x}_N) \end{pmatrix}. $$

This localized approach offers the advantage of inverting a sparse matrix, which mitigates the ill-conditioning issues encountered when dealing with the full matrix generated by the global method.

By solving this sparse system, we obtain the approximate solution $\delta^m$ at the specified points. Once the condition specified in Equation (3.22) is verified, we assign $\psi^{n+1}$ the value of $\delta^m + \psi^{n+1,m}$.

## 4. Numerical experiments

In this section, we present numerical experiments for solving the Richards equation with sink terms due to plant root water uptake. The proposed approach uses the localized RBF meshless method based on the exponential function to solve the governing system (3.21). To assess the effectiveness of the proposed numerical approach, several numerical tests are conducted.

In the first test, two simplified models for root water uptake are used to validate the proposed numerical model by conducting a comparison between the numerical results and analytical solutions. In the second numerical test, the Feddes model [70] is employed to represent the plant root water uptake. The results obtained using the proposed numerical model are compared with those generated by Hydrus software [71], which serves as a benchmark for validation. The last numerical test evaluates the performance of the proposed numerical model to predict the three-dimensional soil moisture distribution profile in the root zone under axi-symmetric irrigation conditions. The numerical tests are performed on a core CPU i7 2.1 GHz computer in a MATLAB 2018a tools.

### 4.1. Unsaturated flow in rooted soils under variable surface flux conditions

In this numerical test, we validate our proposed numerical model by conducting a comparison between the approximate solutions and analytical solutions established in [26]. For this comparison, we use the Gardner model [11] to represent the relative permeability and a simple formulation for the water content [72], given by:

$$ K = K_s \exp(\alpha \psi), \tag{4.1} $$

$$ \theta = \theta_r + (\theta_s - \theta_r) \exp(\alpha \psi). \tag{4.2} $$

Both stepwise and exponential forms [26] are considered for root water uptake [73,74,26], which are respectively given by:

$$ s(z) = R_0 H(z - l_1) = \begin{cases} R_0, & \text{if } l_1 \le z \le l, \\ 0, & \text{if } 0 \le z < l_1. \end{cases} \tag{4.3} $$

$$ s(z) = R_0 \exp[\beta(z - l)], \tag{4.4} $$

where $\alpha [1/L]$ is an empirical parameter related to the macroscopic capillary length of the soil, $R_0 [T^{-1}]$ is the maximum root uptake at the soil surface, $\beta [L^{-1}]$ is a parameter represents the rate of reduction of root uptake, $l [L]$ is the soil depth and $l_1 [L]$ is the maximum root depth.

We perform numerical simulations of soil moisture and pressure head through a rooted soil with the following parameter values $l = 100$ cm, $l_1 = 60$ cm, $\theta_s = 0.45$, $\theta_r = 0.2$, $K_s = 1$ cm $h^{-1}$ and $\beta = 0.04$ $m^{-1}$ [75,26].
<figure><img src="image_4.png" /><figcaption>Fig. 3. Case 1: $\alpha = 0.01 \text{ cm}^{-1}$ and $R_0 = 0.02 \text{ h}^{-1}$. Comparison between approximate and exact solutions for water content and pressure head. Left: with root water uptake. Right: without root water uptake (Solid lines represent the approximate solution and dashed lines represent the exact solution.).</figcaption></figure>

We consider two values of $\alpha$ and $R_0$: $\alpha = 0.01 \text{ cm}^{-1}$ corresponds to $R_0 = 0.02 \text{ h}^{-1}$, and $\alpha = 0.1 \text{ cm}^{-1}$ corresponds to $R_0 = 0.0025 \text{ h}^{-1}$. The following initial and boundary conditions are used:

$$ \left\{ \begin{array}{l} \psi(z, 0) = \psi_0(z), \\ \psi(0, t) = 0, \\ \left[-K \left(\frac{\partial \psi}{\partial z} + 1\right)\right]_{z=L} = q_1(t), \end{array} \right. \qquad (4.5) $$

where $\psi_0$ is the initial pressure head and $q_1$ is the flux at the soil surface, which changes over time. Our numerical test considers both steady-state and time-variant surface fluxes.

First, we apply a steady-state surface flux as the upper boundary condition, considering a constant infiltration flux of $0.9 \text{ cm h}^{-1}$ (i.e. $q_1 = -0.9 \text{ cm h}^{-1}$). For the root water absorption, we use the formulation given by equation (4.3). We present the distribution of water content and pressure head in Figs. 3 and 4, respectively. The simulations are performed using the two values of $\alpha$: $0.01 \text{ cm}^{-1}$ and $0.1 \text{ cm}^{-1}$.

The numerical simulations are conducted over a 50-hour duration. To investigate the impact of root water uptake on soil moisture and pressure head distribution, we compare the time evolution of the water content and pressure head with and without considering root water uptake. The corresponding results are displayed for both cases, with and without root water uptake.

We observe a significant impact of root water uptake, particularly when $\alpha = 0.01 \text{ m}^{-1}$ and $R_0 = 0.02 \text{ h}^{-1}$. The water content encounters significant changes under these conditions. However, for the case with $\alpha = 0.1 \text{ m}^{-1}$ and $R_0 = 0.0025 \text{ h}^{-1}$, due to the relatively low value of the maximum water uptake, the effect on water content is relatively weak.
M. Boujoudar, A. Beljadid and A. Taik

Applied Numerical Mathematics xxx (xxxx) xxx

<figure><img src="image_5.png" /><figcaption>Fig. 4. Case 1: $\alpha = 0.01 \text{ cm}^{-1}$ and $R_0 = 0.02 \text{ h}^{-1}$. Comparison between approximate and exact solutions for water content and pressure head. Left: with root water uptake. Right: without root water uptake (Solid lines represent the approximate solution and dashed lines represent the exact solution).</figcaption></figure>

With the same localized EXP-RBF parameters $N_z = 1001$, $n_s = 3$, $\epsilon = 0.1$, and $\Delta t = 0.01$, we display in Fig. 5 the root mean square error (RMSE) associated to water content as a function of the time step for both first-order backward differentiation formula (BDF1) and BDF2 schemes.

We confirm that the BDF2 scheme consistently achieves a convergence rate of approximately 2, while the BDF1 method exhibits a convergence rate of 1. Furthermore, BDF2 outperforms BDF1 in terms of accuracy.

Secondly, we consider a time-variable surface flux as the upper boundary condition, providing a more realistic representation of real-world conditions influenced by factors such as evaporation, rainfall, and irrigation. We consider that the surface flux decreases exponentially over time, specifically $q_1(t) = q_0 + \delta_r \exp(k_1 t)$, with $\delta_r = -0.8 \text{ cm h}^{-1}$ and $k_1 = -0.1 \text{ h}^{-1}$. Root water uptake is modeled using the exponential formulation (4.4). Fig. 6 illustrates the variation of water content in both time and space for the rooted soils. The left figure corresponds to $\alpha = 0.01 \text{ cm}^{-1}$ and $R_0 = 0.02 \text{ h}^{-1}$, while the right figure corresponds to $\alpha = 0.1 \text{ cm}^{-1}$ and $R_0 = 0.0025 \text{ h}^{-1}$.

Both soils receive an equal water supply from the surface, but their moisture content patterns differ significantly. Soil profile 1, characterized by $R_0 = 0.02 \text{ h}^{-1}$ and $\alpha = 0.01 \text{ cm}^{-1}$, generally exhibits higher moisture levels than soil profile 2 with $R_0 = 0.0025 \text{ h}^{-1}$ and $\alpha = 0.1 \text{ cm}^{-1}$. Although soil profile 1 absorbs more water through its roots than profile soil 2, it benefits from the capillary rise, which facilitates water transfer from the water table to the root zones. This accounts for the noticeable difference in moisture content patterns between the two soils.

Here, we will show the advantage of using BDF2, for temporal discretization in the proposed numerical model, in terms of computational time. To this end, we will analyze the CPU time of the proposed numerical techniques which are based on LRBF approach while using BDF2 or BDF1 for temporal integration. Numerical simulations are performed for the same numerical example until the final time $T = 50$ hours for both temporal schemes. With a time step $\Delta t = 0.1$, BDF2 yields an RMSE of $1.64 \times 10^{-5}$ and requires a CPU time of 53 seconds. For consistency in terms of comparison, we repeat the simulation using BDF1 in order to achieve the same level of accuracy which requires a time step $\Delta t = 0.015$ and leads to a CPU time of 293 seconds. This comparison emphasizes the advantage of using BDF2 in our numerical model. The performance of the proposed numerical model in terms of total mass
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Root mean square error (RMSE) as a function of the time step for BDF1 and BDF2 schemes. Left: $\alpha = 0.1 \text{ cm}^{-1}$ and $R_0 = 0.0025 \text{ h}^{-1}$. Right: $\alpha = 0.01 \text{ cm}^{-1}$ and $R_0 = 0.02 \text{ h}^{-1}$.</figcaption></figure>

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

<figure><img src="image_9.png" /><figcaption>Fig. 6. Time evolution of soil water content. Left: $\alpha = 0.01 \text{ cm}^{-1}$ and $R_0 = 0.02 \text{ h}^{-1}$. Right: $\alpha = 0.1 \text{ cm}^{-1}$ and $R_0 = 0.0025 \text{ h}^{-1}$.</figcaption></figure>

conservation is analyzed under both constant and variable surface flux conditions. We study the evolution of the total mass ratio $MB(t)$ of water in the computational domain. Here, we will present the general methodology which will be applied for both 1D and 2D cases to compute the following mass balance ratio [40,76,53]:

$$MB(t) = \frac{\text{Total additional mass in the domain}}{\text{Total net flux into the domain}}, \qquad (4.6)$$

where the total additional mass in the domain is the difference between the mass measured at any time $t$ and the initial mass in the domain, and the total net flux into the domain is the flux balance integrated over time up to $t$. These quantities can be obtained from spatial and temporal integrations of the balance law (2.1) and are formulated as follows [71]:

$$\text{Total additional mass in the domain} = \int_{\Omega} \theta(\mathbf{x}, t)d\Omega - \int_{\Omega} \theta(\mathbf{x}, 0)d\Omega, \qquad (4.7)$$

$$\text{Total net flux into the domain} = \int_{0}^{t} \int_{\partial\Omega} q(-n) d\Gamma dt - \int_{0}^{t} \int_{\Omega} s d\Omega dt, \qquad (4.8)$$

where,

$$q = -K \nabla(\psi + z), \qquad (4.9)$$

$d\Gamma$ denotes an infinitesimal arc length of the boundary $\partial\Omega$ and as already defined $n$ is the outward unit normal vector of the boundary. The relative error in mass balance, $R_e M(t)$, is defined as:

$$R_e M(t) = |1 - MB(t)|. \qquad (4.10)$$
<figure><img src="image_10.png" /><figcaption>Fig. 7. Mass balance ratio for unsaturated flow in rooted soils. Left: Constant surface flux. Right: Variable surface flux.</figcaption></figure>

<table><caption>Table 1
Plants parameters.</caption><thead><tr><th>Plant</th><th>ψ<sub>1</sub><br>(m)</th><th>ψ<sub>2</sub><br>(m)</th><th>ψ<sub>3,low</sub><br>(m)</th><th>ψ<sub>3,high</sub><br>(m)</th><th>ψ<sub>4</sub><br>(m)</th></tr></thead><tbody><tr><td>Pasture</td><td>-0.1</td><td>-0.25</td><td>-2</td><td>-8</td><td>-80</td></tr><tr><td>Wheat</td><td>0</td><td>-0.01</td><td>-5</td><td>-9</td><td>-160</td></tr></tbody></table>

Fig. 7 illustrates the time evolution of the mass balance ratio for both constant and variable surface flux. Since the total net flux into the domain in the denominator of (4.6) is zero at $t = 0$, the mass balance ratio is computed from $t = 1$ h. We obtain accurate results in terms of conservation and the relative errors at the final time $T = 50$ h are $6.11 \times 10^{-4}$ and $2.09 \times 10^{-3}$, respectively, for constant and variable surface fluxes.

## 4.2. Root water uptake soil profile with groundwater table

In this numerical example, we used Equations (2.2), (2.3) and (2.4) to describe root water uptake. We consider a rooted soil with a depth of $1.2$ m and a root depth of $0.9$ m. The root distribution is assumed to decrease linearly with depth. The bottom of the soil is in hydrostatic equilibrium with the groundwater table at the beginning of the simulation ($\psi(z, t = 0) = -z$). The upper boundary condition is characterized by a non-flow boundary condition due to the absence of precipitation or irrigation. As a result, the plant's transpiration depletes the soil water, which is partially recharged by capillary action from the groundwater.

Two types of plants, wheat and pasture, are considered, with their respective stress function parameters outlined in Table 1. $\psi_{3,low}$ and $\psi_{3,high}$ represent the critical pressure head values, below which roots can no longer efficiently extract water at the maximum rates, $r_{2L}$ and $r_{2H}$, respectively. These potential transpiration rates, $r_{2L}$ and $r_{2H}$ are presently calibrated at 0.1 cm/day and 0.5 cm/day, as reported in [71]. To get the value of $\psi_3$, we use the following interpolation proposed by [71]:

$$ \psi_3 = \begin{cases} \psi_{3,low} + \frac{(h_{3,high} - h_{3,low})}{(r_{2L} - r_{2H})} (r_{2h} - T_p), & r_{2L} < T_p < r_{2H}, \\ \psi_{3,low}, & T_p \le r_{2L}, \\ \psi_{3,high}, & T_p \ge r_{2H}. \end{cases} \qquad (4.11) $$

The hypothetical plants are presumed to exhibit a potential transpiration rate $T_p = 4$ mm/day [77,78]. These simulations were conducted over a period of 50 days. This numerical test is treated in [77,78], where the authors used Hydrus to simulate root water uptake. In this case, the van Genuchten and the van Genuchten-Mualem [25,13] models are used to describe the capillary pressure and the relative permeability, respectively. The soil investigated in this study is identified as loamy soil [77], and its hydraulic properties are given by $\theta_s = 0.430$, $\theta_r = 0.078$, $\alpha_{vg} = 3.6 m^{-1}$, $K_s = 0.2496$ m/day, $n_{vg} = 1.56$ and $m_{vg} = 0.3590$. We selected 1001 points along the z-axis and used a time step of 0.01. The localized EXP-RBF function was configured with the parameters $\epsilon = 0.1$ and $n_s = 3$. Figs. 8 and 9 display the pressure head profiles and root water uptake at different time levels for two plants: pasture (left) and wheat (right), respectively. A comparison is made between the approximate solutions obtained using the proposed numerical method and the results provided by Hydrus-1D. A good agreement is observed across all cases, indicating a good level of consistency. The accuracy of our proposed numerical method in predicting soil moisture, accounting for root water uptake, is reinforced by the obtained root mean square error (RMSE) values.
M. Boujoudar, A. Beljadid and A. Taik

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

<figure><img src="image_12.png" /><figcaption>Fig. 8. Pressure head profiles for two considered plants. Left: Pasture. Right: Wheat (Solid lines represent the approximate solution and dashed lines represent the reference solution using Hydrus-1D).</figcaption></figure>

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

<figure><img src="image_14.png" /><figcaption>Fig. 9. Root water uptake profiles for two considered plants. Left: Pasture. Right: Wheat (Solid lines represent the approximate solution and dashed lines represent the reference solution using Hydrus-1D).</figcaption></figure>

The figures clearly demonstrate that RMSE values obtained from the numerical simulations are remarkably small, particularly when considering the error associated with the estimation of root water uptake. Fig. 10 presents the temporal evolution of the bottom flux $-K(\psi) [\partial\psi/\partial z + 1]_{z=0}$, potential transpiration flux $T_p$, and actual transpiration flux $T_a$, as well as their respective cumulative values is displayed in Fig. 11. $T_a$ is given by [71]:

$$T_a = T_p \int_{\Omega_R} \alpha(\psi)b_r(z)dz. \qquad (4.12)$$

These variations are observed over a period of 50 days for the two plant types under investigation.

Once again, there is strong agreement between the approximate solutions using localized EXP-RBF meshless method and the results given by Hydrus-1D for all cases, demonstrating the effectiveness of the proposed numerical approach.

Fig. 12 shows the evolution of the mass balance ratio. The corresponding relative error in mass balance at $t = 50$ days is $4.34 \times 10^{-3}$ which demonstrates the effectiveness of the proposed numerical model in terms of mass conservation.

### 4.3. Three-dimensional soil moisture distribution profile under axi-symmetric irrigation conditions

This numerical test was used by Šimůnek and Hopmans [77], where a three-dimensional axi-symmetrical profile is considered to simulate soil moisture in root zone. The computational domain can therefore be restricted to a square as shown in Fig. 13, where $r$ denotes the radial direction. The domain's depth and radius are both set at $1$ m. An irrigation rate $q_z = 0.16$ cm/h is applied to the left half of the soil surface. In this case, the initial pressure head is set at $\psi_0 = -4$ m across the entirety of the soil profile, while the potential transpiration rate is set at $0.04$ cm/h. The Feddes model given by Equation (2.4) is used for the root water uptake model.
M. Boujoudar, A. Beljadid and A. Taik

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

Applied Numerical Mathematics xxx (xxxx) xxx

Wheat

<figure><img src="image_16.png" /><figcaption>**Fig. 11.** The cumulative flux of the potential transpiration, actual transpiration and bottom flux. Left: Pasture. Right: Wheat (Solid lines represent the approximate solution and dashed lines represent the reference solution using Hydrus-1D).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 12.** Mass balance ratio in the case of root water uptake model proposed by Feddes et al. [1].</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. 13. Sketch of the computational domain.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 14. The water content evolution at different times.</figcaption></figure>

The root density is assumed uniform in the radial direction and varied linearly with depth, from its maximum value at the soil surface to zeros at the soil depth 0.5 m.

We conduct numerical simulations using the relationships described in equations (2.2)–(2.3), while employing the same parameters for the loamy soil as in the previous test. Furthermore, we use the same values for the soil water stress response parameters as listed in Table 1. A zero flux boundary condition is imposed on the horizontal side of the computational domain to model the absence of flow across this boundary. Additionally, a free drainage condition is assumed at the bottom of the soil, representing the soil ability to freely drain water without any significant impedance or restriction. The boundary conditions are given by:

$$ \left\{ \begin{array}{ll} -K \left( \frac{\partial \psi}{\partial z} + 1 \right) = -q_z, & z = 1, \\ -K \frac{\partial \psi}{\partial r} = 0, & r = 0, r = 1, \\ -K \frac{\partial \psi}{\partial z} = 0, & z = 0. \end{array} \right. \qquad (4.13) $$

In [77], the authors introduce two distinct sink terms to account for root water uptake. The first sink term corresponds to the conventional Feddes model, representing uncompensated root water uptake used in this study. The second sink term is known as compensated root water uptake.

In this numerical test, we use the following parameters: $\varepsilon = 0.2$, $\Delta t = 0.005$, $n_s = 5$, $N_x = 100$, and $N_z = 200$. Fig. 14 illustrates the temporal evolution of water content in the case of uncompensated root water uptake.

The numerical results align with the expected characteristics observed in the uncompensated scenario, as reported by [77].
<figure><img src="image_20.png" /><figcaption>Fig. 15. Time evolution of the potential and actual transpiration.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 16. Mass balance ratio in the case of axi-symmetric irrigation.</figcaption></figure>

In Fig. 15, we display the temporal evolution of the potential transpiration $T_p$ and actual transpiration $T_a$ computed using Equation (4.12). A comparison is performed between the approximate solutions obtained using the proposed numerical method and the results obtained from Hydrus [77]. The results demonstrate a good agreement, confirming the accuracy and reliability of our proposed numerical approach.

Fig. 16 illustrates the time evolution of the mass balance ratio. We obtain good result in terms of conservation and the relative error in mass balance at $T = 7.5$ days is $9.10 \times 10^{-3}$.

## 5. Conclusion

This study introduces efficient numerical techniques for solving the Richards equation with a sink term due to plant root water uptake. The van Genuchten and Feddes models are employed to describe capillary pressure and plant root water absorption, respectively. Our aim is to develop a numerical model that can accurately predict the soil moisture distribution in the root zone and to study the impact of plant root water absorption on soil moisture distribution.

An efficient approach is proposed, combining the localized exponential radial basis function method for space and the second-order backward differentiation formula for temporal discretization. The modified Picard's iteration method is used to linearize the system represented in the mixed form of the Richards equation. The localized RBF methods, which use a unique set of scattered nodes distributed throughout the computational domain and its boundaries, eliminate the need for mesh generation and simplify the computational process. These techniques effectively overcome ill-conditioning problems often encountered in the use of global RBF methods since a sparse matrix for the global system is obtained. This localized approach improved memory usage and computational time. Additionally, the second-order accurate BDF2 scheme ensures reasonable time steps and enhances the accuracy of the method.
M. Boujoudar, A. Beljadid and A. Taik

Applied Numerical Mathematics xxx (xxxx) xxx

The validation of the proposed numerical model for infiltration with plant root water absorption is performed using numerical experiments. The numerical results obtained using the proposed numerical model are compared against analytical and reference solutions.

The numerical model was validated for modeling unsaturated flow in rooted soils under variable surface flux conditions. Accurate results are obtained in comparison to the corresponding analytical solutions within a short CPU time. These results are confirmed by the RMSE values obtained across all considered cases, underlining the accuracy of the proposed numerical model in predicting the soil moisture dynamics in the root zone. Numerical simulations are performed using the proposed numerical techniques for modeling root water uptake soil profile with groundwater table. We obtained good agreement between the numerical simulations and the results generated by the Hydrus software. The numerical simulations of three-dimensional soil moisture distribution profiles under axi-symmetric irrigation conditions, confirmed the accuracy and reliability of the proposed numerical model by comparing its findings with those presented in existing literature. The numerical simulations confirm that the combination of localized EXP-RBF and BDF2 methods leads to efficient numerical model for simulating soil moisture dynamics in the presence of root water uptake which is important in understanding soil-water-plant interactions.

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.

CRediT authorship contribution statement

Mohamed Boujoudar: Writing – original draft, Software, Investigation. Abdelaziz Beljadid: Writing – review & editing, Supervision, Investigation, Conceptualization. Ahmed Taik: Writing – review & editing, Supervision.

Acknowledgements

AB gratefully acknowledges funding from UM6P/OCP Group of Morocco (UM6P-OCP-2019). AB and AT gratefully acknowledge funding for the APRD research program from the Moroccan Ministry of Higher Education, Scientific Research and Innovation and the OCP Foundation (25219).

References

[1] R.A. Feddes, P.J. Kowalik, H. Zaradny, Simulation of Field Water Use and Crop Yield, Pudoc for the Centre for Agricultural Publishing and Documentation, Wageningen, 1978.

[2] R. Feddes, P. Raats, Parameterizing the soil-water-plant root system, in: Unsaturated-zone modeling; Progress, challenges and applications, no. 6, Kluwer Academic Publishers, 2004, pp. 95-141.

[3] F. Tardieu, T. Simonneau, Variability among species of stomatal control under fluctuating soil water status and evaporative demand: modelling isohydric and anisohydric behaviours, J. Exp. Bot. 49 (1998) 419-432.

[4] L.A. Richards, Capillary conduction of liquids through porous mediums, Physics 1 (5) (1931) 318-333, https://doi.org/10.1063/1.1745010.

[5] J. Bear, A. Verruijt, Theory and Applications of Transport in Porous Media, Modeling of Groundwater Flow and Pollution, Reidel, Dordrecht, 1987.

[6] Mehdi Eliassi, Robert J. Glass, On the porous-continuum modeling of gravity-driven fingers in unsaturated materials: Extension of standard theory with a hold-back-pile-up effect, Water Resour. Res. 38 (11) (2002) 16-1, https://doi.org/10.1029/2001WR001131.

[7] Luis Cueto-Felgueroso, Ruben Juanes, A phase field model of unsaturated flow, Water Resour. Res. 45 (2009) W10409, https://doi.org/10.1029/2009WR007945.

[8] Yunwu Xiong, Flow of water in porous media with saturation overshoot: A review, J. Hydrol. 510 (2014) 353-362, https://doi.org/10.1016/j.jhydrol.2013.12.043.

[9] David A. DiCarlo, Stability of gravity-driven multiphase flow in porous media: 40 years of advancements, Water Resour. Res. 49 (2013) 4531-4544, https://doi.org/10.1002/wrcr.20359.

[10] Abdelaziz Beljadid, Luis Cueto-Felgueroso, Ruben Juanes, A continuum model of unstable infiltration in porous media endowed with an entropy function, Adv. Water Resour. 144 (2020) 103684, https://doi.org/10.1016/j.advwatres.2020.103684.

[11] W. Gardner, Some steady-state solutions of the unsaturated moisture flow equation with application to evaporation from a water table, Soil Sci. 85 (4) (1958) 228-232.

[12] R. Brooks, A. Corey, Hydraulic Properties of Porous Media: Hydrology Papers, Colorado State University, Fort Collins, Colorado, 1964.

[13] M.T. Van Genuchten, A closed-form equation for predicting the hydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. J. 44 (5) (1980) 892-898, https://doi.org/10.2136/sssj1980.03615995004400050002x.

[14] W.R. Gardner, Dynamic aspects of water availability to plants, Soil Sci. 89 (2) (1960) 63-73.

[15] F. Molz, I. Remson, A. Fungaroli, R. Drake, Soil moisture availability for transpiration, Water Resour. Res. 4 (6) (1968) 1161-1169, https://doi.org/10.1029/WR004i006p01161.

[16] F.J. Molz, Models of water transport in the soil-plant system: a review, Water Resour. Res. 17 (5) (1981) 1245-1260, https://doi.org/10.1029/WR017i005p01245.

[17] D. Hillel, V. Beek, H. Talpaz, A microscopic-scale model of soil water uptake and salt movement to plant roots, Soil Sci. 120 (5) (1975) 385-399.

[18] J. Hainsworth, L. Aylmore, Water extraction by single plant roots, Soil Sci. Soc. Am. J. 50 (4) (1986) 841-848, https://doi.org/10.2136/sssj1986.03615995005000040003x.

[19] M.T. Van Genuchten, A Numerical Model for Water and Solute Movement in and Below the Root Zone, United States Department of Agriculture Agricultural Research Service US, 1987.

[20] R. Albasha, J.-C. Mailhol, B. Chevron, Compensatory uptake functions in empirical macroscopic root water uptake models-experimental and numerical analysis, Agric. Water Manag. 155 (2015) 22-39, https://doi.org/10.1016/j.agwat.2015.03.010.

[21] F.J. Molz, I. Remson, Extraction term models of soil moisture use by transpiring plants, Water Resour. Res. 6 (5) (1970) 1346-1356, https://doi.org/10.1029/WR006i005p01346.
M. Boujoudar, A. Beljadid and A. Taik

Applied Numerical Mathematics xxx (xxxx) xxx

[22] S. Mathur, S. Rao, Modeling water uptake by plant roots, J. Irrig. Drain. Eng. 125 (3) (1999) 159-165, https://doi.org/10.1061/(ASCE)0733-9437(1999)125:3(159).

[23] P. de Willigen, J.C. van Dam, M. Javaux, M. Heinen, Root water uptake as simulated by three soil water flow models, Vadose Zone J. 11 (3) (2012), https://doi.org/10.2136/vzj2012.0018.

[24] P.J. Kramer, Plant and Soil Water Relationships: A Modern Synthesis, McGraw-Hill Book Company, 1969.

[25] Y. Mualem, A new model for predicting the hydraulic conductivity of unsaturated porous media, Water Resour. Res. 12 (3) (1976) 513-522, https://doi.org/10.1029/WR012i003p00513.

[26] F. Yuan, Z. Lu, Analytical solutions for vertical flow in unsaturated, rooted soils with variable surface fluxes, Vadose Zone J. 4 (4) (2005) 1210-1218, https://doi.org/10.2136/vzj2005.0043.

[27] O. Wilderotter, An adaptive numerical method for the Richards equation with root growth, Plant Soil 251 (2003) 255-267, https://doi.org/10.1023/A:1023031924963.

[28] T. Janz, R. Stonier, Modeling water flow in cropped soils: water uptake by plant roots, Environ. Int. 21 (5) (1995) 705-709, https://doi.org/10.1016/0160-4120(95)00069-W.

[29] G.J. Machado, R.M. Pereira, S. Clain, N. Araújo, S.O. Lopes, A new stabilised scheme for the Richards' equation with evapotranspiration, Groundw. Sustain. Dev. 17 (2022) 100736, https://doi.org/10.1016/j.gsd.2022.100736.

[30] M. Berardi, A. Andrisani, L. Lopez, M. Vurro, A new data assimilation technique based on ensemble Kalman filter and Brownian bridges: an application to Richards' equation, Comput. Phys. Commun. 208 (2016) 43-53, https://doi.org/10.1016/j.cpc.2016.07.025.

[31] H.T. Shahraiyni, B. Ataie-Ashtiani, Mathematical forms and numerical schemes for the solution of unsaturated flow equations, J. Irrig. Drain. Eng. 138 (1) (2012) 63-72, https://doi.org/10.1061/(ASCE)IR.1943-4774.0000377.

[32] Y. Zha, J. Yang, J. Zeng, C.-H.M. Tso, W. Zeng, L. Shi, Review of numerical solution of Richardson-Richards equation for variably saturated flow in soils, Wiley Interdiscip. Rev. Water 6 (5) (2019) e1364, https://doi.org/10.1002/wat2.1364.

[33] Y. Zha, M.C.-H. Tso, L. Shi, J. Yang, Comparison of noniterative algorithms based on different forms of Richards' equation, Environ. Model. Assess. 21 (2016) 357-370, https://doi.org/10.1007/s10666-015-9467-1.

[34] C. Chávez-Negrete, F. Domínguez-Mota, D. Santana-Quinteros, Numerical solution of Richards' equation of water flow by generalized finite differences, Comput. Geotech. 101 (2018) 168-175, https://doi.org/10.1016/j.compgeo.2018.05.003.

[35] W. Cai, J. Wang, K. Wang, Convergence analysis of Crank-Nicolson Galerkin-Galerkin FEMs for miscible displacement in porous media, J. Sci. Comput. 83 (2020) 1-26, https://doi.org/10.1007/s10915-020-01194-0.

[36] V. Baron, Y. Coudière, P. Sochała, Adaptive multistep time discretization and linearization based on a posteriori error estimates for the Richards equation, Appl. Numer. Math. 112 (2017) 104-125, https://doi.org/10.1016/j.apnum.2016.10.005.

[37] B. Cumming, T. Moroney, I. Turner, A mass-conservative control volume-finite element method for solving Richards' equation in heterogeneous porous media, BIT Numer. Math. 51 (2011) 845-864, https://doi.org/10.1007/s10543-011-0335-3.

[38] S. Keita, A. Beljadid, Y. Bourgault, Implicit and semi-implicit second-order time stepping methods for the Richards equation, Adv. Water Resour. 148 (2021) 103841, https://doi.org/10.1016/j.adwatres.2020.103841.

[39] F. List, F.A. Radu, A study on iterative methods for solving Richards' equation, Comput. Geosci. 20 (2) (2016) 341-353, https://doi.org/10.1007/s10596-016-9566-3.

[40] M.A. Celia, E.T. Bouloutas, R.L. Zarba, A general mass-conservative numerical solution for the unsaturated flow equation, Water Resour. Res. 26 (7) (1990) 1483-1496, https://doi.org/10.1029/WR026i007p01483.

[41] F. Lehmann, P. Ackerer, Comparison of iterative methods for improved solutions of the fluid flow equation in partially saturated porous media, Transp. Porous Media 31 (3) (1998) 275-292.

[42] C.K. Lee, X. Liu, S.C. Fan, Local multiquadric approximation for solving boundary value problems, Comput. Mech. 30 (5-6) (2003) 396-409, https://doi.org/10.1007/s00466-003-0416-5.

[43] M. Boujoudar, A. Beljadid, A. Taik, Modelling of unsaturated flow through porous media using meshless methods, in: Proceedings of the Canadian Society of Civil Engineering Annual Conference 2021, Springer, 2021, pp. 565-576.

[44] M. Boujoudar, A. Beljadid, A. Taik, Localized RBF methods for modeling infiltration using the Kirchhoff-transformed Richards equation, Eng. Anal. Bound. Elem. 152 (2023) 259-276, https://doi.org/10.1016/j.enganabound.2023.03.034.

[45] W. Madych, S. Nelson, Bounds on multivariate polynomials and exponential error estimates for multiquadric interpolation, J. Approx. Theory 70 (1) (1992) 94-114, https://doi.org/10.1016/0021-9045(92)90058-V.

[46] A.-D. Cheng, M. Golberg, E. Kansa, G. Zammito, Exponential convergence and H-c multiquadric collocation method for partial differential equations, Numer. Methods Partial Differ. Equ. 19 (5) (2003) 571-594, https://doi.org/10.1002/num.10062, An International Journal.

[47] A.-D. Cheng, Multiquadric and its shape parameter—a numerical investigation of error estimate, condition number, and round-off error by arbitrary precision computation, Eng. Anal. Bound. Elem. 36 (2) (2012) 220-239, https://doi.org/10.1016/j.enganabound.2011.07.008.

[48] G.E. Fasshauer, Meshfree Approximation Methods with MATLAB, vol. 6, World Scientific, 2007.

[49] S. Dharmaraja, An analysis of the TR-BDF2 integration scheme, Ph.D. thesis, Massachusetts Institute of Technology, 2007.

[50] A. Beljadid, A. Mohammadian, M. Charron, C. Girard, Theoretical and numerical analysis of a class of semi-implicit semi-Lagrangian schemes potentially applicable to atmospheric models, Mon. Weather Rev. 142 (12) (2014) 4458-4476.

[51] D. Illiano, I.S. Pop, F.A. Radu, Iterative schemes for surfactant transport in porous media, Comput. Geosci. 25 (2021) 805-822, https://doi.org/10.1007/s10596-020-09949-2.

[52] A. Dogan, L.H. Motz, Saturated-unsaturated 3D groundwater model. I: development, J. Hydrol. Eng. 10 (6) (2005) 492-504, https://doi.org/10.1061/(ASCE)1084-0699(2005)10:6(492).

[53] G. Manzini, S. Ferraris, Mass-conservative finite volume methods on 2-D unstructured grids for the Richards' equation, Adv. Water Resour. 27 (12) (2004) 1199-1215, https://doi.org/10.1016/j.adwatres.2004.08.008.

[54] M. Boujoudar, A. Beljadid, A. Taik, Localized MQ-RBF meshless techniques for modeling unsaturated flow, Eng. Anal. Bound. Elem. 130 (2021) 109-123, https://doi.org/10.1016/j.enganabound.2021.05.011.

[55] G.J. Hoffman, M.T. Van Genuchten, Soil properties and efficient water use: water management for salinity control, in: Limitations to Efficient Water Use in Crop Production, 1983, pp. 73-85.

[56] J. Vrugt, M.T. van Wijk, J. Hopmans, J. Šimunek, One-, two-, and three-dimensional root water uptake functions for transient modeling, Water Resour. Res. 37 (10) (2001) 2457-2470, https://doi.org/10.1029/2000WR000027.

[57] M.R. Kirkland, R. Hills, P. Wierenga, Algorithms for solving Richards' equation for variably saturated soils, Water Resour. Res. 28 (8) (1992) 2049-2058, https://doi.org/10.1029/92WR00082.

[58] K. Huang, B. Mohanty, M.T. Van Genuchten, A new convergence criterion for the modified Picard iteration method to solve the variably saturated flow equation, J. Hydrol. 178 (1) (1996) 69-91, https://doi.org/10.1016/0022-1694(95)02799-8.

[59] B. Šarler, R. Vertnik, Meshfree explicit local radial basis function collocation method for diffusion problems, Comput. Math. Appl. 51 (8) (2006) 1269-1282, https://doi.org/10.1016/j.camwa.2006.04.013.

[60] M. Li, W. Chen, C. Chen, The localized RBFs collocation methods for solving high dimensional PDEs, Eng. Anal. Bound. Elem. 37 (10) (2013) 1300-1304, https://doi.org/10.1016/j.enganabound.2013.06.001.
M. Boujoudar, A. Beljadid and A. Taik

Applied Numerical Mathematics xxx (xxxx) xxx

[61] J.L. Bentley, Multidimensional binary search trees used for associative searching, Commun. ACM 18 (9) (1975) 509-517, https://doi.org/10.1145/361002.361007.

[62] G. Yao, J. Duo, C. Chen, L. Shen, Implicit local radial basis function interpolations based on function values, Appl. Math. Comput. 265 (2015) 91-102, https://doi.org/10.1016/j.amc.2015.04.107.

[63] G.E. Fasshauer, J.G. Zhang, On choosing “optimal” shape parameters for RBF approximation, Numer. Algorithms 45 (1-4) (2007) 345-368, https://doi.org/10.1007/s11075-007-9072-8.

[64] R.L. Hardy, Multiquadric equations of topography and other irregular surfaces, J. Geophys. Res. 76 (8) (1971) 1905-1915, https://doi.org/10.1029/JB076i008p01905.

[65] C. Franke, R. Schaback, Solving partial differential equations by collocation using radial basis functions, Appl. Math. Comput. 93 (1) (1998) 73-82, https://doi.org/10.1016/S0096-3003(97)10104-7.

[66] A. Khoshfetrat, M. Abedini, Numerical modeling of long waves in shallow water using LRBF-DQ and hybrid DQ/LRBF-DQ, Ocean Model. 65 (2013) 1-10, https://doi.org/10.1016/j.ocemod.2013.01.006.

[67] C.A. Micchelli, Interpolation of scattered data: distance matrices and conditionally positive definite functions, in: Approximation Theory and Spline Functions, Springer, 1984, pp. 143-145.

[68] M.T. Musavi, W. Ahmed, K.H. Chan, K.B. Faris, D.M. Hummels, On the training of radial basis function classifiers, Neural Netw. 5 (4) (1992) 595-603, https://doi.org/10.1016/S0893-6080(05)80038-3.

[69] G. Garmanjani, R. Cavoretto, M. Esmaeilbeigi, A RBF partition of unity collocation method based on finite difference for initial-boundary value problems, Comput. Math. Appl. 75 (11) (2018) 4066-4090, https://doi.org/10.1016/j.camwa.2018.03.014.

[70] R.A. Feddes, P. Kowalik, K. Kolinska-Malinka, H. Zaradny, Simulation of field water uptake by plants using a soil water dependent root extraction function, J. Hydrol. 31 (1) (1976) 13-26, https://doi.org/10.1016/0022-1694(76)90017-2.

[71] J. Simunek, M.T. Van Genuchten, M. Sejna, The HYDRUS-1D software package for simulating the one-dimensional movement of water, heat, and multiple solutes in variably-saturated media, University of California-Riverside Research Reports 3 2005, pp. 1-240.

[72] A.W. Warrick, Soil Water Dynamics, Oxford University Press, 2003.

[73] P. Raats, Steady flows of water and salt in uniform soil profiles with plant roots, Soil Sci. Soc. Am. J. 38 (5) (1974) 717-722, https://doi.org/10.2136/sssaj1974.03615995003800050012x.

[74] Y. Rubin, D. Or, Stochastic modeling of unsaturated flow in heterogeneous soils with water uptake by plant roots: the parallel columns model, Water Resour. Res. 29 (3) (1993) 619-631, https://doi.org/10.1029/92WR02292.

[75] R. Srivastava, T.J. Yeh, Analytical solutions for one-dimensional, transient infiltration toward the water table in homogeneous and layered soils, Water Resour. Res. 27 (5) (1991) 753-762, https://doi.org/10.1029/90WR02772.

[76] M. Berardi, F.V. Difonzo, A quadrature-based scheme for numerical solutions to Kirchhoff transformed Richards' equation, J. Comput. Dyn. 9 (2) (2022) 69-84.

[77] J. Šimúnek, J.W. Hopmans, Modeling compensated root water and nutrient uptake, Ecol. Model. 220 (4) (2009) 505-521, https://doi.org/10.1016/j.ecolmodel.2008.11.004.

[78] S.M. Naghedifar, A.N. Ziaei, H. Ansari, Numerical analysis and optimization of triggered furrow irrigation system, Irrig. Sci. 38 (2020) 287-306, https://doi.org/10.1007/s00271-020-00672-5.