# LRBF meshless methods for predicting soil moisture distribution in root zone

Mohamed Boujoudar, Abdelaziz Beljadid, Ahmed Taik

► To cite this version:

Mohamed Boujoudar, Abdelaziz Beljadid, Ahmed Taik. LRBF meshless methods for predicting soil moisture distribution in root zone. Computers & Mathematics with Applications, 2025, 179, pp.178-196. 10.1016/j.camwa.2024.11.028. insu-04848960

HAL Id: insu-04848960

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

Submitted on 15 Jul 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
LRBF meshless methods for predicting soil moisture distribution
in root zone

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

<sup>a</sup> Mohammed VI Polytechnic University, Green city, Morocco

<sup>b</sup> University of Ottawa, Ottawa, Canada

<sup>c</sup> University of Hassan II -Casablanca, Morocco

**Abstract**

In this paper, we first propose a coupled numerical model of unsaturated flow in soils and plant root water uptake. The Richards equation and different formulations are used in the developed numerical model to describe infiltration in root zone and to investigate the impact of the plant root on the distribution of soil moisture. The Kirchhoff transformed Richards equation is used and the Gardner model is considered for capillary pressure. In our approach, we employ a meshless method based on localized radial basis functions (LRBF) to solve the resulting system of equations. The LRBF approach is an accurate and computationally efficient method that does not require mesh generation and is flexible in addressing high-dimensional problems with complex geometries. Furthermore, this method leads to a sparse matrix system, which avoids ill-conditioning issues. We implement the coupled numerical model of infiltration and plant root water uptake for one, two, and three-dimensional soils. Numerical experiments are performed using nontrivial analytical solutions and available experimental data to validate the coupled numerical model. The numerical results demonstrate the performance and ability of the proposed numerical method to predict soil moisture dynamics in root zone.

**Keywords:** Richards equation, Unsaturated water flow, Root water uptake, Kirchhoff transformation, Soil-water-plant interactions, Meshfree methods, Localized Radial Basis Function

**1. Introduction**

Studying water flow dynamics is crucial in modeling hydrological processes. Water up-
take by plant roots is a significant process in subsurface unsaturated flow modeling and has
a significant impact on the evolution of soil moisture and nutrient transport [1]. Predicting
moisture fluxes in unsaturated soils and their interaction with vegetation has practical im-
plications in agriculture, water management, and climate science [2]. Understanding root
water uptake process is needed to estimate crops water requirements and design best water

<footnote>Email address: abdelaziz.beljadid@um6p.ma (Abdelaziz Beljadid)</footnote>
management practices. For instance, optimum irrigation scheduling for efficient water man-
agement of crops requires a proper understanding of soil, water and plant interactions [3, 4].
Plant root water absorption is a complex process that involves various interactions between
soil, plant, and climate [1]. The interactions of these complex processes can be explained
using robust and efficient coupled numerical models.

The Richards equation [5] is the most commonly used model to describe water flow
in unsaturated soils due to the effects of gravity and capillarity. This equation is highly
nonlinear due to the nonlinear constitutive relations which describe the relationship among
capillary pressure, relative permeability, and saturation [6, 7, 8].

Several models have been developed to describe plant root water uptake. Molz [9] classi-
fied these models into two categories. The first is based on a microscopic approach, whereas
the second is based on a macroscopic one. The microscopic approach describes water ex-
traction at an individual root level and considers soil water flow toward the particular root,
predominantly radial [10, 11, 12]. The microscopic models are physically-based and are
efficiently used in various studies [13, 12, 14].

The macroscopic approach [15, 16, 17] considers properties of the overall root system. In
this case, the root water extraction is considered as a volumetric sink term in the Richards
equation. The macroscopic models are typically empirical [17, 18, 19]. Several macroscopic
models have been developed based on different approaches. For instance, Gardner [20]
proposed a macroscopic model describing the dynamic aspects of water availability to plants.
Feddes et al. [17] proposed a relatively simple empirical model to describe root water uptake.
Molz and Remsen [16] proposed a linear model where the extraction rate is assumed to vary
linearly with the soil depth. Prasad [21] proposed another linear model for root water uptake
which depends on the root depth and rate of evapotranspiration. In [21], the author showed
that the linear models give satisfactory results in predicting soil moisture in the root zone.
Li et al. [22] suggested an exponential root water uptake model derived from measured
root length distributions. In comparison with field data, they showed that the exponential
model [22] performs well compared to constant and linear models [17, 21]. Ojha and Rai
[23] proposed a nonlinear root water extraction model and investigated its performance on
irrigation scheduling compared to experimental data [24]. Several studies have favoured the
use of macroscopic models for their simplicity since they don't require information about
the detailed geometry of the root system. For more details, we refer to the previous studies
[16, 9, 25, 26].

Incorporating the plant root water uptake as a sink term in the Richard equation adds
more complexity to its analytical and numerical solutions. Except for some simple cases,
analytical solutions have been developed [27, 28, 29, 30, 31, 32, 33, 34]. For instance, Lomen
and Warrick [29] presented analytical solutions for the Richards equation under steady-state
conditions. Their solutions were based on the Gardner model for capillary pressure, and
they incorporated a sink term that depends on the matrix flux potential. Yuan and Lu [34]
developed exact solutions using the Gardner model [6] for the soil hydraulic conductivity and
different formulations for root water uptake sink term. Broadbridge et al. [35] developed
exact solutions for the Richards equation with a sink term nonlinearly dependent on soil
water content. In their solutions, they used the Gardner model for the capillary pressure.
Most of these analytical solutions assume a simplified model for root water uptake and consider specific initial and boundary conditions. In more practical situations, efficient numerical approaches in terms of computing cost and accuracy are required.

Various numerical approaches have been proposed to solve the Richards equation without a sink term [36, 37, 38, 39, 40]. Nevertheless, the design of appropriate numerical schemes for the Richards equation that includes a sink term remains challenging. Few numerical techniques have been proposed to solve the Richards equation considering root water uptake as a sink term. For instance, Neuman et al. [41] used a finite-element method to solve the soil water flow equation, taking into account water uptake by roots. Feddes et al. [42] proposed an implicit finite-difference method to solve the Richards equation considering plant root water uptake as sink term. Vrugt et al. [43] introduce the Galerkin finite element method to solve the Richards equation with a sink term of water uptake by roots. Rees and Ali [44] used the finite-element method to solve the Richards equation incorporating a sink term for root water extraction. Difonzo et al. [45] proposed a numerical approach based on the shooting method to solve the unsaturated flow equation with root water uptake models. Despite the studies conducted, there remains a need for more advanced numerical methods for computing water flow in unsaturated soils taking into account the absorption of water by the roots.

This study presents a coupled numerical model that accounts for both unsaturated soil flow and plant root water uptake. The proposed model uses the Richards equation and includes different root water uptake formulations as sink terms in the governing equation. The first formulation is proposed by Yoan et al. [34], where the stepwise and exponential functional forms are considered for the root water uptake. The second formulation is proposed by Broadbridge et al. [35], where root water uptake's sink term is assumed nonlinearly dependent on soil water content. In all considered formulations, the Gardner model [6] is employed for the capillary pressure.

In our approach, we employ a meshless method based on localized radial basis function (LRBF) to solve the resulting system of equations. The RBF meshless approach is an accurate and computationally efficient method that eliminates the need for mesh generation and is flexible in addressing high-dimensional problems with complex geometries. RBF meshless methods can be categorized into global [46, 47] and local [48, 49] methods. In the global approach, the collocation is performed globally over the entire computational domain, while the local approach performs collocation locally over a set of influence domains. Although the global method is simple to implement [46, 47], it faces two challenges which are: ill-conditioned resultant matrices and the issue of selecting an appropriate shape parameter for certain radial basis functions [50]. To overcome these drawbacks, LRBF methods are proposed [48, 51, 49]. These local approaches result in a sparse matrix that avoids ill-conditioning issues encountered in the global approach, where a full matrix is produced [49]. LRBF meshless methods have been effective in solving various partial differential equations, such as the Richards equation without the sink term [52, 53, 54, 37, 38]. To our best knowledge, no studies have explored the use of LRBF methods in solving the Richards equation that includes a sink term accounting for root water uptake.

