Computational orders of convergence of iterative methods for Richards’ equation

Nicolae Suciu
(preprint), Convergence order, Finite element (FEM), html, Numerical Analysis, paper, PDEs (numerical)

Abstract

Numerical solutions for flows in partially saturated porous media pose challenges related to the non-linearity and elliptic-parabolic degeneracy of the governing Richards’ equation. Iterative methods are therefore required to manage the complexity of the flow problem. Norms of successive corrections in the iterative procedure form sequences of positive numbers. Definitions of computational orders of convergence and theoretical results for abstract convergent sequences can thus be used to evaluate and compare different iterative methods. We analyze in this frame Newton’s and L-scheme methods for an implicit finite element method (FEM) and the L-scheme for an explicit finite difference method (FDM). We also investigate the effect of the Anderson Acceleration (AA) on both the implicit and the explicit L-schemes. Considering a two-dimensional test problem, we found that the AA halves the number of iterations and renders the convergence of the FEM scheme two times faster. As for the FDM approach, AA does not reduce the number of iterations and even increases the computational effort. Instead, being explicit, the FDM L-scheme without AA is faster and as accurate as the FEM L-scheme with AA.

Authors

Nicolae Suciu
Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy

Florin A. Radu
Center for Modeling of Coupled Subsurface Dynamics, University of Bergen, Bergen, Norway

Jacob S. Stokke
Center for Modeling of Coupled Subsurface Dynamics, University of Bergen, Bergen, Norway

Emil Cătinaș
Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy

Andra Malina
Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy

Keywords

Richards’ equation; flows in partially saturated porous media; Newton’s method; L-scheme method; finite element method (FEM); finite difference method; Anderson Acceleration; computational convergence orders; Q-order; C-order; corrections; errors.

Paper coordinates

N. Suciu, F.A. Radu, J.S. Stokke, E. Cătinaș, A. Malina, Computational orders of convergence of iterative methods for Richards’ equation, in: Sequeira, A., Silvestre, A., Valtchev, S.S., Janela, J. (eds) Numerical Mathematics and Advanced Applications ENUMATH 2023, vol. 2. ENUMATH 2023. Lecture Notes in Computational Science and Engineering, vol. 154 (2025), pp. 379–388. Springer, Cham. https://doi.org/10.1007/978-3-031-86169-7_39

PDF

https://arxiv.org/pdf/2402.00194.pdf

About this paper

Journal

Preprint, Arxiv

Publisher Name
DOI

https://doi.org/10.48550/arXiv.2402.00194

Print ISSN
Online ISSN

google scholar link

[1] Anderson, D.G., Iterative procedures for nonlinear integral equations. JACM 12(4), 547-560, (1965)
[2] Beyer, W. A., Ebanks, B. R., Qualls, C. R., Convergence rates and convergence-order profiles for sequences. Acta Appl. Math. 20, 267–284, (1990)
[3] Catinas¸, E., A survey on the high convergence orders and computational convergence orders of sequences. Appl. Math. Comput. 343, 1–20, (2019)
[4] Catinas¸, E., How many steps still left to x*? SIAM Rev. 63(3), 585–624, (2021)
[5] List, F., Radu, F.A., A study on iterative methods for solving Richards’ equation. Comput. Geosci. 20(2) 341-353, (2016)
[6] Potra, F. A., On Q-order and R-order of convergence. J. Optim. Theory Appl. 63(3), 415–431,(1989)
[7] Radu, F.A., Pop, I.S., Knabner, P., On the convergence of the Newton method for the mixed finite element discretization of a class of degenerate parabolic equations. In: Numerical Mathematics and Advanced Applications 42, 1192-1200, (2006)
[8] Suciu, N., Illiano, D., Prechtel, A., Radu, F.A., Global random walk solvers for fully coupled flow and transport in saturated/unsaturated porous media. Adv. Water Resour. 152, 103935, (2021)
[9] Suciu, N., Radu, F.A., Global random walk solvers for reactive transport and biodegradation processes in heterogeneous porous media. Adv. Water Res., 166, 104268, (2022)
[10] Stokke, J.S., Mitra, K., Storvik, E., Both, J.W., Radu, F.A., An adaptive solution strategy for Richards’ equation. Comput. Math. Appl. 152, 155-167, (2023)

Paper (preprint) in HTML form

2402.00194v1

Computational orders of convergence of iterative methods for Richards' equation

Nicolae Suciu, Florin A. Radu, Jakob S. Stokke, Emil Cătinaş, Andra Malina

Abstract

Numerical solutions for flows in partially saturated porous media pose challenges related to the non-linearity and elliptic-parabolic degeneracy of the governing Richards' equation. Iterative methods are therefore required to manage the complexity of the flow problem. Norms of successive corrections in the iterative procedure form sequences of positive numbers. Definitions of computational orders of convergence and theoretical results for abstract convergent sequences can thus be used to evaluate and compare different iterative methods. We analyze in this frame Newton's and L L LLL-scheme methods for an implicit finite element method (FEM) and the L L LLL-scheme for an explicit finite difference method (FDM). We also investigate the effect of the Anderson Acceleration (AA) on both the implicit and the explicit L L LLL schemes. Considering a two-dimensional test problem, we found that the AA halves the number of iterations and renders the convergence of the FEM scheme two times faster. As for the FDM approach, AA does not reduce the number of iterations and even increases the computational effort. Instead, being explicit, the FDM L L LLL-scheme without AA is faster and as accurate as the FEM L L LLL-scheme with AA.

1 Introduction

