# Efficient Solvers for Nonstandard Models for Flow and Transport in
Unsaturated Porous Media

Davide Illiano<EMAIL_ADDRESS>Jakub Wiktor Both<EMAIL_ADDRESS>Iuliu Sorin Pop<EMAIL_ADDRESS>Florin Adrian Radu<EMAIL_ADDRESS>Department of Mathematics, University of Bergen, Allegaten 41, Bergen, Norway
Faculty of Science, University of Hasselt,Agoralaan Building D, BE 3590
Diepenbeek, Belgium

###### Abstract

We study several iterative methods for fully coupled flow and reactive
transport in porous media. The resulting mathematical model is a coupled,
nonlinear evolution system. The flow model component builds on the Richards
equation, modified to incorporate nonstandard effects like dynamic capillarity
and hysteresis, and a reactive transport equation for the solute. The two
model components are strongly coupled. On one hand, the flow affects the
concentration of the solute; on the other hand, the surface tension is a
function of the solute, which impacts the capillary pressure and,
consequently, the flow. After applying an Euler implicit scheme, we consider a
set of iterative linearization schemes to solve the resulting nonlinear
equations, including both monolithic and two splitting strategies. The latter
include a canonical nonlinear splitting and an alternate linearized splitting,
which appears to be overall faster in terms of numbers of iterations, based on
our numerical studies. The (time discrete) system being nonlinear, we
investigate different linearization methods. We consider the linearly
convergent L-scheme, which converges unconditionally, and the Newton method,
converging quadratically but subject to restrictions on the initial guess.
Whenever hysteresis effects are included, the Newton method fails to converge.
The L-scheme converges; nevertheless, it may require many iterations. This
aspect is improved by using the Anderson acceleration. A thorough comparison
of the different solving strategies is presented in five numerical examples,
implemented in MRST, a toolbox based on MATLAB.

##  1 Introduction

Mathematical models for complex physical phenomena are generally neglecting
several processes, in order to guarantee that the result is sufficiently
simple and to facilitate the numerical simulations. With a particular focus on
porous media applications, in this sense we mention enhanced oil recovery,
diffusion of substances in living tissues, and pollution of underground
aquifers. With the increase of computational power, and the development of
efficient simulation algorithms, mathematical models are improved
continuously, and more and more of the neglected effects are included.

When studying unsaturated flow, the equilibrium capillary pressure plays a
fundamental role. It is typically assumed to be a nonlinear, monotone function
of the water content. Explicit representations have been obtained thanks to
numerous experiments under equilibrium conditions (no flowing phases). Even
though this formulation has been the most commonly used in the last decades,
it has been observed [18, 25, 53, 65], that changes in time of the water
content, thus its time derivatives, do influence the profile of the capillary
pressure. In terms of modeling, this is achieved by including the so-called
dynamic effects [13, 30, 47]. Numerous papers investigate the existence of a
solution for systems including such effects, among them we cite [20, 41, 48].
Furthermore, the problem has been already studied numerically in, e.g., [2, 3,
21].

The hysteresis effect is another phenomenon often neglected. Again,
experiments have revealed that the curve obtained when investigating the
imbibition process, is different from the one observed during the drainage,
[23, 34, 50, 51]. This is sketched in Fig. 1.

Figure 1: Primary hysteresis loop as presented in [50].

In this article, we study unsaturated flow in porous media, modeled by the
Richards equation [14, 32], however including both dynamic and hysteresis
effects. Furthermore, we include a solute component, e.g., a surfactant, in
the wetting phase, which can directly influence the fluid properties ([4,
56]). The study of the transportation of an external components, e.g.,
surfactant, in variably saturated porous media has been already investigated
both numerically [36, 42, 59] and experimentally [33, 38].

Here, we will mainly concentrate on numerical studies, extending the solution
techniques in [36] to include dynamic capillarity and hysteresis We consider
the following model for coupled unsaturated flow and reactive transport

$\begin{split}\partial_{t}\theta(\Psi,c)-\nabla\cdot(K(\theta(\Psi,c))\nabla(\Psi+z))\
=\ \mathbb{S}_{1},\\\ \partial_{t}(\theta(\Psi,c)c)-\nabla\cdot(D\nabla
c-{\bf{u_{w}}}c)+R(c)\ =\ \mathbb{S}_{2}.\\\ \end{split}$ (1)

Here, $\theta(\Psi,c)$ is the water content, expressed as a function of both
the unknown pressure head $\Psi$ and the concentration of the external
component $c$. $K$, a function of the water content $\theta$, is the
conductivity, $z$ the vertical coordinate of $\vec{x}$, pointing against
gravity, $D$ the dispersion/diffusion coefficient,
${\bf{u_{w}}}:=-K(\theta(\Psi,c))\nabla(\Psi+z)$ the water flux, $R(c)$ the
reaction term and $\mathbb{S}_{1},\mathbb{S}_{2}$ the external sink/source
terms involved.

Next to a concentration dependence of $\theta$ and $\Psi$, here we include
also play-type hysteresis and dynamic capillary effects as introduced in [13].
More precisely,

$\Psi\in-
p_{cap}(\theta,c)+\tau(\theta)\partial_{t}\theta+\gamma(\theta)\operatorname{sign}(\partial_{t}\theta),$
(2)

where $p_{cap}$ is the equilibrium capillary pressure, expressed as a function
of $\theta$ and $c$, $\tau(\theta)$ the dynamic effects, and $\gamma(\theta)$
the width of the primary hysteresis loop. Later on, for ease of presentation,
we consider $\gamma$ as a positive constant, $\gamma\in\mathbb{R}_{\geq 0}$.
Note that (2) is a differential inclusion as the sign graph is multi-valued
and defined as follow,

$\begin{split}sign(\xi)=\begin{cases}1&\mbox{for}\ \xi>0,\\\
[-1,1]&\mbox{for}\ \xi=0,\\\ -1&\mbox{for}\ \xi<0.\end{cases}\end{split}$ (3)

The multi-valued graph allows switching between the imbibition and drainage
curves in the play-type hysteresis. For more details on the formulation we
refer to [13].

The primary unknowns of the system are the pressure $\Psi$, the concentration
$c$ and the water content $\theta$. In standard models, also obtained as
special case for $\gamma=\tau(\theta)=0$, $\theta$ is a function of pressure
and concentration. Therefore, (2) is replaced by an algebraic relationship,
which simplifies the model and allows eliminating $\theta$ as an unknown. In
the extended/nonstandard formulation, $\theta$ is an unknown and (2) is
required as additional equation of the model. Initial and boundary conditions
will complete the system.

To avoid working with a graph, we consider the following regularization,

$\Phi(\xi)=\begin{cases}\operatorname{sign}(\xi)&\mbox{if}\qquad|\xi|\geq\delta,\\\
\frac{\xi}{\delta}&\mbox{if}\qquad|\xi|<\delta,\end{cases}$ (4)

where $\delta\in\mathbb{R}^{+}$ is a small parameter. Using this in (2) gives
the regularized system of equations

$\begin{split}\partial_{t}\theta-\nabla\cdot(K(\theta)\nabla(\Psi+z))\ &=\
\mathbb{S}_{1},\\\
\Psi=-p_{cap}(\theta,c)+\tau(\theta)\partial_{t}\theta&+\gamma\Phi(\partial_{t}\theta),\\\
\partial_{t}(\theta c)-\nabla\cdot(D\nabla c-{\bf{u_{w}}}c)+R(c)\ &=\
\mathbb{S}_{2}.\\\ \end{split}$ (5)

From now on, the system (5) will be further investigated. We will discretize
the equations and study different solving algorithms.

######  Remark 1.

An inverse formulation is proposed in [13], obtained by solving (2), as its
regularized counterpart in (5), in terms of $\partial_{t}\theta$. This gives

$\partial_{t}\theta=F(\Psi,\theta,c),$ (6)

for a suitable function $F$. The time derivative in the flow equation can then
be substituted by $F$,

$F(\Psi,\theta,c)-\nabla\cdot\big{(}K(\theta)\nabla(\Psi+z)\big{)}\ =\
\mathbb{S}_{1}.$ (7)

This formulation is used for the mathematical analysis of such models,
[19, 61]. It has been observed, e.g., in [46], that such formulation can
reduce the number of iterations required to solve the system of equations,
compared to the formulation in (2). However, for the particular test cases
investigated here, no remarkable improvements are observed. Thus, for ease of
presentation, we will report the results obtained only for the formulation
given by (5).

We point out that the concentration of the external component directly
influences the capillary pressure. The presence of such a component results in
a non-constant surface tension, which induces a rescaling of the pressures
[35, 63, 64].

To solve the system (5) numerically, one first needs to discretize in time and
space, and then develop solvers for the discretized equations. In this paper,
due to the expected low regularity of the solutions [6] and the desire of
relatively large time steps, we choose to use the backward Euler method for
the time discretization. Certain processes investigated in porous media flow
can take place on time intervals longer than decades, thus the need for large
time steps. Multiple spatial discretization techniques are available, e.g.,
the Galerkin Finite Element Method (_FEM_) [11, 52, 60], Discontinuous
Galerkin Method (_DGM_) [9, 39, 43, 66], the Mixed Finite Element Method
(_MFEM_) [8, 21, 58, 59, 67, 70], the Finite Volume Method (_FVM_)
[27] and the Multi-Point Flux Approximation (_MPFA_)
[1, 10, 12, 40]. We will here concentrate on _FEM_ and _TPFA_ (Two Points Flux
Approximation), a particular case of _MPFA_. In particular, we cite [15, 26,
71] for papers on improved numerical schemes applied to the Richards equation.

Since the equations investigated here are characterized by several nonlinear
quantities, $K(\theta)$, $p_{cap}(\theta,c)$, $\tau(\partial_{t}\theta)$, and
$R(c)$, and the time discretization is not explicit, one needs to solve a
nonlinear system at each time step, requiring a linearization procedure.
Examples of possible linearization schemes are: the Newton method [54], the
modified Picard method [24] and the L-scheme [45, 57]. In this paper, we
investigate the Newton method and the L-scheme. The former is a commonly used
linearization scheme which is quadratically convergent. However, this
convergence is only local and one needs to compute the Jacobian matrix, which
can be expensive. The L-scheme is instead globally (linearly) convergent,
under mild restrictions, and it does not require the computation of any
derivative. The L-scheme is in general slower in terms of numbers of
iterations than the Newton method. Moreover, the linear systems to be solved
within each iteration are better conditioned when compared to the ones given
by the Newton method [36, 45]. Furthermore, the rate of convergence of the
scheme strongly depends on user-defined parameters. Such aspects are
investigated for numerous nonlinear problems, including Richards equation, and
two-phase flow in porous media, in [36, 45, 49, 57, 62]. Finally there
numerous papers proposing improved formulation of the L-scheme, among them we
cite [5, 49].

In this work, we test the L-scheme on more complex problems involving
hysteresis and dynamic effects, and coupled reactive transport and flow.
Furthermore, we investigate a post-processing technique, the Anderson
Acceleration (AA) [7], which can drastically improve linearly convergent
schemes. The acceleration tool requires user-defined parameters. As will be
seen below, choosing the suitable parameters for the AA, significantly relaxes
the choice of the parameters for the L-scheme linearization.

We observe that the system (5) is fully coupled. This is due to the dependence
of the capillary pressure on both $\theta$ and $c$. Therefore, we will
investigate multiple solution algorithms, combining different linearization
schemes and decoupling techniques. Decoupling/splitting the equations may
present multiple advantages such as: an easier implementation, a better
conditioned problem to solve, similar convergence properties but faster
computations. We divide the schemes into three main categories: monolithic
(Mono), nonlinear splitting (NonLinS) and alternate splitting (AltS).
Subsequently, we denote, e.g., by Newton-Mono, the monolithic scheme obtained
by applying the Newton method as linearization. Such schemes have already been
investigated for the standard model in [36].

The paper is organized as follows. In Section 2, we present the linearization
and discretization techniques including monolithic or decoupled solution
approaches. Section 3 presents five different numerical examples, which allow
to compare the efficiency and robustness of the solving algorithms. Section 4
concludes this work with the final remarks.

##  2 Problem formulation, discretization and iterative schemes

In the following, we use the standard notations of functional analysis. The
domain $\Omega\subset\mathbb{R}^{d}$, $d=1,2$ or $3$, is bounded and has a
Lipschitz continuous boundary $\partial\Omega$. The final time is $T>0$, and
the time domain is $(0,T]$. $L^{2}(\Omega)$ denotes the space of real valued,
square integrable functions defined on $\Omega$ and $H^{1}(\Omega)$ its
subspace containing the functions also having weak first derivatives in
$L^{2}(\Omega)$. $H_{0}^{1}(\Omega)$ is the space of functions belonging to
$H^{1}(\Omega)$, having zero trace on the boundary $\partial\Omega$.
Furthermore, we denote by $<\cdot,\cdot>$ the standard $L^{2}(\Omega)$ scalar
product and by $\left\lVert\cdot\right\rVert$ the associated norm.

To numerically solve the system of equations (5), one needs to discretize both
in time and space. We combine the backward Euler method with linear Galerkin
finite elements. Let $N\in\mathbb{N}$ be a strictly positive natural number.
We define the time step size $\Delta t\ =\ T/N$ and $t_{n}\ =\ n\Delta t\
(n\in{1,2,\dots,N})$. Furthermore, let $T_{h}$ be a regular decomposition of
$\Omega$, $\overline{\Omega}=\underset{T\in T_{h}}{\cup}T$, with $h$ denoting
the mesh diameter. The finite element spaces $V_{h}\subset H_{0}^{1}(\Omega)$
and $W_{h}\subset L^{2}(\Omega)$ are defined by

$V_{h}:=\left\\{v_{h}\in H_{0}^{1}(\Omega)\ s.t.\
v_{h|T}\in\mathbb{P}_{1}(T),\ T\in T_{h}\right\\},\qquad
W_{h}:=\left\\{w_{h}\in L^{2}(\Omega)\ s.t.\ w_{h|T}\in\mathbb{P}_{1}(T),\
T\in T_{h}\right\\},$ (8)

where $\mathbb{P}_{1}(T)$ denotes the space of the linear polynomials on $T$.
The fully discrete Galerkin formulation of the system (5) can now be written
as:

Problem Pn: Let $n\geq 1$ be fixed. Assuming that
$\Psi^{n-1}_{h},c^{n-1}_{h}\in V_{h}$ and $\theta^{n-1}_{h}\in W_{h}$ are
given, find $\Psi^{n}_{h},c^{n}_{h}\in V_{h}$ and $\theta^{n}_{h}\in W_{h}$
such that

