# A Numerical Method for a Nonlocal Form of Richards’ Equation Based on
Peridynamic Theory

Marco Berardi  Istituto di Ricerca sulle Acque, Consiglio Nazionale delle
Ricerche, Via F. de Blasio 5, 70132 Bari, Italy<EMAIL_ADDRESS>,  Fabio V. Difonzo  Dipartimento di Matematica, Università degli Studi di
Bari Aldo Moro, Via E. Orabona 4, 70125 Bari, Italy<EMAIL_ADDRESS>and  Sabrina F. Pellegrino  Dipartimento di Management, Finanza e Tecnologia,
Università LUM Giuseppe Degennaro, S.S. 100 Km 18, 70010 Casamassima,Italy
<EMAIL_ADDRESS>
###### Abstract.

Forecasting water content dynamics in heterogeneous porous media has
significant interest in hydrological applications; in particular, the
treatment of infiltration when in presence of cracks and fractures can be
accomplished resorting to peridynamic theory, which allows a proper modeling
of non localities in space. In this framework, we make use of Chebyshev
transform on the diffusive component of the equation and then we integrate
forward in time using an explicit method. We prove that the proposed spectral
numerical scheme provides a solution converging to the unique solution in some
appropriate Sobolev space. We finally exemplify on several different soils,
also considering a sink term representing the root water uptake.

###### Key words and phrases:

Richards’ equation, Peridynamic, Nonlocal Model, Spectral Numerical Method

###### 1991 Mathematics Subject Classification:

65M70, 42B30

##  1\. Introduction