The convergence of the iterative methods for Richards' equation can be proved under more or less restrictive smoothness and boundedness conditions that may not be met in some practical applications. Alternatively, the convergence may be assessed or verified numerically by investigating the behavior of the successive corrections.
This could be mainly useful in the case of L L LLL-schemes, which do not require the evaluation of the derivatives and, as such, can be applied to problems with nonsmooth coefficients and boundary conditions. Although the decay of the correction norms allows the estimation of the convergence order, this may not be enough to prove that a specific convergence behavior holds, mainly when trying to prove the linear convergence of a L L LLL-scheme. General methods developed for convergent sequences of real numbers allow a more in-depth convergence analysis and the evaluation of the relaxation strategies, such as AA, used to improve the iterative methods.
The convergence of an arbitrary sequence of real numbers x s → x ∗ ∈ R x s → x ∗ ∈ R x_(s)rarrx^(**)inRx_{s} \rightarrow x^{*} \in \mathbb{R}xs→x∗∈R is characterized by the behavior of the successive errors e s = | x ∗ − x s | e s = x ∗ − x s e_(s)=|x^(**)-x_(s)|e_{s}=\left|x^{*}-x_{s}\right|es=|x∗−xs|. The sequence { x s } x s {x_(s)}\left\{x_{s}\right\}{xs} converges with the (classical) C C CCC-order p ≥ 1 p ≥ 1 p >= 1p \geq 1p≥1 if
(1) (C) lim s → ∞ e s + 1 ( e s ) p = Q p ∈ { ( 0 , ∞ ) p > 1 ( 0 , 1 ) p = 1 , (1)  (C)  lim s → ∞   e s + 1 e s p = Q p ∈ ( 0 , ∞ ) p > 1 ( 0 , 1 ) p = 1 , {:(1)" (C) "quadlim_(s rarr oo)(e_(s+1))/((e_(s))^(p))=Q_(p)in{[(0","oo),p > 1],[(0","1),p=1","]:}:}\text { (C) } \quad \lim _{s \rightarrow \infty} \frac{e_{s+1}}{\left(e_{s}\right)^{p}}=Q_{p} \in \begin{cases}(0, \infty) & p>1 \tag{1}\\ (0,1) & p=1,\end{cases}(1) (C) lims→∞es+1(es)p=Qp∈{(0,∞)p>1(0,1)p=1,
which implies the (weaker) orders Q Q QQQ and R 1 R 1 R^(1)R^{1}R1 (see [6, 2, 3, 4]),
(3) ( Q ) lim s → ∞ ln ⁡ e s + 1 ln ⁡ e s = p ( R ) lim s → ∞ | ln ⁡ e s | 1 / s = p (3) ( Q ) lim s → ∞   ln ⁡ e s + 1 ln ⁡ e s = p ( R ) lim s → ∞   ln ⁡ e s 1 / s = p {:(3){:[(Q),lim_(s rarr oo)(ln e_(s+1))/(ln e_(s))=p],[(R),lim_(s rarr oo)|ln e_(s)|^(1//s)=p]:}:}\begin{array}{ll} (Q) & \lim _{s \rightarrow \infty} \frac{\ln e_{s+1}}{\ln e_{s}}=p \\ (R) & \lim _{s \rightarrow \infty}\left|\ln e_{s}\right|^{1 / s}=p \tag{3} \end{array}(3)(Q)lims→∞ln⁡es+1ln⁡es=p(R)lims→∞|ln⁡es|1/s=p
If the limit x ∗ x ∗ x^(**)x^{*}x∗ is not known, the errors e s e s e_(s)e_{s}es are replaced by the corrections | x s + 1 − x s | x s + 1 − x s |x_(s+1)-x_(s)|\left|x_{s+1}-x_{s}\right||xs+1−xs| and (1-3) define "computational orders of convergence" C ′ , Q ′ C ′ , Q ′ C^('),Q^(')C^{\prime}, Q^{\prime}C′,Q′, and R ′ R ′ R^(')R^{\prime}R′. If p > 1 p > 1 p > 1p>1p>1, computational and error-based orders of convergence are equivalent, denoted by using curly braces, and they are related by [4]
{ C , C ′ } ⇒ ⇍ { Q , Q ′ } ⇒ ⇍ { R , R ′ } . C , C ′ ⇒ ⇍ Q , Q ′ ⇒ ⇍ R , R ′ . {C,C^(')}=>⇍{Q,Q^(')}=>⇍{R,R^(')}.\left\{C, C^{\prime}\right\} \underset{\nLeftarrow}{\Rightarrow}\left\{Q, Q^{\prime}\right\} \underset{\nLeftarrow}{\Rightarrow}\left\{R, R^{\prime}\right\} .{C,C′}⇒⇍{Q,Q′}⇒⇍{R,R′}.
However, when p = 1 p = 1 p=1p=1p=1, the above equivalences do not hold. In [2] was given an example of a sequence with C C CCC-linear order but no C ′ C ′ C^(')C^{\prime}C′-linear order. Also, in Sect. 4.2 below we present the situation of two sequences where the numerical estimates of (1) indicate C C CCC-sublinear and C ′ C ′ C^(')C^{\prime}C′-linear orders of convergence.
In this framework, we examine the Newton's and L L LLL-scheme methods utilized in implicit FEM approaches [5,10] as well as an explicit FDM L L LLL-scheme [8, 9].

2 The setup of the flow problem

Combining the mass balance equation and Darcy's law, one obtains the Richards' equation for the pressure head ψ ( x , z , t ) ψ ( x , z , t ) psi(x,z,t)\psi(x, z, t)ψ(x,z,t) in variably saturated soils,
(4) ∂ θ ( ψ ) ∂ t − ∇ ⋅ [ K ( θ ( ψ ) ∇ ( ψ + z ) ] = f (4) ∂ θ ( ψ ) ∂ t − ∇ ⋅ [ K ( θ ( ψ ) ∇ ( ψ + z ) ] = f {:(4)(del theta(psi))/(del t)-grad*[K(theta(psi)grad(psi+z)]=f:}\begin{equation*} \frac{\partial \theta(\psi)}{\partial t}-\nabla \cdot[K(\theta(\psi) \nabla(\psi+z)]=f \tag{4} \end{equation*}(4)∂θ(ψ)∂t−∇⋅[K(θ(ψ)∇(ψ+z)]=f
where z z zzz is the height against the gravitational direction, θ θ theta\thetaθ is the water content, K K KKK stands for the hydraulic conductivity of the medium, and f f fff is a source/sink term.
Equation (4) is closed by relationships defining the water content θ ( ψ ) θ ( ψ ) theta(psi)\theta(\psi)θ(ψ) and the hydraulic conductivity K ( θ ( ψ ) ) K ( θ ( ψ ) ) K(theta(psi))K(\theta(\psi))K(θ(ψ)) given by the van Genuchten-Mualem model
(5) Θ ( ψ ) = { ( 1 + ( − α ψ ) n ) − m ψ < 0 1 ψ ≥ 0 , (6) K ( Θ ( ψ ) ) = { K s Θ ( ψ ) 1 2 [ 1 − ( 1 − Θ ( ψ ) 1 m ) m ] 2 , ψ < 0 K s ψ ≥ 0 , (5) Θ ( ψ ) = 1 + ( − α ψ ) n − m ψ < 0 1 ψ ≥ 0 , (6) K ( Θ ( ψ ) ) = K s Θ ( ψ ) 1 2 1 − 1 − Θ ( ψ ) 1 m m 2 , ψ < 0 K s ψ ≥ 0 , {:[(5)Theta(psi)={[(1+(-alpha psi)^(n))^(-m),psi < 0],[1,psi >= 0","]:}],[(6)K(Theta(psi))={[K_(s)Theta(psi)^((1)/(2))[1-(1-Theta(psi)^((1)/(m)))^(m)]^(2)",",psi < 0],[K_(s),psi >= 0","]:}]:}\begin{gather*} \Theta(\psi)= \begin{cases}\left(1+(-\alpha \psi)^{n}\right)^{-m} & \psi<0 \\ 1 & \psi \geq 0,\end{cases} \tag{5}\\ K(\Theta(\psi))= \begin{cases}K_{s} \Theta(\psi)^{\frac{1}{2}}\left[1-\left(1-\Theta(\psi)^{\frac{1}{m}}\right)^{m}\right]^{2}, & \psi<0 \\ K_{s} & \psi \geq 0,\end{cases} \tag{6} \end{gather*}(5)Θ(ψ)={(1+(−αψ)n)−mψ<01ψ≥0,(6)K(Θ(ψ))={KsΘ(ψ)12[1−(1−Θ(ψ)1m)m]2,ψ<0Ksψ≥0,
where Θ = ( θ − θ r ) / ( θ s − θ r ) Θ = θ − θ r / θ s − θ r Theta=(theta-theta_(r))//(theta_(s)-theta_(r))\Theta=\left(\theta-\theta_{r}\right) /\left(\theta_{s}-\theta_{r}\right)Θ=(θ−θr)/(θs−θr) is the normalized water content, θ r θ r theta_(r)\theta_{r}θr and θ s θ s theta_(s)\theta_{s}θs are the residual and the saturated water content, K s K s K_(s)K_{s}Ks is the saturated hydraulic conductivity, and α , n , m = 1 − 1 / n α , n , m = 1 − 1 / n alpha,n,m=1-1//n\alpha, n, m=1-1 / nα,n,m=1−1/n are model parameters depending on the soil type.
To illustrate the use of the computational orders in analyzing and comparing iterative methods for Richards' equation, we solve (4) closed by (5) and (6), with parameters from [10]. We consider the domain Ω = Ω 1 ∪ Ω 2 Ω = Ω 1 ∪ Ω 2 Omega=Omega_(1)uuOmega_(2)\Omega=\Omega_{1} \cup \Omega_{2}Ω=Ω1∪Ω2, where Ω 1 = [ 0 , 1 ] × [ 0 , 1 / 4 ) Ω 1 = [ 0 , 1 ] × [ 0 , 1 / 4 ) Omega_(1)=[0,1]xx[0,1//4)\Omega_{1}=[0,1] \times [0,1 / 4)Ω1=[0,1]×[0,1/4) and Ω 2 = [ 0 , 1 ] × [ 1 / 4 , 1 ] Ω 2 = [ 0 , 1 ] × [ 1 / 4 , 1 ] Omega_(2)=[0,1]xx[1//4,1]\Omega_{2}=[0,1] \times[1 / 4,1]Ω2=[0,1]×[1/4,1] and the time interval [ 0 , T ] [ 0 , T ] [0,T][0, T][0,T], where T = 0.003 T = 0.003 T=0.003T=0.003T=0.003. Sequences of corrections are stored at T / 3 , 2 T / 3 T / 3 , 2 T / 3 T//3,2T//3T / 3,2 T / 3T/3,2T/3, and T T TTT.
As an initial condition, we choose the pressure head
ψ 0 ( x , z ) = { − z + 1 / 4 ( x , z ) ∈ Ω 1 − 3 ( x , z ) ∈ Ω 2 ψ 0 ( x , z ) = − z + 1 / 4      ( x , z ) ∈ Ω 1 − 3      ( x , z ) ∈ Ω 2 psi^(0)(x,z)={[-z+1//4,(x","z)inOmega_(1)],[-3,(x","z)inOmega_(2)]:}\psi^{0}(x, z)= \begin{cases}-z+1 / 4 & (x, z) \in \Omega_{1} \\ -3 & (x, z) \in \Omega_{2}\end{cases}ψ0(x,z)={−z+1/4(x,z)∈Ω1−3(x,z)∈Ω2
On the top boundary, we employ constant Dirichlet conditions, equal to the initial condition at all times. On the remaining boundary no-flow conditions are used. We also consider the following source term
f ( x , z ) = { 0 ( x , z ) ∈ Ω 1 0.006 cos ⁡ ( 4 3 π ( z − 1 ) sin ⁡ ( 2 π x ) ) ( x , z ) ∈ Ω 2 f ( x , z ) = 0      ( x , z ) ∈ Ω 1 0.006 cos ⁡ 4 3 π ( z − 1 ) sin ⁡ ( 2 π x )      ( x , z ) ∈ Ω 2 f(x,z)={[0,(x","z)inOmega_(1)],[0.006 cos((4)/(3)pi(z-1)sin(2pi x)),(x","z)inOmega_(2)]:}f(x, z)= \begin{cases}0 & (x, z) \in \Omega_{1} \\ 0.006 \cos \left(\frac{4}{3} \pi(z-1) \sin (2 \pi x)\right) & (x, z) \in \Omega_{2}\end{cases}f(x,z)={0(x,z)∈Ω10.006cos⁡(43π(z−1)sin⁡(2πx))(x,z)∈Ω2

3 Linearization methods

Newton's method and the L L LLL-scheme employed in the Galerkin FEM are reviewed in Sect. 3.1 and a recently introduced L L LLL-scheme for the FDM is shortly described in Sec. 3.2. Theoretical convergence results are presented in Sect. 3.3.

3.1 Implicit FEM approach

The numerical scheme for fully discrete Richards' equation with backward Euler time discretization and Galerkin FEM reads: given ψ h , k − 1 ψ h , k − 1 psi_(h,k-1)\psi_{h, k-1}ψh,k−1 find ψ h , k ψ h , k psi_(h,k)\psi_{h, k}ψh,k such that
(7) ⟨ θ ( ψ h , k ) − θ ( ψ h , k − 1 ) , v h ⟩ + Δ t ⟨ K ( θ ( ψ h , k ) ) ∇ ( ψ h , k + z ) , ∇ v h ⟩ = Δ t ⟨ f k , v h ⟩ , (7) θ ψ h , k − θ ψ h , k − 1 , v h + Δ t K θ ψ h , k ∇ ψ h , k + z , ∇ v h = Δ t f k , v h , {:(7)(:theta(psi_(h,k))-theta(psi_(h,k-1)),v_(h):)+Delta t(:K(theta(psi_(h,k)))grad(psi_(h,k)+z),gradv_(h):)=Delta t(:f_(k),v_(h):)",":}\begin{equation*} \left\langle\theta\left(\psi_{h, k}\right)-\theta\left(\psi_{h, k-1}\right), v_{h}\right\rangle+\Delta t\left\langle K\left(\theta\left(\psi_{h, k}\right)\right) \nabla\left(\psi_{h, k}+z\right), \nabla v_{h}\right\rangle=\Delta t\left\langle f_{k}, v_{h}\right\rangle, \tag{7} \end{equation*}(7)⟨θ(ψh,k)−θ(ψh,k−1),vh⟩+Δt⟨K(θ(ψh,k))∇(ψh,k+z),∇vh⟩=Δt⟨fk,vh⟩,
∀ v h ∈ V h ∀ v h ∈ V h AAv_(h)inV_(h)\forall v_{h} \in V_{h}∀vh∈Vh, where V h V h V_(h)V_{h}Vh is the Galerkin linear element space, ⟨ ⋅ , ⋅ ⟩ ⟨ ⋅ , ⋅ ⟩ (:*,*:)\langle\cdot, \cdot\rangle⟨⋅,⋅⟩ denotes the inner product, Δ t = T / K , K ∈ N Δ t = T / K , K ∈ N Delta t=T//K,KinN\Delta t=T / \mathcal{K}, \mathcal{K} \in \mathbb{N}Δt=T/K,K∈N, is the time step, t k = k Δ t t k = k Δ t t_(k)=k Delta tt_{k}=k \Delta ttk=kΔt, and k ∈ { 1 , … , K } k ∈ { 1 , … , K } k in{1,dots,K}k \in\{1, \ldots, \mathcal{K}\}k∈{1,…,K}.
Newton's scheme linearizes both the θ ( ⋅ ) θ ( ⋅ ) theta(*)\theta(\cdot)θ(⋅) and K ( ⋅ ) K ( ⋅ ) K(*)K(\cdot)K(⋅) nonlinearities in Eq. (7) as follows: given ψ h , k − 1 , ψ h , k s − 1 ∈ V h ψ h , k − 1 , ψ h , k s − 1 ∈ V h psi_(h,k-1),psi_(h,k)^(s-1)inV_(h)\psi_{h, k-1}, \psi_{h, k}^{s-1} \in V_{h}ψh,k−1,ψh,ks−1∈Vh find ψ h , k s ∈ V h ψ h , k s ∈ V h psi_(h,k)^(s)inV_(h)\psi_{h, k}^{s} \in V_{h}ψh,ks∈Vh such that
⟨ θ ′ ( ψ h , k s − 1 ) ( ψ h , k s − ψ h , k s − 1 ) , v h ⟩ + Δ t ⟨ K ′ ( θ ( ψ h , k s − 1 ) ) ∇ ( ψ h , k s + z ) ( ψ h , k s − ψ h , k s − 1 ) , ∇ v h ⟩ = Δ t ⟨ f k , v h ⟩ − ⟨ θ ( ψ h , k s − 1 ) − θ ( ψ h , k − 1 ) , v h ⟩ − Δ t ⟨ K ( θ ( ψ h , k s − 1 ) ) ∇ ( ψ h , k s + z ) , ∇ v h ⟩ . θ ′ ψ h , k s − 1 ψ h , k s − ψ h , k s − 1 , v h + Δ t K ′ θ ψ h , k s − 1 ∇ ψ h , k s + z ψ h , k s − ψ h , k s − 1 , ∇ v h = Δ t f k , v h − θ ψ h , k s − 1 − θ ψ h , k − 1 , v h − Δ t K θ ψ h , k s − 1 ∇ ψ h , k s + z , ∇ v h . {:[(:theta^(')(psi_(h,k)^(s-1))(psi_(h,k)^(s)-psi_(h,k)^(s-1)),v_(h):)+Delta t(:K^(')(theta(psi_(h,k)^(s-1)))grad(psi_(h,k)^(s)+z)(psi_(h,k)^(s)-psi_(h,k)^(s-1)),gradv_(h):)],[=Delta t(:f_(k),v_(h):)-(:theta(psi_(h,k)^(s-1))-theta(psi_(h,k-1)),v_(h):)-Delta t(:K(theta(psi_(h,k)^(s-1)))grad(psi_(h,k)^(s)+z),gradv_(h):).]:}\begin{aligned} & \left\langle\theta^{\prime}\left(\psi_{h, k}^{s-1}\right)\left(\psi_{h, k}^{s}-\psi_{h, k}^{s-1}\right), v_{h}\right\rangle+\Delta t\left\langle K^{\prime}\left(\theta\left(\psi_{h, k}^{s-1}\right)\right) \nabla\left(\psi_{h, k}^{s}+z\right)\left(\psi_{h, k}^{s}-\psi_{h, k}^{s-1}\right), \nabla v_{h}\right\rangle \\ & =\Delta t\left\langle f_{k}, v_{h}\right\rangle-\left\langle\theta\left(\psi_{h, k}^{s-1}\right)-\theta\left(\psi_{h, k-1}\right), v_{h}\right\rangle-\Delta t\left\langle K\left(\theta\left(\psi_{h, k}^{s-1}\right)\right) \nabla\left(\psi_{h, k}^{s}+z\right), \nabla v_{h}\right\rangle . \end{aligned}⟨θ′(ψh,ks−1)(ψh,ks−ψh,ks−1),vh⟩+Δt⟨K′(θ(ψh,ks−1))∇(ψh,ks+z)(ψh,ks−ψh,ks−1),∇vh⟩=Δt⟨fk,vh⟩−⟨θ(ψh,ks−1)−θ(ψh,k−1),vh⟩−Δt⟨K(θ(ψh,ks−1))∇(ψh,ks+z),∇vh⟩.
The L L LLL-scheme only solves the non-linearity of θ ( ⋅ ) θ ( ⋅ ) theta(*)\theta(\cdot)θ(⋅) by replacing the derivative θ ′ θ ′ theta^(')\theta^{\prime}θ′ with a positive constant L L LLL. The linearization problem reads: given ψ h , k − 1 , ψ h , k s − 1 ∈ V h ψ h , k − 1 , ψ h , k s − 1 ∈ V h psi_(h,k-1),psi_(h,k)^(s-1)inV_(h)\psi_{h, k-1}, \psi_{h, k}^{s-1} \in V_{h}ψh,k−1,ψh,ks−1∈Vh find ψ h , k s ∈ V h ψ h , k s ∈ V h psi_(h,k)^(s)inV_(h)\psi_{h, k}^{s} \in V_{h}ψh,ks∈Vh such that
L ⟨ ψ h , k s − ψ h , k s − 1 , v h ⟩ + Δ t ⟨ K ( θ ( ψ h , k s − 1 ) ) ∇ ( ψ h , k s + z ) , ∇ v h ⟩ = Δ t ⟨ f k , v h ⟩ − ⟨ θ ( ψ h , k s − 1 ) − θ ( ψ h , k − 1 ) , v h ⟩ L ψ h , k s − ψ h , k s − 1 , v h + Δ t K θ ψ h , k s − 1 ∇ ψ h , k s + z , ∇ v h = Δ t f k , v h − θ ψ h , k s − 1 − θ ψ h , k − 1 , v h {:[L(:psi_(h,k)^(s)-psi_(h,k)^(s-1),v_(h):)+Delta t(:K(theta(psi_(h,k)^(s-1)))grad(psi_(h,k)^(s)+z),gradv_(h):)],[quad=Delta t(:f_(k),v_(h):)-(:theta(psi_(h,k)^(s-1))-theta(psi_(h,k-1)),v_(h):)]:}\begin{aligned} & L\left\langle\psi_{h, k}^{s}-\psi_{h, k}^{s-1}, v_{h}\right\rangle+\Delta t\left\langle K\left(\theta\left(\psi_{h, k}^{s-1}\right)\right) \nabla\left(\psi_{h, k}^{s}+z\right), \nabla v_{h}\right\rangle \\ & \quad=\Delta t\left\langle f_{k}, v_{h}\right\rangle-\left\langle\theta\left(\psi_{h, k}^{s-1}\right)-\theta\left(\psi_{h, k-1}\right), v_{h}\right\rangle \end{aligned}L⟨ψh,ks−ψh,ks−1,vh⟩+Δt⟨K(θ(ψh,ks−1))∇(ψh,ks+z),∇vh⟩=Δt⟨fk,vh⟩−⟨θ(ψh,ks−1)−θ(ψh,k−1),vh⟩
In both schemes, the iterations start with the solution at the last time step, ψ h , k s = 0 = ψ h , k − 1 ψ h , k s = 0 = ψ h , k − 1 psi_(h,k)^(s=0)=psi_(h,k-1)\psi_{h, k}^{s=0}= \psi_{h, k-1}ψh,ks=0=ψh,k−1 and are performed until the corrections to the vector ψ k s ψ k s psi_(k)^(s)\boldsymbol{\psi}_{k}^{s}ψks (with components ψ h , k s ψ h , k s psi_(h,k)^(s)\psi_{h, k}^{s}ψh,ks in a basis of V h V h V_(h)V_{h}Vh ) reach the accuracy ‖ ψ k s − ψ k s − 1 ‖ ≤ ε ψ k s − ψ k s − 1 ≤ ε ||psi_(k)^(s)-psi_(k)^(s-1)|| <= epsi\left\|\boldsymbol{\psi}_{k}^{s}-\boldsymbol{\psi}_{k}^{s-1}\right\| \leq \varepsilon‖ψks−ψks−1‖≤ε for an appropriate norm and some given tolerance ε ε epsi\varepsilonε. We notice that while Newton's method is 2nd order but only locally convergent, the L L LLL-scheme is a 1st order globally convergent method [5].
Different relaxation strategies can improve the efficiency of Newton's and L L LLL schemes. In particular, with the AA the current approximation is given by a linear combination of previous approximations with coefficients that minimize the weighted average of the successive corrections [1].

3.2 Explicit FDM approach

Explicit L L LLL-schemes were introduced in [8]. For brevity and clarity, we consider the simpler case of the one-dimensional Richards' equation,
(8) ∂ θ ( ψ ) ∂ t − ∂ ∂ z [ K ( θ ( ψ ) ) ∂ ∂ z ( ψ + z ) ] = f , (8) ∂ θ ( ψ ) ∂ t − ∂ ∂ z K ( θ ( ψ ) ) ∂ ∂ z ( ψ + z ) = f , {:(8)(del theta(psi))/(del t)-(del)/(del z)[K(theta(psi))(del)/(del z)(psi+z)]=f",":}\begin{equation*} \frac{\partial \theta(\psi)}{\partial t}-\frac{\partial}{\partial z}\left[K(\theta(\psi)) \frac{\partial}{\partial z}(\psi+z)\right]=f, \tag{8} \end{equation*}(8)∂θ(ψ)∂t−∂∂z[K(θ(ψ))∂∂z(ψ+z)]=f,
where z ∈ [ 0 , Z ] , t ∈ [ 0 , T ] z ∈ [ 0 , Z ] , t ∈ [ 0 , T ] z in[0,Z],t in[0,T]z \in[0, Z], t \in[0, T]z∈[0,Z],t∈[0,T]. The solution of (8) at times t = k Δ t , t ≤ T , k ∈ N t = k Δ t , t ≤ T , k ∈ N t=k Delta t,t <= T,k inNt=k \Delta t, t \leq T, k \in \mathbb{N}t=kΔt,t≤T,k∈N, and positions z = i Δ z , i = 1 , … , Z / Δ z z = i Δ z , i = 1 , … , Z / Δ z z=i Delta z,i=1,dots,Z//Delta zz=i \Delta z, i=1, \ldots, Z / \Delta zz=iΔz,i=1,…,Z/Δz is first approximated by an implicit FDM scheme
with backward discretization in time,
θ ( ψ i , k ) − θ ( ψ i , k − 1 ) = Δ t Δ z 2 { [ K ( ψ i + 1 / 2 , k ) ( ψ i + 1 , k − ψ i , k ) − K ( ψ i − 1 / 2 , k ) ( ψ i , k − ψ i − 1 , k ) ] (9) + ( K ( ψ i + 1 / 2 , k ) − K ( ψ i − 1 / 2 , k ) ) Δ z } + Δ t f θ ψ i , k − θ ψ i , k − 1 = Δ t Δ z 2 K ψ i + 1 / 2 , k ψ i + 1 , k − ψ i , k − K ψ i − 1 / 2 , k ψ i , k − ψ i − 1 , k (9) + K ψ i + 1 / 2 , k − K ψ i − 1 / 2 , k Δ z + Δ t f {:[theta(psi_(i,k))-theta(psi_(i,k-1))=(Delta t)/(Deltaz^(2)){[K(psi_(i+1//2,k))(psi_(i+1,k)-psi_(i,k))-K(psi_(i-1//2,k))(psi_(i,k)-psi_(i-1,k))]:}],[(9){:+(K(psi_(i+1//2,k))-K(psi_(i-1//2,k)))Delta z}+Delta tf]:}\begin{align*} \theta\left(\psi_{i, k}\right)-\theta\left(\psi_{i, k-1}\right) & =\frac{\Delta t}{\Delta z^{2}}\left\{\left[K\left(\psi_{i+1 / 2, k}\right)\left(\psi_{i+1, k}-\psi_{i, k}\right)-K\left(\psi_{i-1 / 2, k}\right)\left(\psi_{i, k}-\psi_{i-1, k}\right)\right]\right. \\ & \left.+\left(K\left(\psi_{i+1 / 2, k}\right)-K\left(\psi_{i-1 / 2, k}\right)\right) \Delta z\right\}+\Delta t f \tag{9} \end{align*}θ(ψi,k)−θ(ψi,k−1)=ΔtΔz2{[K(ψi+1/2,k)(ψi+1,k−ψi,k)−K(ψi−1/2,k)(ψi,k−ψi−1,k)](9)+(K(ψi+1/2,k)−K(ψi−1/2,k))Δz}+Δtf
Next, we assign to ψ i , k ψ i , k psi_(i,k)\psi_{i, k}ψi,k the iteration count s s sss, we add the stabilization term L ( ψ i , k s + 1 − ψ i , k s L ψ i , k s + 1 − ψ i , k s L(psi_(i,k)^(s+1)-:}psi_(i,k)^(s)L\left(\psi_{i, k}^{s+1}-\right. \psi_{i, k}^{s}L(ψi,ks+1−ψi,ks ) to the left-hand side of (9) and perform iterations s = 1 , 2 , … s = 1 , 2 , … s=1,2,dotss=1,2, \ldotss=1,2,… starting with ψ i , k 1 = ψ i , k − 1 ψ i , k 1 = ψ i , k − 1 psi_(i,k)^(1)=psi_(i,k-1)\psi_{i, k}^{1}=\psi_{i, k-1}ψi,k1=ψi,k−1 until the corrections to the solution ψ k s = ( ψ i , k s , … , ψ Z / Δ z , k s ) ψ k s = ψ i , k s , … , ψ Z / Δ z , k s psi_(k)^(s)=(psi_(i,k)^(s),dots,psi_(Z//Delta z,k)^(s))\psi_{k}^{s}=\left(\psi_{i, k}^{s}, \ldots, \psi_{Z / \Delta z, k}^{s}\right)ψks=(ψi,ks,…,ψZ/Δz,ks) fall under a given tolerance, ‖ ψ k s + 1 − ψ k s ‖ ≤ ε ψ k s + 1 − ψ k s ≤ ε ||psi_(k)^(s+1)-psi_(k)^(s)|| <= epsi\left\|\boldsymbol{\psi}_{k}^{s+1}-\boldsymbol{\psi}_{k}^{s}\right\| \leq \varepsilon‖ψks+1−ψks‖≤ε. The result is the following explicit L L LLL-scheme:
(10) ψ i , k s + 1 = [ 1 − ( r i + 1 / 2 , k s + r i − 1 / 2 , k s ) ] ψ i , k s + r i + 1 / 2 , k s ψ i + 1 , k s + r i − 1 / 2 , k s ψ i − 1 , k s + f s (10) ψ i , k s + 1 = 1 − r i + 1 / 2 , k s + r i − 1 / 2 , k s ψ i , k s + r i + 1 / 2 , k s ψ i + 1 , k s + r i − 1 / 2 , k s ψ i − 1 , k s + f s {:(10)psi_(i,k)^(s+1)=[1-(r_(i+1//2,k)^(s)+r_(i-1//2,k)^(s))]psi_(i,k)^(s)+r_(i+1//2,k)^(s)psi_(i+1,k)^(s)+r_(i-1//2,k)^(s)psi_(i-1,k)^(s)+f^(s):}\begin{equation*} \psi_{i, k}^{s+1}=\left[1-\left(r_{i+1 / 2, k}^{s}+r_{i-1 / 2, k}^{s}\right)\right] \psi_{i, k}^{s}+r_{i+1 / 2, k}^{s} \psi_{i+1, k}^{s}+r_{i-1 / 2, k}^{s} \psi_{i-1, k}^{s}+f^{s} \tag{10} \end{equation*}(10)ψi,ks+1=[1−(ri+1/2,ks+ri−1/2,ks)]ψi,ks+ri+1/2,ksψi+1,ks+ri−1/2,ksψi−1,ks+fs
where
r i ± 1 / 2 , k s = K ( ψ i ± 1 / 2 , k s ) Δ t / ( L Δ z 2 ) , r i ± 1 / 2 , k s ≤ 1 / 2 f s = Δ t f / L + ( r i + 1 / 2 , k s − r i − 1 / 2 , k s ) Δ z − [ θ ( ψ i , k s ) − θ ( ψ i , k − 1 ) ] / L r i ± 1 / 2 , k s = K ψ i ± 1 / 2 , k s Δ t / L Δ z 2 , r i ± 1 / 2 , k s ≤ 1 / 2 f s = Δ t f / L + r i + 1 / 2 , k s − r i − 1 / 2 , k s Δ z − θ ψ i , k s − θ ψ i , k − 1 / L {:[r_(i+-1//2,k)^(s)=K(psi_(i+-1//2,k)^(s))Delta t//(L Deltaz^(2))","quadr_(i+-1//2,k)^(s) <= 1//2],[f^(s)=Delta tf//L+(r_(i+1//2,k)^(s)-r_(i-1//2,k)^(s))Delta z-[theta(psi_(i,k)^(s))-theta(psi_(i,k-1))]//L]:}\begin{aligned} & r_{i \pm 1 / 2, k}^{s}=K\left(\psi_{i \pm 1 / 2, k}^{s}\right) \Delta t /\left(L \Delta z^{2}\right), \quad r_{i \pm 1 / 2, k}^{s} \leq 1 / 2 \\ & f^{s}=\Delta t f / L+\left(r_{i+1 / 2, k}^{s}-r_{i-1 / 2, k}^{s}\right) \Delta z-\left[\theta\left(\psi_{i, k}^{s}\right)-\theta\left(\psi_{i, k-1}\right)\right] / L \end{aligned}ri±1/2,ks=K(ψi±1/2,ks)Δt/(LΔz2),ri±1/2,ks≤1/2fs=Δtf/L+(ri+1/2,ks−ri−1/2,ks)Δz−[θ(ψi,ks)−θ(ψi,k−1)]/L

3.3 Theoretical convergence results

The FEM approaches benefit from well-established theoretical results. For instance, the quadratic convergence of the Newton's method can be proved by using the Kirchhoff transformation and a regularization step [7]. The linear convergence of the L L LLL-scheme is ensured by the monotonically increasing θ θ theta\thetaθ, Lipschitz continuity of θ θ theta\thetaθ and K K KKK, bounded ∇ ψ h , k ∇ ψ h , k gradpsi_(h,k)\nabla \psi_{h, k}∇ψh,k, and the condition L ≥ L θ = sup ψ { θ ′ ( ψ ) } L ≥ L θ = sup ψ   θ ′ ( ψ ) L >= L_(theta)=s u p_(psi){theta^(')(psi)}L \geq L_{\theta}=\sup _{\psi}\left\{\theta^{\prime}(\psi)\right\}L≥Lθ=supψ{θ′(ψ)} [5].
Similar results for the explicit FDM L L LLL-scheme are not yet available. However, we can look for convergence conditions by investigating the evolution of the difference e s = ψ s − ψ e s = ψ s − ψ e^(s)=psi^(s)-psie^{s}=\psi^{s}-\psies=ψs−ψ of the solutions of (10) and (9) at fixed position i i iii and time k k kkk. After denoting by Δ t Δ z 2 D ( ψ s ) Δ t Δ z 2 D ψ s (Delta t)/(Deltaz^(2))D(psi^(s))\frac{\Delta t}{\Delta z^{2}} \mathcal{D}\left(\psi^{s}\right)ΔtΔz2D(ψs) the r.h.s. of (9) and subtracting (9) from (10), we have
L ( e s + 1 − e s ) + θ ( ψ s ) − θ ( ψ ) = Δ t Δ z 2 [ D ( ψ s ) − D ( ψ ) ] L e s + 1 − e s + θ ψ s − θ ( ψ ) = Δ t Δ z 2 D ψ s − D ( ψ ) L(e^(s+1)-e^(s))+theta(psi^(s))-theta(psi)=(Delta t)/(Deltaz^(2))[D(psi^(s))-D(psi)]L\left(e^{s+1}-e^{s}\right)+\theta\left(\psi^{s}\right)-\theta(\psi)=\frac{\Delta t}{\Delta z^{2}}\left[\mathcal{D}\left(\psi^{s}\right)-\mathcal{D}(\psi)\right]L(es+1−es)+θ(ψs)−θ(ψ)=ΔtΔz2[D(ψs)−D(ψ)]
Assuming (A1) θ ′ ( ⋅ ) > 0 θ ′ ( ⋅ ) > 0 theta^(')(*) > 0\theta^{\prime}(\cdot)>0θ′(⋅)>0 and (A2) | D ( ψ s ) − D ( ψ ) | ≤ L D | ψ s − ψ | = L D | e s | D ψ s − D ( ψ ) ≤ L D ψ s − ψ = L D e s |D(psi^(s))-D(psi)| <= L_(D)|psi^(s)-psi|=L_(D)|e^(s)|\left|\mathcal{D}\left(\psi^{s}\right)-\mathcal{D}(\psi)\right| \leq L_{\mathcal{D}}\left|\psi^{s}-\psi\right|=L_{\mathcal{D}}\left|e^{s}\right||D(ψs)−D(ψ)|≤LD|ψs−ψ|=LD|es| (strictly monotonic θ θ theta\thetaθ and Lipschitz continuous D D D\mathcal{D}D ) one obtains
| e s + 1 | ≤ | e s | | 1 − θ ′ ( ψ ) L + L D Δ t L Δ z 2 | e s + 1 ≤ e s 1 − θ ′ ( ψ ) L + L D Δ t L Δ z 2 |e^(s+1)| <= |e^(s)||1-(theta^(')(psi))/(L)+(L_(D)Delta t)/(L Deltaz^(2))|\left|e^{s+1}\right| \leq\left|e^{s}\right|\left|1-\frac{\theta^{\prime}(\psi)}{L}+\frac{L_{\mathcal{D}} \Delta t}{L \Delta z^{2}}\right||es+1|≤|es||1−θ′(ψ)L+LDΔtLΔz2|
Assuming the condition inf ψ { θ ′ ( ψ ) } > L D Δ t / ( Δ z 2 ) inf ψ   θ ′ ( ψ ) > L D Δ t / Δ z 2 i n f_(psi){theta^(')(psi)} > L_(D)Delta t//(Deltaz^(2))\inf _{\psi}\left\{\theta^{\prime}(\psi)\right\}>L_{\mathcal{D}} \Delta t /\left(\Delta z^{2}\right)infψ{θ′(ψ)}>LDΔt/(Δz2), the above relation becomes a contraction. This implies | e s + 1 | < | e s | e s + 1 < e s |e^(s+1)| < |e^(s)|\left|e^{s+1}\right|<\left|e^{s}\right||es+1|<|es|, which in turn implies ‖ ψ k s + 1 − ψ k ‖ < ‖ ψ k s − ψ k ‖ ‖ ψ k s + 1 − ψ k ‖ < ‖ ψ k s − ψ k ‖ ||psi_(k)^(s+1)-psi_(k)|| < ||psi_(k)^(s)-psi_(k)||\| \boldsymbol{\psi}_{k}^{s+1}- \boldsymbol{\psi}_{k}\|<\| \boldsymbol{\psi}_{k}^{s}-\boldsymbol{\psi}_{k} \|‖ψks+1−ψk‖<‖ψks−ψk‖. The condition above has to be compatible with the condition
K ( ψ s ) Δ t / ( L Δ z 2 ) ≤ 1 / 2 K ψ s Δ t / L Δ z 2 ≤ 1 / 2 K(psi^(s))Delta t//(L Deltaz^(2)) <= 1//2K\left(\psi^{s}\right) \Delta t /\left(L \Delta z^{2}\right) \leq 1 / 2K(ψs)Δt/(LΔz2)≤1/2 which constrains the L L LLL-scheme (10). These are necessary conditions for the C C CCC-linear convergence without being however sufficient conditions, i.e. they do not imply the existence of the limit (1).

4 Numerical results

Computational orders of convergence are assessed by investigating the behavior of the sequence { x s } x s {x_(s)}\left\{x_{s}\right\}{xs} of successive corrections x s = ‖ ψ s − ψ s − 1 ‖ x s = ψ s − ψ s − 1 x_(s)=||psi^(s)-psi^(s-1)||x_{s}=\left\|\psi^{s}-\psi^{s-1}\right\|xs=‖ψs−ψs−1‖, which are used to assess the accuracy of the approximation in the general case when the exact solution of the problem is not available. If the linearization method converges x ∗ = lim s → ∞ x s = 0 x ∗ = lim s → ∞   x s = 0 x^(**)=lim_(s rarr oo)x_(s)=0x^{*}=\lim _{s \rightarrow \infty} x_{s}=0x∗=lims→∞xs=0 and e s = | x ∗ − x s | = x s e s = x ∗ − x s = x s e_(s)=|x^(**)-x_(s)|=x_(s)e_{s}=\left|x^{*}-x_{s}\right|=x_{s}es=|x∗−xs|=xs. The error-based orders C , Q C , Q C,QC, QC,Q, and R R RRR of the sequence { x s } x s {x_(s)}\left\{x_{s}\right\}{xs} computed according to (1-3), if they exist 2 2 ^(2){ }^{2}2, are in fact computational orders C ′ , Q ′ C ′ , Q ′ C^('),Q^(')C^{\prime}, Q^{\prime}C′,Q′, and R ′ R ′ R^(')R^{\prime}R′ for the sequence { ψ s } ψ s {psi^(s)}\left\{\psi^{s}\right\}{ψs}. To proceed, we use the definitions (2) and (3) to assess a convergence order p p ppp, which is further used in (1) to check whether the classical C C CCC convergence of the sequence { x s } x s {x_(s)}\left\{x_{s}\right\}{xs} takes place. The linearization methods presented in Sect. 3 are tested in Sects. 4.1 and 4.2 on the benchmark problem formulated in Sec. 2 with a fixed tolerance ε = 10 − 7 ε = 10 − 7 epsi=10^(-7)\varepsilon=10^{-7}ε=10−7. In Sects. 4.2 we also assess error-based orders of convergence for a test problem with known solution modeling variably saturated flow fully coupled with surfactant transport.

4.1 FEM Newton's method and L L L\boldsymbol{L}L-scheme

Newton's method converges after only a few iterations (Fig. 1). The estimate p ( s ) = ln ⁡ e s + 1 / ln ⁡ e s ≈ 2 p ( s ) = ln ⁡ e s + 1 / ln ⁡ e s ≈ 2 p(s)=ln e_(s+1)//ln e_(s)~~2p(s)=\ln e_{s+1} / \ln e_{s} \approx 2p(s)=ln⁡es+1/ln⁡es≈2 shown in Fig. 2 suggests the Q Q QQQ-quadratic convergence of the sequence { x s } x s {x_(s)}\left\{x_{s}\right\}{xs}, according to definition (2). The sequence of successive approximations { ψ s } ψ s {psi^(s)}\left\{\psi^{s}\right\}{ψs} of the Newton's method is then Q ′ Q ′ Q^(')Q^{\prime}Q′-quadratic convergent.
The estimate p ( s ) = | ln ⁡ e s | 1 / s ≈ 2 p ( s ) = ln ⁡ e s 1 / s ≈ 2 p(s)=|ln e_(s)|^(1//s)~~2p(s)=\left|\ln e_{s}\right|^{1 / s} \approx 2p(s)=|ln⁡es|1/s≈2 presented in Fig. 3 also indicates the R ′ R ′ R^(')R^{\prime}R′ quadratic convergence according to definition (3), sustaining the hypothesis of the Q ′ Q ′ Q^(')Q^{\prime}Q′-quadratic order. The finite values Q 2 ( s ) = e s + 1 / e s 2 Q 2 ( s ) = e s + 1 / e s 2 Q_(2)(s)=e_(s+1)//e_(s)^(2)Q_{2}(s)=e_{s+1} / e_{s}^{2}Q2(s)=es+1/es2 shown in Fig. 4 are consistent with the definition (1) and we may expect the C ′ C ′ C^(')C^{\prime}C′-quadratic convergence. Provided that the equivalence of computational and error-based orders for p > 1 p > 1 p > 1p>1p>1 in R n R n R^(n)\mathbb{R}^{n}Rn [4] holds in the normed space of the solutions, the sequence { ψ s } ψ s {psi^(s)}\left\{\psi^{s}\right\}{ψs} also possesses C , Q C , Q C,QC, QC,Q, and R R RRR quadratic orders. Longer sequences necessary to validate these results could only be obtained by doing multi-precision computations, because the e s e s e_(s)e_{s}es values from Fig. 1 are already close to the machine epsilon in double precision.
Similarly, using (2) and (3) we first verified the Q ′ Q ′ Q^(')Q^{\prime}Q′ - and R ′ R ′ R^(')R^{\prime}R′-convergence and we identified the order p = 1 p = 1 p=1p=1p=1 for the sequence { ψ s } ψ s {psi^(s)}\left\{\psi^{s}\right\}{ψs} of solutions provided by the FEM L L LLL-scheme. The behavior of the quotient Q 1 ( s ) = e s + 1 / e s Q 1 ( s ) = e s + 1 / e s Q_(1)(s)=e_(s+1)//e_(s)Q_{1}(s)=e_{s+1} / e_{s}Q1(s)=es+1/es shown in Fig. 5 provides
Fig. 1 FEM: Norms of successive corrections provided by Newton's method.
Fig. 3 FEM: Estimation of R R RRR-order 2 for Newton's method.
Fig. 2 FEM: Estimation of Q Q QQQ-order p p ppp for the convergence of Newton's method.
Fig. 4 FEM: Estimation of the quotient Q 2 Q 2 Q_(2)Q_{2}Q2 from definition (1) for Newton's method.
a strong indication of the C ′ C ′ C^(')C^{\prime}C′-linear convergence according to the definition (1) with p = 1 p = 1 p=1p=1p=1. The scheme benefits from the use of AA which, for the setup of Sect. 2, halves the number of iterations and renders the convergence two times faster (Fig. 6).
Fig. 5 FEM: Estimation of the quotient Q 1 Q 1 Q_(1)Q_{1}Q1 from definition (1) for L L LLL-scheme ( L = 0.15 L = 0.15 L=0.15L=0.15L=0.15 ).
Fig. 6 FEM: Estimation of the quotient Q 1 Q 1 Q_(1)Q_{1}Q1 from definition (1) for L L LLL-scheme with AA.

4.2 FDM L L LLL-scheme

The number of iterations for the explicit FDM L L LLL-scheme (Fig. 7) is one order of magnitude larger than for the implicit FEM L L LLL-scheme without AA. However, the computing time is generally one order of magnitude smaller than for implicit schemes in solving similar problems (see [8]). The Q ′ Q ′ Q^(')Q^{\prime}Q′ - and R ′ R ′ R^(')R^{\prime}R′-linear convergence and the order p = 1 p = 1 p=1p=1p=1 are assessed in the same way as for the FEM L L LLL-scheme. Estimates of Q 1 Q 1 Q_(1)Q_{1}Q1 shown in Fig. 8 indicate the C ′ C ′ C^(')C^{\prime}C′-linear convergence. Applying AA to the explicit L L LLL-scheme only results in increasing the computing time, without any significant improvement of the explicit FDM L L LLL-scheme.
Fig. 7 FDM: Norms of successive corrections provided by the L L LLL-scheme ( L = 0.5 L = 0.5 L=0.5L=0.5L=0.5 ).
Fig. 8 FDM: Estimation of the quotient Q 1 Q 1 Q_(1)Q_{1}Q1 from definition (1) for L L LLL-scheme.
Error-based orders of convergence for linearly convergent L L LLL-schemes also can be inferred when the exact solution of the problem is known. As an example, we consider a test problem for fully coupled flow and surfactant transport with the manufactured exact one-dimensional solution ψ m ( z , t ) = − t z ( 1 − z ) + z / 4 ψ m ( z , t ) = − t z ( 1 − z ) + z / 4 psi_(m)(z,t)=-tz(1-z)+z//4\psi_{m}(z, t)=-t z(1-z)+z / 4ψm(z,t)=−tz(1−z)+z/4, for the pressure, and c m ( z , t ) = t z ( 1 − z ) + 1 c m ( z , t ) = t z ( 1 − z ) + 1 c_(m)(z,t)=tz(1-z)+1c_{m}(z, t)=t z(1-z)+1cm(z,t)=tz(1−z)+1, for the concentration, with θ ( ψ ) = 1 / ( 1 − ψ − c / 10 ) θ ( ψ ) = 1 / ( 1 − ψ − c / 10 ) theta(psi)=1//(1-psi-c//10)\theta(\psi)=1 /(1-\psi-c / 10)θ(ψ)=1/(1−ψ−c/10) and K ( θ ( ψ ) ) = ψ 2 K ( θ ( ψ ) ) = ψ 2 K(theta(psi))=psi^(2)K(\theta(\psi))=\psi^{2}K(θ(ψ))=ψ2 (the one-dimensional version of the problem considerd in [8, Sect. 5.2.1]). The problem is solved by utilising the FDM L L LLL-scheme for Richards' equation coupled with an explicit global random walk L L LLL-scheme for the advection-diffusion equation. The parameter L = 100 L = 100 L=100L=100L=100 is the same for both schemes and the stopping criterion uses errors instead of corrections 3 , max ( ‖ ψ k s − ψ m , k ‖ , ‖ c k s − c m , k ‖ ) ≤ 10 − 3 3 , max ψ k s − ψ m , k , c k s − c m , k ≤ 10 − 3 ^(3),max(||psi_(k)^(s)-psi_(m,k)||,||c_(k)^(s)-c_(m,k)||) <= 10^(-3){ }^{3}, \max \left(\left\|\psi_{k}^{s}-\psi_{m, k}\right\|,\left\|c_{k}^{s}-c_{m, k}\right\|\right) \leq 10^{-3}3,max(‖ψks−ψm,k‖,‖cks−cm,k‖)≤10−3.
The quotient Q 1 ( s ) Q 1 ( s ) Q_(1)(s)Q_{1}(s)Q1(s) computed with errors of the sequence { ψ s } ψ s {psi^(s)}\left\{\psi^{s}\right\}{ψs} shows a tendency towards the C C CCC-sublinear convergence, i.e., Q 1 → 1 Q 1 → 1 Q_(1)rarr1Q_{1} \rightarrow 1Q1→1 (Fig. 9). The evolution of the corresponding quotient of corrections, with Q 1 < 1 Q 1 < 1 Q_(1) < 1Q_{1}<1Q1<1, indicates the C ′ C ′ C^(')C^{\prime}C′-linear convergence (Fig. 10). The sequence of concentrations, { c s } c s {c^(s)}\left\{c^{s}\right\}{cs}, behaves similarly, i.e. C C CCC-sublinear convergence of the errors and C ′ C ′ C^(')C^{\prime}C′-linear convergence of the corrections (Figs. 11 and12).
Fig. 9 FDM: Quotient Q 1 Q 1 Q_(1)Q_{1}Q1 for L L LLL-scheme estimated with errors ‖ ψ k s − ψ m , k ‖ ψ k s − ψ m , k ||psi_(k)^(s)-psi_(m,k)||\left\|\psi_{k}^{s}-\psi_{m, k}\right\|‖ψks−ψm,k‖.
Fig. 11 FDM: Quotient Q 1 Q 1 Q_(1)Q_{1}Q1 for L L LLL-scheme estimated with errors ‖ c k s − c m , k ‖ c k s − c m , k ||c_(k)^(s)-c_(m,k)||\left\|c_{k}^{s}-c_{m, k}\right\|‖cks−cm,k‖.
Fig. 10 FDM: Quotient Q 1 Q 1 Q_(1)Q_{1}Q1 for L L LLL-scheme estimated with corrections ‖ ψ k s − ψ k s − 1 ‖ ψ k s − ψ k s − 1 ||psi_(k)^(s)-psi_(k)^(s-1)||\left\|\psi_{k}^{s}-\psi_{k}^{s-1}\right\|‖ψks−ψks−1‖.
Fig. 12 FDM: Quotient Q 1 Q 1 Q_(1)Q_{1}Q1 for L L LLL-scheme estimated with corrections ‖ c k s − c k s − 1 ‖ c k s − c k s − 1 ||c_(k)^(s)-c_(k)^(s-1)||\left\|c_{k}^{s}-c_{k}^{s-1}\right\|‖cks−cks−1‖.
Thanks to the manufactured solution, we also assessed the estimated order of convergence (EOC) of the numerical FDM solution towards the exact solution ( ψ m , c m ψ m , c m psi_(m),c_(m)\psi_{m}, c_{m}ψm,cm ), by successively halving Δ z Δ z Delta z\Delta zΔz, with max ( ‖ ψ k s − ψ k s − 1 ‖ , ‖ c k s − c k s − 1 ‖ ) ≤ 10 − 6 max ψ k s − ψ k s − 1 , c k s − c k s − 1 ≤ 10 − 6 max(||psi_(k)^(s)-psi_(k)^(s-1)||,||c_(k)^(s)-c_(k)^(s-1)||) <= 10^(-6)\max \left(\left\|\psi_{k}^{s}-\psi_{k}^{s-1}\right\|,\left\|c_{k}^{s}-c_{k}^{s-1}\right\|\right) \leq 10^{-6}max(‖ψks−ψks−1‖,‖cks−cks−1‖)≤10−6 as stopping criterion for the iterations of the coupled L L LLL-schemes (Table 1).
Δ z Δ z Delta z\Delta zΔz ‖ ψ − ψ m ‖ ψ − ψ m ||psi-psi_(m)||\left\|\psi-\psi_{m}\right\|‖ψ−ψm‖ EOC ‖ c − c m ‖ c − c m ||c-c_(m)||\left\|c-c_{m}\right\|‖c−cm‖ EOC
1.00 e − 1 1.00 e − 1 1.00e-11.00 \mathrm{e}-11.00e−1 4.59 e − 02 4.59 e − 02 4.59e-024.59 \mathrm{e}-024.59e−02 - 7.64 e − 03 7.64 e − 03 7.64e-037.64 \mathrm{e}-037.64e−03 -
5.00 e − 1 5.00 e − 1 5.00e-15.00 \mathrm{e}-15.00e−1 1.14 e − 02 1.14 e − 02 1.14e-021.14 \mathrm{e}-021.14e−02 2.00 1.95 e − 03 1.95 e − 03 1.95e-031.95 \mathrm{e}-031.95e−03 1.97
2.50 e − 2 2.50 e − 2 2.50e-22.50 \mathrm{e}-22.50e−2 2.95 e − 03 2.95 e − 03 2.95e-032.95 \mathrm{e}-032.95e−03 1.95 4.38 e − 04 4.38 e − 04 4.38e-044.38 \mathrm{e}-044.38e−04 2.16
1.25 e − 2 1.25 e − 2 1.25e-21.25 \mathrm{e}-21.25e−2 9.97 e − 04 9.97 e − 04 9.97e-049.97 \mathrm{e}-049.97e−04 1.57 8.99 e − 05 8.99 e − 05 8.99e-058.99 \mathrm{e}-058.99e−05 2.28
Delta z ||psi-psi_(m)|| EOC ||c-c_(m)|| EOC 1.00e-1 4.59e-02 - 7.64e-03 - 5.00e-1 1.14e-02 2.00 1.95e-03 1.97 2.50e-2 2.95e-03 1.95 4.38e-04 2.16 1.25e-2 9.97e-04 1.57 8.99e-05 2.28| $\Delta z$ | $\left\\|\psi-\psi_{m}\right\\|$ | EOC | $\left\\|c-c_{m}\right\\|$ | EOC | | :---: | :---: | :---: | :---: | :---: | | $1.00 \mathrm{e}-1$ | $4.59 \mathrm{e}-02$ | - | $7.64 \mathrm{e}-03$ | - | | $5.00 \mathrm{e}-1$ | $1.14 \mathrm{e}-02$ | 2.00 | $1.95 \mathrm{e}-03$ | 1.97 | | $2.50 \mathrm{e}-2$ | $2.95 \mathrm{e}-03$ | 1.95 | $4.38 \mathrm{e}-04$ | 2.16 | | $1.25 \mathrm{e}-2$ | $9.97 \mathrm{e}-04$ | 1.57 | $8.99 \mathrm{e}-05$ | 2.28 |
Table 1 Orders of convergence of the numerical FDM solution.

5 Discussion

The convergence of the linearization methods for Richards' equation can be investigated with methods developed for sequences of real numbers. To do that, we follow the recipe: first, use the definition of the Q Q QQQ - and R R RRR-order convergence to identify the convergence order p p ppp, then, check whether the sequence { x s } = { ‖ ψ s − ψ s − 1 ‖ } x s = ψ s − ψ s − 1 {x_(s)}={||psi^(s)-psi^(s-1)||}\left\{x_{s}\right\}=\left\{\left\|\psi^{s}-\psi^{s-1}\right\|\right\}{xs}={‖ψs−ψs−1‖} converges with the classical C C CCC-order p p ppp.
We have to keep in mind that the error-based convergence orders C , Q C , Q C,QC, QC,Q, and R R RRR for the sequence of positive numbers { x s } x s {x_(s)}\left\{x_{s}\right\}{xs} correspond to the computational orders C ′ C ′ C^(')C^{\prime}C′, Q ′ Q ′ Q^(')Q^{\prime}Q′, and R ′ R ′ R^(')R^{\prime}R′ of the sequence of solutions { ψ s } ψ s {psi^(s)}\left\{\psi^{s}\right\}{ψs} provided by the linearization scheme. Error-based orders for the linearization scheme only can be inferred when exact (e.g., manufactured) solutions of the flow problem are available.
The numerical analysis of the iterative methods for nonlinear systems of equations involved in modeling processes in porous media presented in the literature aim at showing the convergence to zero of the error of the scheme. Numerically inferred computational or error-based orders of convergence complete the characterization of the convergence and may be useful in several ways: to verify the theoretically predicted convergence and to check the convergence in cases where the conditions of the convergence theorems are hard to meet, to evaluate the efficiency of the relaxation strategies, or to enable comparisons of the convergence speed for different linearization procedures.

References

  1. Anderson, D.G.: Iterative procedures for nonlinear integral equations. JACM 12(4), 547-560, (1965)
  2. Beyer, W. A., Ebanks, B. R., Qualls, C. R.: Convergence rates and convergence-order profiles for sequences. Acta Appl. Math. 20, 267-284, (1990)
  3. Cătinaş, E.: A survey on the high convergence orders and computational convergence orders of sequences. Appl. Math. Comput. 343, 1-20, (2019)
  4. Cătinaș, E.: How many steps still left to x ∗ x ∗ x^(**)x^{*}x∗ ? SIAM Rev. 63(3), 585-624, (2021)
  5. List, F., Radu, F.A.: A study on iterative methods for solving Richards' equation. Comput. Geosci. 20(2) 341-353, (2016)
  6. Potra, F. A.: On Q-order and R-order of convergence. J. Optim. Theory Appl. 63(3), 415-431, (1989)
  7. Radu, F.A., Pop, I.S., Knabner, P.: On the convergence of the Newton method for the mixed finite element discretization of a class of degenerate parabolic equations. In: Numerical Mathematics and Advanced Applications 42, 1192-1200, (2006)
  8. Suciu, N., Illiano, D., Prechtel, A., Radu, F.A.: Global random walk solvers for fully coupled flow and transport in saturated/unsaturated porous media. Adv. Water Resour. 152, 103935, (2021)
  9. Suciu, N., Radu, F.A.: Global random walk solvers for reactive transport and biodegradation processes in heterogeneous porous media. Adv. Water Res., 166, 104268, (2022)
  10. Stokke, J.S., Mitra, K., Storvik, E., Both, J.W., Radu, F.A.: An adaptive solution strategy for Richards' equation. Comput. Math. Appl. 152, 155-167, (2023)
  1. Nicolae Suciu ⋅ ⋅ *\cdot⋅ Emil Cătinaş ⋅ ⋅ *\cdot⋅ Andra Malina
    Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy, Cluj-Napoca, Romania, e-mail: nsuciu@ictp.acad.ro, ecatinas@ictp.acad.ro, andra.malina@ictp.acad.ro
    Florin A. Radu • Jakob S. Stokke
    Center for Modeling of Coupled Subsurface Dynamics, University of Bergen, Bergen, Norway, e-mail: florin.radu@uib.no,jakob.stokke@uib.no
  2. 1 1 ^(1){ }^{1}1 These are usually denoted by Q L Q L Q_(L)Q_{L}QL and R L R L R_(L)R_{L}RL to distinguish from other types of Q Q QQQ and R R RRR orders. Since only the definitions (2) and (3) will be used in the following, we disregard the subscript L L LLL.
  3. 2 2 ^(2){ }^{2}2 In a multidimensional setting the C C CCC orders of convergence are norm-dependent [3]. However, the numerical examples presented below indicate the existence of the limit from the definition (1).
  4. 3 3 ^(3){ }^{3}3 We note that once this stopping criterion for errors is satisfied, the corresponding corrections for both { ψ s } ψ s {psi^(s)}\left\{\psi^{s}\right\}{ψs} and { c s } c s {c^(s)}\left\{c^{s}\right\}{cs} are of the order 10 − 8 10 − 8 10^(-8)10^{-8}10−8.
2025

Related Posts