$\begin{split}<\theta^{n}_{h}-\theta^{n-1}_{h},v_{1,h}>+\Delta
t<K(\theta^{n}_{h})(\nabla\Psi^{n}_{h}+{{\bf{e_{z}}}}),\nabla v_{1,h}>\ =\
\Delta t<\mathbb{S}_{1},v_{1,h}>\\\ \Delta t<\Psi^{n}_{h},w_{1,h}>+\Delta
t<p_{cap}(\theta^{n}_{h},c^{n}_{h}),w_{1,h}>-<\tau(\theta^{n}_{h})(\theta^{n}_{h}-\theta^{n-1}_{h}),w_{1,h}>\
=\ \Delta t\gamma<\Phi\left(\frac{\theta^{n}_{h}-\theta^{n-1}_{h}}{\Delta
t}\right),w_{1,h}>\\\
<\theta^{n}_{h}(c^{n}_{h}-c^{n-1}_{h})+c^{n}_{h}(\theta^{n}_{h}-\theta^{n-1}_{h}),v_{2,h}>+\Delta
t<D\nabla\Psi^{n}_{h}+{\bf{u^{n}_{w}}}c^{n}_{h},\nabla v_{2,h}>+\Delta
t<R(c^{n}_{h}),v_{2,h}>\ =\ \Delta t<\mathbb{S}_{2},v_{2,h}>\end{split}$

(9)

holds for all $v_{1,h},v_{2,h}\in V_{h}$ and for all $w_{1,h}\in W_{h}$. We
denote by ${{\bf{e_{z}}}}$ the unit vector in the direction opposite to
gravity.

Observe that choosing the space $H^{1}_{0}(\Omega)$ implies that homogeneous
boundary conditions have been adopted for the pressure and the concentration.
However, this choice is made for the ease of presentation, the extension to
other boundary conditions being possible without major complications. We also
mention that, for $n=1$, we use the approximation in $V_{h}$ of the initial
water content and concentration, respectively $\theta^{0}_{h}$ and
$c^{0}_{h}$.

In the following, we investigate different iterative schemes for solving
Problem Pn. These schemes are based on the ones discussed in [36], extending
them, not only to the case of dynamic capillary pressure $(\tau(\theta)\neq
0)$ [37], but also to the case of hysteresis. Among the numerous papers
investigating numerically the effects of hysteresis and dynamic capillarity
pressure, we cite [55, 72]. As mentioned, we compare monolithic (Mono) and
splitting (NonLinS and AltS) solvers, combined with two different
linearization schemes, the Newton method and the L-scheme. Furthermore, the
Anderson acceleration [7] will be taken into account to speed up the linearly
convergent L-scheme.

###  2.1 Solving algorithms

In what follows, when solving (9) iteratively, the index $n$ will always refer
to the time step level, whereas $j$ will denote the iteration index. As a
rule, the iterations will start with the solution at the last time step,
$t_{n-1}$, for example $\Psi^{n,1}=\Psi^{n-1}$. As mentioned, this choice is
not required for L-type schemes but it is a natural one.

In a compact form Problem Pn can be seen as the system

$\begin{cases}F_{1}(\Psi_{h}^{n},\theta_{h}^{n})&=0,\\\
F_{2}(\Psi_{h}^{n},\theta_{h}^{n},c_{h}^{n})&=0,\\\
F_{3}(\Psi_{h}^{n},\theta_{h}^{n},c_{h}^{n})&=0,\end{cases}$ (10)

with $F_{1},F_{2}$ resulting from the flow equations and $F_{3}$ from the
transport. In the following we will indicate with $F^{lin}$, the linearized
formulation of $F$ obtained by either the Newton method or the L-scheme.
Finally, we can proceed to present monolithic and splitting solvers.

In the monolithic approach one solves the three equations of the system (10)
at once. Formally, one iteration is:
Find $\Psi_{h}^{n,j+1}$, $\theta_{h}^{n,j+1}$ and $c_{h}^{n,j+1}$ such that

$\begin{cases}F^{lin}_{1}(\Psi_{h}^{n,j+1},\theta_{h}^{n,j+1})&=0,\\\
F^{lin}_{2}(\Psi_{h}^{n,j+1},\theta_{h}^{n,j+1},c_{h}^{n,j+1})&=0,\\\
F^{lin}_{3}(\Psi_{h}^{n,j+1},\theta_{h}^{n,j+1},c_{h}^{n,j+1})&=0,\end{cases}$
(11)

where $F_{i}^{Lin}$ is the linearization of $F_{i}$, $i\in\\{1,2,3\\}$.
Depending on which linearization technique is used, we refer to the Newton-
monolithic scheme (Newton-Mono) or monolithic-$L$-scheme (LS-Mono). These two
schemes will be presented in details below. Fig. 2 displays the sketched
version of the monolithic solver.

Figure 2: The monolithic approach.

In the iterative splitting approach, the flow and the transport equations are
solved subsequently, iterating between them. We will distinguish between two
primary splitting schemes: the nonlinear splitting (NonLinS) and the alternate
linearized splitting (AltLinS), illustrated in Figure 3 and Figure 4,
respectively. Such schemes have already been studied, in the case of a
standard flow model in [36].
In the nonlinear splitting, one iteration step is:
Find first $\Psi_{h}^{n,j+1}$, $\theta_{h}^{n,j+1}$ such that

$\begin{split}\begin{cases}F_{1}(\Psi_{h}^{n,j+1},\theta_{h}^{n,j+1})&=0,\\\
F_{2}(\Psi_{h}^{n,j+1},\theta_{h}^{n,j+1},c_{h}^{n,j})&=0,\end{cases}\\\
\end{split}$ (12)

and then find $c_{h}^{n,j+1}$ such that

$\begin{split}F_{3}(\Psi_{h}^{n,j+1},\theta_{h}^{n,j+1},c_{h}^{n,j+1})=0.\end{split}$
(13)

The two flow equations are solved at once. Each of the nonlinear systems (12)
and (13) is solved until some convergence criterion is met. Once the pressure
and water content are obtained, $\Psi_{h}^{n,j+1}$ and $\theta_{h}^{n,j+1}$,
are then used in the transport equation (13) to compute $c_{h}^{n,j+1}$. The
resulting $F_{1}$, $F_{2}$ and $F_{3}$, being nonlinear, are linearized using
the Newton method or the $L$-scheme.

In contrast, the alternate linearized splitting (AltLinS) schemes perform only
one linearization step per iteration, see Figure 4. One iteration in the
alternate splitting scheme can be written as:
Find $\Psi_{h}^{n,j+1},\theta_{h}^{n,j+1}$ such that

$\begin{split}\begin{cases}F^{lin}_{1}(\Psi_{h}^{n,j+1},\theta_{h}^{n,j+1})&=0,\\\
F^{lin}_{2}(\Psi_{h}^{n,j+1},\theta_{h}^{n,j+1},c_{h}^{n,j})&=0,\end{cases}\\\
\end{split}$ (14)

and then $c_{h}^{n,j+1}$ such that

$\begin{split}F^{lin}_{3}(\Psi_{h}^{n,j+1},\theta_{h}^{n,j+1},c_{h}^{n,j+1})=0.\end{split}$
(15)

Again, depending on which linearization is used, we refer to alternate
splitting Newton (AltS-Newton) or alternate splitting $L$-scheme (AltS-LS).
Both schemes will be presented in detail below.

In the following sections we will illustrate, in the details, the different
schemes here investigated.

Figure 3: The nonlinear splitting approach. Figure 4: The alternate splitting
approach.

####  2.1.1 The monolithic Newton method (Newton-Mono)

The standard monolithic Newton method applied to (9) reads as:

Problem P-Newton-Mono: Let $j>1$ be fixed. Let
$\Psi^{n-1}_{h},\Psi^{n,j}_{h},c_{h}^{n-1},c_{h}^{n,j}\in V_{h}$, and
$\theta^{n-1}_{h},\theta^{n,j}_{h}\in W_{h}$ be given, find
$\Psi^{n,j+1}_{h},c_{h}^{n,j+1}\in V_{h}$, and $\theta^{n,j+1}_{h}\in W_{h}$
such that

$\begin{split}&<\theta^{n,j+1}_{h}-\theta^{n-1}_{h},v_{1,h}>+\Delta
t<K(\theta^{n,j}_{h})(\nabla(\Psi^{n,j+1}_{h})+{\bf{e_{z}}}),\nabla
v_{1,h}>\\\ +\Delta
t&<\partial_{\theta}K(\theta^{n,j}_{h})(\nabla(\Psi^{n,j}_{h})+{\bf{e_{z}}})(\theta^{n,j+1}_{h}-\theta^{n,j}_{h}),\nabla
v_{1,h}>=\Delta t<\mathbb{S}_{1},v_{1,h}>\end{split}$ (16)
$\begin{split}\Delta t<\Psi^{n,j+1}_{h},&w_{1,h}>=-\Delta
t<p_{cap}(\theta^{n,j}_{h},c^{n,j}_{h}),w_{1,h}>-\Delta
t<\partial_{\theta}p_{cap}(\theta^{n,j}_{h},c^{n,j}_{h})(\theta^{n,j+1}_{h}-\theta^{n,j}_{h}),w_{1,h}>\\\
&-\Delta
t<\partial_{c}p_{cap}(\theta^{n,j}_{h},c^{n,j}_{h})(c^{n,j+1}_{h}-c^{n,j}_{h}),w_{1,h}>+<\tau(\theta^{n,j}_{h})(\theta^{n,j+1}_{h}-\theta^{n-1}_{h}),w_{1,h}>\\\
&+<\partial_{\theta}\tau(\theta^{n,j}_{h})(\theta^{n,j}_{h}-\theta^{n-1}_{h})(\theta^{n,j+1}_{h}-\theta^{n,j}_{h}),w_{1,h}>+\Delta
t\gamma<\Phi\left(\frac{\theta^{n,j}_{h}-\theta^{n-1}_{h}}{\Delta
t}\right),w_{1,h}>\\\ \end{split}$ (17)

and

$\begin{split}<\theta^{n,j}_{h}&(c^{n,j+1}_{h}-c^{n-1}_{h})+c^{n,j}_{h}(\theta^{n,j+1}_{h}-\theta^{n-1}_{h}),v_{2,h}>+\Delta
t<D\nabla c^{n,j+1}_{h}+{\bf{u^{n,j}_{w}}}c^{n,j+1}_{h},\nabla v_{2,h}>\\\
&+\Delta t<R(c^{n,j}_{h}),v_{2,h}>+\Delta
t<\partial_{c}R(c^{n,j}_{h})(c^{n,j+1}_{h}-c^{n,j}_{h})>=\Delta
t<\mathbb{S}_{2},v_{2,h}>\end{split}$ (18)

hold true for all $v_{1,h},v_{2,h}\in V_{h}$, and for all $w_{1,h}\in W_{h}$.
By $\partial_{\theta}$ we denote the partial derivative with respect to the
water content $\theta$, and by $\partial_{c}$ the partial derivative with
respect to the concentration $c$, and
${\bf{u^{n,j}_{w}}}:=-K(\theta^{n,j}_{h})\nabla(\Psi^{n,j}_{h}+{\bf{e_{z}}})$.

####  2.1.2 The monolithic $L$-scheme (LS-Mono)

The monolithic $L$-scheme for solving (9) reads:
Problem P-LS-Mono: Let $j>1$ be fixed. Let
$\Psi^{n-1}_{h},\Psi^{n,j}_{h},c^{n-1}_{h},c_{h}^{n,j}\in V_{h}$, and
$\theta^{n-1}_{h},\theta^{n,j}_{h}\in W_{h}$ be given, find
$\Psi^{n,j+1}_{h},c^{n,j+1}_{h}\in V_{h}$, and $\theta^{n,j+1}_{h}\in W_{h}$
such that

$<\theta^{n,j+1}_{h}-\theta^{n-1}_{h},v_{1,h}>+\Delta
t<K(\theta^{n,j}_{h})(\nabla(\Psi^{n,j+1}_{h})+{\bf{e_{z}}}),\nabla
v_{1,h}>+L_{1}<\Psi^{n,j+1}_{h}-\Psi^{n,j}_{h},v_{1,h}>=\Delta
t<\mathbb{S}_{1},v_{1,h}>$ (19) $\begin{split}\Delta
t<\Psi^{n,j+1}_{h},w_{1,h}>&=-\Delta
t<p_{cap}(\theta^{n,j}_{h},c^{n,j}_{h}),w_{1,h}>+<\tau(\theta^{n,j}_{h})(\theta^{n,j+1}_{h}-\theta^{n-1}_{h}),w_{1,h}>\\\
&+\Delta t\gamma<\Phi\left(\frac{\theta^{n,j}_{h}-\theta^{n-1}_{h}}{\Delta
t}\right),w_{1,h}>+L_{2}<(\theta^{n,j+1}_{h}-\theta^{n,j}_{h}),w_{1,h}>\\\
\end{split}$ (20)

and

$\begin{split}<\theta^{n,j}_{h}&(c^{n,j+1}_{h}-c^{n-1}_{h})+c^{n,j}_{h}(\theta^{n,j+1}_{h}-\theta^{n-1}_{h}),v_{2,h}>+\Delta
t<D\nabla c^{n,j+1}_{h}+{\bf{u^{n,j}_{w}}}c^{n,j+1}_{h},\nabla v_{2,h}>\\\
&+\Delta
t<R(c^{n,j}_{h}),v_{2,h}>+L_{3}<c^{n,j+1}_{h}-c^{n,j}_{h},v_{2,h}>=\Delta
t<\mathbb{S}_{2},v_{2,h}>\end{split}$ (21)

hold true for all $v_{1,h},v_{2,h}\in V_{h}$, and for all $w_{1,h}\in W_{h}$.
$L_{1},L_{2}$ and $L_{3}$ are three positive, user-defined parameters on which
only mild conditions are imposed. We refer to [36, 45, 57] for the analysis of
the numerical schemes which have inspired the ones presented here. Often, one
needs to properly tune these parameters to obtain a robust and relatively fast
solver.

####  2.1.3 The nonlinear splitting approach (NonLinS)

The nonlinear splitting approach for solving (9) reads:

Problem P-NonLinS: Let $j>1$ be fixed. Let
$\Psi^{n-1}_{h},\Psi^{n,j}_{h},c^{n-1}_{h},c^{n,j}_{h}\in V_{h}$ and
$\theta^{n-1}_{h},\theta^{n,j}_{h}\in W_{h}$ be given, find
$\Psi^{n,j+1}_{h}\in V_{h}$, and $\theta^{n,j+1}_{h}\in W_{h}$ such that

$<\theta^{n,j+1}_{h}-\theta^{n-1}_{h},v_{1,h}>+\Delta
t<K(\theta^{n,j+1}_{h})(\nabla(\Psi^{n,j+1}_{h})+{\bf{e_{z}}}),\nabla
v_{1,h}>=\Delta t<\mathbb{S}_{1},v_{1,h}>$ (22) $\begin{split}\Delta
t<\Psi^{n,j+1}_{h},w_{1,h}>=-\Delta
t<&p_{cap}(\theta^{n,j+1}_{h},c^{n,j}_{h}),w_{1,h}>+<\tau(\theta^{n,j+1}_{h})(\theta^{n,j+1}_{h}-\theta^{n-1}_{h}),w_{1,h}>\\\
&+\Delta t\gamma<\Phi\left(\frac{\theta^{n,j}_{h}-\theta^{n-1}_{h}}{\Delta
t}\right),w_{1,h}>\end{split}$ (23)

holds true for all $v_{1,h}\in V_{h}$ and for all $w_{1,h}\in W_{h}$.

Then let $\Psi^{n-1}_{h},\Psi^{n,j}_{h},c^{n-1}_{h},c^{n,j}_{h}\in V_{h}$ and
$\theta^{n-1}_{h},\theta^{n,j}_{h}\in W_{h}$ be given, $\Psi^{n,j+1}_{h}\in
V_{h}$ and $\theta^{n,j+1}_{h}\in W_{h}$ are obtained from the equations
above, find $c^{n,j+1}_{h}\in V_{h}$ such that