Environmental protection and related sustainability management policies demand
a thorough understanding of complex coupling between hydrology, soil sciences,
ecology, agronomy, atmospheric sciences, calling for deeper mathematical
modeling and numerical methods able to deal with the multiphysics processes
involved in these environmental phenomena. In particular, flow processes in
unsaturated media have to be studied for a better understanding of the whole
water cycle; a correct managing of irrigation needs relies, for instance, on
robust numerical solvers for unsaturated flows with root water uptake (see for
instance, [21, 44]), or it is the basis for forecasting contaminant transport
in the vadose zone (see for instance [50]). Classical local advection-
diffusion equations in porous media often fail to describe accurately such
complex phenomena.
The idea of incorporating non-local behaviors in standard unsaturated flow
models is gaining interest in recent times. Besides non-localities in space,
which are the focus of this paper, also non-local effects in time can be
considered, that generally account for memory terms in the advection-diffusion
equations: in some pioneering works in the early ’60s [43]) it had been
already noticed that diffusivity depends not only on water content, but also
explicitly on time, and this argument has been then extended also to hydraulic
conductivity (see [23]; later on a model, in which derivative of water content
on time is fractional, has been first proposed in [41] and then generalized in
[28]. A memory component has been observed also when modeling water stress in
the root water uptake: the experimental evidence of such ”ecological memory”
of plant roots has been noticed, for instance, in [48, 13] and has been
recently formalized in [6].

When dealing with spatial discontinuities or significant heterogeneities,
classical local formulations of flow and transport phenomena present severe
limitations; for instance, in some cases, standard unsaturated flow models can
not forecast correctly water dynamics; as reported in [39], when modeling fast
infiltration processes (for instance infiltration after a heavy rainfall
event), ”first arrival time at the groundwater […] are often underpredicted”
because of preferential flow paths. These preferential flows can be ascribed
to non-equilibrium of water pressure at a local scale. As a matter of facts,
there is an experimental evidence that pore structures in natural soils
dynamically change due to alternating swelling and shrinkage processes (see
for instance [17]): this phenomenon can be described by a dual permeability
approach, by which the bulk porous medium consists of two dynamic interacting
pore domains: (i) the fracture (from shrinkage) pore domain and (ii) the
aggregate (interparticles plus structural pores), respectively (see [18]): in
practice, two different unsaturated flow equations are considered in each part
of this domain. Analogously, in the context of solute transport, the solute
exchange between mobile and immobile water has been modeled by a delay term in
[36], and, in a computational framework, this approach has been implemented in
[38].

More in detail, multirate mass transfer is modeled assuming advection-
diffusion on the fast mobile continuum and only diffusion in the slow immobile
continuum: after solving analytically the diffusion model, the consequent fast
domain model results non-local in time ([14]). In this dual-continuum
framework, the pioneering work [39] shows that the linearization of the
nonlinear diffusion equation, governing capillary flow in the slow continuum,
ensures a good description of the averaged water content dynamics in the slow
domain: therefore, they derive a non-local Richards’ equation in the mobile
domain, endowed with a memory kernel encoding mass transfer dynamics of the
slow domains.
From the viewpoint of applications, in this context, dessication cracks impact
the efficiency of irrigation and provokes a fast leakage of nutrients and
contaminants below the root zone into the groundwater. Even worst is the
effect of such cracks into an earth dam, since it may lead to the failure of
the dam itself.
On the other hand, several efforts have been accomplished towards a
comprehensive modeling and efficient numerical solvers for such nonlocal
problems. In [49] the coupling of peridynamic formulation of chemical
transport with water flow is proposed in the unsaturated context, and an
implicit numerical solver is implemented, and tested over different case
studies, in order to show the ability of the model to recognize
discontinuities and heterogeneities, including stationary cracks, propagating
cracks, and randomly distributed permeable and impermeable inclusions. In [19]
authors discuss how a single continuum model can properly catch the
contributions from all the flow paths only if the control volume (i.e. the
computational cell) is much larger than the longest connections between pores:
therefore a non-local model is proposed therein, showing that if the longest
connections are much smaller than the size of the control volume, these models
converge to Darcy’s law. A significant work has be presented in [40], in which
the peridynamic theory is employed for simulating hydraulic fracture
propagation in saturated porous media, and it is coupled with FEM for
forecasting fluid flow therein. In this paper we aim at introducing a tailored
numerical method for the corresponding peridynamic model of Richards’ equation
describing the unsaturated flow; for the sake of clarity we should say that
peridynamic theory was introduced by Silling in [46] as a nonlocal version of
elasticity theory, for modeling long-range interactions occurring in real
materials, ruling several phenomena like fractures, instabilities and cracks.
In general, peridynamic models consist of an integro-differential equation not
involving spatial derivatives and describe the motion of a material body
subjected to external loading conditions. The theory prescribes the existence
of a domain influence, called horizon, which represents a measure of the
nonlocality of the model and defines the range of interactions between
material particles.

In this framework, the remaining of the paper is structured as follows. In
Section 2 an introduction to nonlocal framework and a peridynamic formulation
of Richards’ equation is given, with all the necessary assumptions to justify
the current setting. Then, in Section 3, we propose a numerical method to
integrate forward in time a semi-discretized version of the equation,
leveraging spectral theory and Chebyshev transform properties to prove
convergence results of the discretized solution to the exact one. The
implementation of Chebyshev collocation method provides a good accuracy and
does not require to impose periodic boundary conditions. Finally, in Section 4
we exemplify on different soils with several type of Dirichlet boundary
conditions to support our findings.

###  1.1. A short overview on Richards’ equation

It is well known that Richards’ equation is a mass conservation law in terms
of the volumetric moisture content $\theta$ and of the soil matric head
$h_{m}$ defined on some compact domain $\Omega\subseteq\mathbb{R}^{3}$,
coupled with the Buckingham-Darcy’s law for the description on the flux:

$\displaystyle\frac{\partial\theta}{\partial t}(\mathbf{x},t)$
$\displaystyle=-\nabla
q(\mathbf{x},t)+S(\mathbf{x},\theta),\quad\mathbf{x}\in\Omega$ $\displaystyle
q(\mathbf{x},t)$ $\displaystyle=-K(h_{m})\nabla(h_{m}-z),$

where $z$ is the elevation component of the space variable $\mathbf{x}$,
$\theta$ represents the volumetric water content, $K$ is the so called
hydraulic conductivity and $S(\mathbf{x},\theta)$ is a source or sink term
describing, for instance, the root water uptake. Thus, Richards’ equation
reads as

(1.1) $\frac{\partial\theta}{\partial
t}(\mathbf{x},t)-\nabla\left(K(h_{m})\nabla(h_{m}-z)\right)=S(\mathbf{x},\theta),\quad\mathbf{x}\in\Omega,\,t\in[0,T],$

endowed with suitable initial and boundary conditions.

With the hypothesis that air pressure in the pores is constant, Richards’
equation assumes that matric head at a given location is in equilibrium and
that there exists a bijective function relating $\theta$ with $h_{m}$, called
water retention curve (see [39]), which is generally defined according to
empirical functions. Moreover, for Richards’ equation to be well posed, $K$
must be smooth enough to guarantee existence and uniqueness of solutions, also
in case of heterogeneous soils with smooth boundary (see [8] and references
therein). In particular, hereafter and through the whole paper, $K$ and
$h_{m}$ will be assumed to be locally Lipschitz on their respective domains.
However, in case of desiccation cracks or anisotropic soils could affect well-
posedness of Richards’ equation (1.1) and prevent existence of any solution.
An alternative approach has been proposed in [26], where theory of elasticity
for solid mechanics has been applied to unsaturated, heterogeneous,
anisotropic soils. In this case, though, the flow density function depends on
the position, matric head or moisture content, instead of the relative
distance and relative displacement [46, 47].
The numerous numerical issues arising when solving Richards’ equation in a
computational framework rely mainly in its nature of highly nonlinear
degenerate elliptic parabolic PDE. Here we just mention some significant
references for the main numerical problems arising in Richards’ equation. For
instance, since implicit methods are generally used for time integration, the
arising nonlinear problems have been studied with different methods, such as
Newton’s (e.g. [15, 10], Picard ([16]), L-Scheme or its variants ([42, 37]).
Even richer is the literature on spatial discretization techniques, for which
we refer to [2, 30, 35, 29] and references therein. As regards numerical
integration over layered discontinuous geological formations, a domain
decomposition approach is followed in [45], while a transversal method of
lines is adopted in [5, 9].
In this paper, we are looking at the 1D version of (1.1) equipped with initial
and Dirichlet boundary conditions, in which diffusion evolves exclusively
along the depth, so that $\Omega=[0,Z]$ for some $Z>0$, and the forcing term
$S$ only depends on $z\in[0,Z]$. Thus, one considers

(1.2) $\frac{\partial\theta}{\partial t}(z,t)-\frac{\partial}{\partial
z}\left(K(h_{m})\frac{\partial}{\partial z}(h_{m}-z)\right)=S(z),\quad
z\in[0,Z],\,t\in[0,T].$

##  2\. Peridynamic Model: assumptions and derivation

Let us consider a compact domain $\Omega$ with smooth boundary and let us
define

(2.1)
$B_{\delta}(z)\mathrel{\mathop{:}}=\\{z^{\prime}\in\Omega\,:\,\|z^{\prime}-z\|\leq\delta\\},$

the _horizon_ of $z$ of radius $\delta>0$. We assume that moisture dynamics at
$z$ is only affected by pairwise interaction with $z^{\prime}\in
B_{\delta}(z)$; points outside the horizon of $z$ do not contribute to any
dynamics therein.

The model is built on the concept of peripipes. Given any $z\in\Omega$, we
assume that for each $z^{\prime}\in B_{\delta}(z)$ there exists a fictitious
pipe, called _peripipe_ , connecting every $z$ to $z^{\prime}$. We assume that
the following requirements hold for any peripipe (see [25]):

  1. (1)

Moisture is stored at the endpoints $z,x^{\prime}$ of a peripipe, and zero
moisture content is located along a peripipe;

  2. (2)

moisture flows in the direction of the peripipe and no transversal flux
crosses its boundaries

  3. (3)

a peripipe is purely resistive, it has zero reactance and its response is
proportional to $H(z)-H(z^{\prime})$;

  4. (4)

a peripipe has uniform conductivity;

  5. (5)

peripipe conductivity is function of medium conductivity at its endpoints;

  6. (6)

the length of a peripipe is $\|z-z^{\prime}\|$;

  7. (7)

peripipe response may also depend on the its length.

Following [26] and requirements above, we assume that the rate of volumetric
moisture flow from a point $z^{\prime}$ to a point $z$ per unit volume of $z$
and per unit volume of $z^{\prime}$ is given by

(2.2) $J(z,z^{\prime},t)=C(z,z^{\prime})(H(z^{\prime},t)-H(z,t)),$

where $C(z,z^{\prime})$ is the peridynamic hydraulic conductance density and
$H(z,t)$ is the total hydraulic potential, defined as

$H(z,t)=h_{m}(z,t)+z.$

Hereafter, for the sake of readability, we omit time dependence, unless
required by the context.

The peripipe conductance depends on the peridynamic hydraulic conductivity
$\kappa(z^{\prime},z)$, which is an intrinsic material property (related to
the classical hydraulic conductivity $K$), in the following way:

(2.3) $C(z,z^{\prime})=\frac{\kappa(z^{\prime},z)}{\|z-z^{\prime}\|},$

where

(2.4) $\kappa(z,z^{\prime})\mathrel{\mathop{:}}=K\varphi(z-z^{\prime}).$

The function $\varphi(z-z^{\prime})$ is the so-called _influence function_ ,
representing a convolution kernel relating the horizon (2.1) with the nature
of boundary conditions assigned to (1.2). The shape of such an even function
and the way to select it turns out to be crucial, as we will see in Section
2.1.

Therefore, the changes of moisture stored at $z$ and at $z^{\prime}$, mediated
by the peripipe $zz^{\prime}$, are given by

$\displaystyle\Delta V_{m}(z,z^{\prime})$
$\displaystyle=\kappa(z^{\prime},z)\frac{H(z^{\prime})-H(z)}{\|z-z^{\prime}\|}\mathrm{d}V_{z^{\prime}}\mathrm{d}V_{z},$
$\displaystyle\Delta V_{m}(z^{\prime},z)$
$\displaystyle=\kappa(z,z^{\prime})\frac{H(z)-H(z^{\prime})}{\|z^{\prime}-z\|}\mathrm{d}V_{z}\mathrm{d}V_{z}^{\prime}.$

As an immediate consequence it must hold
$\kappa(z,z^{\prime})=\kappa(z^{\prime},z)$.
In case of inhomogeneous soils in unsaturated regime, above relations could be
leveraged to define a peridynamic conductivity density by setting

(2.5)
$\kappa(z,z^{\prime})\mathrel{\mathop{:}}=\frac{\kappa(z)+\kappa(z^{\prime})}{2},$

where $\kappa(z)\equiv\kappa(z,0)$, as proposed in [49, 26].
Now, since the change over time of volumetric moisture content due to
$z^{\prime}$ at time $t$, on the account of (2.2), is given by

$\frac{\partial\theta}{\partial t}(z|z^{\prime},t)=J(z,z^{\prime}),$

from which

$\frac{\partial\theta}{\partial
t}(z,t)=\int_{B_{\delta}(z)}J(z,z^{\prime})\,\mathrm{d}V_{z^{\prime}}+S(z).$

Thus, using (2.3) and with peridynamic conductivity given by (2.5), our model
(1.2), endowed with Dirichlet boundary conditions, reads

(2.6a) $\displaystyle\frac{\partial\theta}{\partial t}$
$\displaystyle=\int_{B_{\delta}(z)}\frac{\varphi(z^{\prime}-z)}{\|z^{\prime}-z\|}\frac{K(z)+K(z^{\prime})}{2}[H(z^{\prime})-H(z)]\,\mathrm{d}V_{z^{\prime}}+S(z),$
(2.6b) $\displaystyle\theta(z,0)$ $\displaystyle=\theta^{0}(z),\quad
z\in[0,Z],$ (2.6c) $\displaystyle\theta(0,t)$
$\displaystyle=\theta_{0}(t),\quad t\in[0,T],$ (2.6d)
$\displaystyle\theta(Z,t)$ $\displaystyle=\theta_{Z}(t),\quad t\in[0,T].$

###  2.1. Selection of the influence function

Usually (see, e.g. [49, 11, 26]) $\varphi(z)$ in (2.4) represents a
convolution kernel, which can be chosen as a uniform influence function

$\varphi(z)\mathrel{\mathop{:}}=\begin{cases}\frac{2}{\delta},&\|z\|\leq\delta,\\\
0,&\|z\|>\delta,\end{cases}$

or as a linear influence function

$\varphi(z)\mathrel{\mathop{:}}=\begin{cases}1-\frac{\|z\|}{\delta},&\|z\|\leq\delta,\\\
0,&\|z\|>\delta.\end{cases}$

However, since such kernels would suggest the model to weigh more those cells
where they are nonzero, and since our boundary conditions would typically be
of Dirichlet type, we propose to consider a _distributed influence function_
(see Figure 1), concentrated on the domain boundary, of the form

(2.7)
$\varphi(z)\mathrel{\mathop{:}}=\begin{cases}\frac{\|z\|-1+\delta}{\delta},&\|z\|\geq
1-\delta,\\\ 0,&\|z\|<1-\delta.\end{cases}$

$-1$$-0.8$$-0.6$$-0.4$$-0.2$$0$$0.2$$0.4$$0.6$$0.8$$1$$0$$0.2$$0.4$$0.6$$0.8$$1$$z$$\varphi(z)$
Figure 1. Distributed influence function defined in (1) with $\delta=0.15$.

In so doing, we are suggesting the model to averaging out not just what
happens in the middle of the dynamics, but rather the behavior around each
point of the spatial domain. In all our experiments, presented in Section 4,
uniform and linear influence functions do not make our proposed numerical
method converge, resulting in instabilities and blow-ups after a relatively
small amount of time integration; on the other hand, and as presented below,
using (2.7) guarantees stability and convergence, plus a reasonable shape of
the numerical solutions.

##  3\. Numerical Method

The nonlocal Richards’ equation (2.6) can be discretized in space by using
Chebyshev polynomials. This approach is typically used when the integral
operator can be expressed in terms of convolution products [34, 33, 32].
Moreover, the choice of such kind of polynomials allows us to overcome the
limitation of imposing periodic boundary condition, which is necessary when
dealing with Fourier trigonometric polynomials.

The proposed technique consists in looking for an approximation of
$\theta(x,t)$ in the form of a finite linear combination of Chebyshev
polynomials of the first kind. To do so, we can assume the spatial domain to
be $[-1,1]$, as we can benefit of the orthogonality properties of the
polynomials. However, a more general interval can be used as spatial domain by
linearity. Moreover, for time integration we use the explicit Euler method, as
in [49].

In this section we briefly make a review on Chebyshev polynomials, then we
derive the semidiscrete model of (2.6) and finally prove the convergence of
the proposed method.

###  3.1. Basic overview on Chebyshev polynomials

Chebyshev polynomials of the first kind, $T_{k}(z)$ are defined by

$T_{k}(z)=\cos\left(k\arccos(z)\right),\qquad z\in[-1,1],\,\,k\in\mathbb{N},$

and are orthogonal with respect to the weight function
$w(z)\mathrel{\mathop{:}}=\left(\sqrt{1-z^{2}}\right)^{-1}$.

These polynomials are commonly used in the context of spectral approximation
because they satisfy an interpolation property: given an integer $N$, any
sufficiently smooth function $u$ defined on $[-1,1]$ can be expanded as an
$(N+1)$-term linear combination of polynomials $u^{N}$ given by

(3.1) $u^{N}(z)\mathrel{\mathop{:}}=\sum_{k=0}^{N}\widetilde{u}_{k}\
T_{k}(z),$

where $\widetilde{u}_{k}$ are the coefficients of the expansion and
approximate the Chebyshev coefficients

$\widehat{u}_{k}=\frac{2}{\pi
c_{k}}\int_{-1}^{1}u(x)T_{k}(z)w(x)\,\mathrm{d}z,$

with

$c_{k}=\begin{cases}2&\quad k=0\\\ 1&\quad k\neq 0.\end{cases}$

The explicit expression of $\widetilde{u}_{k}$ depends on the choice of the
grid points used to discretize $[-1,1]$. In particular, if we choose the
Gauss-Lobatto collocation points

(3.2) $z_{h}\mathrel{\mathop{:}}=\cos\left(\frac{h\pi}{N}\right),\qquad
h=0,\dots,N,$

the expression of $\widetilde{u}_{k}$ is

(3.3) $\widetilde{u}_{k}=\frac{1}{\gamma_{k}}\sum_{h=0}^{N}u(x_{h})\
T_{k}(z_{h})w_{h},$

where $\gamma_{k}$ is a normalization constant defined by

(3.4) $\gamma_{k}\mathrel{\mathop{:}}=\begin{cases}\pi\quad&k=0,N\\\
\frac{\pi}{2}\quad&k=1,\dots,N-1\end{cases}$

and

(3.5) $w_{h}\mathrel{\mathop{:}}=\begin{cases}\frac{\pi}{2N}\quad&h=0,N\\\
\frac{\pi}{N}\quad&h=1,\dots,N-1.\end{cases}$

Equation (3.1) represents the inverse discrete Chebyshev transform, while the
coefficients $\widetilde{u}_{k}$ in (3.3) correspond to the discrete Chebyshev
transform. They can be efficiently computed using the Fast Fourier Transform.
Additionally, they fulfill the same properties of the Fourier transform. In
particular, we can rewrite a convolution product in the physic space as a
multiplication of the Chebyshev transform of each factor in the frequency
space.

The following result shows the rate of convergence of the Chebyshev
approximation.

######  Theorem 3.1 (see [12]).

For any $0\leq\mu\leq 2$ and $u\in L^{2}([-1,1])$, there exists a positive
constant $C$ independent on $N$, such that

(3.6)
$\left\|u-u^{N}\right\|_{L^{2}([-1,1])}\leq\frac{C}{N^{2-\mu}}\left\|u\right\|_{L^{2}([-1,1])}.$

In the next section, to lighten the notation, we denote the Chebyshev
transform by $\mathcal{F}$ and the inverse Chebyshev transform by
$\mathcal{F}^{-1}$.

###  3.2. Chebyshev semi-discrete collocation method for the nonlocal
Richards’ equation

In what follows, we develop a spectral approximation of (2.6) by using the
Chebyshev transform. We fix $N>0$ and assume $\Omega=[-1,1]$. We can
discretize the spatial domain by the Guass-Lobatto points $z_{h}$,
$h=0,\dots,N$ defined in (3.2).

If we set

$\displaystyle\Lambda(z)$ $\displaystyle\mathrel{\mathop{:}}=K(z)H(z),$
$\displaystyle\overline{\varphi}(z)$
$\displaystyle\mathrel{\mathop{:}}=\frac{\varphi(z)}{\left\|z\right\|},$
$\displaystyle\beta$
$\displaystyle\mathrel{\mathop{:}}=\int_{-1}^{1}\overline{\varphi}(z)\,\mathrm{d}z,$

then, since from distributed influence function definition (2.7) it follows
that $\beta=2\left(1+\frac{1-\delta}{\delta}\ln(1-\delta)\right)$, we can
rewrite model (2.6a) as

(3.7) $\begin{split}\frac{\partial\theta}{\partial
t}&=\int_{B_{\delta}(z)}\frac{\varphi(z^{\prime}-z)}{\|z^{\prime}-z\|}\
\frac{K(z^{\prime})+K(z)}{2}\
[H(z^{\prime})-H(z)]\,\mathrm{d}V_{z^{\prime}}+S(z),\\\
&=\frac{1}{2}\left[\left(\overline{\varphi}\ast\Lambda\right)(z)+K(z)\left(\overline{\varphi}\ast
H\right)(z)-H(z)\left(\overline{\varphi}\ast
K\right)(z)-\beta\Lambda(z)\right]+S(z).\end{split}$

Thus, the right hand side of (3.7) can be computed by means of the finite
discrete Chebyshev transform. Indeed, we have

(3.8) $\displaystyle\left(\overline{\varphi}\ast\Lambda\right)(z)$
$\displaystyle=\mathcal{F}^{-1}\left(\mathcal{F}\left(\overline{\varphi}\right)\mathcal{F}\left(\Lambda\right)\right)(z),$
(3.9) $\displaystyle\left(\overline{\varphi}\ast H\right)(z)$
$\displaystyle=\mathcal{F}^{-1}\left(\mathcal{F}\left(\overline{\varphi}\right)\mathcal{F}\left(H\right)\right)(z),$
(3.10) $\displaystyle\left(\overline{\varphi}\ast K\right)(z)$
$\displaystyle=\mathcal{F}^{-1}\left(\mathcal{F}\left(\overline{\varphi}\right)\mathcal{F}\left(K\right)\right)(z).$

So, at each collocation point $z_{h}$, the semi-discretization of the model
reads

(3.11) $\begin{split}\frac{\partial\theta}{\partial
t}(z_{h},t)&=\frac{1}{2}\left(\mathcal{F}^{-1}\left(\mathcal{F}\left(\overline{\varphi}\right)\mathcal{F}\left(\Lambda\right)\right)(z_{h})+K(z_{h})\
\mathcal{F}^{-1}\left(\mathcal{F}\left(\overline{\varphi}\right)\mathcal{F}\left(H\right)\right)(z_{h})\right)\\\
&\quad-\frac{1}{2}\left(H(z_{h})\
\mathcal{F}^{-1}\left(\mathcal{F}\left(\overline{\varphi}\right)\mathcal{F}\left(K\right)\right)(z_{h})+\beta\Lambda(z_{h})\right)+S(z_{h})\end{split}$

The function $\Lambda$ is defined as the product between the conductivity $K$
and the hydraulic potential $H$: therefore, to compute its Chebyshev
transform, we first need to compute a product. The computational cost to
obtain this term could be efficiently reduced by observing that the Chebyshev
coefficients of $\Lambda$ can be obtained from the Chebyshev coefficients of
$H$ and $K$.

Indeed, the following result holds (see [3]).

######  Theorem 3.2.

Let $N\in\mathbb{N}$. If $H$ and $K$ are approximated by a finite series of
Chebyshev polynomials $H^{N}$ and $K^{N}$, respectively, given by

$H^{N}(z)=\sum_{j=0}^{N}\widetilde{H}_{k}T_{k}(z),\quad
K^{N}(z)=\sum_{j=0}^{N}\widetilde{K}_{j}T_{j}(z),$

then the product $\Lambda(z)=H(z)K(z)$ can be approximated by the following
$2N+1$ combination of Chebyshev polynomials

$\Lambda^{N}(z)=\sum_{j=0}^{2N}\widetilde{\Lambda}_{j}T_{j}(z),$

where the coefficients $\widetilde{\Lambda}_{j}$ are given by

$2\widetilde{\Lambda}_{j}=\begin{cases}2\widetilde{H}_{0}\widetilde{K}_{0}+\sum_{\ell=1}^{N}\widetilde{H}_{\ell}\widetilde{K}_{\ell},\quad&j=0\\\
\sum_{\ell=0}^{j}\widetilde{H}_{j-\ell}\widetilde{K}_{\ell}+\sum_{\ell=0}^{N-j}\widetilde{H}_{j+\ell}\widetilde{K}_{\ell}+\sum_{\ell=j}^{N}\widetilde{H}_{\ell-j}\widetilde{K}_{\ell},\quad&j=1,\dots,N\\\
\sum_{\ell=j-N}^{N}\widetilde{H}_{j-\ell}\widetilde{K}_{\ell},\quad&j=N+1,\dots,2N.\end{cases}$

The application of Theorem 3.2 implies that the first term in the right hand
side of (3.11) is discretized by $2N+1$ mesh points. Therefore, to maintain
the consistency of the scheme, the discretization of the remaining terms on
the right hand side of (3.11) is accomplished by considering $2N+1$ Gauss-
Lobatto collocation points.

###  3.3. Convergence of the semi-discrete scheme

We prove the convergence of the spectral semi-discrete method in a suitable
weighted Hilbert space. Throughout this section, $C$ denotes a generic
constant.

We consider the weighted Lebesgue space

$L^{2}_{w}([-1,1])=\left\\{u\in
L^{2}\,:\,\int_{-1}^{1}u^{2}(z)w(z)\mathrm{d}z<+\infty\right\\}$

equipped with the inner product and the norm respectively

$(u,v)_{w}=\int_{-1}^{1}u(z)v(z)w(z)\,\mathrm{d}z,\qquad\left\|u\right\|^{2}_{w}=(u,u)_{w},$

where $w(z)=\left(\sqrt{1-z^{2}}\right)^{-1}$.

For any $s\geq 0$, we set

$H^{s}_{w}([-1,1])=\left\\{u\in L^{2}_{w}([-1,1])\ |\
\left\|u\right\|_{s,w}<+\infty\right\\},$

where

$\left\|u\right\|_{s,w}^{2}=\sum_{|\alpha\leq
s|}\left\|D^{\alpha}u\right\|_{w}^{2}.$

Let $S_{N}$ be the space of Chebyshev polynomials of degree $N$,

$S_{N}\mathrel{\mathop{:}}=\text{span}\left\\{T_{h}(x)\ |\ 0\leq h\leq
N\right\\}\subset L^{2}_{w}([-1,1]),$

and $P_{N}:L^{2}_{w}([-1,1])\to S_{N}$ be an orthogonal projection operator

$P_{N}u(x)\mathrel{\mathop{:}}=\sum_{h=0}^{N}\widehat{u}_{h}T_{h}(x)w_{h},$

for $w_{h}$ defined in (3.5), such that for any $u\in L^{2}_{w}([-1,1])$, the
following equality holds

(3.12) $(u-P_{N}u,\varphi)_{w}=0,\quad\text{for every $\varphi\in S_{N}$}.$

The projection operator $P_{N}$ commutes with derivatives in the
distributional sense:

$\partial_{t}^{q}P_{N}u=P_{N}\partial_{t}^{q}u,\quad q\in\mathbb{N},q\geq 1,$

where, as usual, $\partial_{t}\mathrel{\mathop{:}}=\frac{\partial}{\partial
t}$.
Letting $s\geq 1$, we denote by
$X_{s}\mathrel{\mathop{:}}=\mathcal{C}^{0}\left(0,T;H^{s}_{w}([-1,1])\right)$
the space of all continuous functions in the weighted Sobolev space
$H^{s}_{w}([-1,1])$, with norm

$\left\|u\right\|_{X_{s}}^{2}\mathrel{\mathop{:}}=\max_{t\in[0,T]}\left\|u(\cdot,t)\right\|_{s,w}^{2},$

for any $T>0$. We denote by $\mathcal{L}$ the nonlocal integral operator of
(2.6), namely

(3.13)
$\mathcal{L}\left(\theta\right)\mathrel{\mathop{:}}=\int_{B_{\delta}(z)}\frac{\varphi(z^{\prime}-z)}{\|z^{\prime}-z\|}\frac{K(z)+K(z^{\prime})}{2}[H(z^{\prime})-H(z)]\,\mathrm{d}V_{z^{\prime}}.$

Then, the semi-discrete spectral scheme for (2.6) can be rewritten as

(3.14) $\displaystyle\frac{\partial\theta^{N}}{\partial
t}=P_{N}\mathcal{L}(\theta^{N})+P_{N}S(z),$ (3.15)
$\displaystyle\theta^{N}(z,0)=P_{N}\theta^{0}(z),$

where $\theta^{N}(z,t)\in S_{N}$ for every $0\leq t\leq T$.

To obtain the convergence of the semi-discrete scheme, we need of the
following lemma.

######  Lemma 3.3 ([12, Theorem 3.1]).

For any real $0\leq\mu\leq s$, there exists a positive constant $C$ such that

(3.16)
$\left\|u-P_{N}u\right\|_{H^{\mu}_{s,w}{\mu}([-1,1])}\leq\frac{C}{N^{s-\mu}}\left\|\theta\right\|_{H^{s}_{w}([-1,1])},\quad\text{for
every $\theta\in H^{s}_{w}([-1,1])$}.$

Recalling that $K$ and $h_{m}$ are locally Lipschitz in their respective
domains, we can prove the following theorem.

######  Theorem 3.4.

Let $s\geq 1$ and $\theta(z,t)\in X_{s}$ be the solution to the initial-
boundary-valued problem (2.6) and $\theta^{N}(z,t)$ be the solution to the
semi-discrete scheme (3.14)-(3.15). Then, there exists a positive constant
$C$, independent on $N$, such that

(3.17) $\left\|\theta-\theta^{N}\right\|_{X_{1}}\leq
C(T)\left(\frac{1}{N}\right)^{s-1}\left\|u\right\|_{X_{s}},$

for any initial data $\theta^{0}\in H^{s}_{w}([-1,1])$ and for any $T>0$.

###### Proof.

Let $s\geq 1$. Using the triangular inequality, we have

(3.18) $\left\|\theta-\theta^{N}\right\|_{X_{1}}\leq\left\|\theta-
P_{N}\theta\right\|_{X_{1}}+\left\|P_{N}\theta-\theta^{N}\right\|_{X_{1}}.$

Lemma 3.3 implies

$\left\|(\theta-
P_{N}\theta)(\cdot,t)\right\|_{H^{1}_{w}([-1,1])}\leq\frac{C}{N^{s-1}}\left\|\theta(\cdot,t)\right\|_{H^{s}_{w}([-1,1])}.$

Therefore,

(3.19) $\left\|\theta-
P_{N}\theta\right\|_{X_{1}}\leq\frac{C}{N^{s-1}}\left\|\theta\right\|_{X_{s}}.$

Subtracting (3.14) from (2.6) and taking the weighted inner product with
$P_{N}\theta-\theta^{N}\in S_{N}$, we have

(3.20)
$\begin{split}0=&\underbrace{\int_{-1}^{1}\left(\partial_{t}\theta(z,t)-\partial_{t}\theta^{N}(z,t)\right)\left(P_{N}\theta(z,t)-\theta^{N}(z,t)\right)w(z)\,\mathrm{d}z}_{{}=:I_{1}}\\\
&-\underbrace{\int_{-1}^{1}\left(\mathcal{L}(\theta(z,t))-P_{N}\mathcal{L}(\theta^{N}(z,t))\right)\left(P_{N}\theta(z,t)-\theta^{N}(z,t)\right)w(z)\,\mathrm{d}z}_{{}=:I_{2}}\\\
&-\underbrace{\int_{-1}^{1}\left(S(z,t)-P_{N}S(z,t)\right)\left(P_{N}\theta(z,t)-\theta^{N}(z,t)\right)w(z)\,\mathrm{d}z}_{{}=:I_{3}}.\end{split}$

The orthogonal condition (3.12) implies that

$\int_{-1}^{1}\left(\partial_{t}\theta(z,t)-P_{N}\partial_{t}\theta(z,t)\right)\left(P_{N}\theta(z,t)-\theta^{N}(z,t)\right)w(z)\,\mathrm{d}z=0,$

and

$\int_{-1}^{1}\left(S(z)-P_{N}S(z)\right)\left(P_{N}\theta(z,t)-\theta^{N}(z,t)\right)w(z)\,\mathrm{d}z=0.$

Thus,

(3.21)
$\begin{split}I_{1}&=\int_{-1}^{1}\left(\partial_{t}\theta(z,t)-P_{N}\partial_{t}\theta(z,t)\right)\left(P_{N}\theta(z,t)-\theta^{N}(z,t)\right)w(z)\,\mathrm{d}z\\\
&+\int_{-1}^{1}\left(P_{N}\partial_{t}\theta(z,t)-\theta^{N}_{t}(z,t)\right)\left(P_{N}\theta(z,t)-\theta^{N}(z,t)\right)w(z)\,\mathrm{d}z\\\
&=\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\left\|(P_{N}\theta-\theta^{N})(\cdot,t)\right\|^{2}_{H^{1}_{w}([-1,1])}.\end{split}$

Since it is straightforward to see that $I_{3}=0$, we can focus on $I_{2}$.
From (3.12) it follows that

$\int_{-1}^{1}\left(\mathcal{L}(\theta^{N}(z,t))-P_{N}\mathcal{L}(\theta^{N}(z,t))\right)\left(P_{N}\theta(z,t)-\theta^{N}(z,t)\right)w(z)\,\mathrm{d}z=0,$

and so, due to the locally Lipschitzianity of $K$ and $h_{m}$ and from Cauchy
and triangular inequalities, we obtain

(3.22)
$\begin{split}I_{2}&=\int_{-1}^{1}\left(\mathcal{L}(\theta(z,t))-\mathcal{L}(\theta^{N}(z,t))\right)\left(P_{N}\theta(z,t)-\theta^{N}(z,t)\right)w(z)\,\mathrm{d}z\\\
&\leq
C\left\|\left(\theta-\theta^{N}\right)(\cdot,t)\right\|^{2}_{H^{1}_{w}([-1,1])}+C\left\|(P_{N}\theta-\theta^{N})(\cdot,t)\right\|^{2}_{H^{1}_{w}([-1,1])}\\\
&\leq 2C\left\|\left(\theta-
P_{N}\theta\right)(\cdot,t)\right\|^{2}_{H^{1}_{w}([-1,1])}+3C\left\|(P_{N}\theta-\theta^{N})(\cdot,t)\right\|^{2}_{H^{1}_{w}([-1,1])}.\end{split}$

Plugging (3.21) and (3.22) in (3.20), we have

(3.23)
$\frac{\mathrm{d}}{\mathrm{d}t}\left\|(P_{N}\theta-\theta^{N})(\cdot,t)\right\|^{2}_{H^{1}_{w}([-1,1])}\leq
4C\left\|(\theta-
P_{N}\theta)(\cdot,t)\right\|^{2}_{H^{1}_{w}([-1,1])}+6C\left\|(P_{N}\theta-\theta^{N})(\cdot,t)\right\|^{2}_{H^{1}_{w}([-1,1])}.$

Since
$\left\|(P_{N}\theta-\theta^{N})(\cdot,0)\right\|_{H^{1}_{w}([-1,1])}=0$,
Lemma 3.3 and Gronwall’s inequality imply that

$\displaystyle\left\|(P_{N}\theta-\theta^{N})(\cdot,t)\right\|^{2}_{H^{1}([-1,1])}$
$\displaystyle\leq 4C\int_{0}^{t}e^{6C(t-\tau)}\left\|(\theta-
P_{N}\theta)(\cdot,\tau)\right\|^{2}_{H^{1}_{w}([-1,1])}\,\mathrm{d}\tau$
$\displaystyle\leq\frac{C(T)}{N^{2s-2}}\int_{0}^{t}\left\|\theta(\cdot,\tau)\right\|^{2}_{H^{1}_{w}([-1,1])}\,\mathrm{d}\tau.$

Hence,

(3.24)
$\left\|P_{N}\theta-\theta^{N}\right\|_{X_{1}}\leq\frac{C(T)}{N^{s-1}}\left\|\theta\right\|_{X_{s}}.$

Finally, using (3.19) and (3.24) in (3.18), we complete the proof. ∎

##  4\. Numerical Simulations

In this section we test our proposed method on different soils, possibly with
a sink forcing term, representing the water uptake due to root systems. Van
Genuchten - Mualem constitutive relations are considered in the following
numerical simulations.

######  Example 4.1.

Drawing from [4], we consider a soil with the following parameters:

$\theta_{r}=0.075,\,\theta_{S}=0.287,\,\alpha=0.036,\,n=1.56,\,K_{S}=0.94e-3\,\textrm{cm/s}.$

Figure 2. Numerical solution relative to Example 4.1.

We added a sink term $S=-700\,\textrm{s}^{-1}$ and parameter $\delta=0.15$ in
(2.1). We set our initial and boundary conditions as follows

$\displaystyle\theta(0,t)$
$\displaystyle=0.2234\left(1-\frac{t}{T}\right)+0.181\frac{t}{T},\,\,t\in[0,T],$
$\displaystyle\theta(Z,t)$
$\displaystyle=0.1368\left(1-\frac{t}{T}\right)+0.1174\frac{t}{T},\,\,t\in[0,T],$
$\displaystyle\theta(z,0)$
$\displaystyle=0.2234-\left(1-\frac{z}{Z}\right)\frac{0.0848}{2},\,\,z\in[0,Z].$

Here $Z=30$ cm, $T=60$ s; moreover, we used $\Delta t=0.06$ s and $\Delta
x=0.3$ cm. Results are shown in Figure 2.

######  Example 4.2.

As already considered by [24], we select a Glendale clay loam, characterized
by the following parameters

$\theta_{r}=0.1060,\,\theta_{S}=0.4686,\,\alpha=0.0104,\,n=1.3954,\,K_{S}=1.5162e-4\,\textrm{cm/s}.$

Figure 3. Numerical solution relative to Example 4.2.

We put neither sink nor source on this simulation. Peridynamic parameter is
$\delta=0.15$ in (2.1). Our boundary conditions are constant with values

$\displaystyle\theta(0,t)$ $\displaystyle=0.2,\,\,t\in[0,T],$
$\displaystyle\theta(Z,t)$ $\displaystyle=0.3,\,\,t\in[0,T],$

while initial condition follows a nonlinear profile of the form

$\theta(z,0)=-0.05z^{3}+0.25,\,\,z\in[0,Z].$

We select $Z=70$ cm, $T=2400$ s; moreover, we used $\Delta t=2.4$ s and
$\Delta x=0.3$ cm. The resulting water content profiles are shown in Figure 3.

######  Example 4.3.

As in [24], we consider a Berino loamy fine sand, with parameters

$\theta_{r}=0.0286,\,\theta_{S}=0.3658,\,\alpha=0.0280,\,n=2.2390,\,K_{S}=0.0063\,\textrm{cm/s}.$

Figure 4. Numerical solution relative to Example 4.3.

We added a sink term $S=-100\,\textrm{s}^{-1}$ and parameter $\delta=0.15$ in
(2.1). We set our initial and boundary conditions as follows

$\displaystyle\theta(0,t)$
$\displaystyle=0.3\left(1-\frac{t}{T}\right)+0.29\frac{t}{T},\,\,t\in[0,T],$
$\displaystyle\theta(Z,t)$ $\displaystyle=0.2,\,\,t\in[0,T],$

while initial condition has a nonlinear profile

$\theta(z,0)=0.05z^{3}+0.25,\,\,z\in[0,Z].$

We select $Z=70$ cm, $T=600$ s; moreover, we used $\Delta t=0.06$ s and
$\Delta x=0.3$ cm. Results are shown in Figure 4.

##  5\. Conclusions and future works

Starting from an appropriate and physically based rewriting of Richards’
equation using non-locality theory, we propose to compute its numerical
solution using a semi-discretized time forward scheme based on Chebyshev
transform. We prove that such approach converges in suitable Sobolev spaces,
providing a theoretical background for further extensions of the present work
to higher dimensional domains. We also propose a new kind of convolutional
kernel, or influence function, in order to manage the peridynamic behavior of
the proposed model. Such an influence function distributes its effect on the
domain so to correctly catch boundary conditions. In fact, we experienced
major benefits from this choice, as numerical convergence turns out to be
robust with respect to time, and compared to results coming from classical
choices of influence functions. To testify our theoretical analysis, we have
performed several experiments, on a wide range on soils, in MATLAB. We have
considered different Dirichlet boundary conditions and linear and non-linear
initial conditions and show that, with suitable discretization step-sizes, our
method is reliable and accurate.

The present work suggests several possible directions for future and already
ongoing research studies. For instance, it would be of interest applying
Eulerian-Lagrangian methods (e.g. [1]) in the proposed peridynamic framework
for Richards’ equation, while the idea of introducing a basic control approach
on the boundary conditions, as in [7], could be adapted as well. Another
development would consider non-local terms in time, so to resort to numerical
solvers coming from specific tools in fractional differential calculus (see
[22, 20]), or, even more promisingly, by spectral methods in 2D (see [33, 27,
31]).

## Acknowledgments

All authors are member of the INdAM Research group GNCS. MB and FVD also
acknowledge the partial support of the 2022 project “Modelli di evoluzione non
locali: analisi, trattamento numerico e algoritmi” funded by GNCS-INdAM. SFP
acknowledges the partial support of “Finanziamento giovani ricercatori 2022”
funded by GNCS-INdAM. MB acknowledges the partial support of the CNR project
“MENTOR”. FVD has been supported by REFIN Project, grant number 812E4967
funded by Regione Puglia. SFP has been supported by REFIN Project, grant
number D1AB726C funded by Regione Puglia.

## References

  * [1] E. Abreu, R. De la cruz, J.C. Juajibioy, and W. Lambert.  Lagrangian-Eulerian Approach for Nonlocal Conservation Laws.  Journal of Dynamics and Differential Equations, 2022.
  * [2] T. Arbogast, M.F. Wheeler, and N.Y. Zhang.  A Nonlinear Mixed Finite Element Method for a Degenerate Parabolic Equation Arising in Flow in Porous Media.  SIAM Journal on Numerical Analysis, 33(4):1669––1687, 1996.
  * [3] Günter Baszenski and Manfred Tasche.  Fast polynomial multiplication and convolutions related to the discrete cosine transform.  Linear Algebra and Its Applications, 252(1-3):1 – 25, 1997.
  * [4] M. Berardi, F. Difonzo, F. Notarnicola, and M. Vurro.  A transversal method of lines for the numerical modeling of vertical infiltration into the vadose zone.  Applied Numerical Mathematics, 135:264 – 275, 2019.
  * [5] M. Berardi, F. Difonzo, M. Vurro, and L. Lopez.  The 1D Richards’ equation in two layered soils: a Filippov approach to treat discontinuities.  Advances in Water Resources, 115:264–272, may 2018.
  * [6] M. Berardi and G. Girardi.  Richards’ equation with non-local root water uptake modeling plant water deficit.  submitted, December 2022.
  * [7] Marco Berardi, Marcello D’Abbicco, Giovanni Girardi, and Michele Vurro.  Optimizing water consumption in Richards’ equation framework with step-wise root water uptake: a simplified model.  Transport in Porous Media, 142:469–498, 2022.
  * [8] Marco Berardi and Fabio V. Difonzo.  Strong solutions for Richards’ equation with Cauchy conditions and constant pressure gradient.  Environmental Fluid Mechanics, 20(1):165–174, Feb 2020.
  * [9] Marco Berardi, Fabio Vito Difonzo, and Luciano Lopez.  A mixed MoL-TMoL for the numerical solution of the 2D Richards’ equation in layered soils.  Computers & Mathematics with Applications, 79:1990–2001, 2020\.
  * [10] Luca Bergamaschi and Mario Putti.  Mixed finite elements and Newton-type linearizations for the solution of Richards’equation.  International Journal for Numerical methods in Engineering, 45:1025–1046, 1999.
  * [11] Florin Bobaru and Monchai Duangpanya.  A peridynamic formulation for transient heat conduction in bodies with evolving discontinuities.  Journal of Computational Physics, 231(7):2764–2785, 2012.
  * [12] C. Canuto and A. Quarteroni.  Approximation results for orthogonal polynomials in Sobolev spaces.  Math. Comp., 38:67–86, 1982.
  * [13] Andrea Carminati.  A Model of Root Water Uptake Coupled with Rhizosphere Dynamics.  Vadose Zone Journal, 11, 2012.
  * [14] Jesús Carrera, Xavier Sánchez-Vila, Inmaculada Benet, Agustín Medina, Germán Galarza, and Jordi Guimerà.  On matrix diffusion: formulations, solution methods and qualitative effects.  Hydrogeology Journal, 6:178–190, 1998.
  * [15] V. Casulli and P. Zanolli.  A Nested Newton-Type Algorithm for Finite Volume Methods Solving Richards’ Equation in Mixed Form.  SIAM Journal on Scientific Computing, 32(4):2255–2273, 2010.
  * [16] Michael A. Celia, Efthimios T. Bouloutas, and Rebecca L. Zarba.  A general mass-conservative numerical solution for the unsaturated flow equation.  Water Resources Research, 26(7):1483–1496, 1990.
  * [17] Antonio Coppola, Alessandro Comegna, Giovanna Dragonetti, Horst H. Gerke, and Angelo Basile.  Simulated Preferential Water Flow and Solute Transport in Shrinking Soils.  Vadose Zone Journal, 14(9):vzj2015.02.0021, 2015.
  * [18] Antonio Coppola, Horst H. Gerke, Alessandro Comegna, Angelo Basile, and Vincenzo Comegna.  Dual-permeability model for flow in shrinking soil with dominant horizontal deformation.  Water Resources Research, 48(8), 2012.
  * [19] Amir H. Delgoshaie, Daniel W. Meyer, Patrick Jenny, and Hamdi A. Tchelepi.  Non-local formulation for multiscale flow in porous media.  Journal of Hydrology, 531:649–654, 2015.
  * [20] F. V. Difonzo and R. Garrappa.  A numerical procedure for fractional-time-space differential equations with the spectral fractional Laplacian, 2022.  accepted on Springer INdAM Series.
  * [21] Fabio V. Difonzo, Costantino Masciopinto, Michele Vurro, and Marco Berardi.  Shooting the numerical solution of moisture flow equation with root uptake: a Python tool.  Water Resources Management, 35:2553–2567, 2021.
  * [22] Roberto Garrappa and Marina Popolizio.  Generalized exponential time differencing methods for fractional order problems.  Computers & Mathematics with Applications, 62(3):876–890, 2011\.  Special Issue on Advances in Fractional Differential Equations II.
  * [23] Ivan A. Guerrini and D. Swartzendruber.  Soil Water Diffusivity as Explicitly Dependent on Both Time and Water Content.  Soil Science Society of America Journal, 56(2):335–340, 1992.
  * [24] R. G. Hills, I. Porro, D. B. Hudson, and P. J. Wierenga.  Modeling one-dimensional infiltration into very dry soils: 1. Model development and evaluation.  Water Resources Research, 25(6):1259–1269, 1989.
  * [25] Rami Jabakhanji.  Peridynamic Modeling of Coupled Mechanical Deformations and Transient Flow in Unsaturated Soils.  PhD thesis, Purdue University, 2013.  https://docs.lib.purdue.edu/open_access_dissertations/147.
  * [26] Rami Jabakhanji and Rabi Mohtar.  A peridynamic model of flow in porous media.  Advances in Water Resources, 78, 02 2015.
  * [27] S. Jafarzadeh, L. Wang, A. Larios, and F. Bobaru.  A fast convolution-based method for peridynamic transient diffusion in arbitrary domains.  Computer Methods in Applied Mechanics and Engineering, 375:113633, 2021.
  * [28] M. L. Kavvas, A. Ercan, and J. Polsinelli.  Governing equations of transient soil water flow and soil water flux in multi-dimensional fractional anisotropic media and fractional time.  Hydrology and Earth System Sciences, 21(3):1547–1557, 2017.
  * [29] C.E. Kees, M.W. Farthing, and C.N. Dawson.  Locally conservative, stabilized finite element methods for variably saturated flow.  Computer Methods in Applied Mechanics and Engineering, 197(51):4610–4625, 2008.
  * [30] H. Li, M.W. Farthing, and C.T. Miller.  Adaptive local discontinuous Galerkin approximation to Richards’ equation.  Advances in Water Resources, 30(9):1883–1901, 2007.
  * [31] L. Lopez and S. F. Pellegrino.  A space-time discretization of a nonlinear peridynamic model on a 2D lamina.  Computers and Mathematics with Applications, 116:161–175, 2022\.
  * [32] L. Lopez and S.F. Pellegrino.  Computation of Eigenvalues for Nonlocal Models by Spectral Methods.  Journal of Peridynamics and Nonlocal Modeling. In Press. https://doi.org/10.1007/s42102-021-00069-8, 2021.
  * [33] L. Lopez and S.F. Pellegrino.  A fast-convolution based space–time Chebyshev spectral method for peridynamic models.  Advances in Continuous and Discrete Models. In press, 2022.
  * [34] L. Lopez and S.F. Pellegrino.  A non-periodic Chebyshev spectral method avoiding penalization techniques for a class of nonlinear peridynamic models.  International Journal for Numerical Methods in Engineering, 123(20):4859–4876, 2022.
  * [35] Gianmarco Manzini and Stefano Ferraris.  Mass-conservative finite volume methods on 2-D unstructured grids for the Richards’ equation.  Advances in Water Resources, 27(12):1199 – 1215, 2004.
  * [36] Costantino Masciopinto and Giuseppe Passarella.  Mass-transfer impact on solute mobility in porous media: A new mobile-immobile model.  Journal of Contaminant Hydrology, 215:21–28, 2018.
  * [37] K. Mitra and I.S. Pop.  A modified L-scheme to solve nonlinear diffusion problems.  Computers & Mathematics with Applications, 2018.
  * [38] Federico Municchi, Nicodemo Di Pasquale, Marco Dentz, and Matteo Icardi.  Heterogeneous Multi-Rate mass transfer models in OpenFOAM®.  Computer Physics Communications, 261:107763, 2021.
  * [39] Insa Neuweiler, Daniel Erdal, and Marco Dentz.  A Non-Local Richards Equation to Model Unsaturated Flow in Highly Heterogeneous Media under Nonequilibrium Pressure Conditions.  Vadose Zone Journal, 11(3), 08 2012.
  * [40] Tao Ni, Francesco Pesavento, Mirco Zaccariotto, Ugo Galvanetto, Qi-Zhi Zhu, and Bernhard A. Schrefler.  Hybrid FEM and peridynamic simulation of hydraulic fracture propagation in saturated porous media.  Computer Methods in Applied Mechanics and Engineering, 366:113101, 2020.
  * [41] Yakov Pachepsky, Dennis Timlin, and Walter Rawls.  Generalized Richards’ equation to simulate water transport in unsaturated soils.  Journal of Hydrology, 272(1):3–13, 2003.
  * [42] I.S. Pop, F. Radu, and P. Knabner.  Mixed finite elements for the Richards’ equation: linearization procedure.  Journal of Computational and Applied Mathematics, 168(1–2):365 – 373, 2004.
  * [43] Stephen L. Rawlins and Walter H. Gardner.  A Test of the Validity of the Diffusion Equation for Unsaturated Flow of Soil Water.  Soil Science Society of America Journal, 27(5):507–511, 1963.
  * [44] M. I. Romashchenko, V. O. Bohaienko, T. V. Matiash, V. P. Kovalchuk, and A. V. Krucheniuk.  Numerical simulation of irrigation scheduling using fractional Richards equation.  Irrigation Science, 39:385–396, 2021.
  * [45] David Seus, Koondanibha Mitra, Iuliu Sorin Pop, Florin Adrian Radu, and Christian Rohde.  A linear domain decomposition method for partially saturated flow in porous media.  Computer Methods in Applied Mechanics and Engineering, 333:331 – 355, 2018.
  * [46] S.A. Silling.  Reformulation of elasticity theory for discontinuities and long-range forces.  Journal of the Mechanics and Physics of Solids, 48(1):175–209, 2000\.
  * [47] S.A. Silling and R.B. Lehoucq.  Peridynamic Theory of Solid Mechanics.  In Hassan Aref and Erik van der Giessen, editors, Advances in Applied Mechanics, volume 44 of Advances in Applied Mechanics, pages 73–168. Elsevier, 2010.
  * [48] Xun Wu, Qiang Zuo, Jianchu Shi, Lichun Wang, Xuzhang Xue, and Alon Ben-Gal.  Introducing water stress hysteresis to the Feddes empirical macroscopic root water uptake model.  Agricultural Water Management, 240:106293, 2020.
  * [49] Huaxiang Yan, Majid Sedighi, and Andrey P. Jivkov.  Peridynamics modelling of coupled water flow and chemical transport in unsaturated porous media.  Journal of Hydrology, 591:125648, 2020.
  * [50] Jiří Šimůnek and Martinus Th. van Genuchten.  Contaminant Transport in the Unsaturated Zone: Theory and Modeling.  In J.H. Cushman and D.M Tartakovsky, editors, The Handbook of Groundwater Engineering, chapter 8, pages 266–290. CRC Press, Boca Raton, 2016\.