In this paper, we propose LRBF meshless method [49] to solve the Richards equation with
the presence of a sink term of water uptake by plant roots. To deal with the nonlinearity of the governing equation, we employed the Kirchhoff transformation. To linearize the nonlinear problem resulting from the Kirchhoff transformation technique, Picard's iterations are used. The numerical results are validated and compared with nontrivial analytical solutions and available experimental data.

The structure of the paper is organized as follows: Section 2 introduces the coupled mathematical model of infiltration and plant root water uptake. Section 3 presents the proposed LRBF meshless method to solve the governing system. Numerical experiments for modeling soil moisture distribution in the root zone are performed in Section 4. Conclusions are provided in Section 5.

## 2. Mathematical model

### 2.1. The governing equations

The Richards equation [5], which governs water flow in the unsaturated zone while accounting for root water uptake with a sink term, is given by:

$$ \frac{\partial \theta}{\partial t} - \nabla \cdot (D(\theta) \nabla \theta) - \frac{\partial K(\theta)}{\partial z} = -R(\mathbf{x}, \theta), \quad \mathbf{x} \in \Omega, 0 < t \le T, \qquad (2.1) $$

where $\theta$ [$L^3/L^3$] is the volumetric water content, $D = Kdh/d\theta$ [$L^2/T$] is the soil-water diffusivity, $h$ [$L$] is the pressure head, $K$ [$L/T$] is the unsaturated hydraulic conductivity which is given by $K = K_s k_r$ where $K_s$ [$L/T$] is the saturated hydraulic conductivity and $k_r$ [–] is the water relative permeability, $R$ [$T^{-1}$] is the root water uptake, the coordinate vector $\mathbf{x} = (x, y, z)^T$ consists of the lateral directions $x$ [$L$] and $y$ [$L$], and the vertical direction $z$ [$L$], which is positive upward, $\Omega$ is an open bounded set of $\mathbb{R}^d$, where $d$ represents the dimension of the computational domain and $T$ is a fixed time.

The Richards equation is highly nonlinear due to the nonlinear constitutive relations which describe the relationship among capillary pressure, relative permeability and saturation [6, 8, 7]. The inclusion of the root-water uptake term in the Richards equation makes it further complicated in terms of numerical resolution. In this study, we simplify our analysis by considering a homogeneous soil. We adopt the assumption made by Broadbridge et al. [35], wherein the derivative of $K$ with respect to $\theta$ is expressed as follows:

$$ \frac{dK}{d\theta} = \alpha D, \qquad (2.2) $$

where $\alpha$ [$1/L$] is a pore size parameter, then the Richards equation becomes:

$$ \frac{\partial \theta}{\partial t} - \nabla \cdot (D(\theta) \nabla \theta) - \alpha D \frac{\partial \theta}{\partial z} = -R(\mathbf{x}, \theta). \qquad (2.3) $$

### 2.2. Root water uptake models

The sink term $R$ represents the rate of water uptake by plants, expressed as the volume of water removed per unit of time from a unit volume of soil. Different mathematical
models are developed to describe plant root water uptake, which are classified into two approaches: microscopic [10, 11, 12] and macroscopic [55, 56, 22, 16] approaches. Here, different macroscopic models are considered sink terms in the Richards equation to study the impact of plant root absorption on soil moisture distribution.

### 2.2.1. Model I

We present the first model used in [34] where the stepwise and exponential functional forms are considered for the root water uptake. In this case, the Gardner model [6] for the soil hydraulic conductivity function and a simple formulation [57] for the volumetric water content are used:

$$K = K_s \exp(\alpha h), \qquad (2.4)$$

$$\theta = \theta_r + (\theta_s - \theta_r) \exp(\alpha h), \qquad (2.5)$$

where $\theta_r$ [$L/L$] and $\theta_s$ [$L/L$] are the residual and saturated moisture content respectively. We consider an unsaturated soil zone stretching from the soil surface ($z = L$) to the water table ($z = 0$). As shown in Figure 1, $L_1$ denotes the maximum root depth [$L$], $q_1$ is time-dependent flux applied at the soil surface [$L/T$] and $ET$ is the evapotranspiration through the root zone [$L/T$]. Assuming $R_0$ represents the maximum uptake at the land surface [$T^{-1}$],

<figure><img src="image_2.png" /><figcaption>Figure 1: Schematic of vadose zone.</figcaption></figure>

the stepwise form of the root water uptake may be written as [34]:

$$R(z) = R_0 \delta(z - L_1), \qquad (2.6)$$

where $\delta(z - L_1)$ is the Heaviside function given by:

$$\delta(z - L_1) = \begin{cases} 1, & \text{if } L_1 \le z \le L, \\ 0, & \text{if } 0 \le z < L_1. \end{cases} \qquad (2.7)$$

Root water extraction can be represented using the exponential form [27, 30, 34] given by:

$$R(z) = R_0 \exp(\beta(z - L)), \qquad (2.8)$$

where $\beta$ [$L^{-1}$] is a parameter represents the rate of reduction of root uptake.
### 2.2.2. Model II

The second model is proposed by Broadbridge et al. [35]. In this case, water uptake by plant roots is assumed nonlinearly dependent on soil water content. The soil diffusivity, unsaturated hydraulic conductivity, and root water uptake rate are expressed as follows [35]:

$$ egin{cases} D = \frac{1}{\alpha^2 t_s} \frac{m}{(e^m - 1)} e^{m\Theta}, \\ K = K_s \frac{e^{m\Theta} - 1}{(e^m - 1)}, \\ R = \left(\frac{\theta_s}{t_s}\right) \left(\frac{-k}{(e^m - 1)}(e^{\Theta m} - 1) - \frac{A}{m}(1 - e^{-\Theta m})\right), \end{cases} \quad (2.9) $$

where $\Theta = (\theta - \theta_r)/(\theta_s - \theta_r)$ [-] is the normalized water content, $t_s = \theta_s/(\alpha K_s)$ [T] is the gravity time scale, $A$ [-], $m$ [-] and $k$ [-] are constant parameters depend on plant root water uptake rate and verify [35]:

$$ |k| < \frac{-A}{m} e^{-m}(1 - e^{-m}), \quad A < 0 \text{ and } m > 0. \quad (2.10) $$

At high moisture contents, $R$ is very high and close to a maximum value called the potential extraction rate $R_s$ [$T^{-1}$]. As soil moisture decreases and is near the wilting point, $R$ decreases as well and becomes very low and close to zero.

### 2.3. Kirchhoff-transformed governing equation

We use the Kirchhoff transformation to circumvent some difficulties encountered in solving Richards' equation due to its high nonlinearity. This transformation is given by:

$$ \mu = \int_{0}^{\theta} D(\theta) d\theta. \quad (2.11) $$

According to models presented before, this integral can be analytically calculated, which results in the following:

$$ \mu = \begin{cases} \frac{K_s}{(\theta_s - \theta_r)} \theta, & \text{Model I,} \\ \left(\frac{\theta_s - \theta_r}{\alpha^2 t_s}\right) \left(\frac{e^{m\Theta} - 1}{e^m - 1}\right), & \text{Model II.} \end{cases} \quad (2.12) $$

Using this transformation, Equation (2.3) can be simplified to:

$$ \frac{1}{D(\mu)} \frac{\partial \mu}{\partial t} - \nabla^2 \mu - \alpha \frac{\partial \mu}{\partial z} = -R(\mu). \quad (2.13) $$

The boundary conditions are transformed according to the Kirchhoff variable: For Dirichlet boundary conditions, we assume $\theta = \theta_d$, which leads to $\mu = \mu(\theta_d)$. In the case of Neumann boundary conditions, two cases are considered. First, when $-K \frac{\partial h}{\partial x^i} = q_x$, it
follows that $-\frac{\partial \mu}{\partial x^i} = q_x$. Second, when $-K \left[ \frac{\partial h}{\partial z} + 1 \right] = q_z$, it implies that $-\left[ \frac{\partial \mu}{\partial z} + \alpha \mu \right] = q_z$, where, $\theta_d$, $q_x$ and $q_z$ are given functions associated with the boundary conditions. The variable $i$ can take the values of 1, 2, which correspond to the horizontal spatial dimensions $(x^1, x^2) = (x, y)$. The high nonlinearity of the governing equation is reduced using the Kirchhoff transformation, which has the advantage of improving the convergence properties of the proposed numerical method [38, 40, 53]. Once the approximate values of $\mu$ are determined by solving Equation (2.13), we can deduce the estimated values of $\theta$ based on Equation (2.12):