$\begin{split}<\theta^{n,j+1}_{h}(c^{n,j+1}_{h}-c^{n-1}_{h})+c^{n,j}_{h}&(\theta^{n,j+1}_{h}-\theta^{n-1}_{h}),v_{2,h}>+\Delta
t<D\nabla c^{n,j+1}_{h}+{\bf{u^{n,j+1}_{w}}}c^{n,j+1}_{h},\nabla v_{2,h}>\\\
&+\Delta t<R(c^{n,j+1}_{h}),v_{2,h}>=\Delta
t<\mathbb{S}_{2},v_{2,h}>\end{split}$ (24)

holds true for all $v_{2,h}\in V_{h}$. The water flux is given by
${\bf{u^{n,j+1}_{w}}}:=-K(\theta^{n,j+1}_{h})\nabla(\Psi^{n,j+1}_{h}+{\bf{e_{z}}})$.

Observe that (22)–(23) and (24) are nonlinear. To approximate their respective
solutions, one can employ, e.g., the Newton method (NonLinS-Newton) or the
$L$-scheme (NonLinS-LS).

####  2.1.4 The alternate splitting Newton method (Newton-AltLinS)

Applied to (9), the alternate splitting Newton method reads:

Problem P-Newton-AltLinS: Let $j>1$ be fixed. Let
$\Psi^{n-1}_{h},\Psi^{n,j}_{h},c^{n-1}_{h},c^{n,j}_{h}\in V_{h}$ and
$\theta^{n-1}_{h},\theta^{n,j}_{h}\in W_{h}$ be given, find
$\Psi^{n,j+1}_{h}\in V_{h}$, and $\theta^{n,j+1}_{h}\in W_{h}$ such that

$\begin{split}&<\theta^{n,j+1}_{h}-\theta^{n-1}_{h},v_{1,h}>+\Delta
t<K(\theta^{n,j}_{h})(\nabla(\Psi^{n,j+1}_{h})+{\bf{e_{z}}}),\nabla
v_{1,h}>\\\ +\Delta
t&<\partial_{\theta}K(\theta^{n,j}_{h})(\nabla(\Psi^{n,j}_{h})+{\bf{e_{z}}})(\theta^{n,j+1}_{h}-\theta^{n,j}_{h}),\nabla
v_{1,h}>=\Delta t<\mathbb{S}_{1},v_{1,h}>\end{split}$ (25)
$\begin{split}\Delta t<\Psi^{n,j+1}_{h},w_{1,h}>=-\Delta
t<p_{cap}&(\theta^{n,j}_{h},c^{n,j}_{h}),w_{1,h}>-\Delta
t<\partial_{\theta}p_{cap}(\theta^{n,j}_{h},c^{n,j}_{h})(\theta^{n,j+1}_{h}-\theta^{n,j}_{h}),w_{1,h}>\\\
+<\tau(\theta^{n,j}_{h})(\theta^{n,j+1}_{h}-&\theta^{n-1}_{h}),w_{1,h}>+<\partial_{\theta}\tau(\theta^{n,j}_{h})(\theta^{n,j}_{h}-\theta^{n-1}_{h})(\theta^{n,j+1}_{h}-\theta^{n,j}_{h}),w_{1,h}>\\\
&+\Delta t\gamma<\Phi\left(\frac{\theta^{n,j}_{h}-\theta^{n-1}_{h}}{\Delta
t}\right),w_{1,h}>\end{split}$ (26)

hold true for all $v_{1,h}\in V_{h}$ and $w_{1,h}\in W_{h}$.

Then, with given $\Psi^{n-1}_{h},\Psi^{n,j}_{h},c^{n-1}_{h},c^{n,j}_{h}\in
V_{h}$ and $\theta^{n-1}_{h},\theta^{n,j}_{h}\in W_{h}$, $\Psi^{n,j+1}_{h}\in
V_{h}$ and $\theta^{n,j+1}_{h}\in W_{h}$ are obtained from the equations
above, find $c^{n,j+1}_{h}\in V_{h}$ such that

$\begin{split}<\theta^{n,j+1}_{h}&(c^{n,j+1}_{h}-c^{n-1}_{h})+c^{n,j}_{h}(\theta^{n,j+1}_{h}-\theta^{n-1}_{h}),v_{2,h}>+\Delta
t<D\nabla c^{n,j+1}_{h}+{\bf{u^{n,j+1}_{w}}}c^{n,j+1}_{h},\nabla v_{2,h}>\\\
&+\Delta t<R(c^{n,j}_{h}),v_{2,h}>+\Delta
t<\partial_{c}R(c^{n,j}_{h})(c^{n,j+1}_{h}-c^{n,j}_{h}),v_{2,h}>=\Delta
t<\mathbb{S}_{2},v_{2,h}>\end{split}$ (27)

hold true for all $v_{2,h}\in V_{h}$.

####  2.1.5 The alternate splitting $L$-scheme (LS-AltLinS)

The alternate splitting $L$-scheme for solving (9) is:

Problem P-LS-AltLinS: Let $j>1$ be fixed. Let
$\Psi^{n-1}_{h},\Psi^{n,j}_{h},c^{n-1}_{h},c^{n,j}_{h}\in V_{h}$ and
$\theta_{h}^{n-1},\theta^{n,j}_{h}\in W_{h}$ be given, find
$\Psi^{n,j+1}_{h}\in V_{h}$, and $\theta^{n,j+1}_{h}\in W_{h}$ such that

$\begin{split}<\theta^{n,j+1}_{h}-\theta^{n-1}_{h},v_{1,h}>+\Delta
t<K(\theta^{n,j}_{h})(\nabla(\Psi^{n,j+1}_{h})+{\bf{e_{z}}}),\nabla
v_{1,h}>+L_{1}<\Psi^{n,j+1}_{h}-\Psi^{n,j}_{h},v_{1,h}>=\Delta
t<\mathbb{S}_{1},v_{1,h}>\end{split}$ (28) $\begin{split}\Delta
t<\Psi^{n,j+1}_{h},&w_{1,h}>=-\Delta
t<p_{cap}(\theta^{n,j}_{h},c^{n,j}_{h}),w_{1,h}>+<\tau(\theta^{n,j}_{h})(\theta^{n,j+1}_{h}-\theta^{n-1}_{h}),w_{1,h}>\\\
&+L_{2}<\theta^{n,j+1}_{h}-\theta^{n,j}_{h},w_{1,h}>+\Delta
t\gamma<\Phi\left(\frac{\theta^{n,j}_{h}-\theta^{n-1}_{h}}{\Delta
t}\right),w_{1,h}>\end{split}$ (29)

holds true for all $v_{1,h}\in V_{h}$ and $w_{1,h}\in W_{h}$.

Then, with given $\Psi^{n-1}_{h},\Psi^{n,j}_{h},c^{n-1}_{h},c^{n,j}_{h}\in
V_{h}$ and $\theta^{n-1}_{h},\theta^{n,j}_{h}\in W_{h}$, and
$\Psi^{n,j+1}_{h}\in V_{h}$ and $\theta^{n,j+1}_{h}\in W_{h}$ from the
equations above. We find $c^{n,j+1}_{h}\in V_{h}$ such that

$\begin{split}<\theta^{n,j+1}_{h}&(c^{n,j+1}_{h}-c^{n-1}_{h})+c^{n,j}_{h}(\theta^{n,j+1}_{h}-\theta^{n-1}_{h}),v_{2,h}>+\Delta
t<D\nabla c^{n,j+1}_{h}+{\bf{u^{n,j+1}_{w}}}c^{n,j+1}_{h},\nabla v_{2,h}>\\\
&+\Delta
t<R(c^{n,j}_{h}),v_{2,h}>+L_{3}<c^{n,j+1}_{h}-c^{n,j}_{h},v_{2,h}>=\Delta
t<\mathbb{S}_{2},v_{2,h}>\end{split}$ (30)

hold true for all $v_{2,h}\in V_{h}$.

######  Remark 2.

There exist multiple improved formulations of both the Newton method and
L-scheme. We refer, among others, to the trust region techniques [69], and the
modified L-scheme in [49].

######  Remark 3.

(Stopping criterion) For all schemes (monolithic or splitting), the iterations
are stopped when,

$\left\lVert\Psi_{h}^{n,j+1}-\Psi_{h}^{n,j}\right\rVert_{\infty}\leq\epsilon_{1},\qquad\left\lVert\theta_{h}^{n,j+1}-\theta_{h}^{n,j}\right\rVert_{\infty}\leq\epsilon_{2}\qquad\text{and}\qquad\left\lVert
c_{h}^{n,j+1}-c_{h}^{n,j}\right\rVert_{\infty}\leq\epsilon_{3},$

where by $\left\lVert\cdot\right\rVert_{\infty}$ we mean the
$L^{\infty}(\Omega)$ norm. Later on, for ease of presentation, we consider
$\epsilon_{1}=\epsilon_{2}=\epsilon_{3}=\epsilon$. The parameter $\epsilon$
will be defined in the numerical section.

###  2.2 Anderson acceleration

Although the L-scheme is robust and converges under mild restrictions, the
convergence rate depends strongly on the linearization parameters. We refer to
[45, 57, 62] for the analysis in case of standard Richards equation, and to
[39] for the nonstandard model. Tuning the parameters to obtain optimal
results in terms of numbers of iterations and thus of computational times, can
be tedious and time-consuming. The Anderson Acceleration (AA) is a powerful
post-processing tool which can drastically reduce the numbers of iterations
required by linearly convergent schemes, such as the L-scheme here
investigated. In addiction, it reduces the need for finding close to optimal
linearization parameters.

D. G. Anderson introduced the acceleration tool in 1965 [7], and since then it
has been investigated in multiple works, to name a few [16, 28, 68]. We recall
here the definition of AA, presented in [68], formulated for a general fixed
point problem, of the form: given $g:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$,
solve the system $x=g(x)$.

Algorithm 1 Classical Fixed-Point iteration

1:Given $x_{0}$

2:for $k=0,1,$… until convergence do

3: $x_{k+1}=g(x_{k})$

4:end for

Opposed to utilize only the last iteration $x_{k}$, in the AA the new
approximation is a linear combination of previously computed ones, see
Algorithm 2. In the following, we denote by AA$(m)$ the Anderson acceleration
where $m+1$ previously computed iterates are taken into account. With this,
AA$(0)$ is the non-accelerated formulation.

Algorithm 2 Anderson Acceleration AA$(m)$

1:Given $x_{0}$

2:for $k=1,2$… until convergence do

3: Set $m_{k}=min\\{m,k-1\\}$

4: Define the matrix $F_{k}=(f_{k-m_{k}-1},\cdots,f_{k-1})$, where
$f_{i}=g(x_{i})-x_{i}$

5: Find $\alpha\in\mathbb{R}^{m_{k}+1}$ that solves
$\min\limits_{\alpha=(\alpha_{0},\cdots,\alpha_{m_{k}})^{T}}\left\lVert
F_{k}\alpha\right\rVert\ \text{s.t.}\ \sum_{i=0}^{m_{k}}\alpha_{i}=1.$

6: Define $x_{k}:=\sum_{i=0}^{m_{k}}\alpha_{i}g(x_{k-m_{k}+i-1})$

7:end for

As revealed in the test cases below, this technique can drastically reduce the
number of iterations required by the L-scheme.

The original formulation presented in [7] allows a for more general step,

$x_{k}:=\beta_{k}\sum_{i=0}^{m_{k}}\alpha_{i}g(x_{k-m_{k}+i-1})+(1-\beta_{k})\sum_{i=0}^{m_{k}}\alpha_{i}x_{k-m_{k}+i-1},$

for a user-defined tuning parameter $\beta_{k}\in(0,1]$. We considered the
simplified formulation, obtained with $\beta_{k}=1$, because no improvements
have been observed in the numerical results when using the extended one.

We remark that large values for the depth $m$ can result in an instability of
the solution algorithm. When implementing the Anderson acceleration, one has
to tune this parameter properly. A small $m$ could produce only a small
reduction in the numbers of iterations; too large $m$ could result in a non-
converging algorithm [29].

######  Remark 4.

The definition of the nonlinear splitting solvers allows for different ways to
apply the AA. We study three different loops: the coupling one and the
linearizing ones, one for each set of equations. We apply the Anderson
acceleration to each of them. Two different parameters, $m$ and $m_{lin}$, are
defined. The former is used for the AA on the coupling loop, the latter for
the implementation on the linearization ones. The same $m_{lin}$ will be used
for the loop regarding the flow equations and for the one regarding the
transport.

##  3 Numerical examples

In the following, we consider four numerical examples with increasing
complexity, based on a manufactured solution, and an example in which the
boundary conditions drive the flow but no manufactured solution is given. The
first four will differ in the different values for $\gamma$, $\delta$ and
$\tau(\theta)$ taken into account. We have implemented the models and solving
schemes in MRST, a toolbox based on Matlab for the simulations of flow in
porous media [44]. We use the two point flux approximation, one of the most
common spatial discretization techniques. We remark that the linearization
schemes and solving algorithms do not depend on the particular choice of the
spatial discretization, so one may apply these solvers to other methods as
well, without any difficulty.

The domain is the unit square $\Omega$ and the final time taken into
consideration is $T=3$. The simulations are performed on regular meshes,
consisting of squares with sides $dx=1/10,\ 1/20,$ and $1/40$. The time steps
are $\Delta t=T/25,\ T/50$ and $T/100$. The $L$ parameters, used in the
L-scheme formulations, are $L_{1}=L_{2}=L_{3}=0.1$, if not specified
otherwise. We took into consideration different values, but the aforementioned
choice seems to produce a robust algorithm which required fewest iterations to
achieve the convergence. For the ease of the presentation, we set the three
$L$ parameters equal to each other; one could define different values for each
parameter, investigating even further the linearization of each equation. We
avoided this due to the application of the AA. We will observe that the
schemes can be drastically accelerated, even though the $L$ parameters are not
optimal.

The condition numbers, for the stiffness matrices resulting from the different
solving algorithms are computed using the $L^{1}$ norm, and we here report the
averaged values over the full simulation. A minus sign $(-)$, in the tables
reporting on iterations and condition numbers, implies that the method failed
to converge for the particular combination of the time step and mesh size. The
tolerance $\epsilon$ used in the stopping criterion presented in Remark 3 is
$\epsilon=1e-6$. We always report the total numbers of iterations required by
the full simulation, not the average number required by each time step. For
the splitting solver, we present, separately, the condition numbers of both
flow and transport equations. Furthermore, for the nonlinear splitting, the
iterations are divided in two, the ones required by the flow equations and the
ones for the transport. Finally, the condition numbers reported are obtained
by averaging over the full simulations.

We apply the Anderson acceleration to each solving algorithm, always reporting
the depths $m$ and $m_{lin}$ used. Once more, $m_{lin}$ is the Anderson
parameter used for the acceleration of the linearization loops regarding the
flow and transport equations in the nonlinear splitting solvers.

Inspired by [46], the first four examples are constructed in such a way that
the following is an exact solution:

$\theta_{m}(x,y,t)=\begin{cases}1-\frac{1}{2}\cos((t_{1}(x,y)-t)^{2})&\mbox{if}\qquad
t<t_{1}(x,y),\\\ \frac{1}{2}&\mbox{if}\qquad t_{1}(x,y)\leq t\leq
t_{2}(x,y),\\\ 1-\frac{1}{2}\cos((t-t_{2}(x,y))^{2})\qquad&\mbox{if}\qquad
t>t_{2}(x,y),\end{cases}$ (31)
$\Psi_{m}(x,y,t)=\begin{cases}-p_{cap}(\theta_{m})+\tau(\theta_{m})\partial_{t}\theta_{m}-\gamma&\mbox{if}\qquad\partial_{t}\theta_{m}<-\delta,\\\
-p_{cap}(\theta_{m})+\tau(\theta_{m})\partial_{t}\theta_{m}+\frac{\gamma}{\delta}\partial_{t}\theta_{m}&\mbox{if}\qquad-\delta\leq\partial_{t}\theta_{m}\leq\delta,\\\
-p_{cap}(\theta_{m})+\tau(\theta_{m})\partial_{t}\theta_{m}+\gamma&\mbox{if}\qquad\partial_{t}\theta_{m}>\delta,\end{cases}$
(32) $c_{m}(x,y,t)=x(x-1)y(y-1)t,$ (33)

where $t_{1}(x,y)=xy$, $t_{2}(x,y)=xy+2$. Once the manufactured water content
is defined, one obtains the pressure by simply using the second equation in
(5). The capillary pressure is expressed as
$p_{cap}(\theta,c)=1-\theta^{2}-0.1\ c^{3}$ and the conductivity as
$K(\theta)=1+\theta^{2}$. Even though such a formulation may appear non-
realistic, we are mainly interested in the nonlinearities. Furthermore, a
nonlinear reaction term, $R(c)=c/(c+1)$, is taken into account in the
transport equation and the diffusion/dispersion coefficient $D$ is set equal
to $1$.

Given the analytical expressions above, we can easily define the initial
conditions, the Dirichlet boundary conditions on the unit square and compute
the source terms $\mathbb{S}_{1}$ and $\mathbb{S}_{2}$ such that $\Psi_{m}$,
$\theta_{m}$ and $c_{m}$ are solutions of the system. In particular, the
initial concentration and water content are

$\displaystyle c(x,y,0)$ $\displaystyle=0\qquad$
$\displaystyle\mbox{on}\qquad\Omega,$ $\displaystyle\theta(x,y,0)$
$\displaystyle=1-\frac{1}{2}\cos\big{(}t_{1}(x,y)^{2}\big{)}\qquad$
$\displaystyle\mbox{on}\qquad\Omega.$

We impose a zero concentration $c(x,y,t)=0$ on the boundary of the domain. The
remaining boundary conditions, concerning the pressure, are time dependent.
One needs to compute $t_{1}$ and $t_{2}$ on every side of the unit square.
Once the time intervals given by $t_{1}$ and $t_{2}$ are obtained, the
pressure can easily be imposed. For example, on the left side $x=0$, thus
$t_{1}=0$ and $t_{2}=2$. The water content $\theta$ becomes

$\theta(0,y,t)=\theta_{left}(t)=\begin{cases}\frac{1}{2}&\mbox{if}\qquad
0<t<2,\\\ 1-\frac{1}{2}\cos((t-2)^{2})\qquad&\mbox{if}\qquad t\geq
2,\end{cases}$

and the resulting pressure boundary condition is

$\Psi_{left}(0,y,t)=\begin{cases}-p_{cap}(\theta_{left})+\tau(\theta_{left})\partial_{t}\theta_{left}-\gamma&\mbox{if}\qquad\partial_{t}\theta_{left}<-\delta,\\\
-p_{cap}(\theta_{left})+\tau(\theta_{left})\partial_{t}\theta_{left}+\frac{\gamma}{\delta}\partial_{t}\theta_{left}&\mbox{if}\qquad-\delta\leq\partial_{t}\theta_{left}\leq\delta,\\\
-p_{cap}(\theta_{left})+\tau(\theta_{left})\partial_{t}\theta_{left}+\gamma&\mbox{if}\qquad\partial_{t}\theta_{left}>\delta.\end{cases}$

Analogously, one can compute the pressure boundary conditions on the remaining
sides.

###  3.1 Example 1, $\gamma=0$, $\tau(\theta)=1$

In the first example we impose $\gamma=0$, thus, the hysteresis effects are
neglected but we include a dynamic effect by considering a constant
$\tau(\theta)=1$. We compare the different algorithms presented in Section
2.1, reporting in the Tables 1 and 2 the total numbers of iterations required
by each algorithm, and the condition numbers of the systems associated with
each scheme. In the former, we investigate a fixed time step size, $\Delta
t=T/25$, in the latter a fixed mesh size, $dx=1/10$. As expected, a finer mesh
results in worse conditioned systems, while smaller time steps give better
conditioned ones. Moreover, the total number of iterations is increasing as we
reduce the time step; smaller $\Delta t$ implies more time steps and thus more
iterations.

The schemes based on the L-scheme appear to be better conditioned than those
based on the Newton method. The result is coherent with the theory [36, 45,
49, 57, 62]. One could even improve the condition numbers by using larger $L$
parameters. However, larger values would have also increased the total numbers
of iterations.

The alternate splitting schemes are converging much faster than the nonlinear
ones. It is also interesting to observe that the numbers of iterations,
required by the alternate splitting schemes, are comparable with the ones
associated with the monolithic solvers. In [36], we observed similar results
when solving the models without hysteresis and dynamic effects.

We notice also some reduction in the number of iterations required by the
L-schemes thanks to the Anderson acceleration. The results obtained for the
non-accelerated L-schemes ($m=0$) are already optimal in terms of numbers of
iterations; thus, the improvement can only be minimal. We report the total
number of iterations for the full simulation, but mention that on average, for
each time step, the L-scheme requires only five or six iterations. This is
already a remarkable result, achieved thanks to the optimal L parameters.
Furthermore, the Newton solvers have resulted in being slower when combined
with the AA. This is coherent with the theory where it has been observed that
quadratically convergent schemes, cannot be improved and the resulting
accelerated solvers appear slower [28].

| Monolithic |  |  | NonLinS |  |  |  | AltLinS |  |
---|---|---|---|---|---|---|---|---|---|---
| Newton |  |  | Newton |  |  |  | Newton |  |
dx | # iterations | condition # | # iterations | cond. # Flow |  | cond. # Transport | # iterations | cond. # Flow |  | cond. # Transport
1/10 | 65 | 4.91e+02 | 56 - 50 | 1.82e+02 |  | 1.60e+02 | 66 | 1.83e+02 |  | 1.60e+02
1/20 | 69 | 2.22e+03 | 57 - 50 | 7.61e+02 |  | 6.33e+02 | 66 | 7.60e+02 |  | 6.3359e+02
1/40 | 70 | 1.11e+04 | 58 - 50 | 3.29e+03 |  | 2.52e+03 | 66 | 3.42e+03 |  | 2.5165e+03
| Newton | (AA m = 1) |  | Newton (AA $m=m_{lin}$ = 1) |  |  |  | Newton (AA m = 1) |  |
dx | # iterations | condition # | # iterations | cond. # Flow |  | cond. # Transport | # iterations | cond. # Flow |  | cond. # Transport
1/10 | 93 | 4.76e+02 | 65 - 50 | 1.82e+02 |  | 1.60e+02 | 74 | 1.82e+02 |  | 1.60e+02
1/20 | 98 | 2.10e+03 | 66 - 50 | 7.54e+02 |  | 6.33e+02 | 74 | 7.59e+02 |  | 6.34e+02
1/40 | 100 | 9.92e+03 | 69 - 50 | 3.28e+03 |  | 2.52e+03 | 76 | 3.36e+03 |  | 2.52e+03
| L-scheme |  |  | L-scheme |  |  |  | L-scheme |  |
dx | # iterations | condition # | # iterations | cond. # Flow |  | cond. # Transport | # iterations | cond. # Flow |  | cond. # Transport
1/10 | 134 | 4.15e+02 | 117 - 116 | 1.43e+02 |  | 1.36e+02 | 140 | 1.62e+02 |  | 1.3717e+02
1/20 | 140 | 1.82e+03 | 119 - 115 | 6.95e+02 |  | 5.41e+02 | 144 | 7.37e+02 |  | 5.41e+02
1/40 | 146 | 8.52e+03 | 128 - 116 | 3.12e+03 |  | 2.15e+03 | 150 | 3.18e+03 |  | 2.15 e03
| L-scheme | (AA m = 1) |  | L-scheme (AA $m=m_{lin}$ = 1) |  |  |  | L-scheme (AA m= 1) |  |
dx | # iterations | condition # | # iterations | cond. # Flow |  | cond. # Transport | # iterations | cond. # Flow |  | cond. # Transport
1/10 | 127 | 4.14e+02 | 107 - 100 | 1.9690e+02 |  | 1.3681e+02 | 136 | 1.44e+02 |  | 1.33e+02
1/20 | 129 | 1.85e+03 | 111 - 100 | 7.9408e+02 |  | 5.4122e+02 | 142 | 7.13e+02 |  | 5.31e+02
1/40 | 130 | 8.80e+03 | 118 - 100 | 3.3661e+03 |  | 2.1522e+03 | 146 | 3.01e+03 |  | 2.02e+03

Table 1: Example 1: Total number of iterations and condition numbers for fixed
$\Delta t=T/25$, and different $dx$. Here, $L_{1}=L_{2}=L_{3}=0.1$ and
$m=m_{lin}=1$.

| Monolithic |  |  | NonLinS |  |  |  | AltLinS |  |
---|---|---|---|---|---|---|---|---|---|---
| Newton |  |  | Newton |  |  |  | Newton |  |
$\Delta t$ | # iterations | condition # | # iterations | cond. # Flow |  | cond. # Transport | # iterations | cond. # Flow |  | cond. # Transport
T/25 | 65 | 4.91e+02 | 56 - 50 | 1.82e+02 |  | 1.60e+02 | 66 | 1.83e+02 |  | 1.60e+02
T/50 | 103 | 2.75e+02 | 98 - 50 | 1.23e+02 |  | 8.86e+01 | 99 | 1.24e+02 |  | 8.89e+01
T/100 | 186 | 1.97e+02 | 172 - 50 | 8.69e+01 |  | 4.73e+01 | 172 | 8.74e+01 |  | 4.74e+01
| Newton | (AA m = 1) |  | Newton (AA $m=m_{lin}$ = 1) |  |  |  | Newton (AA m = 1) |  |
$\Delta t$ | # iterations | condition # | # iterations | cond. # Flow |  | cond. # Transport | # iterations | cond. # Flow |  | cond. # Transport
T/25 | 93 | 4.76e+02 | 65 - 50 | 1.82e+02 |  | 1.60e+02 | 74 | 1.82e+02 |  | 1.60e+02
T/50 | 137 | 2.75e+02 | 114 - 50 | 1.23e+02 |  | 8.86e+01 | 114 | 1.23e+02 |  | 8.89e+01
T/100 | 244 | 1.96e+02 | 201 - 100 | 8.59e+01 |  | 4.73e+01 | 200 | 8.69e+01 |  | 4.74e+01
| L-scheme |  |  | L-scheme |  |  |  | L-scheme |  |
$\Delta t$ | # iterations | condition # | # iterations | cond. # Flow |  | cond. # Transport | # iterations | cond. # Flow |  | cond. # Transport
T/25 | 134 | 4.16e+02 | 117 - 116 | 1.43e+02 |  | 1.36e+02 | 140 | 1.62e+02 |  | 1.3717e+02
T/50 | 219 | 2.39e+02 | 182 - 200 | 1.12e+02 |  | 7.51e+01 | 218 | 1.2460e+02 |  | 7.5203e+01
T/100 | 425 | 1.64e+02 | 346 - 400 | 8.35e+01 |  | 3.98e+01 | 438 | 8.2955e+01 |  | 3.9923e+01
| L-scheme | (AA m = 1) |  | L-scheme (AA $m=m_{lin}$ = 1) |  |  |  | L-scheme (AA m = 1) |  |
$\Delta t$ | # iterations | condition # | # iterations | cond. # Flow |  | cond. # Transport | # iterations | cond. # Flow |  | cond. # Transport
T/25 | 127 | 4.14e+02 | 107 - 100 | 1.97e+02 |  | 1.37e+02 | 136 | 1.44e+02 |  | 1.33e+02
T/50 | 217 | 2.36e+02 | 192 - 160 | 1.25e+02 |  | 7.51e+01 | 238 | 1.14e+02 |  | 7.50e+01
T/100 | 387 | 1.66e+02 | 347 - 300 | 8.33e+01 |  | 3.98e+01 | 432 | 8.21e+01 |  | 3.91e+01

Table 2: Example 1: Total number of iterations and condition numbers for fixed
$dx=1/10$, and different $\Delta t$. Here $L_{1}=L_{2}=L_{3}=0.1$,
$m=m_{lin}=1$.

In Table 3 we present the numerical errors and the estimated order of
convergence of the spatial discretization based on the successively refined
meshes investigated. Given the manufactured solution $\Psi_{m}$, we compute
the numerical error $e_{\Psi}=\left\lVert\Psi_{m}-\Psi_{n}\right\rVert$, where
$\Psi_{n}$ is the numerical pressure computed. Similarly, we can define
$e_{\theta}$ and $e_{c}$. Furthermore, $e_{\Psi,1}$ is the numerical error
obtained for the mesh size $dx=1/10$ and $\Delta t=T/25$, $e_{\Psi,2}$ for
$dx=1/20$ and $\Delta t=T/50$, and finally $e_{\Psi,3}$ for $dx=1/40$ and
$\Delta t=T/100$.
$EOC=\log\left(\frac{e_{\Psi,i}}{e_{\Psi,i+1}}\right)/\log(2)$ is the
estimated order of convergence. These results are independent from the solving
algorithm taken into account, only the discretization approach plays a role.
In this case we use a TPFA which is known to have a order of convergence equal
to 1, as also reported here in the table.

| $e_{1}$ | EOC | $e_{2}$ | EOC | $e_{3}$
---|---|---|---|---|---
$\Psi$ | 4.61e-02 | 0.98 | 2.33e-02 | 1.00 | 1.16e-02
$\theta$ | 1.31e-02 | 0.97 | 6.71e-03 | 0.99 | 3.38e-03
$c$ | 6.24e-03 | 1.53 | 2.15e-03 | 1.32 | 8.60e-04
Table 3: Example 1: Numerical error and estimated order of convergence (EOC)
of the discretization method.

In Table 4 we tested different values of $m$ and L. We can observe that, for
large L, the Monolithic L-Scheme, here investigated, requires more iterations
than for smaller parameters. If many iterations are required to achieve the
convergence at each time step, one can take in consideration larger $m$
values. For $L=0.1$, the optimal choice, in terms of numbers of iterations, is
$m=1$; larger values result in slightly slower schemes. For the largest $L$
tested, $L=2.3$, the optimal choice is $m=2$. Such L value corresponds to the
theoretical L, $\max\limits_{\theta}\\{\frac{\partial
p_{cap}}{\partial\theta}\\}\approx 2.3$. To ensure the monotone convergence of
the scheme it has been proved, that the parameter chosen must be larger than
the Lipschitz constant of the nonlinearity, in this case the capillary
pressure as a function of $\theta$ (see [45, 57, 62]). We set $L_{1}$, $L_{2}$
and $L_{3}$ equal to the theoretical L computed for the capillary pressure.

We can conclude that it is possible to obtain significant improvements by
investigating the AA and thus finding the appropriate depth $m$. In this work,
we have focused more on individuating the optimal L parameters, as refining
the AA can be done more easily. The depths used are, in fact, small natural
numbers.

L | $m=0$ | $m=1$ | $m=2$ | $m=3$ | $m=5$
---|---|---|---|---|---
.1 | 146 | 130 | 146 | 179 | 215
.5 | 247 | 161 | 152 | 172 | 200
1 | 411 | 202 | 168 | 165 | 195
2 | 711 | 290 | 199 | 206 | 217
2.3 | 810 | 317 | 207 | 215 | 227
Table 4: Example 1: Comparison of number of iterations for different $m$ and L
parameters for L-Mono.
Here, $dx=1/40$, $\Delta t=T/25$ and $L_{1}=L_{2}=L_{3}=L$.

Finally, we investigate the order of convergence of the linearization schemes.
In Figure 5, we plot the residuals of pressure, water content and
concentration, obtained at the final time step for the finest mesh size,
$dx=1/40$ and $\Delta t=T/25$. We can deduce the rates of convergence of the
different linearization schemes. The L-schemes appear to be, as expected,
linearly convergent in term of numbers of iterations. The AA improves the
results only slightly, as already observed in Tables 1 and 2. This is
justified by the fact that the L parameters ($L_{1},L_{2},L_{3}$) chosen here,
appear to be optimal. The Newton methods are instead quadratically convergent.
More precise results are observable in Table 5. Here we present the exact
order of convergence of the different schemes. Given the residual of each
unknown ($res_{\Psi},res_{\theta},res_{c}$), at each iteration $j$, we compute
the order of convergence as follow:
$ORD_{j}=\big{(}\log(res_{j+1}/res_{j})\big{)}/\big{(}\log(res_{j}/res_{j-1})\big{)}$.
For a fixed time step, we can average the orders obtained over the number of
iterations required to achieve the convergence. We report below the values
obtained by investigating the final time step; similar results have been
observed for previous time steps.

(a) Pressure residuals.

(b) Water content residuals.

(c) Concentration residuals.

Figure 5: Example 1: Residuals of each unknown at the final time step, for the
different schemes.
Here, $L_{1}=L_{2}=L_{3}=0.1$, $m=m_{lin}=1$, $dx=1/40$, $\Delta t=T/25$. | LS-Mono | LS-Mono And. | LS-NonLinS | LS-NonLinS And. | LS-AltLinS | LS-AltLinS And.
---|---|---|---|---|---|---
$\Psi$ | 1.00 | 1.26 | 1.13 | 1.46 | 1.13 | 2.07
$c$ | 0.96 | 1.55 | 1.01 | 1.77 | 1.01 | 1.30
$\theta$ | 1.01 | 1.40 | 0.89 | 1.01 | 0.94 | 1.32
| New.-Mono | New.-Mono And. | New.-NonLinS | New.-NonLinS And. | New.-AltLinS | New.-AltLinS And.
$\Psi$ | 2.03 | 1.61 | 1.94 | 1.86 | 2.04 | 2.03
$c$ | 2.41 | 2.01 | 1.61 | 1.41 | 1.85 | 1.46
$\theta$ | 1.97 | 0.57 | 1.81 | 0.97 | 1.61 | 0.95
Table 5: Example 1: Estimated order of convergence for the different
linearization schemes.

######  Remark 5.

Since the Anderson acceleration with small depth is a cheap post-processing
step, reducing the number of iterations directly reduces the CPU time almost
proportionally.

###  3.2 Example 2, $\gamma=0$, $\tau(\theta)=1+\theta^{2}$

In the second example, the setup of the first is extended by adopting a
nonlinear $\tau$, precisely $\tau(\theta)=1+\theta^{2}$. In Tables 6 and 7, we
present the condition numbers and the required iteration counts associated
with each solving algorithm.

The introduction of a nonlinear $\tau(\theta)$ increases the numbers of
iterations required by each solver. The L-scheme is linearly convergent while
the Newton method is quadratically convergent. Furthermore, the conclusions
from Example 1 concerning the AA remain the same. In particular, we observe
only small reductions in the numbers of iterations required by the L-schemes.
Once more, this is justified by the optimal choice of the L parameters. We can
observe that, for each time step, only a few iterations are required; thus, no
further acceleration is expected.

| Monolithic |  |  | NonLinS |  |  |  | AltLinS |  |
---|---|---|---|---|---|---|---|---|---|---
| Newton |  |  | Newton |  |  |  | Newton |  |
dx | # iterations | condition # | # iterations | cond. # Flow |  | cond. # Transport | # iterations | cond. # Flow |  | cond. # Transport
1/10 | 67 | 481.17 | 59 - 50 | 161.39 |  | 161.00 | 64 | 152.62 |  | 161.76
1/20 | 69 | 2.20e+03 | 57 - 50 | 656.26 |  | 636.33 | 66 | 644.98 |  | 644.99
1/40 | 70 | 1.10e+04 | 59 - 50 | 2.85e+03 |  | 2.53e+03 | 68 | 2.93e+03 |  | 2.53e+03
| Newton | (Anderson m = 1) |  | Newton (AA $m=m_{lin}$ = 1) |  |  |  | Newton (Anderson m = 1) |  |
dx | # iterations | condition # | # iterations | cond. # Flow |  | cond. # Transport | # iterations | cond. # Flow |  | cond. # Transport
1/10 | 94 | 467.04 | 69 - 50 | 160.36 |  | 161.00 | 76 | 152.45 |  | 161.94
1/20 | 98 | 2.07e+03 | 68 - 50 | 636.14 |  | 636.33 | 76 | 642.25 |  | 638.40
1/40 | 100 | 9.85e+03 | 71 - 50 | 2.80e+03 |  | 2.53e+03 | 80 | 2.87e+03 |  | 2.53e+03
| L-scheme |  |  | L-scheme |  |  |  | L-scheme |  |
dx | # iterations | condition # | # iterations | cond. # Flow |  | cond. # Transport | # iterations | cond. # Flow |  | cond. # Transport
1/10 | 136 | 403.49 | 112 - 117 | 158.33 |  | 137.81 | 130 | 158.05 |  | 138.50
1/20 | 139 | 1.77e+03 | 120 - 116 | 644.85 |  | 544.34 | 136 | 651.07 |  | 545.85
1/40 | 144 | 8.34e+03 | 125 - 116 | 2.73e+03 |  | 2.16e+03 | 142 | 2.81e+03 |  | 2.16e+03
| L-scheme | (Anderson m = 1) |  | L-scheme (AA $m=m_{lin}$ = 1) |  |  |  | L-scheme (AA m = 1) |  |
dx | # iterations | condition # | # iterations | cond. # Flow |  | cond. # Transport | # iterations | cond. # Flow |  | cond. # Transport
1/10 | 129 | 401.61 | 106 - 100 | 158.87 |  | 137.74 | 132 | 156.89 |  | 138.21
1/20 | 131 | 1.78e+03 | 113 - 100 | 644.63 |  | 543.86 | 134 | 647.35 |  | 545.18
1/40 | 137 | 8.44e+03 | 117 - 100 | 2.76e+03 |  | 2.16e+03 | 142 | 2.78e+03 |  | 2.16e+03

Table 6: Example 2: Total number of iterations and condition numbers for fixed
$\Delta t=T/25$, and different $dx$. Here, $L_{1}=L_{2}=L_{3}=0.1$ and
$m=m_{lin}=1$.

| Monolithic |  |  | NonLinS |  |  |  | AltLinS |  |
---|---|---|---|---|---|---|---|---|---|---
| Newton |  |  | Newton |  |  |  | Newton |  |
$\Delta t$ | # iterations | condition # | # iterations | cond. # Flow |  | cond. # Transport | # iterations | cond. # Flow |  | cond. # Transport
T/25 | 67 | 481.17 | 59 - 50 | 161.39 |  | 161.00 | 64 | 152.62 |  | 161.76
T/50 | 107 | 256.90 | 102 - 50 | 95.40 |  | 89.26 | 102 | 95.08 |  | 89.80
T/100 | 193 | 159.33 | 184 - 100 | 73.42 |  | 47.70 | 184 | 72.37 |  | 47.91
| Newton | (AA m = 1) |  | Newton (AA $m=m_{lin}$ = 1) |  |  |  | Newton (AA m = 1) |  |
$\Delta t$ | # iterations | condition # | # iterations | cond. # Flow |  | cond. # Transport | # iterations | cond. # Flow |  | cond. # Transport
T/25 | 94 | 467.04 | 69 - 50 | 160.36 |  | 161.00 | 76 | 152.45 |  | 161.94
T/50 | 142 | 254.28 | 117 - 50 | 94.64 |  | 89.26 | 118 | 95.71 |  | 89.88
T/100 | 259 | 159.96 | 207 - 100 | 73.05 |  | 47.70 | 208 | 72.78 |  | 47.95
| L-scheme |  |  | L-scheme |  |  |  | L-scheme |  |
$\Delta t$ | # iterations | condition # | # iterations | cond. # Flow |  | cond. # Transport | # iterations | cond. # Flow |  | cond. # Transport
T/25 | 136 | 403.49 | 112 - 117 | 157.86 |  | 137.81 | 130 | 158.05 |  | 138.50
T/50 | 220 | 211.60 | 185 - 201 | 98.63 |  | 75.62 | 220 | 97.80 |  | 75.86
T/100 | 433 | 125.49 | 356 - 400 | 65.18 |  | 40.16 | 446 | 65.24 |  | 40.28
| L-scheme | (AA m = 1) |  | L-scheme (AA $m=m_{lin}$ = 1) |  |  |  | L-scheme (AA m = 1) |  |
$\Delta t$ | # iterations | condition # | # iterations | cond. # Flow |  | cond. # Transport | # iterations | cond. # Flow |  | cond. # Transport
T/25 | 129 | 401.61 | 106 - 100 | 158.87 |  | 137.74 | 132 | 156.89 |  | 138.21
T/50 | 214 | 212.13 | 192 - 161 | 97.72 |  | 75.69 | 240 | 97.63 |  | 75.66
T/100 | 388 | 126.15 | 355 - 302 | 64.73 |  | 40.16 | 440 | 64.45 |  | 40.17

Table 7: Example 2: Total number of iterations and condition numbers for fixed
$dx=1/10$, and different $\Delta t$. Here $L_{1}=L_{2}=L_{3}=0.1$,
$m=m_{lin}=1$.

As for the results presented in Tables 6 and 7, the numerical errors and EOC
reported in Table 8 are similar to the ones from the first example.

| $e_{1}$ | EOC | $e_{2}$ | EOC | $e_{3}$
---|---|---|---|---|---
$\Psi$ | 0.0759 | 0.9688 | 0.0388 | 0.9913 | 0.0195
$\theta$ | 0.0140 | 0.9556 | 0.0072 | 0.9830 | 0.0036
$c$ | 0.0084 | 1.5198 | 0.0029 | 1.3037 | 0.0012
Table 8: Example 2: Numerical error and estimated order of convergence (EOC)
of the discretization method.

In Figure 6, we report the residuals of the pressure, water content and
concentration at the final time step. The L-schemes are linearly convergent,
and applying the AA does not result in significant improvements. The
convergence rates and number of iterations remain the same. Also for the
Newton solvers, since they are quadratically convergent, the AA cannot improve
this aspect. Table 9 presents the precise order of convergence of the
different linearization schemes.

(a) Pressure residuals.

(b) Water content residuals.

(c) Concentration residuals.

Figure 6: Example 2: Residuals of each unknown at the final time step, for the different schemes. Here, $L_{1}=L_{2}=L_{3}=0.1$, $m=m_{lin}=1$, $dx=1/40$, $\Delta t=T/25$.  | LS-Mono | LS-Mono And. | LS-NonLinS | LS-NonLinS And. | LS-AltLinS | LS-AltLinS And.
---|---|---|---|---|---|---
$\Psi$ | 1.07 | 1.40 | 1.11 | 1.36 | 1.14 | 1.24
$c$ | 0.99 | 1.23 | 0.98 | 1.45 | 0.96 | 1.10
$\theta$ | 1.03 | 1.15 | 0.97 | 1.25 | 0.93 | 0.98
| New.-Mono | New.-Mono And. | New.-NonLinS | New.-NonLinS And. | New.-AltLinS | New.-AltLinS And.
$\Psi$ | 1.61 | 1.58 | 1.97 | 1.69 | 2.15 | 2.14
$c$ | 2.68 | 1.47 | 1.98 | 1.54 | 1.98 | 1.83
$\theta$ | 1.99 | 1.61 | 1.89 | 1.70 | 2.17 | 1.95
Table 9: Example 2: Estimated order of convergence for the different
linearization schemes.

###  3.3 Example 3, $\gamma=1$, $\delta=5e-3$, $\tau(\theta)=0$

With the same manufactured solutions, we now consider the case without dynamic
effects ($\tau(\theta)=0$), but include hysteresis by choosing $\gamma=1$ and
$\delta=5e-3$.

From the results in Tables 10 and 11, we notice that the Newton method, in all
its formulations, fails to converge. In Table 11, smaller time steps are
taken, but no improvements are observable. A further reduction of the time
step could have resulted in converging Newton solvers but the total numbers of
iterations for the full simulation would have been larger than the ones
required by the L-schemes on fewer but larger time steps. In contrast, the
L-schemes are more robust and, even though requiring a higher number of
iterations than previously, they converge. We take $L_{1}=L_{2}=L_{3}=L=1$,
which appear to be the optimal choice in terms of numbers of iterations.

The AA improves the convergence of the L-schemes. This is the first example of
this study in which the results obtained thanks to the AA are improved
substantially. This is due to the presence of the hysteresis, requiring a
large $L$ for the overall convergence, and thus the total numbers of
iterations is larger. On average, the monolithic L-scheme solver requires
circa 18 iterations per time step. For $m=1$, the AA reduces the iterations by
circa $50\%$. Different $m$ values have been tested but none of the ones
investigated lead to the convergence of the Newton schemes. On all tests,
Newton has failed to converge, whereas the L-schemes converged and the AA
yields further improvement.