$$ \theta = \begin{cases} \frac{(\theta_s - \theta_r)}{K_s} \mu, & \text{Model I,} \\ \theta_r + \frac{(\theta_s - \theta_r)}{m} \log \left( 1 + \alpha^2 t_s \left( \frac{e^m - 1}{\theta_s - \theta_r} \right) \mu \right), & \text{Model II.} \end{cases} \quad (2.14) $$

Evenly, one can deduce the approximate solutions of $h$ and $\Theta$ using the constitutive relations given by Equations (2.5) and (2.9).

### 3. Material and methodology

In this section, we present the numerical approach proposed to solve Equation (2.13). Specifically, we employ LRBF meshless techniques to solve the Richards equation, which includes a sink term to account for the root water uptake. For temporal discretization, we utilize an implicit scheme based on the backward differentiation formula (BDF). This choice ensures stability when solving Richards' equation and permits the use of reasonable time steps, which enhance the computational efficiency. To linearize the system, Picard's iterations are utilized.

#### 3.1. Linearization techniques

The backward Euler method of Equation (2.13) is given by:

$$ \frac{1}{D^{n+1}} \frac{\mu^{n+1} - \mu^n}{\Delta t} - \nabla^2 \mu^{n+1} - \alpha \frac{\partial \mu^{n+1}}{\partial z} = -R^{n+1}, \quad (3.1) $$

where $n \ge 0$, $\Delta t$ is the time step, $\mu^{n+1}$, $D^{n+1}$ and $R^{n+1}$ are the estimate values of $\mu$, $D$ and $R$ at time $t^{n+1} = (n+1)\Delta t$, respectively. Starting from a given initial condition $\mu^0 = \mu(\theta_{\{t=0\}})$, Equation (3.1) is linearized using the Picard iteration scheme [36], which is given by:

$$ \frac{1}{D^{n+1,m}} \frac{\mu^{n+1,m+1} - \mu^n}{\Delta t} - \nabla^2 \mu^{n+1,m+1} - \alpha \frac{\partial \mu^{n+1,m+1}}{\partial z} = -R^{n+1,m}, \quad (3.2) $$

where $m$ denotes iteration levels and $\mu^{n+1,0} = \mu^n$ is the initial guess of $\mu$. At $(m+1)^{th}$ iteration, Equation (3.1) is solved for $\mu^{n+1,m+1}$ using the values of $R^{n+1,m}$ and $D^{n+1,m}$ obtained
from the previous iteration. Subject to boundary and initial conditions, Equation (3.2) can be written as:

$$ \begin{cases} \mathcal{L}^m \mu^{n+1,m+1}(\boldsymbol{x}) = f^{n,m}(\boldsymbol{x}), & \boldsymbol{x} \in \Omega, \\ \mathcal{B} \mu^{n+1,m+1}(\boldsymbol{x}) = g(\boldsymbol{x}), & \boldsymbol{x} \in \partial\Omega, \\ \mu^0(\boldsymbol{x}) = \mu_0(\boldsymbol{x}), & \boldsymbol{x} \in \Omega, \end{cases} \qquad (3.3) $$

where $\mathcal{L}^m$ and $f^{n,m}$ are given by:

$$ \mathcal{L}^m = \left( \frac{1}{D^{n+1,m}} \frac{1}{\Delta t} \cdot \nabla^2 \cdot \alpha \frac{\partial}{\partial z} \right), \qquad (3.4) $$

$$ f^{n,m} = \left( \frac{\mu^n}{D^{n+1,m}} \frac{1}{\Delta t} - R^{n+1,m} \right). \qquad (3.5) $$

$\mathcal{L}^m$ denotes the differential operator, which is linearized during each Picard iteration ($m+1$), $\mathcal{B}^m$ is a linear border operator, which can be either Dirichlet or Neumann. At each time level $n$ and iteration ($m+1$), LRBF meshless method [48] is used to solve Equation (3.3) until the stop condition is satisfied, which is determined by:

$$ |\delta^{n+1,m+1}| = |\mu^{n+1,m+1} - \mu^{n+1,m}| < TOL, \qquad (3.6) $$

where $TOL$ is the error tolerance.

## 3.2. Meshless method based on LRBF

Let $\{\boldsymbol{x}_j\}_{j=1}^{N_i}$ be a set of $N_i$ distinct interior collocation points located in $\Omega$, and $\{\boldsymbol{x}_j\}_{j=N_i+1}^N$ be the boundary points, where $N$ is the total number of collocation points in $\Omega \cup \partial\Omega$. For each point $\boldsymbol{x}_s \in \Omega$, $s = 1, 2, \dots, N$, a localized influence domain $\Omega^{[s]}$ is created using the kd-algorithm [58]. It contains $n_s$ nearest neighbors interpolation points $\{\boldsymbol{x}_j^{[s]}\}_{j=1}^{n_s}$ to $\boldsymbol{x}_s$. Figure 2 shows three examples of local influence domains including $n_s = 3, 5$ and $9$. According to

<figure><img src="image_3.png" /><figcaption>Figure 2: Schematic of influence domains $\Omega^{[s]}$ for $n_s = 3, 5$, and $9$.</figcaption></figure>

previous research [49, 59, 40], the number of neighboring points $n_s$ needed for an accurate
scheme may vary depending on the dimension of the computational domain. For instance, Lee et al. [48] recommend to select 2 dim(Ω)+1 nodes in each influence domain. The approximate solution of $\mu_{[s]}^{m+1,n+1}$ can be expressed as a linear combination of RBFs within each influence domain:

$$ \mu_{[s]}^{n+1,m+1}(\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.7) $$

where $\{\lambda_k^{[s]^{n+1,m+1}}\}_{k=1}^{n_s}$ are unknown coefficients and $\varphi$ is a RBF. In our case, we employ the exponential RBF defined as:

$$ \varphi(r_k) = \exp(-\varepsilon^2 r_k^2), \quad (3.8) $$

where $r_k = \|\mathbf{x}_s - \mathbf{x}_k^{[s]}\|$ indicates the distance between $\mathbf{x}_s$ and $\mathbf{x}_k^{[s]}$, and $\varepsilon > 0$ is an arbitrary shape parameter. The exponential RBF is chosen due to its proven positive definiteness, as demonstrated in various studies [60, 61, 62]. This characteristic is crucial to guarantee that the resulting matrix will be non-singular [60, 61, 62, 63, 64]. The choice of an appropriate shape parameter is crucial for the accuracy of RBF meshless methods [65]. Various studies have proposed optimal choices for the shape parameter [66, 67]. In this study, the use of LRBF methods allows us to overcome the issue of selecting the optimal shape parameter, as shown in previous studies [48, 68], they are less sensitive to the choice of shape parameter compared to the global RBF methods. From Equation (3.7), we obtain the matrix form of the solution $\mu_{[s]}^{n+1,m+1}$:

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