| Monolithic |  |  | NonLinS |  |  |  | AltLinS |  |
---|---|---|---|---|---|---|---|---|---|---
| Newton |  |  | Newton |  |  |  | Newton |  |
dx | # iterations | condition # | # iterations | cond. # Flow |  | cond. # Transport | # iterations | cond. # Flow |  | cond. # Transport
1/10 | - | - | - | - |  | - | - | - |  | -
1/20 | - | - | - | - |  | - | - | - |  | -
1/40 | - | - | - | - |  | - | - | - |  | -
| L-scheme |  |  | L-scheme |  |  |  | L-scheme |  |
dx | # iterations | condition # | # iterations | cond. # Flow |  | cond. # Transport | # iterations | cond. # Flow |  | cond. # Transport
1/10 | 448 | 409.16 | 210 - 441 | 484.83 |  | 69.22 | 450 | 361.76 |  | 69.34
1/20 | 456 | 1.65e+03 | 266 - 439 | 1.98e+03 |  | 259.35 | 452 | 1.46e+03 |  | 260.02
1/40 | 468 | 6.62e+03 | 276 - 438 | 7.74e+03 |  | 996.36 | 460 | 5.88e+03 |  | 999.13
| L-scheme | (AA m = 2) |  | L-scheme (AA m = 2, $m_{lin}$ = 5) |  |  |  | L-scheme (AA m = 1) |  |
dx | # iterations | condition # | # iterations | cond. # Flow |  | cond. # Transport | # iterations | cond. # Flow |  | cond. # Transport
1/10 | 226 | 468.28 | 179 - 150 | 497.56 |  | 70.10 | 328 | 450.28 |  | 71.81
1/20 | 278 | 1.97e+03 | 187 - 141 | 2.03e+03 |  | 261.40 | 408 | 2.09e+03 |  | 269.72
1/40 | 303 | 8.24e+03 | - | - |  | - | 378 | 7.76e+03 |  | 967.29

Table 10: Example 3: Total number of iterations and condition numbers for
fixed $\Delta t=T/25$, and different $dx$. Here, $L_{1}=L_{2}=L_{3}=1$,
different $m$ and $m_{lin}$ are taken into account.

| Monolithic |  |  | NonLinS |  |  |  | AltLinS |  |
---|---|---|---|---|---|---|---|---|---|---
| Newton |  |  | Newton |  |  |  | Newton |  |
dx | # iterations | condition # | # iterations | cond. # Flow |  | cond. # Transport | # iterations | cond. # Flow |  | cond. # Transport
T/25 | - | - | - | - |  | - | - | - |  | -
T/50 | - | - | - | - |  | - | - | - |  | -
T/100 | - | - | - | - |  | - | - | - |  | -
| L-scheme |  |  | L-scheme |  |  |  | L-scheme |  |
|  |  |  |  | cond. # |  |  |  | cond. # |
$\Delta t$ | # iterations | condition # | # iterations | Flow |  | Transport | # iterations | Flow |  | Transport
T/25 | 448 | 409.16 | 210 - 441 | 484.83 |  | 69.22 | 450 | 361.76 |  | 69.34
T/50 | 836 | 363.15 | 513 - 846 | 422.16 |  | 35.96 | 838 | 332.02 |  | 36.03
T/100 | 1787 | 395.55 | 1261 - 1597 | 428.24 |  | 18.91 | 1764 | 363.40 |  | 19.07
| L-scheme | (AA m = 1) |  | L-scheme (AA m = 2, $m_{lin}$ = 5) |  |  |  | L-scheme (AA m = 1) |  |
dx | # iterations | condition # | # iterations | cond. # Flow |  | cond. # Transport | # iterations | cond. # Flow |  | cond. # Transport
T/25 | 226 | 468.28 | 179 - 150 | 497.56 |  | 70.10 | 328 | 450.28 |  | 71.81
T/50 | 533 | 410.93 | 346 - 316 | 504.00 |  | 36.85 | 664 | 442.54 |  | 37.43
T/100 | 1217 | 467.07 | 861 - 944 ($m_{lin}$ = 1) | 491.24 |  | 19.15 | 1842 | 473.32 |  | 19.75

Table 11: Example 3: Total number of iterations and condition numbers for
fixed $dx=1/10$, and different $\Delta t$. Here, $L_{1}=L_{2}=L_{3}=1$,
different $m$ and $m_{lin}$ are taken into account

Once more, we report the numerical errors and the estimated orders of
convergence associated with the discretization technique here implemented
(TPFA). In Table 12, we present the values obtained for the monolithic
L-scheme. The EOC depends only on the discretization technique, not the
linearization scheme or solving algorithm.

| $e_{1}$ | EOC | $e_{2}$ | EOC | $e_{3}$
---|---|---|---|---|---
$\Psi$ | 0.0308 | 1.0345 | 0.0150 | 0.9281 | 0.0085
$\theta$ | 0.0350 | 1.2985 | 0.0142 | 1.0997 | 0.0066
$c$ | 0.0060 | 1.3202 | 0.0024 | 1.3000 | 0.0010
Table 12: Example 3: Numerical error and estimated order of convergence (EOC)
of the discretization method.

In Figure 7, we report the residuals of pressure, water content and
concentration, at the final time step. The differences between the accelerated
and non-accelerated schemes seem to be minimal at the final time step but we
observe in Tables 10 and 11 that the total improvements are actually
substantial. The precise orders of convergence for the different solving
algorithms are reported in Table 13. The non-accelerated L-schemes have an
order of convergence equal to one, while the accelerated ones have slightly
larger values. No result was reported for the Newton schemes due to the lack
of convergence.

(a) Pressure residuals.

(b) Water content residuals.

(c) Concentration residuals.

Figure 7: Example 3: residuals of each unknown at the final time step, for the
different schemes.
Here, $L_{1}=L_{2}=L_{3}=1$, $dx=1/40$, $\Delta t=T/25$ and $m\neq m_{lin}$. | LS-Mono | LS-Mono And. | LS-NonLinS | LS-NonLinS And. | LS-AltLinS | LS-AltLinS And.
---|---|---|---|---|---|---
$\Psi$ | 1.14 | 1.33 | 1.00 | 1.19 | 1.00 | 1.81
$c$ | 1.00 | 1.32 | 1.00 | 1.51 | 0.99 | 1.52
$\theta$ | 1.00 | 1.29 | 1.00 | 1.39 | 1.00 | 1.16
Table 13: Example 3: Estimated order of convergence for the different
linearization schemes.

###  3.4 Example 4, $\gamma=1$, $\delta=5e-3$, $\tau(\theta)=1+\theta^{2}$

Finally, we study a problem which includes both hysteresis and dynamic
capillary effects. We choose $\delta=5e-3$, $\gamma=1$ and
$\tau(\theta)=1+\theta^{2}$. As for the previous examples, we report the total
numbers of iterations required by each algorithm, the condition numbers
associated with the linearized equations, the EOC of the discretization
technique and the residual for each unknown.

In Tables 14 and 15, we present the total number of iterations required by
each algorithm, and the condition numbers associated with each system. As in
the previous example, the Newton method fails to converge, while the L-scheme
based solvers present no difficulties. The L parameters are all set equal to
$0.1$. This leads to a faster convergence, when compared to the previous
example, where larger values have been required for robustness. This is
explained by the fact that, since the dynamic effects are introduced
($\tau>0$), the solution is more regular [19, 47].

We have tested different values of $m$ on the Newton methods, but none ensured
the convergence of the schemes. As in the previous test cases, we have
investigated smaller time steps, but the Newton solvers has still failed to
converge.

Regarding the results obtained thanks to the AA, we can notice some
improvements which are smaller than the ones observed for the previous test
cases. Once more, this is due to the optimal choice of the $L$ parameters,
ensuring that the L-scheme converges, on average, in 5 iterations per time
step. Therefore further improvements are not expected. Note that, compared to
the first example (Table 4), larger $L$ values are used leading to larger
numbers of iterations. This explains why the AA with proper parameters $m$
have improved the convergence behaviour of the L-scheme there.

| Monolithic |  |  | NonLinS |  |  |  | AltLinS |  |
---|---|---|---|---|---|---|---|---|---|---
| Newton |  |  | Newton |  |  |  | Newton |  |
|  |  |  |  | cond. # |  |  |  | cond. # |
dx | # iterations | condition # | # iterations | Flow |  | Transport | # iterations | Flow |  | Transport
1/10 | - | - | - | - |  | - | - | - |  | -
1/20 | - | - | - | - |  | - | - | - |  | -
1/40 | - | - | - | - |  | - | - | - |  | -
| L-scheme |  |  | L-scheme |  |  |  | L-scheme |  |
|  |  |  |  | cond. # |  |  |  | cond. # |
dx | # iterations | condition # | # iterations | Flow |  | Transport | # iterations | Flow |  | Transport
1/10 | 152 | 290.13 | 128 - 122 | 208.09 |  | 162.46 | 251 | 206.24 |  | 174.77
1/20 | 160 | 768.54 | 137 - 121 | 558.01 |  | 621.77 | 259 | 486.23 |  | 598.01
1/40 | 165 | 3.04e+03 | 141 - 120 | 2.10e+03 |  | 2.35e+03 | 328 | 2.14e+03 |  | 2.4415e+03
| L-scheme | (AA m = 1) |  | L-scheme (AA m = 2, $m_{lin}$ = 3) |  |  |  | L-scheme (AA m = 2) |  |
|  |  |  |  | cond. # |  |  |  | cond. # |
dx | # iterations | condition # | # iterations | Flow |  | Transport | # iterations | Flow |  | Transport
1/10 | 139 | 288.40 | 112 - 85 | 212.29 |  | 165.24 | 152 | 198.15 |  | 166.47
1/20 | 144 | 752.41 | 117 - 89 | 550.05 |  | 630.43 | 162 | 522.01 |  | 636.00
1/40 | 149 | 3.04e+03 | 127 - 88 | 2.19e+03 |  | 2.39e+03 | 166 | 2.06e+03 |  | 2.40e+03

Table 14: Example 4: Total number of iterations and condition numbers for
fixed $\Delta t=T/25$, and different $dx$. Here, $L_{1}=L_{2}=L_{3}=0.1$ and
$m\neq m_{lin}$.

| Monolithic |  |  | NonLinS |  |  |  | AltLinS |  |
---|---|---|---|---|---|---|---|---|---|---
| Newton |  |  | Newton |  |  |  | Newton |  |
|  |  |  |  | cond. # |  |  |  | cond. # |
$\Delta t$ | # iterations | condition # | # iterations | Flow |  | Transport | # iterations | Flow |  | Transport
T/25 | - | - | - | - |  | - | - | - |  | -
T/50 | - | - | - | - |  | - | - | - |  | -
T/100 | - | - | - | - |  | - | - | - |  | -
| L-scheme |  |  | L-scheme |  |  |  | L-scheme |  |
|  |  |  |  | cond. # |  |  |  | cond. # |
$\Delta t$ | # iterations | condition # | # iterations | Flow |  | Transport | # iterations | Flow |  | Transport
T/25 | 152 | 290.13 | 128 - 122 | 208.09 |  | 162.46 | 251 | 206.24 |  | 174.77
T/50 | 248 | 310.93 | 201 - 225 | 264.09 |  | 89.64 | 424 | 263.91 |  | 97.09
T/100 | 508 | 415.80 | 403- 405 | 403.31 |  | 47.53 | 768 | 403.81 |  | 52.72
| L-scheme | (AA m = 1) |  | L-scheme (AA m = 2, $m_{lin}$ = 3) |  |  |  | L-scheme (AA m = 2) |  |
|  |  |  |  | cond. # |  |  |  | cond. # |
$\Delta t$ | # iterations | condition # | # iterations | Flow |  | Transport | # iterations | Flow |  | Transport
T/25 | 139 | 288.40 | 112 - 85 | 212.29 |  | 165.24 | 152 | 198.15 |  | 166.47
T/50 | 233 | 312.95 | 195 - 167 | 267.27 |  | 90.55 | 250 | 260.41 |  | 91.42
T/100 | 448 | 416.05 | 358 - 308 | 404.34 |  | 47.65 | 506 | 403.30 |  | 48.38

Table 15: Example 4: Total number of iterations and condition numbers for
fixed $dx=1/10$, and different $\Delta t$. Here $L_{1}=L_{2}=L_{3}=0.1$ and
$m\neq m_{lin}$

The numerical errors and the estimated orders of convergence of the
discretization technique (TPFA), presented in Table 16, are consistent with
the ones from the previous test cases.

| $e_{1}$ | EOC | $e_{2}$ | EOC | $e_{3}$
---|---|---|---|---|---
$\Psi$ | 0.0759 | 0.9558 | 0.0391 | 0.8837 | 0.0212
$\theta$ | 0.0138 | 0.9115 | 0.0073 | 0.9463 | 0.0038
$c$ | 0.0101 | 1.2531 | 0.0042 | 1.2655 | 0.0018
Table 16: Example 4: Numerical error and estimated order of convergence (EOC)
of the discretization method.

Finally, regarding the order of convergence of the different solving
algorithms, Figure 8 presents the residuals for each unknown, and Table 17 the
precise orders computed averaging over iterations at the final time step.

(a) Pressure residuals.

(b) Water content residuals.

(c) Concentration residuals.

Figure 8: Example 4: Residuals of each unknown at the final time step, for the
different schemes.
Here, $L_{1}=L_{2}=L_{3}=0.1$, $dx=1/40$, $\Delta t=T/25$ and different $m\neq m_{lin}$. | LS-Mono | LS-Mono And. | LS-NonLinS | LS-NonLinS And. | LS-AltLinS | LS-AltLinS And.
---|---|---|---|---|---|---
$\Psi$ | 1.00 | 1.29 | 0.99 | 1.27 | 0.99 | 1.30
$c$ | 1.29 | 1.29 | 1.02 | 1.15 | 1.07 | 1.11
$\theta$ | 1.00 | 1.16 | 1.00 | 1.35 | 1.00 | 1.20
Table 17: Example 4: Estimated order of convergence for the different
linearization schemes.

###  3.5 Physical example

As final numerical study, we investigate a test case that involves realistic
parameters, but without having a manufactured solution. The flow will be
driven by the boundary conditions. The domain $\Omega$ is the unit square and
the final time is $T=4$. The capillary pressure and conductivity expressions
are given by the van Genuchten formulation [31],
$K(\theta)=\theta_{e}^{l}\Big{(}1-\big{(}1-\theta_{e}^{1/M}\big{)}^{M}\Big{)}^{2}$
and
$p_{cap}(\theta,c)=\big{(}1-b\ln(c/a+1)\big{)}^{-1}(-\theta^{-1/M})^{1-M}$,
where $\theta_{e}=(\theta-\theta_{r})/(\theta_{s}-\theta_{r})$ is the
effective water content, $\theta_{s}=0.9$, $\theta_{r}=0.005$, $M=2$,
$l=0.31$, $a=0.04$ and $b=0.47$. Furthermore, we take
$\tau(\theta)=1+\theta^{2}$, and the hysteresis effects are included by
setting $\gamma=1$ and $\delta=5e-3$, as in Example 4.

Dirichlet boundary conditions are imposed at the left side of the unit square

$\begin{split}\Psi|_{x=0}&=1+\begin{cases}0.5t\qquad&\text{if}\ t<1,\\\
0.5\qquad&\text{if}\ 1\leq t<2,\\\ 0.5(3-t)\qquad&\text{if}\ 2\leq t<3,\\\
-0.4\qquad&\text{if}\ 3\leq t\leq 4,\\\ \end{cases}\\\
c|_{x=0}&=2,\end{split}$

whereas, at the remaining sides, we consider homogeneous Neumann boundary
conditions. The discontinuity in time $t=3$ makes solving the problem
numerically even more complex. The initial conditions are

$\theta^{0}=x,\qquad\text{and}\qquad c^{0}=1.$

All L parameters are set to $0.5$. We have tested different values, but,
$L=0.5$ seems to give the best results in terms of numbers of iterations.
Furthermore, the results may be improved even further by choosing different
values for each parameter, $L_{1},L_{2}$ and $L_{3}$, but for ease of
presentation this has been omitted here.

In Tables 18 and 19, we report the total numbers of iterations and condition
numbers associated to each algorithm. We observe that, due to the higher
nonlinearities of the conductivity $K$ and capillary pressure $p_{cap}$
involved, the results are different compared to the ones presented for the
previous examples.

Again, the Newton solvers have failed to converge and the systems associated
with the linearized equations are badly conditioned. Considering smaller time
steps did not resolve this.

The L-schemes on the other hand converge, but require high numbers of
iterations. We observe a significant improvement thanks to the AA. The
performance of the monolithic solver is for example drastically improved, the
iterations required are reduced by circa 50%. In case of finer meshes, one
needs to use a larger L, precisely $L=1$. The AA can introduce some
instabilities, and thus a larger L may be required. Clearly, this leads to an
increase in the number of iterations. Such results are still better than the
one obtained for smaller L without acceleration. Similar observations can also
be made for the splitting solvers. Even though larger L parameters may be
required, the accelerated schemes perform better then the non-accelerated
ones. The nonlinear splitting seems to be less stable than the alternate
linear one. For a coarse mesh, we could use a large $m$, $m=5$, resulting in
an extreme reduction in the numbers of iterations. For finer meshes, $m$ had
to be set equal to 1, otherwise the schemes did not converge, and a larger L
parameter was required. The alternate linearized splitting seems to be more
stable, as the L-scheme linearization acts as stabilization term for both the
nonlinearities and the decoupling. Unfortunately, for the first time, the
results are slower than for the nonlinear splitting. The main differences are
observable in the second table, Table 19. The mesh size is fixed, $dx=1/10$,
and the nonlinear solvers converge even for $m=5$. The large AA parameter $m$
ensures a strong reduction in the numbers of iterations. The alternate
linearized splitting converges only for $m=1$, anyhow the improvements are
remarkable.

| Monolithic |  |  | NonLinS |  |  |  | AltLinS |  |
---|---|---|---|---|---|---|---|---|---|---
| Newton |  |  | Newton |  |  |  | Newton |  |
|  |  |  |  | cond. # |  |  |  | cond. # |
dx | # iterations | condition # | # iterations | Flow |  | Transport | # iterations | Flow |  | Transport
1/10 | - | - | - | - |  | - | - | - |  | -
1/20 | - | - | - | - |  | - | - | - |  | -
1/40 | - | - | - | - |  | - | - | - |  | -
| L-scheme |  |  | L-scheme |  |  |  | L-scheme |  |
|  |  |  |  | cond. # |  |  |  | cond. # |
dx | # iterations | condition # | # iterations | Flow |  | Transport | # iterations | Flow |  | Transport
1/10 | 899 | 702.16 | 1522 - 235 | 790.37 |  | 155.18 | 1490 | 623.78 |  | 155.10
1/20 | 930 | 3.58e+03 | 1515 - 240 | 3.32e+03 |  | 615.72 | 1428 | 2.54e+03 |  | 596.60
1/40 | 941 | 1.52e+04 | 1680 - 243 | 1.44e+04 |  | 2.46e+03 | 1548 | 1.24e+04 |  | 2.45e+03
| L-scheme | (AA m = 1) |  | L-scheme (AA m = 5, $m_{lin}$ = 1) |  |  |  | L-scheme (AA m = 1) |  |
|  |  |  |  | cond. # |  |  |  | cond. # |
dx | # iterations | condition # | # iterations | Flow |  | Transport | # iterations | Flow |  | Transport
1/10 | 480 | 1.37e+03 | 210 - 96 | 811.29 |  | 155.16 | 795 | 774.00 |  | 155.18
1/20 | 541 ($L=1$) | 5.27e+03 | 758 - 170 ($m=1,L=1$) | 2.76e+03 |  | 395.74 | 532 | 2.72e+03 |  | 468.04
1/40 | 603 ($L=1$) | 2.29e+04 | 1214 - 261 ($m=1,L=2$) | 8.84e+03 |  | 911.31 | 1798 ($L_{2}=2$) | 1.89e+04 |  | 2.46e+03

Table 18: Example 5: Total number of iterations and condition numbers for
fixed $\Delta t=T/25$, and different $dx$. Here, $L_{1}=L_{2}=L_{3}=0.5$ and
$m\neq m_{lin}$.

| Monolithic |  |  | NonLinS |  |  |  | AltLinS |  |
---|---|---|---|---|---|---|---|---|---|---
| Newton |  |  | Newton |  |  |  | Newton |  |
|  |  |  |  | cond. # |  |  |  | cond. # |
dx | # iterations | condition # | # iterations | Flow |  | Transport | # iterations | Flow |  | Transport
T/25 | - | - | - | - |  | - | - | - |  | -
T/50 | - | - | - | - |  | - | - | - |  | -
T/100 | - | - | - | - |  | - | - | - |  | -
| L-scheme |  |  | L-scheme |  |  |  | L-scheme |  |
|  |  |  |  | cond. # |  |  |  | cond. # |
$\Delta t$ | # iterations | condition # | # iterations | Flow |  | Transport | # iterations | Flow |  | Transport
T/25 | 899 | 702.16 | 1522 - 235 | 790.37 |  | 155.18 | 1490 | 623.78 |  | 155.10
T/50 | 2892 | 908.39 | 5218 - 450 | 747.79 |  | 84.75 | 4942 | 533.58 |  | 84.83
T/100 | 9809 | 707.12 | 18261 - 849 | 692.77 |  | 46.72 | 17406 | 485.71 |  | 46.84
| L-scheme | (AA m =1) |  | L-scheme (AA m = 5, $m_{lin}$ = 1) |  |  |  | L-scheme (AA m = 1) |  |
|  |  |  |  | cond. # |  |  |  | cond. # |
$\Delta t$ | # iterations | condition # | # iterations | Flow |  | Transport | # iterations | Flow |  | Transport
T/25 | 480 | 1.37e+03 | 210 - 96 | 811.29 |  | 155.16 | 795 | 774.00 |  | 155.18
T/50 | 897 | 1.65e+03 | 517 - 214 | 618.31 |  | 53.28 | 1784 | 775.76 |  | 84.77
T/100 | 2106 | 1.04e+03 | 1354 - 424 | 576.54 |  | 28.62 | 4928 | 1.05e+03 |  | 46.73

Table 19: Example 5: Total number of iterations and condition numbers for
fixed $dx=1/20$, and different $\Delta t$. Here, $L_{1}=L_{2}=L_{3}=0.5$ and
$m\neq m_{lin}$.

It is also interesting to notice that, for decreasing time steps, the number
of L iterations per time step is increasing. This is coherent with the theory.
With the AA, this is mitigated; the average number of iterations remains more
stable.

In Table 20, we report the different $m$ and $m_{lin}$ values investigated for
the AA. We observe that for the L-scheme Mono solver, the optimal choice, in
terms of numbers of iterations is $m=1$. For the nonlinear splitting solver,
only one value of $m_{lin}$ has been taken in consideration for the coupling
loop, precisely $m_{lin}=1$. This is justified by the fact that the majority
of the iterations have taken place in the inside loops, the nonlinearities of
the equations are playing a larger role than the coupling aspect.

| # it. m=0 | # it. m=1 | # it. m=2 | # it. m=5
---|---|---|---|---
Mono L-scheme | 899 | 480 | 775 | –
NonLinS L-scheme | 1522 - 325 ($m_{lin}$=0) | 322 - 118 ($m_{lin}$=1) | 232 - 96 ($m_{lin}$=1) | 210 - 96 ($m_{lin}$=1)
AltLinS L-scheme | 1490 | 795 | – | –
Newton-Mono | – | – | – | –
Table 20: Example 5: Numbers of iterations associated to different $m$ and
$m_{lin}$ values, $dx=1/10$ $dt=T/25$.

(a) Pressure residuals.

(b) Water content residuals.

(c) Concentration residuals.

Figure 9: Example 5: Residuals of each unknown at the final time step,
monolithic L-scheme. Here, $L_{1}=L_{2}=L_{3}=0.1$, different $m$ are tested,
$dx=1/10$, and $\Delta t=T/25$.

(a) Pressure residuals.

(b) Water content residuals.

(c) Concentration residuals.

Figure 10: Example 5: Residuals of each unknown at the final time step,
AltLinS L-scheme. Here, $L_{1}=L_{2}=L_{3}=0.1$, different $m$ are tested,
$m_{lin}=1$, $dx=1/40$, and $\Delta t=T/25$.

(a) Pressure residuals.

(b) Water content residuals.

(c) Concentration residuals.

Figure 11: Example 5: Residuals of each unknown at the final time step,
NonLinS L-scheme. Here, $L_{1}=L_{2}=L_{3}=0.1$, different $m$ are tested ,
$dx=1/40$, and $\Delta t=T/25$.

Finally, in the Figures 9, 10 and 11 we report the residuals at the final time
step, for each unknown and each algorithm, and for different values of $m$.
The results are coherent with the ones presented in the Tables 18 and 19. The
monolithic solver shows a clear improvement thanks to the AA, for both $m=1$
and $m=2$. For the alternate linearized splitting solver, the AA does not seem
to produce any improvement, as the rates of convergence seem to be worse. The
result does not directly contradict the ones presented in the Tables 18 and
19; it simply states that, at the final time step, the AA does not produce any
improvement. On the full simulation, we could observe a clear reduction in the
numbers of iterations. For the nonlinear splitting, we can observe some
improvements but once more they are not as evident as for the monolithic
solvers. We already observed that both splitting solvers became unstable once
the AA was applied, either requiring a larger L or a smaller $m$.

Table 21 presents the precise rates of convergence of the different
linearization schemes. Once more we can observe as the AA improves the rates
of convergence of the solvers based on the L-scheme.

| LS-Mono | LS-Mono And. | LS-NonLinS | LS-NonLinS And. | LS-AltLinS | LS-AltLinS And.
---|---|---|---|---|---|---
$\Psi$ | 0.94 | 1.30 | 0.98 | 1.87 | 0.98 | 1.21
$c$ | 1.01 | 1.78 | 0.95 | 1.42 | 0.93 | 1.80
$\theta$ | 0.95 | 1.34 | 0.99 | 1.35 | 0.95 | 1.02
Table 21: Example 5: Order of convergence of the linearization schemes.

##  4 Conclusions

We consider models for flow and reactive transport in a porous medium. Next,
to account for the influence of the solute concentration on the flow
parameters, we incorporate effects like dynamic capillary pressure and
hysteresis. The problem results being fully coupled.

For solving the time discrete equations (9), obtained after applying the Euler
implicit scheme, we investigate different approaches: a monolithic solution
algorithm and two splitting ones. Furthermore, for solving the nonlinear
problem, two linearizations are studied: the Newton method and the L-scheme.
The latter appears to be more stable than the former, which is more commonly
implemented.

Finally, we have studied the effects of the Anderson acceleration. We observed
that its implementation is particularly simple and can result in significant
improvements. There were cases in which the differences between the
accelerated and non-accelerated schemes were minimal, but due to its
simplicity and the possibility of the great reduction in the numbers of
iterations, we think it should always be tested. Particularly, one can either
invest time in finding the optimal L parameters or the best depth $m$ for
which the AA results in the fastest scheme. Often, finding the most suitable
$m$ is simpler, and it can results in impressive improvements.

######  Acknowledgments.

The research of D. Illiano was funded by VISTA, a collaboration between the
NorwegianAcademy of Science and Letters and Equinor, project number 6367,
project name: adaptive model and solver simulation of enhanced oil recovery.
The research of J.W. Both was supported by the Research Council of Norway
Project 250223, as well as the FracFlow project funded by Equinor through
Akademiaavtalen. The research of I.S. Pop was supported by the Research
Foundation-Flanders (FWO), Belgium through the Odysseus programme (project
G0G1316N) and Equinor through the Akademia grant.