where $\mu_{[s]}^{n+1,m+1} = [\mu_{[s]}^{n+1,m+1}(\mathbf{x}_1^{[s]}), \mu_{[s]}^{n+1,m+1}(\mathbf{x}_2^{[s]}), \dots, \mu_{[s]}^{n+1,m+1}(\mathbf{x}_{n_s}^{[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}^{[s]})]^T$ and $\varphi^{[s]} = [\varphi(\|\mathbf{x}_i^{[s]} - \mathbf{x}_j^{[s]}\|)]_{1 \le i,j \le n_s}$

is an $n_s \times n_s$ real symmetric coefficient matrix. The unknown coefficients $\lambda^{[s]^{n+1,m+1}}$ can be obtained as follow:

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

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

$$ \begin{aligned} \mathcal{L}^m \mu_{[s]}^{n+1,m+1}(\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} \mu_{[s]}^{n+1,m+1} = \Lambda_{[s]}^m \mu_{[s]}^{n+1,m+1}, \end{aligned} \quad (3.11) $$

where $\Psi^m = \mathcal{L}^m \varphi$, $\vartheta_{[s]}^m = [\Psi(\|\mathbf{x}_s - \mathbf{x}_1^{[s]}\|), \dots, \Psi(\|\mathbf{x}_s - \mathbf{x}_{n_s}^{[s]}\|)]$ and $\Lambda_{[s]}^m = \vartheta_{[s]}^m (\varphi^{[s]})^{-1}$. In order to reformulate Equation (3.11) in terms of the global $\mu^{n+1,m+1}$ instead of the local $\mu_{[s]}^{n+1,m+1}$, $\Lambda^m$ is considered as the expansion of $\Lambda_{[s]}^m$ obtained by padding the local vector with zeros at the proper positions. It follows that:

$$ \mathcal{L}^m \mu_{[s]}^{n+1,m+1}(\mathbf{x}_s) = \Lambda^m \mu^{n+1,m+1}, \quad (3.12) $$
where $\mu^{n+1,m+1} = [\mu^{n+1,m+1}(\mathbf{x}_1), \mu^{n+1,m+1}(\mathbf{x}_2), \dots, \mu^{n+1,m+1}(\mathbf{x}_N)]^T$. In the same way, we apply the linear operator $\mathcal{B}$ to Equation (3.7) for $\mathbf{x}_s \in \partial\Omega$:

$$ \begin{aligned} \mathcal{B}\mu_{[s]}^{m+1,n+1}(\mathbf{x}_s) &= \sum_{k=1}^{n_s} \lambda_k^{[s]^{m+1,n+1}} \mathcal{B}\varphi(\|\mathbf{x}_s - \mathbf{x}_k^{[s]}\|) = (\mathcal{B}\varphi^{[s]})\lambda_k^{[s]^{m+1,n+1}} \\ &= (\mathcal{B}\varphi^{[s]})(\varphi^{[s]})^{-1}\mu_{[s]}^{n+1,m+1} = \sigma^{[s]}\mu_{[s]}^{n+1,m+1} = \sigma\mu^{n+1,m+1}, \end{aligned} \quad (3.13) $$

where $\sigma^{[s]}$ is defined as $(\mathcal{B}\varphi^{[s]})(\varphi^{[s]})^{-1}$, and $\sigma$ is the expansion of $\sigma^{[s]}$ that is obtained by adding zeros at the appropriate positions. Combining Equations (3.12) and (3.13) into Equation (3.3), we obtain the following system:

$$ \begin{aligned} \mathcal{L}^m \mu^{n+1,m+1}(\mathbf{x}_s) &= \Lambda^m(\mathbf{x}_s)\mu^{n+1,m+1} = f^{m,n}(\mathbf{x}_s), \\ \mathcal{B}\mu^{n+1,m+1}(\mathbf{x}_s) &= \sigma(\mathbf{x}_s)\mu^{n+1,m+1} = g(\mathbf{x}_s). \end{aligned} \quad (3.14) $$

As a result, we obtain a sparse system:

$$ \left( \begin{array}{c} \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{array} \right) \left( \begin{array}{c} \mu^{n+1,m+1}(\mathbf{x}_1) \\ \mu^{n+1,m+1}(\mathbf{x}_2) \\ \vdots \\ \mu^{n+1,m+1}(\mathbf{x}_{N_i}) \\ \mu^{n+1,m+1}(\mathbf{x}_{N_i+1}) \\ \vdots \\ \mu^{n+1,m+1}(\mathbf{x}_N) \end{array} \right) = \left( \begin{array}{c} f^{n,m}(\mathbf{x}_1) \\ f^{n,m}(\mathbf{x}_2) \\ \vdots \\ f^{n,m}(\mathbf{x}_{N_i}) \\ g(\mathbf{x}_{N_i+1}) \\ \vdots \\ g(\mathbf{x}_N) \end{array} \right). \quad (3.15) $$

This localized approach leads to inverting a sparse matrix, which avoids ill-conditioning problems that occur in the full matrix generated using the global method. By solving this sparse system, we get the approximate solution $\mu^{n+1,m+1}$ at all given points. Once the condition given by Equation (3.6) is satisfied, we assign $\mu^{n+1} = \mu^{n+1,m+1}$.

## 4. Numerical simulations

This section presents numerical experiments to solve the Richards equation with a sink term, which represents the plant root water uptake. Our approach is based on LRBF meshless method based on the exponential function to solve the governing system (3.5). The first numerical test investigates soil water content dynamics during evaporation process. In the second example, root water uptake is described using model I. The third numerical test focuses on a 2D irrigation furrows system and aims to explore the influence of plant roots uptake on soil moisture distribution. The final numerical test applies the proposed numerical model to a 3D irrigation system from a circular source, evaluating its ability to predict 3D soil moisture distribution in the root zone. The numerical results are validated in comparison with non trivial analytical solutions and available experimental data. Note that, in the numerical tests where $z$ is considered positive downward, the same equation (2.1) is considered by replacing $z$ with $-z$.
### 4.1. Test 1: Soil water dynamics during evaporation process

In this numerical test, we consider a laboratory evaporation experience which is carried out by Teng et al. [69, 70]. A kind of soil in Japan named K-7 sand was selected in this experience to investigate soil water content dynamics during evaporation. The soil was filled into a cylinder with a height of 20 cm and a diameter of 10 cm and then was wetted to saturation. Three soil samples were subjected to different environmental conditions which include relative humidity, wind speed and temperature. These conditions were controlled in three cases which are summarized in [70]. Five water content probes were injected into the cylinder at depths of 1, 5, 10, 15 and 19 cm, respectively. The hydraulic parameters of the considered soil are $\alpha = 4.8 \ m^{-1}$, $K_s = 3.9 \times 10^{-6} \ m/s$, $\theta_r = 0$ and $\theta_s = 0.4$. The initial and boundary conditions corresponding to this experiment are:

$$ \begin{cases} \Theta(z, 0) = 1, \\ \Theta(0, t) = \exp(-\beta t), \end{cases} \qquad (4.1) $$

where $\beta$ is a positive constant given by $\beta = 4Db^2$ [70] and $b \le \alpha/4$ is a fitting parameter. The values of $b$ are 0.413, 0.489 and 0.37 for cases 1, 2, and 3, respectively. The exact solution associated to this numerical test is given by [70]:

$$ \Theta(z,t) = \frac{1}{2} \left[ \operatorname{erfc} \left( \frac{D\alpha t - z}{2\sqrt{Dt}} \right) - e^{\alpha z} \operatorname{erfc} \left( \frac{D\alpha t + z}{2\sqrt{Dt}} \right) \right] + \\ \frac{1}{2} e^{\left( \frac{\alpha z}{2} - 4b^2 Dt \right)} \left[ e^{2\alpha z} \operatorname{erfc} \left( \frac{z}{2\sqrt{Dt}} + 2c\sqrt{Dt} \right) + e^{-2\alpha z} \operatorname{erfc} \left( \frac{z}{2\sqrt{Dt}} - 2c\sqrt{Dt} \right) \right], \quad (4.2) $$

where erfc is the complementary error function and $c = \sqrt{(\alpha/4)^2 - b^2}$. We display in Figure 3 (left) a comparison of water content distributions at different times of measured, approximate and exact solutions.

We can observe that the drying rate is higher for the case where the $b$ value is large. It is evident that a higher drying rate results in a lower moisture content profile, whereas a slower drying rate persists for a more extended period. For instance, in the second case where $b = 0.489$, the water content takes nearly 200 h to reach 10%. In contrast, the third case requires a minimum of 300 h. The evaporation rates under the three cases are displayed in Figure 3 (right). The computed evaporation rate $E(t)$ by neglecting the gravity is given by [70]:

$$ -E(t) = -D \frac{\partial \theta}{\partial z} \Big|_{z=0} \qquad (4.3) $$

Scatter points, dashed and solid lines present the experimental profile, the exact and approximate solutions, respectively. A good agreement is observed in comparison to the exact solution. Similarly, a satisfactory correspondence with measured data is obtained, confirming the efficacy of the proposed numerical method in accurately predicting evaporation. For each case, both experimental and numerical results show two distinct evaporation phases. The constant rate stage persists for approximately 90, 60, and 180 hours in cases 1, 2, and 3, respectively. Subsequently, a falling rate phase occurs, followed by the initiation of a residual
<figure><img src="image_4.png" /></figure>

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

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

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

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

<figure><img src="image_9.png" /><figcaption>Figure 3: Left: measured and computed water content profile for three cases. Right: measured and computed evaporative rate. Symbols present the experimental profile, dashed lines are for exact solution and the solid lines are for approximate solution.</figcaption></figure>

falling-rate phase at 250, 200, and 370 hours, for cases 1, 2, and 3 respectively. The results are obtained using $\varepsilon = 0.4$, $n_s = 3$, $N_z = 200$ and $\Delta t = 0.01$. The root mean squared errors (RMSE) for water content in the three cases, obtained by comparing the approximate and exact solutions at $T = 600 \text{ h}$, $350 \text{ h}$, and $900 \text{ h}$, are $3.8 \times 10^{-8}$, $5.31 \times 10^{-8}$, and $1.98 \times 10^{-8}$ for cases 1, 2, and 3, respectively. These results show the efficacy of the LRBF technique
in providing accurate numerical solutions. The numerical results confirm that the proposed numerical model is able to predict the water content distribution during the evaporation process.

## 4.2. Test 2: Unsaturated flow in rooted soils under variable surface flux conditions

Here, we take into account the presence of the root water uptake in our model, and we use the stepwise (2.6) and exponential (2.7) forms as sink terms in the governing equation. We conduct numerical experiments of the distribution of water content, pressure head, and water flux through a rooted soil with a depth of $L = 100$ cm and a root depth of 40 cm, which corresponds to $L_1 = 60$ cm. We consider the following hydraulic parameters: $\theta_s = 0.45$, $\theta_r = 0.2$, and $K_s = 1$ cm h$^{-1}$ [71, 34]. The rate of reduction in root uptake $\beta$ is set at 0.04 m$^{-1}$, and the maximum water uptake $R_0$ is set at 0.02 h$^{-1}$ for $\alpha = 0.01$ cm$^{-1}$ and 0.0025 h$^{-1}$ for $\alpha = 0.1$ cm$^{-1}$, respectively. The following boundary and initial conditions are used:

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

where $h_0$ is the initial pressure head and $q_1$ is the time-dependent flux at the soil surface. Both steady-state and transient surface fluxes are considered in this numerical test.

We start with steady-state surface flux for the upper boundary condition where we consider a constant infiltration flux of $q_1 = -0.9$ cm h$^{-1}$. In this case the stepwise (2.4) formulation is used for root water uptake. We display in Figures 4 and 5 the distribution of the water content and pressure head for $\alpha = 0.01$ cm$^{-1}$ and $\alpha = 0.1$ cm$^{-1}$ respectively.

The numerical simulations are conducted over a duration of 50 hours. To investigate the impact of water uptake by plan roots on the distribution of soil moisture and pressure head, the time evolution of $h$ and $\theta$ is presented with (right) and without (left) root water uptake.

We can observe that the impact of the root water uptake is significant, particularly when $\alpha = 0.01$ m$^{-1}$ and $R_0 = 0.02$ h$^{-1}$. However, the change is relatively weak for the water content in the case of $\alpha = 0.1$ m$^{-1}$ and $R_0 = 0.0025$ h$^{-1}$ due to the low value of the maximum water uptake considered. In the first case, the time required to reach a steady state is nearly 50 hours in the absence of root water uptake and 30 hours in the presence of root water uptake. However, in the second case, the time required to attain a stable state is quite similar, at about 50 hours.

Now, we use a time-dependent surface flux for the upper boundary condition which is more realistic, because the soil surface conditions vary over time due to factors such as evaporation, rainfall and irrigation. In this case, the flux at the upper boundary is assumed to decrease exponentially as a function of time, represented by $q_1(t) = q_0 + \delta \exp(k_1 t)$, where $\delta = -0.8$ cm h$^{-1}$ and $k_1 = -0.1$ h$^{-1}$. The exponential formulation (2.8) is used for root water uptake. The evolution of water content in time and space for the rooted soils is shown in Figure 6. The left figure illustrates the rooted soil with $\alpha = 0.01$ cm$^{-1}$ and $R_0 = 0.02$ h$^{-1}$. The right figure corresponds to the case where $\alpha = 0.1$ cm$^{-1}$ and
<figure><img src="image_10.png" /><figcaption>Figure 4: Case 1: $\alpha = 0.01 \text{ cm}^{-1}$ and $R_0 = 0.02 \text{ h}^{-1}$. Comparison of the water content and pressure head results of the approximate and exact solutions. Left: with root water uptake. Right: without root water uptake.</figcaption></figure>

$R_0 = 0.0025 \text{ h}^{-1}$. Both soils are supplied with an equal amount of water from the surface, but they exhibit distinct moisture content patterns. The soil profile 1 characterized by $R_0 = 0.02 \text{ h}^{-1}$ and $\alpha = 0.01 \text{ m}^{-1}$ is on average wetter than the soil 2 with $R_0 = 0.0025 \text{ h}^{-1}$ and $\alpha = 0.1 \text{ m}^{-1}$ even though the volume of water absorbed by the roots is greater than
<figure><img src="image_11.png" /><figcaption>Figure 5: Case 2: $\alpha = 0.1 \text{ cm}^{-1}$ and $R_0 = 0.0025 \text{ h}^{-1}$. Comparison of the water content and pressure head results of the approximate and exact solutions. Left: with root water uptake. Right: without root water uptake.</figcaption></figure>

that received from the soil surface. This difference can be attributed to the soil profile 1 favoring capillary rise, which facilitates water transfer from the water table into the root zones.

To assess the transient water flow responding to surface flux changes, we perform nu-
<figure><img src="image_12.png" /><figcaption>Figure 6: Time evolution of soil water content.</figcaption></figure>

merical simulations of the water flux $q_2(t)$ at the interface between the root zone and subsoil ($z = L_1$), as well as the flow $q_3(t)$ at the water table ($z = 0$). Figure 7 displays the flows $q_2(t)$ and $q_3(t)$ for both constant (left) and varying surface flux (right).

<figure><img src="image_13.png" /><figcaption>Figure 7: Time evolution of water flows at the interface between the root zone and subsoil $q_2(t)$ and at the water table $q_3(t)$. Left: constant surface flux. Right: varying surface flux.</figcaption></figure>

As the time reaches 50 hours, $q_2$ and $q_3$ approach 0.5 cm h⁻¹ for the constant surface flux and near 0.6 cm h⁻¹ for the transient surface flux. The results are obtained using $ε = 0.2$, $n_s = 5$, $N_z = 1001$ and $Δt = 0.005$. In all the studied cases, a good agreement is observed between the approximate solutions and the exact solutions which are developed in [34]. The root mean square errors (RMSE) of the water content for these cases are shown in Figure 8. The error values are very small for both constant and varying surface fluxes which highlight the accuracy of the proposed LRBF method. This efficiency is particularly important in cases where the surface fluxes may change over time or under different environmental conditions, as it ensures that the proposed method can provide reliable results under a variety of scenarios.
<figure><img src="image_14.png" /><figcaption>Figure 8: RMSE errors of the water content. Left: constant surface flux. Right: varying surface flux.</figcaption></figure>

### 4.3. Test 3: 2D flow from periodic irrigation furrows

In this numerical test, we consider surface irrigation furrows with an irrigation rate of $F_0 = 4 \text{ cm/h}$. The furrows have a width of $2x_0$ and are periodically placed apart with a spatial period of $2l > 2x_0$, as illustrated in Figure 9 (left). The prescribed uniform vertical

<figure><img src="image_15.png" /><figcaption>Figure 9: Schematic of the periodic irrigation furrows (left). The computational domain (right).</figcaption></figure>

flux through the furrows $V_0(t) = F_0 e^{(A/t_s)t}$ ensures water infiltration into the soil at a specified rate, while the zero flux through the surface between furrows, which is presumed to be protected by mulch, maintains soil moisture by preventing water movement. These boundary conditions at the soil surface exhibit symmetry due to the periodic placement of the irrigation furrows and the uniform behavior of water infiltration and conservation between these furrows. Owing to the symmetry present in this problem, we can simplify our study by restricting the computational domain to a rectangular region $[0, l]$ as shown in Figure 9 (right). The horizontal water movement must be absent along the lines of symmetry at $x = 0$ and $x = l$, maintaining zero flux across these planes.

We use Model II (2.9) for plant root water uptake, where the sink term depends on soil depth and water content. This dependency increases the nonlinearity of the system, resulting in a more complex numerical test compared to previous tests. We select two distinct plants with their respective parameters [35]: Plant 1 with $m = 3$, $A = -0.0021$, and $k = -3.31 \times 10^{-5}$, and plant 2 with $m = 5$, $A = -0.007$, and $k = -3.18 \times 10^{-8}$.
We select the Brindabella silty clay loam soil [72, 73], characterized by the parameters $\theta_r = 0$, $\theta_s = 0.485$, $K_s = 12 \text{ cm h}^{-1}$ and $\alpha = 0.142 \text{ cm}^{-1}$ to simulate water content dynamics, taking into account root water extraction. The chosen values for the parameters $l$ and $x_0$ are $4ls$ and $l_s$ respectively, where $l_s = 7 \text{ cm}$ [35].

The boundary conditions associated to this numerical test can be expressed as follows [35]:

$$ egin{cases} 
-K \left( \frac{\partial h}{\partial z} - 1 \right) = V_0(t), & z = 0 \text{ and } 0 \le x \le x_0, \\ 
-K \left( \frac{\partial h}{\partial z} - 1 \right) = 0, & z = 0 \text{ and } x_0 < x < l, \\ 
-K \frac{\partial h}{\partial x} = 0, & x = 0, x = l, \\ 
-K \left( \frac{\partial h}{\partial z} - 1 \right) = e^{(A/t_{st})} \left( \frac{x_0 F_0}{l} \right) e^{-\alpha(\sqrt{1-4k-1}z/2)}, & z = L. 
\end{cases} 
\quad (4.5) $$

The analytical solution in this case is given by [35]:

$$ \mu = \left( \frac{\theta_s}{\alpha^2 t_s} \right) e^{(A/t_{st})} \Phi(x, z), \quad (4.6) $$

where, $\Phi$ is given by:

$$ \Phi(x, z) = \sum_{j=1}^{\infty} 2A_j \cos \left( \frac{j\pi x}{l} \right) \frac{e^{-\alpha(\sqrt{1+4(j\pi l_s/l)^2 - 4k-1}z/2)}}{1 + \sqrt{1 + 4(j\pi l_s/l)^2 - 4k}}, \quad (4.7) $$

where $A_0 = F_0 x_0 / (K_s l)$ and $A_j = 2F_0 / (K_s j \pi) \sin(j\pi x_0 / l)$. According to Equation (2.14), $\Theta$ can be expressed as follows:

$$ \Theta = \frac{1}{m} \log \left( 1 + \alpha^2 t_s \left( \frac{e^m - 1}{\theta_s - \theta_r} \right) \mu \right). \quad (4.8) $$

In our numerical experiments, we set $\varepsilon = 0.5$, $\Delta t = 0.001$, $n_s = 5$, $N_x = 1000$ and $N_z = 2000$. The evolution of $\Theta$ near the surface ($z^* = z/l_s = 4$) and deeper in the soil ($z^* = z/l_s = 8.5$), associated with plant 1 and plant 2 respectively, is presented in Figures 9 and 10. These figures illustrate the changes in soil moisture levels at different soil depths as influenced by the distinct root water uptake characteristics of the two plants.

We can notice that the soil associated with plant 2 is, on average, more humid than that of plant 1, which shows that the plant corresponding to soil 1 absorbs water at a higher rate than the plant corresponding to soil 2. This observation is confirmed by the values of $\Theta$, which are greater for plant 2 than for plant 1. This difference can be attributed to the distinct values of $k$, $m$, and $A$ chosen for each plant, emphasizing the effect of the plant root system on the soil moisture distribution. For both plant 1 and plant 2, a comparison of the contour plots of $\Theta$ between the approximate and exact solutions is displayed in Figures 9
<figure><img src="image_16.png" /><figcaption>Figure 10: Evolution of Θ near the surface (z* = 4) (a), and deep in the soil (z* = 8.5) (b): Approximate solution (left), exact solution (middle), and comparison contours (right).</figcaption></figure>

<table><caption>Table 1: RMSE errors between approximate and exact solutions.</caption><thead><tr><th>Plant</th><th>N<sub>z</sub></th><th>N<sub>x</sub></th><th>RMSE</th></tr></thead><tbody><tr><td rowspan="3">1</td><td>2000</td><td>1000</td><td>6.35 &times; 10<sup>-5</sup></td></tr><tr><td>3000</td><td>1500</td><td>4.23 &times; 10<sup>-5</sup></td></tr><tr><td>4000</td><td>2000</td><td>3.18 &times; 10<sup>-5</sup></td></tr><tr><td rowspan="3">2</td><td>2000</td><td>1000</td><td>3.73 &times; 10<sup>-5</sup></td></tr><tr><td>3000</td><td>1500</td><td>2.49 &times; 10<sup>-5</sup></td></tr><tr><td>4000</td><td>2000</td><td>1.87 &times; 10<sup>-5</sup></td></tr></tbody></table>

and 10, respectively. The results show a good correspondence between the two solutions, validating the reliability of the model used in this study.

Table 1 presents the RMSE errors between the approximate and exact solutions of Θ for both plant 1 and plant 2. We fix the values of Δt, ε, and n<sub>s</sub>, and then vary the number of nodes distributed in the computational domain. The error values, which are small for
<figure><img src="image_17.png" /><figcaption>Figure 11: Evolution of Θ near the surface ($z^* = 1.5$) (a), and deep in the soil ($z^* = 9$) (b): Approximate solution (left), exact solution (middle), and comparison contours (right).</figcaption></figure>

both plants, diminish as the number of nodes increases. This confirms the convergence and accuracy of the proposed LRBF method in solving the governing equation. The results demonstrate the effectiveness of our approach in simulating soil moisture dynamics while taking into account root water uptake.

### 4.4. Test 4: 3D irrigation from a circular source

In this three-dimensional numerical test, we extend the previous analysis on furrow irrigation to assess the effectiveness of the proposed numerical model in predicting 3D soil moisture distribution in the root zone under cylindrical coordinates. We consider a cylindrical region with a length of $L = 56$ cm and a radius of $R = 7$ cm [35], and assume that the irrigation source is a circular region of radius $r_0$ located at the soil surface, as illustrated in Figure 12. At the top of the cylinder, a prescribed flux boundary $V_0(t)$ is imposed at the circular region of radius $r_0$, and no-flow boundary condition is applied outside of this circular region. Notably, at the lateral sides of the cylinder (i.e., at $r = R$), no-flow boundary
<figure><img src="image_18.png" /><figcaption>Figure 12: The 3D computational domain.</figcaption></figure>

condition is imposed. The boundary conditions described above can be written as follows:

$$ 
\left\{ 
\begin{array}{l} 
-K \left( \frac{\partial h}{\partial z} - 1 \right) = V_0(t), \quad z = 0 \text{ and } r \le r_0, \\ 
-K \left( \frac{\partial h}{\partial z} - 1 \right) = 0, \quad z = 0 \text{ and } r > r_0, \\ 
-K \frac{\partial h}{\partial r} = 0, \quad r = R. 
\end{array} 
\right. 
\qquad (4.9) 
$$

We consider model II for plant root water uptake. The analytical solution associated with this numerical test in case of axisymmetric flow is given by [35]:

$$ \mu = \left( \frac{\theta_s}{\alpha^2 t_s} \right) e^{A/t_s t} \Phi(r), \qquad (4.10) $$

where

$$ \Phi = \frac{\alpha r_0 F_0}{K_s} \int_0^\infty \frac{2J_1(\alpha v r_0) J_0(\alpha v r)}{1 + \sqrt{1 + 4(v^2 - k)}} e^{-\alpha(\sqrt{1+4(v^2-k)}-1)z/2} dv, \qquad (4.11) $$

where $J_0$ and $J_1$ are Bessel functions of the first kind, for orders 0 and 1, respectively. $\Theta$ is calculated using Equation (4.8). With the same parameters of soil and plants taken in the previous test, we show in Figures 13 and 14 the evolution of saturation for $r_0^* = r_0/l_s = 1/5$ and 1/2 for both plants 1 and 2, respectively.

As in the 2D results, we observe that the soil with plant 2 remains more saturated than that with plant 1, which indicates the significant impact of the corresponding root water uptake parameters on soil moisture distribution. The impact of the radius $r_0$ of the circular irrigation source is also evident from the numerical results. As $r_0$ increases, the soil saturation increases accordingly, which is expected since a larger area is irrigated. This effect is more pronounced near the soil surface, where the irrigation is applied. However,
<figure><img src="image_19.png" /><figcaption>Figure 13: Approximate solution of the evolution of $\Theta$ corresponding to plant 1 for $r_0^* = 1/5$ and $r_0^* = 1/2$.</figcaption></figure>

the impact of $r_0$ on the soil moisture distribution decreases with depth as the soil becomes less affected by the irrigation source.

Table 2 presents the computed errors for $r_0^* = 1/2$ with varying numbers of collocation points, indicating that the accuracy of the proposed method increases with the increasing of the number of points. This highlights the effectiveness of the numerical model in accurately predicting 3D soil moisture dynamics in the root zone. Moreover, the results demonstrate the capability of the proposed method in handling the complex geometry of the problem and the ability to capture the dynamics of soil moisture distribution under cylindrical coordinates.

The proposed model offers a robust and efficient numerical approach for simulating soil moisture dynamics in root zone, which can be used to inform irrigation management strategies and improve crop yields. The efficiency of the LRBF method makes it a valuable tool for investigating the impact of different irrigation scenarios and soil properties on soil moisture
<figure><img src="image_20.png" /><figcaption>Figure 14: Approximate solution of the evolution of $\Theta$ corresponding to plant 2 for $r_0^* = 1/5$ and $r_0^* = 1/2$.</figcaption></figure>

<table><caption>Table 2: RMSE errors between approximate and exact solutions.</caption><thead><tr><th>Plant</th><th>N<sub>z</sub></th><th>N<sub>r</sub></th><th>Δt</th><th>ε</th><th>n<sub>s</sub></th><th>RMSE</th></tr></thead><tbody><tr><td rowspan="2">1</td><td>500</td><td>235</td><td>0.01</td><td>0.5</td><td>7</td><td>4.94 × 10<sup>-4</sup></td></tr><tr><td>1000</td><td>125</td><td>0.01</td><td>0.5</td><td>7</td><td>1.62 × 10<sup>-4</sup></td></tr><tr><td rowspan="2">2</td><td>500</td><td>62</td><td>0.01</td><td>0.3</td><td>7</td><td>5.03 × 10<sup>-4</sup></td></tr><tr><td>1000</td><td>125</td><td>0.01</td><td>0.5</td><td>7</td><td>5.83 × 10<sup>-5</sup></td></tr></tbody></table>

distribution.
## 5. Conclusion

In this paper, we focus on modeling soil moisture distribution in the root zone and introducing computational techniques for efficiently solving the coupled numerical model of infiltration in soils and plant root water uptake. The proposed model is based on the Richards equation and different formulations for root water extraction to study the impact of plant root water uptake on the soil moisture distribution. The Gardner model is used for the capillary pressure. Incorporating an implicit sink term of root water uptake in the system increases the complexity in terms of numerical resolution. The Kirchhoff transformation of the governing equation is employed to reduce the nonlinearity of the system. Our numerical approach is based on LRBF meshless method to solve the coupled numerical model. This technique provides significant computational advantages, including reduced computational cost and accurate numerical solutions since it does not require mesh generation. Our proposed coupled model is validated through various numerical tests, including simulations of soil moisture distribution during evaporation and root water uptake processes in one-, two-, and three-dimensional cases. The reliability of our approach is validated against both experimental data and non-trivial exact solutions. The numerical results show the capability of the proposed method to effectively predict the dynamics of unsaturated flow through soils under evaporation and plant root water absorption. Our study offers a robust numerical framework for studying the impact of root water uptake on soil moisture dynamics, which is important for agriculture and ecosystem management. The proposed numerical techniques could be further extended to incorporate more complex root water uptake models, which can provide a promising avenue for future research.

## Acknowledgments

AB gratefully acknowledges funding from UM6P-OCP. AT and AB gratefully acknowledge funding for the APRD research program from the Moroccan Ministry of Higher Education, Scientific Research and Innovation and the OCP Foundation.

## References

[1] J. Vimuunek, J. W. Hopmans, Modeling compensated root water and nutrient uptake, Ecological modelling 220 (4) (2009) 505–521.

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

[3] S. R. Green, M. Kirkham, B. E. Clothier, Root uptake and transpiration: From measurements and models to sustainable irrigation, agricultural water management 86 (1-2) (2006) 165–176.

[4] R. Kumar, M. Jat, V. Shankar, Methods to estimate irrigated reference crop evapotranspiration-a review, Water Science and Technology 66 (3) (2012) 525–535.

[5] L. A. Richards, Capillary conduction of liquids through porous mediums, Physics 1 (5) (1931) 318–333.

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

[7] R. Brooks, T. Corey, Hydraulics of porous media, Hydrology Papers, Colorado State University 24 (1964) 37.
[8] M. T. Van Genuchten, A closed-form equation for predicting the hydraulic conductivity of unsaturated soils, Soil science society of America journal 44 (5) (1980) 892-898.

[9] F. J. Molz, Models of water transport in the soil-plant system: A review, Water resources research 17 (5) (1981) 1245-1260.

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

[11] F. Molz, I. Remson, A. Fungaroli, R. Drake, Soil moisture availability for transpiration, Water resources research 4 (6) (1968) 1161-1169.

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

[13] F. J. Molz, Interaction of water uptake and root distribution 1, Agronomy Journal 63 (4) (1971) 608-610.

[14] J. Hainsworth, L. Aylmore, Water extraction by single plant roots, Soil Science Society of America Journal 50 (4) (1986) 841-848.

[15] G. Ogata, L. Richards, W. R. Gardner, Transpiration of alfalfa determined from soil water content changes, soil Science 89 (4) (1960) 179-182.

[16] F. J. Molz, I. Remson, Extraction term models of soil moisture use by transpiring plants, Water Resources Research 6 (5) (1970) 1346-1356.

[17] R. A. Feddes, P. J. Kowalik, H. Zaradny, Simulation of field water use and crop yield / R. A. Feddes, P. J. Kowalik, and H. Zaradny., Simulation monographs, Pudoc for the Centre for Agricultural Publishing and Documentation, Wageningen, 1978.

[18] M. v. Genuchten, et al., A numerical model for water and solute movement in and below the root zone, Res. Rep. N 121. US Salinity laboratory, USDA, ARS. California, Riverside. 221.

[19] R. Albasha, J.-C. Mailhol, B. Cheviron, Compensatory uptake functions in empirical macroscopic root water uptake models-experimental and numerical analysis, Agricultural Water Management 155 (2015) 22-39.

[20] W. Gardner, Relation of root distribution to water uptake and availability 1, Agronomy Journal 56 (1) (1964) 41-45.

[21] R. Prasad, A linear root water uptake model, Journal of Hydrology 99 (3-4) (1988) 297-306.

[22] K. Li, J. Boisvert, R. D. Jong, An exponential root-water-uptake model, Canadian Journal of soil science 79 (2) (1999) 333-343.

[23] C. S. P. Ojha, A. K. Rai, Nonlinear root-water uptake model, Journal of irrigation and drainage engineering 122 (4) (1996) 198-202.

[24] L. J. Erie, O. F. French, K. Harris, Consumptive use of water by crops in Arizona, College of Agriculture, University of Arizona (Tucson, AZ), 1965.

[25] S. Mathur, S. Rao, Modeling water uptake by plant roots, Journal of Irrigation and Drainage Engineering 125 (3) (1999) 159-165.

[26] P. de Willigen, J. C. van Dam, M. Javaux, M. Heinen, Root water uptake as simulated by three soil water flow models, Vadose Zone Journal 11 (3) (2012) vzj2012-0018.

[27] P. Raats, Steady flows of water and salt in uniform soil profiles with plant roots, Soil Science Society of America Journal 38 (5) (1974) 717-722.

[28] P. Raats, Analytical solutions of a simplified flow equation, Transactions of the ASAE 19 (4) (1976) 683-0689.

[29] D. Lomen, A. Warrick, Solution of the one-dimensional linear moisture flow equation with implicit water extraction functions, Soil Science Society of America Journal 40 (3) (1976) 342-344.

[30] Y. Rubin, D. Or, Stochastic modeling of unsaturated flow in heterogeneous soils with water uptake by plant roots: The parallel columns model, Water Resources Research 29 (3) (1993) 619-631.

[31] P. Broadbridge, The forced burgers equation, plant roots and schrödinger's eigenfunctions, Journal of engineering mathematics 36 (1) (1999) 25-39.

[32] H. Basha, Multidimensional linearized nonsteady infiltration toward a shallow water table, Water Resources Research 36 (9) (2000) 2567-2573.

[33] G. Schoups, J. Hopmans, Analytical model for vadose zone solute transport with root water and solute
uptake, Vadose Zone Journal 1 (1) (2002) 158–171.

[34] F. Yuan, Z. Lu, Analytical solutions for vertical flow in unsaturated, rooted soils with variable surface fluxes, Vadose Zone Journal 4 (4) (2005) 1210–1218.

[35] P. Broadbridge, E. Daly, J. Goard, Exact solutions of the Richards equation with nonlinear plant-root extraction, Water Resources Research 53 (11) (2017) 9679–9691.

[36] M. A. Celia, E. T. Bouloutas, R. L. Zarba, A general mass-conservative numerical solution for the unsaturated flow equation, Water resources research 26 (7) (1990) 1483–1496.

[37] M. Boujoudar, A. Beljadid, A. Taik, Localized MQ-RBF meshless techniques for modeling unsaturated flow, Engineering Analysis with Boundary Elements 130 (2021) 109–123.

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

[39] S. Keita, A. Beljadid, Y. Bourgault, Implicit and semi-implicit second-order time stepping methods for the Richards equation, Advances in Water Resources 148 (2021) 103841.

[40] M. Boujoudar, A. Beljadid, A. Taik, Localized RBF methods for modeling infiltration using the Kirchhoff-transformed Richards equation, Engineering Analysis with Boundary Elements 152 (2023) 259–276.

[41] S. P. Neuman, R. A. Feddes, E. Bresler, Finite element analysis of two-dimensional flow in soils considering water uptake by roots: I. theory, Soil Science Society of America Journal 39 (2) (1975) 224–230.

[42] 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, Journal of Hydrology 31 (1-2) (1976) 13–26.

[43] J. Vrugt, M. T. van Wijk, J. Hopmans, J. Šimunek, One-, two-, and three-dimensional root water uptake functions for transient modeling, Water Resources Research 37 (10) (2001) 2457–2470.

[44] S. W. Rees, N. Ali, Seasonal water uptake near trees: A numerical and experimental study, Geomechanics and Geoengineering: An International Journal 1 (2) (2006) 129–138.

[45] F. V. Difonzo, C. Masciopinto, M. Vurro, M. Berardi, Shooting the numerical solution of moisture flow equation with root water uptake models: A python tool, Water Resources Management 35 (8) (2021) 2553–2567.

[46] E. J. Kansa, Multiquadrics—a scattered data approximation scheme with applications to computational fluid-dynamics—ii solutions to parabolic, hyperbolic and elliptic partial differential equations, Computers & mathematics with applications 19 (8-9) (1990) 147–161.

[47] E. J. Kansa, Multiquadrics—a scattered data approximation scheme with applications to computational fluid-dynamics—i surface approximations and partial derivative estimates, Computers & Mathematics with applications 19 (8-9) (1990) 127–145.

[48] C. K. Lee, X. Liu, S. C. Fan, Local multiquadric approximation for solving boundary value problems, Computational Mechanics 30 (5-6) (2003) 396–409.

[49] M. Li, W. Chen, C. Chen, The localized RBFs collocation methods for solving high dimensional PDEs, Engineering Analysis with Boundary Elements 37 (10) (2013) 1300–1304.

[50] U. Sundin, Global radial basis function collocation methods for pdes, Ph.D. thesis, Uppsala University (2020).

[51] B. Šarler, From global to local radial basis function collocation method for transport phenomena, in: Advances in meshfree techniques, Springer, 2007, pp. 257–282.

[52] D. Stevens, H. Power, M. Lees, H. Morvan, A meshless solution technique for the solution of 3d unsaturated zone problems, based on local hermitian interpolation with radial basis functions, Transport in porous media 79 (2) (2009) 149–169.

[53] D. Stevens, H. Power, A scalable and implicit meshless RBF method for the 3d unsteady nonlinear Richards equation with single and multi-zone domains, International journal for numerical methods in engineering 85 (2) (2011) 135–163.

[54] E. Ben-Ahmed, M. Sadik, M. Wakrim, Radial basis function partition of unity method for modelling water flow in porous media, Computers & Mathematics with Applications 75 (8) (2018) 2925–2941.

[55] J. Hoogland, R. A. Feddes, C. Belmans, Root water uptake model depending on soil water pressure head and maximum extraction rate, in: III International Symposium on Water supply and Irrigation
in the open and under Protected Cultivation 119, 1981, pp. 123–136.

[56] R. A. Feddes, Simulation of field water use and crop yield, in: Simulation of plant growth and crop production, Pudoc, 1982, pp. 194–209.

[57] A. W. Warrick, Soil water dynamics, Oxford University Press, 2003.

[58] J. L. Bentley, Multidimensional binary search trees used for associative searching, Communications of the ACM 18 (9) (1975) 509–517.

[59] M. Hamaidi, A. Naji, A. Charafi, Space-time localized radial basis function collocation method for solving parabolic and hyperbolic equations, Engineering Analysis with Boundary Elements 67 (2016) 152–163.

[60] 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. doi:10.1007/978-94-009-6466-2_7.

[61] A.-D. Cheng, Multiquadric and its shape parameter—A numerical investigation of error estimate, condition number, and round-off error by arbitrary precision computation, Engineering analysis with boundary elements 36 (2) (2012) 220–239. doi:https://doi.org/10.1016/j.enganabound.2011.07.008.

[62] G. E. Fasshauer, Meshfree approximation methods with MATLAB, Vol. 6, World Scientific, 2007.

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

[64] G. Garmanjani, R. Cavoretto, M. Esmaeilbeigi, A RBF partition of unity collocation method based on finite difference for initial-boundary value problems, Computers & Mathematics with Applications 75 (11) (2018) 4066–4090. doi:https://doi.org/10.1016/j.camwa.2018.03.014.

[65] G. E. Fasshauer, J. G. Zhang, On choosing “optimal” shape parameters for RBF approximation, Numerical Algorithms 45 (1-4) (2007) 345–368.

[66] R. L. Hardy, Multiquadric equations of topography and other irregular surfaces, Journal of geophysical research 76 (8) (1971) 1905–1915.

[67] C. Franke, R. Schaback, Solving partial differential equations by collocation using radial basis functions, Applied Mathematics and Computation 93 (1) (1998) 73–82.

[68] A. Khoshfetrat, M. Abedini, Numerical modeling of long waves in shallow water using IRBF-DQ and hybrid DQ/LRBF-DQ, Ocean Modelling 65 (2013) 1–10.

[69] J. Teng, N. Yasufuku, Q. Liu, K. Omine, A climate control apparatus and its application to evaluate evaporation process, in: Proceedings of 8th international symposium on lowland technology, Bali, Indonesia, 2012, pp. 161–166.

[70] J. Teng, N. Yasufuku, Q. Liu, S. Liu, Analytical solution for soil water redistribution during evaporation process, Water science and technology 68 (12) (2013) 2545–2551.

[71] R. Srivastava, T. J. Yeh, Analytical solutions for one-dimensional, transient infiltration toward the water table in homogeneous and layered soils, Water Resources Research 27 (5) (1991) 753–762.

[72] K. Perroux, D. Smiles, I. White, Water movement in uniform soils during constant-flux infiltration, Soil Science Society of America Journal 45 (2) (1981) 237–240.

[73] P. Broadbridge, I. White, Constant rate rainfall infiltration: A versatile nonlinear model: 1. analytic solution, Water Resources Research 24 (1) (1988) 145–154.