## References

  * Aavatsmark- [2001] I. Aavatsmark, _An introduction to multipoint flux approximations for quadrilateral grids, Computational Geosciences Volume 6, Issue 3-4, Pages 405-432_ , 2002\.
  * Abreu et al.- [217] E. Abreu, J. Vieira, _Computing numerical solutions of the pseudo-parabolic Buckley Leverett equation with dynamic capillary pressure, Mathematics and Computers in Simulation, Volume 137, Pages 29-48_ , 2017.
  * Abreu et al.- [2020] E. Abreu, P. Ferraz, J. Vieira, _Numerical resolution of a pseudo-parabolic Buckley-Leverett model with gravity and dynamic capillary pressure in heterogeneous porous media, Journal of computational Physics, Volume 411_ , 2020.
  * Agosti et al.- [2015] A. Agosti, L. Formaggia, A. Scotti, _Analysis of a model for precipitation and dissolution coupled with a Darcy flux, Journal of Mathematical Analysis and Applications, Volume 431, Issue 2, Pages 752-781_ , 2015.
  * Albuja et al.- [2021] G. Albuja, A. I. Avila, _A family of new globally convergent linearization schemes for solving Richards’ equation, Applied Numerical Mathematics, Volume 159, Pages 281-296_ , 2021.
  * Alt et al.- [1983] W. Alt, H. Luckhaus, _Quasilinear elliptic-parabolic differential equations, Mathematische Zeitschrift, Volume 183, Issue 3, Pages 311-341_ , 1983\.
  * Anderson- [1965] D. G. Anderson, _Iterative Procedures for Nonlinear Integral Equations, Journal of the ACM, Volume 12, Issue 4, Pages 547-560_ , 1965.
  * Arbogast- [1996] T. Arbogast, M. F. Wheeler, _A nonlinear mixed finite element method for a degenerate parabolic equation arising in flow in porous media, SIAM Journal on Numerical Analysis, Volume 33, Issue 4, Pages 1669-1687_ , 1996\.
  * Arnold et al.- [2006] D. N. Arnold, F. Brezzi, B. Cockburn, L. D. Marini, _Unified Analysis of Discontinuous Galerkin Methods for Elliptic Problems, SIAM Journal on Numerical Analysis, Volume 39, Issue 5, Pages 1749-1779_ , 2006.
  * Arraras et al.- [2020] A. Arraras, F. J. Gaspar, L. Portero, C. Rodrigo, _Multigrid solvers for multipoint flux approximations of the Darcy problem on rough quadrilateral grids, Computational Geosciences, https://doi.org/10.1007/s10596-020-09979-w,_ 2020\.
  * Barrett et al.- [1997] J. W. Barrett, P. Knabner, _Finite element approximation of the transport of reactive solutes in porous media. Part 1: error estimates for nonequilibrium adsorption processes, SIAM Journal on Numerical Analysis, Volume 34, Issue 1, Pages 201-227_ , 1997\.
  * Bause et al.- [2010] M. Bause, J. Hoffmann, P. Knabner, _First-order convergence of multi-point flux approximation on triangular grids and comparison with mixed finite element methods, Numerische Mathematik, Volume 116, Issue 1, Pages 1-29_ , 2010\.
  * Beliaev et al.- [2001] A.Y. Beliaev, S.M. Hassanizadeh, _A Theoretical Model of Hysteresis and Dynamic Effects in the Capillary Relation for Two-phase Flow in Porous Media. Transport in Porous Media, Volume 43, Issue 3, Pages 487-510_ , 2001.
  * Berardi et al.- [2018] M. Berardi, F. Difonzo, M. Vurro, L. Lopez, _The 1D Richards’ equation in two layered soils: a Filippov approach to treat discontinuities, Advances in Water Resources, Volume 115, Pages 264-272_ , 2018\.
  * Berardi et al.- [2020] M. Berardi, F. Difonzo, L. Lopez, _A mixed MoL-TMoL for the numerical solution of the 2D Richards’ equation in layered soils, Computers & Mathematics with Applications Volume 79, Issue 7, Pages 1990-2001_, 2020.
  * Both et al.- [2019] J. W. Both, K. Kumar, J. M. Nordbotten, F. A. Radu, _Anderson accelerated fixed-stress splitting schemes for consolidation of unsaturated porous media, Computers & Mathematics with Applications, Volume 77, Issue 6, Pages 1479-1502_, 2019.
  * Bottero et al.- [2011] S. Bottero, S. M. Hassanizadeh, P. J. Kleingeld, T. J. Heimovaara, _Nonequilibrium capillarity effects in two-phase flow through porous media at different scales, Water Resources Research, Volume 47, Issue 10_ , 2011.
  * Camps-Roach et al.- [2010] G. Camps-Roach, D. M. O’Carroll, T. A. Newton, T. Sakaki, T. H. Illangasekare, _Experimental investigation of dynamic effects in capillary pressure: Grain size dependency and upscaling, Water Resources Research, Volume 46, Issue 8,_ 2010\.
  * Cao et al.- [2015] X. Cao, I. S. Pop, _Two-phase porous media flows with dynamic capillary effects and hysteresis: Uniqueness of weak solutions, Computers $\&$ Mathematics with Applications, Volume 69, Issue 7, Pages 688-695_, 2015.
  * Cao et al.- [2016] X. Cao, I. S. Pop, _Degenerate two-phase porous media flow model with dynamic capillarity, Journal of Differential Equations Volume 260, Issue 3, Pages 2418-2456,_ 2016\.
  * [21] X. Cao, K. Mitra, _Error estimates for a mixed finite element discretization of a two-phase porous media flow model with dynamic capillarity, Journal of Computational and Applied Mathematics, Volume 353, Pages 164 - 178_ , 2019.
  * [22] X. Cao, S.F. Nemadjieu, I.S. Pop, _Convergence of an MPFA finite volume scheme for two phase porous media flow with dynamic capillarity, IMA Journal of Numerical Analysis, Volume 39, Issue 1, Pages 512–544,_ 2019\.
  * O’Carroll et al.- [2005] D. M. O’Carroll, T. J. Phelan, L. M. Abriola, _Exploring dynamic effects in capillary pressure in multistep outflow experiments, Water Resources Research, Volume 41, Issue 11_ , 2005.
  * Celia et al.- [1990] M. Celia, E. Bouloutas, R. L. Zarba, _A General Mass-Conservative Numerical Solution for the Unsaturated Flow Equation, Advances in Water Resources, Volume 26, Issue 7, Pages 1483-1496_ , 1990\.
  * DiCarlo- [2004] D. A. DiCarlo, _Experimental measurements of saturation overshoot on infiltration, Water Resources Research, Volume 40, Issue 4_ , 2004.
  * Dolejsi et al. [2019] V. Dolejsi, M. Kuraz, P. Solin, _Adaptive higher-order space-time discontinuous Galerkin method for the computer simulation of variably-saturated porous media flows, Applied Mathematical Modelling, Volume 72, Pages 276-305_ , 2019.
  * Eymard et al.- [1999] R. Eymard, M. Gutnic, D. Hilhorst, _The finite volume method for Richards equation, Computational Geosciences, Volume 3, Issue 3-4, Pages 256-294_ , 1999.
  * Evans et al.- [2020] C. Evans, S. Pollock, L. G. Rebholz, M. Xiao, _A proof that Anderson acceleration improves the convergence rate in linearly converging fixed point methods (but not in those converging quadratically), SIAM Journal on Numerical Analysis, Volume 58, Issue 1, Pages 788-810_ , 2020\.
  * Fang et al.- [2008] H. Fang, Y. Saad, _Two classes of multisecant methods for nonlinear acceleration, Numerical Linear Algebra with Application, Volume 16, Issue 3, Pages 197 - 221_ , 2008.
  * Gray et al.- [1998] W. G. Gray, S. M. Hassanizadeh, _Macroscale continuum mechanics for multiphase porous-media flow including phases, interfaces, common lines and common points, Volume 21, Issue 4, Pages 261 - 281_ , 1998\.
  * van Genuchten- [1980] M. van Genuchten, _A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils, Soil Science Society of America Journal, Volume 44, Issue 5, Pages 892-898_ , 1980\.
  * Helmig- [1997] R. Helmig, _two-phase flow and transport processes in the subsurface: a contribution to the modeling of hydrosystems, Springer-Verlag_ , 1997\.
  * Henry et al.- [1999] E. J. Henry, J. E. Smith, A. W. Warrick, _Solubility effects on surfactant-induced unsaturated flow through porous media, Journal of Hydrology, Volume 223, Issues 3-4, Pages 164-174_ , 1999\.
  * Hoa et al.- [1977] N. T. Hoa, R. Gaudu, C. Thirriot, _Influence of the hysteresis effect on transient flows in saturated-unsaturated porous media, Water Resources Research, Volume 13, Issue 6, Pages 992-996_ , 1977.
  * Husseini- [2015] D. Husseini, _Effects of Anions acids on Surface Tension of Water, Undergraduate Research at JMU Scholarly Commons_ , 2015\.
  * [36] D. Illiano, I. S. Pop, F. A. Radu, _Iterative schemes for surfactant transport in porous media, Computational Geosciences, https://doi.org/10.1007/s10596-020-09949-2_ , 2020\.
  * [37] D. Illiano, I. S. Pop, F. A. Radu, _An efficient numerical scheme for fully coupled flow and reactive transport in variably saturated porous media including dynamic capillary effects, accepted author manuscript, Numerical Mathematics and Advanced Applications ENUMATH 2019, Lecture Notes in Computational Science and Engineering 139, https://doi.org10.1007/978-3-030-55874-1 $\\_$55_, 2020\.
  * Karagunduz et al.- [2015] A. Karagunduz, M. H. Young, K. D. Pennell, _Influence of surfactants on unsaturated water flow and solute transport, Water Resources Research, Volume 51, Issue 4, Pages 1977-1988_ , 2015\.
  * Karpinski et al.- [2017] S. Karpinski, I. S. Pop, F. A. Radu, _Analysis of a linearization scheme for an interior penalty discontinuous Galerkin method for two-phase flow in porous media with dynamic capillarity effects, International Journal for Numerical Methods in Engineering, Volume 112, Issue 6, Pages 553-577_ , 2017.
  * Klausen et al.- [2008] R. A. Klausen, F. A. Radu, G. T. Eigestad, _Convergence of MPFA on triangulations and for Richards’ equation, International Journal for Numerical Methods in Fluids, Volume 58, Issue 12, Pages 1327-1351_ , 2008.
  * Koch et al.- [2012] J. Koch, A. Ratz, B. Schweizer, _Two-phase flow equations with a dynamic capillary pressure, European Journal of Applied Mathematics, Volume 24, Issue 1_ , 2012.
  * Knabner et al.- [2003] P. Knabner, S. Bitterlich, R. I. Teran, A. Prechtel, E. Schneid, _Influence of Surfactants on Spreading of Contaminants and Soil Remediation, Jager W., Krebs HJ. (eds) Mathematics — Key Technology for the Future. Springer, https://doi.org/10.1007/978-3-642-55753-8-12_ , 2003\.
  * Li et al.- [2007] H. Li, M. W. Farthing, C. N. Dawson, C. T. Miller, _Local discontinuous Galerkin approximations to Richards’ equation, Advances in Water Resources, Volume 30, Issue 3, Pages 555-575_ , 2007.
  * Lie- [2016] K.-A. Lie, _An Introduction to Reservoir Simulation Using MATLAB: User guide for the Matlab Reservoir Simulation Toolbox (MRST), SINTEF ICT_ , 2016\.
  * List et al.- [2016] F. List, F. A. Radu, _A study on iterative methods for solving Richards’ equation, Computational Geoscience, Volume 20, Issue 2, Pages 341-353_ , 2016\.
  * Lunowa et al.- [2020] S. B. Lunowa, I. S. Pop, B. Koren, _Linearized domain decomposition methods for two-phase porous media flow models involving dynamic capillarity and hysteresis, Computer Methods in Applied Mechanics and Engineering, Volume 372, https://doi.org/10.1016/j.cma.2020.113364_ , 2020.
  * Mikelic- [2010] A. Mikelic, _A global existence result for the equations describing unsaturated flow in porous media with dynamic capillary pressure, Journal of Differential Equations, Volume 248, Issue 6, Pages 1561-1577_ , 2010.
  * Milisic - [2018] J.-P. Milisic, _The unsaturated flow in porous media with dynamic capillary pressure, Journal of Differential Equations, Volume 264, Issue 9, Pages 5629-5658_ , 2018.
  * Mitra et al.- [2019] K. Mitra, I. S. Pop, _A modified L-scheme to solve nonlinear diffusion problems, Computers $\&$ Mathematics with Applications, Volume 77, Issue 6, Pages 1722-1738_, 2019\.
  * McClure et al.- [2018] J. McClure, R. T. Armstrong, M. Berrill, S. Schluter, S. Berg, W. G. Gray, C. T. Miller, _Geometric state function for two-fluid flow in porous media, Physical Review Fluids, Volume 3, Issue 8_ , 2018.
  * Morrow et al.- [2005] N. R. Morrow, C. C. Harris, _Exploring dynamic effects in capillary pressure in multistep outflow experiments, Water Resources Research, Volume 41, Issue 11_ , 2005
  * Nocheto et al.- [1988] R. Nochetto, C. Verdi, _Approximation of degenerate parabolic problems using numerical integration, SIAM Journal on Numerical Analysis, Volume 25, Issue 4, Pages 784-814_ , 1988.
  * Oung et al.- [2005] O, Oung, S. M. Hassanizadeh, A. Bezuijen, _Two-phase flow experiments in a geocentrifuge and the significance of dynamic capillary pressure effect, Journal of Porous Media, Volume 8, Issue 3, Pages 247-257_ , 2005.
  * Paniconi et al.- [1994] C. Paniconi, M. Putti, _A comparison of Picard and Newton iteration in the numerical solution of multidimensional variably saturated flow problems, Water Resources Research, Volume 30, Issue 12, Pages 3357-3374_ , 1994.
  * Peszynska et al.- [2008] M. Peszynska, S. Y. Yi, _Numerical methods for unsaturated flow with dynamic capillary pressure in heterogeneous porous media, International journal of Numerical Analysis and Modeling, Volume 5, Pages 126-149_ , 2008\.
  * Prechtel et al.- [2002] A. Prechtel, P. Knabner, _Accurate and efficient simulation of coupled water flow and nonlinear reactive transport in the saturated and vadose zone - application to surfactant enhanced and intrinsic bioremediation, International Journal of Water Resources Development, Volume 47, Pages 687-694_ , 2002\.
  * Pop et al.- [2004] I. S. Pop, F. A. Radu, P. Knabner, _Mixed finite elements for the Richards’ equation: linearization procedure, Journal of Computational and Applied Mathematics, Volume 168, Issue 1, Pages 365-373_ , 2004\.
  * Radu et al.- [2010] F. A. Radu, I. S. Pop, S. Attinger, _Analysis of an Euler implicit, mixed finite element scheme for reactive solute transport in porous media, Numerical Methods for Partial Differential Equations, Volume 26, Issue 2, Pages 320-344_ , 2010\.
  * Radu et al.- [2013] F. A. Radu, A. Muntean, I. S. Pop, N. Suciu, O. Kolditz, _A mixed finite element discretization scheme for a concrete carbonation model with concentration-dependent porosity, Journal of Computational and Applied Mathematics Volume 246, Pages 74-85,_ 2013\.
  * Russell et al.- [1983] T. F. Russell, M. F. Wheeler, _Finite element and finite difference methods for continuous flows in porous media, SIAM, Pages 35-106_ , 1983\.
  * Schweizer- [2012] B. Schweizer, _The Richards equation with hysteresis and degenerate capillary pressure, Journal of Differential Equations, Volume 252, Issue 10, Pages 5594-5612_ , 2012.
  * Slodicka- [2002] M. Slodicka, _A robust and efficient linearization scheme for doubly nonlinear and degenerate parabolic problems arising in flow in porous media, SIAM Journal on Numerical Analysis, Volume 23, Issue 5, Pages 1593-1614_ , 2002\.
  * Smith et al.- [1994] J. E. Smith, R. W. Gillham, _The effect of concentration-dependent surface tension on the flow of water and transport of dissolved organic compounds: A pressure head-based formulation and numerical model, Water Resources Research, Volume 31, Issue 3, Pages 343-354_ , 1994\.
  * Smith et al.- [1999] J. Smith, R. Gillham, _Effects of solute concentration-dependent surface tension on unsaturated flow: Laboratory sand column experiments, Water Resource Research, Volume 35, Issue 4, Pages 973-982_ , 1999\.
  * Stauffer- [1978] F. Stauffer, _Time dependence of the relations between capillary pressure, water content and conductivity during drainage of porous media, IAHR symposium on scale effects in porous media, Volume 29_ , 1978\.
  * Sun et al.- [2005] S. Sun, M. F. Wheeler, _Discontinuous Galerkin methods for coupled flow and reactive transport problems, Applied Numerical Mathematics Volume 52, Issues 2–3, Pages 273-298_ , 2005.
  * Vohralik- [2007] M. Vohralik, _A posteriori error estimates for lowest-order mixed finite element discretizations of convection-diffusion-reaction equations, SIAM Journal on Numerical Analysis, Volume 45, Issue 4, Pages 1570-1599_ , 2007\.
  * Walker et al.- [2011] H. F. Walker, P. Ni, _Anderson Acceleration for Fixed-Point Iterations, SIAM Journal on Numerical Analysis, Volume 49, Issue 4, Pages 1715 - 1735_ , 2011\.
  * Wang et al.- [2013] X. Wang, H. A. Tchelepi, _Trust-region based solver for nonlinear transport in heterogeneous porous media, Journal of Computational Physics, Volume 253, Pages 114-137_ , 2013.
  * Woodward et al.- [2000] C. S. Woodward, C. N. Dawson, _Analysis of expanded mixed finite element methods for a nonlinear parabolic equation modeling flow into variably saturated porous media, SIAM Journal on Numerical Analysis, Volume 37, Issue 3, Pages 701-724_ , 2000\.
  * Zha et al. [2019] Y. Zha, J. Yang, J. Zeng, C.-H. M. Tso, W. Zeng, L. Shi, _Review of numerical solution of Richardson-Richards equation for variably saturated flow in soils, WIREs Water, Volume 6, Issue 5_ , 2019.
  * Zhang et al.- [2017] H. Zhang, P. A. Zegeling, _A Numerical Study of Two-Phase Flow Models with Dynamic Capillary Pressure and Hysteresis, Transport in Porous Media, Volume 116, Issue 2, Pages 825-846_ , 2017\.