
==== Front
J Sci Comput
J Sci Comput
Journal of Scientific Computing
0885-7474
1573-7691
Springer US New York

2650
10.1007/s10915-024-02650-x
Article
Error Estimates and Adaptivity of the Space-Time Discontinuous Galerkin Method for Solving the Richards Equation
http://orcid.org/0000-0001-6356-934X
Dolejší Vít dolejsi@karlin.mff.cuni.cz

1
Shin Hyun-Geun 1
Vlasák Miloslav 2
1 https://ror.org/024d6js02 grid.4491.8 0000 0004 1937 116X Faculty of Mathematics and Physics, Charles University, Sokolovská 83, Prague, Czech Republic
2 grid.6652.7 0000000121738213 Faculty of Civil Engineering, Czech Technical University, Thakurova 7, Prague 6, 166 29 Czech Republic
20 8 2024
20 8 2024
2024
101 1 1115 1 2024
28 5 2024
28 7 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by/4.0/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
We present a higher-order space-time adaptive method for the numerical solution of the Richards equation that describes a flow motion through variably saturated media. The discretization is based on the space-time discontinuous Galerkin method, which provides high stability and accuracy and can naturally handle varying meshes. We derive reliable and efficient a posteriori error estimates in the residual-based norm. The estimates use well-balanced spatial and temporal flux reconstructions which are constructed locally over space-time elements or space-time patches. The accuracy of the estimates is verified by numerical experiments. Moreover, we develop the hp-adaptive method and demonstrate its efficiency and usefulness on a practically relevant example.

Keywords

Space-time discontinuous Galerkin method
Richards equation
A posteriori error estimate
hp-mesh adaptation
Mathematics Subject Classification

65M60
65M15
65M50
http://dx.doi.org/10.13039/501100001824 Grantová Agentura České Republiky 20-01074S Dolejší Vít http://dx.doi.org/10.13039/100007397 Univerzita Karlova v Praze 28122 Shin Hyun-Geun http://dx.doi.org/10.13039/100011199 FP7 Ideas: European Research Council CZ.02.1.01/0.0/0.0./16 019/0000778 Vlasák Miloslav Charles UniversityOpen access publishing supported by the National Technical Library in Prague.

issue-copyright-statement© Springer Science+Business Media, LLC, part of Springer Nature 2024
==== Body
pmcIntroduction

Fluid flows in variably saturated porous media are usually described by the Richards equation [33], which is expressed in the form1 ∂tϑ(ψ)-∇·(K(θ(ψ))(∇ψ+∇z))=g,

where ∂t denotes the derivative with respect to time, ∇· and ∇ are the divergence and gradient operators, respectively, ψ is the sought pressure head (= normalized pressure), z is the vertical coordinate, θ is the water content function, K is the hydraulic conductivity tensor and g is the source term. In addition, the active pore volume ϑ is related to θ by the following relation2 ϑ(ψ):=θ(ψ)+Ssθs∫-∞ψθ(s)ds,

where Ss,θs≥0 are material parameters. The hydraulic conductivity satisfies K(ψ)=KsKr(ψ), where Ks is the saturated conductivity tensor, and Kr∈[0,1] is the relative saturation. The functions θ and Kr are given by constitutive relations, e.g., by van Genuchten’s law [27] and by Mualem’s law [31], respectively.

The Richards equation belongs to the nonlinear parabolic problems, and it can degenerate, in particular K→0 or dϑdψ→0. Due to the degeneracy, the numerical solution is challenging, and various techniques have been developed for its solution in the last decades, see [25] for a survey.

In [14], we presented the adaptive space-time discontinuous Galerkin (STDG) method for the numerical solution of (1). This technique is based on a piecewise polynomial discontinuous approximation with respect to both the spatial and temporal coordinates. The resulting scheme is sufficiently stable, provides high accuracy, and is suitable for the hp-mesh adaptation. This is an important property, since the weak solution of the Richards equation is (only) piecewise regular and exhibits singularities along the material interfaces and the unsaturated/saturated zone (when ψ≈0). Therefore, an adaptive method that allows different meshes at different time levels, can achieve an accurate approximation with a relatively small number of degrees of freedom.

The numerical experiments presented in [14] showed the potential of the adaptive STDG method. However, the mesh adaptation used is based on interpolation error estimates that do not guarantee an upper error bound. The aim of this work is to overcome this bottleneck, derive a posteriori error estimates, and use them in the hp-mesh adaptation framework.

A posteriori error estimates for the numerical solution of the Richards equation have been treated in many papers for different numerical methods. We mention the finite volume framework with multistep time discretization in [5], the mixed finite element method in [6], the two-point finite volume discretization in [8], the lowest-order discretization on polytopal meshes in [38], finite element techniques in [30] and the references cited therein.

Guaranteed error estimates without unknown constants are usually obtained by measuring the error in a dual norm of the residual. Introducing reconstructed fluxes from the space H1(div,Ω), the upper bound can then be obtained directly. In [18], we developed this approach to the higher-order STDG method for nonlinear parabolic problems, where the temporal discontinuities were treated by temporal flux reconstructions considering the time jumps.

In this paper, we extend the approach [18] to the Richards equation (1). Although the definition of the temporal and spatial flux reconstructions as well as the derivation of the upper bounds is straightforward, the proof of the lower bound (efficiency) is rather tricky since the term θ(ψ) in the time derivative is not a polynomial function for a polynomial ψ. In contrary to [18], the proof of efficiency requires the additional oscillatory data terms. We construct spatial fluxes by solving local Neumann problems defined on space-time patches that generalize the approach from [22]. Moreover, we provide numerical experiments verifying derived error estimates. Compared to [18], the resulting effectivity indices are much closer to one. This is the first novelty of this paper.

Secondly, we deal with the errors arising due to iterative solution of nonlinear algebraic systems. We introduce a cheap stopping criterion for iterative solvers and justify it by numerical experiments. Thirdly, we introduce a space-time adaptive algorithm that employs the anisotropic hp-mesh adaptation technique [15]. The algorithm admits local adaptation of size and shape of mesh elements and the local adaptation of degrees of polynomial approximation with respect to space. However, the size of the time step can vary globally, and the degree of polynomial approximation with respect to time is fixed. Using the equidistribution principle, the algorithm gives an approximate solution with the error estimate under the given tolerance. The performance of the adaptive algorithm is demonstrated numerically, including a practically relevant example.

The rest of the paper is organized as follows. In Sect. 2, we introduce the problem considered, its STDG discretization is briefly described in Sect. 3. The main theoretical results are derived in Sect. 4, where the upper and lower bounds are proved. Two possible spatial reconstructions are discussed in Sect. 5 together with the stopping criteria of iterative solvers. The numerical verification of the error estimates is given in Sect. 6. Furthermore, we present the resulting hp-mesh adaptation algorithm in Sect. 7 and demonstrate its performance by numerical examples. Finally, we conclude with some remarks in Sect. 8.

Problem Formulation

Let Ω⊂Rd (d=2,3) be the domain occupied by a porous medium and T>0 the physical time to be reached. For simplicity, we assume that Ω is polygonal. By Γ:=∂Ω, we denote the boundary of Ω which consists of two disjoint parts: the Dirichlet boundary ΓD and the Neumann boundary ΓN. We write the Richards equation (1) in a different form, which is more suitable for the analysis. We seek a function u=u(x,t):Ω×(0,T)→R, which represents a hydraulic head (with the physical unit L). The quantity u is related to the pressure head ψ by u=ψ+z. The Richards equation (1) reads3 ∂tϑ(u)-∇·(K(u)∇u)=ginΩ×(0,T)u=uDonΓD×(0,T)K(u)∇u·n=gNonΓN×(0,T),u(x,0)=u0inΩ,

where g:Ω×(0,T)→R represents a source term if g is positive or a sink term if g is negative, ϑ:R→R denotes the dimensionless active pore volume, and K:R→Rd×d is the hydraulic conductivity with the physical unit L·T-1 (L = length, T = time). Moreover, uD is a trace of a function u∗∈L2(0,T;H1(Ω)) on ΓD×(0,T), gN∈L2(0,T;L2(ΓN)) and u0∈L2(Ω). We note that with respect to (1), we should write ϑ=ϑ(u-z) and K=K(θ(u-z)), however, we skip this notation for simplicity. We assume that the function ϑ(u) is non-negative, non-decreasing and Lipschitz continuous. Moreover, the tensor K(u) is symmetric, positively semi-definite, and Lipschitz continuous.

In order to introduce the weak solution, we set H(div,Ω)={v∈L2(Ω)d:∇·v∈L2(Ω)} and define the spaces4 X=L2(0,T,H1(Ω)),V={v∈X:v|ΓD=0},Y={v∈X:ϑ′(v)∈L2(0,T,L2(Ω))},Y0={v∈Y:v(0)=u0},

where ϑ′(u)=∂tϑ(u)=dϑdu∂tu denotes the time derivative (in the weak sense). Obviously, if v∈Y then ϑ(v)∈C([0,T],L2(Ω)). In order to shorten the notation, we set the physical flux5 σ(u,∇u):=K(u)∇u,u∈X.

Definition 1

We say that u∈Y is the weak solution of (3) if u-u∗∈V and6 ∫0Tϑ′(u),vΩ+σ(u,∇u),∇vΩ-g,vΩ-(gN,v)ΓNdt=0∀v∈V,

where (u,v)Ω:=∫Ωuvdx and (u,v)ΓN:=∫ΓNuvdS.

The existence and uniqueness of the Richards equation is studied in [2], see also the later works [3, 28].

Space-time discretization

We briefly describe the discretization of (6) by the space-time discontinuous Galerkin (STDG) method, for more details, see [13, 14]. Let 0=t0<t1<…<tr=T be a partition of the time interval (0, T) and set Im=(tm-1,tm) and τm=tm-tm-1. For each m=0,⋯,r, we consider a simplicial mesh Thm covering Ω¯. For simplicity, we assume that Thm, m=0,⋯,r are conforming, i.e., neighbouring elements share an entire edge or face. However, this assumption can be relaxed by the technique from [12].

For each element K∈Thm, we denote by ∂K its boundary, nK its unit outer normal and hK=diam(K) its diameter. In order to shorten the notation, we write ∂KN:=∂K∩ΓN. By the generic symbol γ, we denote an edge (d=2) or a face (d=3) of K∈Thm and hγ denotes its diameter. In the following, we speak only about edges, but we mean faces for d=3. We assume thatThm, m=0,⋯,r are shape regular, i.e., hK/ρK≤C for all K∈Th, where ρK is the radius of the largest d-dimensional ball inscribed in K and constant C does not depend on Thm for h∈(0,h0), m=0,⋯,r.

Thm, m=0,⋯,r are locally quasi-uniform, i.e., hK≤ChK′ for any pair of two neighbouring elements K and K′, where the constant C does not depend on h∈(0,h0), m=0,⋯,r.

Let pK≥1 be an integer denoting the degree of polynomial approximation on K∈Thm, m=0,⋯,r and PpK(K) be the corresponding space of polynomial functions on K. Let7 Shp,m={v∈L2(Ω):v|K∈PpK(K),K∈Thm},m=0,⋯,r

denote the spaces of discontinuous piecewise polynomial functions on Thm with possibly varying polynomial approximation degrees. Furthermore, we consider the space of space-time discontinuous piecewise polynomial functions8 Shpτq={v∈L2(Ω×(0,T)):v|Im∈Pq(Im,Shp,m),m=1,⋯,r},

where q≥0 denotes the time polynomial approximation degree and Pq(Im,Shp,m) is the Bochner space, i.e., v∈Pq(Im,Shp,m) can be written as v(x,t)=∑j=0qtjvj(x), vj∈Shp,m, j=0,⋯,q.

For v∈Shpτq, we define the one-sided limits and time jumps by9 v+m=limt→tm+v(t),m=0,…,r-1,v-m=limt→tm-v(t),m=1,…,r,{v}m=v+m-v-m,m=1,…,r-1,v-0=ϑ(u0),{v}0=v+0-ϑ(u0),

where u0 is the initial condition. In the following, we use the notation10 (u,v)M=∫Muvdx,(u,v)M,m=∫M×Imuvdxdt,m=1,⋯,r,

where M is either element K∈Thm or its (part of) boundary ∂K. The corresponding norms are denoted by ·M and ·M,m, respectively. By ∑K,m=∑m=1r∑K∈Thm, we denote the sum over all space-time elements K×Im, where K∈Thm and m=1,…,r.

Moreover, we define the jumps and mean values of v∈Shp,m on edges γ⊂∂K,K∈Thm by11 [v]=(v(+)-v(-))nKforγ∈Ω,(v(+)-uD)nKforγ⊂ΓD,0forγ⊂ΓN,v=(v(+)+v(-))/2forγ∈Ω,v(+)forγ⊂ΓD,0forγ⊂ΓN,

where v(+) and v(-) denote the traces of v on ∂K from interior and exterior of K, respectively, and uD comes from the Dirichlet boundary condition. For vector-valued v∈[Shp,m]d, we set [v]=(v(+)-v(-))·nK for γ∈Ω and similarly for boundary edges.

For each space-time element K×Im, K∈Thm, m=1,…,r, we define the forms12 aK,m(u,v):=(K(u)∇u,∇v)K,m-(g,v)K,m-(gN,v)∂KN,m,AK,m(u,v):=(K(u)∇u,∇v)K,m-(K(u)∇u·nK-α[u]·nK,v)∂K\ΓN,m+(β-12)(K(u)[u],∇v)∂K\Γ,m+(2β-1)(K(u)[u],∇v)∂K∩ΓD,m-(g,v)K,m-(gN,v)∂KN,m,

where α>0 is a sufficiently large penalization parameter (α∼pK2/hK) and β∈{0,12,1} corresponds to the choice of the variants of the interior penalty discretization (SIPG with β=0, IIPG with β=1/2 and NIPG with β=1), see, e.g., [13, Chapter 2].

We introduce the space-time discontinuous Galerkin discretization of (3).

Definition 2

The function uhτ∈Shpτq is called the approximate solution of (6) obtained by the space-time discontinuous Galerkin method (STDGM), if13 ∑K,mBK,m(uhτ,v)=0∀v∈Shpτq,

where14 BK,m(u,v):=(ϑ′(u),v)K,m+AK,mu,v+({ϑ(u)}m-1,v+m-1)K

with form AK,m given by (12) and {·} defined by (9).

Remark 1

We note that uhτ is discontinuous with respect to time at tm,m=1,⋯,r-1. The solution between Im-1 and Im is stuck together by the “time-penalty” term ({ϑ(u)}m-1,v+m-1)K which also makes sense for u and v belonging to different finite element spaces.

Finally, we derive some identities that will be used later. Let Fhm denote the set of all interior edges γ⊄Γ of mesh Thm and FDm the set of boundary edges of Thm lying on ΓD. Then, the identity15 ∑K∈Thm(w,znK)∂K\ΓN,m=∑γ∈Fhm(w,[z])γ,m+([w],z)γ,m+∑γ∈FDm(w·nK,z)γ,m

holds for a piecewise smooth vector-valued function w and a piecewise smooth scalar function z.

Using identity (15) and the following obvious formulas valid for interior edges K(u)∇u=K(u)∇u, α[u]=α[u], [K(u)∇u]=0, [α[u]]=0, we gain16 ∑K∈Thm(K(u)∇u·nK,v)∂K\ΓN,m=∑γ∈Fhm(K(u)∇u,[v])γ,m+∑γ∈FDm(K(u)∇u·nK,v)γ,m,∑K∈Thm(α[u]·nK,v)∂K\ΓN,m=∑γ∈Fhm(α[u],[v])γ,m+∑γ∈FDm(α[u]·nK,v)γ,m,∑K∈Thm(K(u)[u],∇v)∂K\Γ,m=∑K∈Thm([u],K(u)∇v)∂K\Γ,m=2∑γ∈Fhm([u],K(u)∇v)γ,m,∑K∈Thm(K(u)[u],∇v)∂K∩ΓD,m=∑γ∈FDm([u],K(u)∇v)γ,m.

Consequently, from (12) and (16), we obtain the identity17 ∑K∈ThmAK,m(u,v)=∑K∈Thm(K(u)∇u,∇v)K,m-∑γ∈Fhm(K(u)∇u,[v])γ,m+(2β-1)∑γ∈Fhm([u],K(u)∇v)γ,m-∑γ∈FDm(K(u)∇u·nK,v)γ,m+(2β-1)∑γ∈FDm([u],K(u)∇v)γ,m+∑γ∈Fhm(α[u],[v])γ,m+∑γ∈FDm(α[u]·nK,v)γ,m-(g,v)Ω,m-(gN,v)ΓN,m.

A Posteriori Error Analysis

Error Measures

In order to proceed to the derivation of error estimators, we define the spaces of piecewise continuous functions with respect to time by18 Yτ={v∈X:ϑ′(v)|Im∈L2(Im,L2(Ω))},Vτ={v∈Yτ:v|ΓD×(0,T)=0}.

Obviously, Y0⊂Y⊂Yτ⊂X and Shpτq⊂Yτ. Moreover, we have the following result.

Lemma 1

Let u∈Y0 be the weak solution of (6). Then it satisfies19 ∑K,mbK,m(u,v)=0∀v∈Vτ,

where20 bK,m(u,v):=(ϑ′(u),v)K,m+aK,m(u,v)+({ϑ(u)}m-1,v+m-1)K

with aK,m given by (12) and the time jump {·}m-1 defined by (9). Moreover, there exists a unique solution u∈Yτ such that u-u∗∈Vτ and satisfies (19).

Proof

The proof follows directly by comparing formulas (19)–(20) with (6) and the fact that ({ϑ(u)}m-1,v+m-1)K=0 for u∈Y0. For the proof of uniqueness, we employ the fact that C0∞(Ω) is dense in L2(Ω), i.e., there exists a sequence {vε}⊂C0∞(Ω) for any v∈L2(Ω) such that ‖vε-v‖→0 as ε→0, cf. [34, Theorem 3.14]. We apply v=vs,ε1(x)vt,ε2(t) in (19), where the spatial component vs,ε1∈{v∈H1(Ω):v|ΓD=0} tends to {ϑ(u)}m-1 as ε1→0 and the time component vt,ε2 is given as 0 outside the interval (tm-1,tm-1+ε2) and vt,ε2=1-(t-tm-1)/ε2 on (tm-1,tm-1+ε2), i.e., vt,ε2(t) tends to 0 as ε2→0. Therefore, all the terms containing time integrals in (19) tend to 0 when ε2 tends to 0. Since v+m-1=vs,ε1, the remaining jump term tends to ‖{ϑ(u)}m-1‖2 as ε1 tends to 0. From this it follows that {ϑ(u)}m-1=0. Then it is possible to see that any solution of (19) satisfies the original weak formulation (6). Since the weak problem (6) has a unique solution, cf. [2], the extended problem (19) has a unique solution as well. □

In virtue of [11, § 2.3.1], we define a parameter dK,m associated with the space-time element K×Im, K∈Thm, m=1,⋯,r. The parameter dK,m represents a user-dependent weight, typically with physical units (TL)1/2 so that the error measure has the same physical unit as the energy norm. In this paper, we use two choices 21a dK,m:=hK-2K(uh)m,∞+τm-2Tdϑdu(uh)m,∞-1/2,

21b dK,m:=hK2K(uh)m,∞-1+τm2/Tdϑdu(uh)m,∞-11/2.

where ·m,∞:=·L∞(Ω×Im). We note that the following error analysis is independent of the choice of dK,m. Moreover, we define the norm in the space Vτ (cf. (18)) by22 vVτ2=∑K,mvVK,m2,vVK,m2=dK,m-2hK2∇vK,m2+τm2v′K,m2.

In virtue of (19), we introduce the error measure as a dual norm of the residual23 R(uhτ)=sup0≠v∈Vτ∑K,mbK,m(uhτ,v)vVτ,

where bK,m is given by (20). The residual R(v) represents a natural error measure for u-v∈Vτ, cf. [11, Remark 2.3]. In Sect. 4, we estimate R(uhτ) for uhτ being the solution of (13).

Since the approximate solution uhτ belongs to the space of discontinuous function Shpτq⊄Vτ, we introduce the second building block measuring the nonconformity of the solution in spatial variables. Therefore, similarly to [18], we define the form24 J(v)=∑K,mJK,m(v),JK,m(v)=dK,m2τm-1hK-2CK,m,K,α[v]∂K,m2,

where CK,m,K,α=α2+K(uhτ)L∞(K×Im)2. The scaling factors are chosen such that J(v)1/2 has the same physical unit as R(uhτ).

We note that J(v) measures also the violation of the Dirichlet boundary condition since J(v) contains the term v-uD∂K∩ΓD,m, cf. (11).

The final error measure is then defined by25 E(uhτ):=R(uhτ)2+J(uhτ)1/2,

where R(uhτ) is given by (23) and J(uhτ) by (24).

Lemma 2

The error measure E(uhτ)=0 if and only if uhτ=u is the weak solution given by (6).

Proof

Obviously, if uhτ=u, then J(uhτ)=0 and R(uhτ)=0 due to (19). On the other hand, if J(uhτ)=0, then uhτ∈Yτ and uhτ-u∗∈Vτ. Moreover, R(uhτ)=0 and the uniqueness of (19) imply that uhτ is the weak solution (6). □

Temporal and Spatial Flux Reconstructions

Similarly as in [18], we define a temporal reconstruction Rhτ=Rhτ(x,t) as a continuous function with respect to time that mimics ϑ(uhτ), uhτ∈Shpτq. Let rm∈Pq+1(Im) be the right Radau polynomial on Im, i.e., rm(tm-1)=1 and rm(tm)=0, and rm is orthogonal to Pq-1(Im) with respect to the L2(Im) inner product. Then we set26 Rhτ(x,t):=ϑ(uhτ(x,t))-{ϑ(uhτ)}m-1(x)rm(t),x∈Ω,t∈Im,

where {·} is given by (9). The temporal flux reconstruction Rhτ(x,t) is continuous in time, namely Rhτ∈H1(0,T,L2(Ω)) and it satisfies the initial condition due to27 Rhτ(·,0)=ϑ(uhτ(·,0))-{ϑ(uhτ)}0(·)r1(0)=ϑ(uhτ(·,0))-(ϑ(uhτ(·,0))-ϑ(u0(·))=ϑ(u0(·)).

Moreover, by the integration by parts and the properties rm(tm-1)=1, rm(tm)=0, we obtain28 ((Rhτ-ϑ(uhτ))′,v)K,m=-(rm′{ϑ(uhτ)}m-1,v)K,m=(rm{ϑ(uhτ)}m-1,v′)K,m-rm(tm)({ϑ(uhτ)}m-1,v-m)K+rm(tm-1)({ϑ(uhτ)}m-1,v+m-1)K=(rm{ϑ(uhτ)}m-1,v′)K,m+({ϑ(uhτ)}m-1,v+m-1)K,v∈Vτ,

which together with definition (26) implies29 ((Rhτ-ϑ(uhτ))′,v)K,m-({ϑ(uhτ)}m-1,v+m-1)K=-(Rhτ-ϑ(uhτ),v′)K,m,v∈Vτ.

Finally, using the orthogonality of rm to Pq-1(Im), we obtain from (28), the formula30 (Rhτ-ϑ(uhτ))′,vm,K={ϑ(uhτ)}m-1,v+m-1K∀v∈Pq(Im,L2(K)).

Consequently, if uhτ is the approximate solution given by (13), then it satisfies31 ((Rhτ)′,v)K,m=(ϑ′(uhτ),v)K,m+({ϑ(uhτ)}m-1,v+m-1)K=-AK,m(uhτ,v)∀v∈Pq(Im,PpK(K)).

Obviously, the reconstruction Rhτ is local and explicit, so its computation is fast and easy to implement.

The spatial flux reconstruction needs to define a function σhτ∈L2(0,T,H(div,Ω)) which mimics the flux σ(uhτ,∇uhτ)=K(uhτ)∇uhτ, cf. (5). In particular, σhτ|K×Im∈Pq(Im,RTNp(K)) where32 RTNp(K)=Pp(K)d+xPp(K),K∈Th,m=1,⋯,r

is the Raviart-Thomas-Nedelec finite elements, cf. [7] for more details. We assume that the reconstructed flux σhτ has to be equilibrated with the temporal flux Rhτ33 (∇·σhτ,v)K,m=((Rhτ)′-g,v)K,m∀v∈Pq(Im,PpK(K)),K∈Thm,

and with the Neumann boundary condition34 (σhτ·n,v)γ,m=(gN,v)γ,m∀v∈Pq(Im,PpK(γ))∀γ⊂∂KN,K∈Thm.

In Sect. 5 we present two possible constructions of σhτ including the choice of the spatial polynomial degree p in (32).

Auxiliary Results

In the forthcoming numerical analysis, we need several technical tools. We will employ the scaled space-time Poincarè inequality, cf. [11, Lemma 2.2]: Let φK,m∈P0(K×Im) be the L2-orthogonal projection of φ∈H1(K×Im) onto a constant in each space-time element K×Im, K∈Thm, m=0,⋯,r. Then,35 φ-φK,mK,m≤CPhK2∇φK,m2+τm2φ′K,m21/2=CPdK,mφVK,m,

where CP is the Poincarè constant equal to 1/π for simplicial elements and the last equality follows from (22).

Moreover, we introduce the space-time trace inequality

Lemma 3

Let φγ,m∈P0(γ×Im) be the L2-orthogonal projection of φ∈H1(K×Im) onto a constant on each γ×Im, where γ⊂∂K is an edge of K∈Thm. Then there exists a constant CT>0 such that36 φ-φγ,mγ×Im≤CTmax(1,hγ-1/2)dK,mφVK,m,

where CT=max(cT,CP), CP is from (35) and cT>0 is the constant from the (space) trace inequality.

Proof

The proof is straightforward, we present it for completeness. Let φ∈H1(K×Im) and, for all t∈Im, set φ~(t):=|γ|-1∫γφ(x,t)dS. Observing that (φ-φ~) and (φ~-φγ,m) are L2(γ×Im)-orthogonal, we have37 φ-φγ,mγ×Im2=φ-φ~γ×Im2+φ~-φγ,mγ×Im2.

Using the standard trace inequality (e.g., [21, Lemma 3.32]), we have38 φ(·,t)-φ~(t)γ≤cThγ1/2∇φK∀t∈Im,

where cT>0 is a constant whose values can be set relatively precisely, see the discussion in [37, Section 4.6]. Hence, integrating the square of (38) over Im and using the fact that hγ≤hK, γ⊂hK, we estimate the first term on the right-hand side of (37) as39 φ-φ~γ×Im2≤cT2hγ∇φK×Im2≤cT2hγ-1hK2∇φK×Im2.

Using the fact that φγ,m=τm-1∫Imφ~(t)dt, the one-dimensional Poincarè inequality on In and the Cauchy–Schwarz inequality yield40 φ~-φγ,mγ×Im2=|γ|∫Im|φ~-φγ,m|2(t)dt≤|γ|CP2τm2∫Im|ddtφ~(t)|2dt=CP2τm2|γ|∫Im∫γ∂tφ(x,t)dx2dt≤CP2τm2∫Im∫γ|∂tφ|2dxdt=CP2τm2∂tφγ×Im2.

Collecting bounds (37), (39), (40) and the definition of the norm (22) yields (36). □

Reliability

We presented the upper bound of R(uhτ), cf. (23).

Theorem 1

Let u∈Y be the weak solution of (6) and uhτ∈Shpτq be the approximate solution given by (13). Let Rhτ∈H1(0,T,L2(Ω)) be the temporal reconstruction given by (26) and σhτ∈L2(0,T,H(div,Ω)) be the spatial reconstruction satisfying (33). Then41 R(uhτ)2≤η2:=∑K,mηK,m2,ηK,m:=CPηR,K,m+(ηS,K,m2+ηT,K,m2)1/2+CTηN,K,m,

where CP is the constant from Poincarè inequality (35), CT is the constant from the trace inequality (36) and the estimators ηR,K,m, ηS,K,m, ηT,K,m, and ηN,K,m are given by 42a ηR,K,m:=dK,m(Rhτ)′-∇·σhτ-gK,m,

42b ηS,K,m:=dK,mhKσhτ-σ(uhτ,∇uhτ)K,m,

42c ηT,K,m:=dK,mτmRhτ-ϑ(uhτ)K,m,

42d ηN,K,m:=∑γ⊂∂KNmax(1,hγ-1/2)dK,mσhτ·n-gN∂KN,m.

The proof of Theorem 1 can be found in [19] for the case of the homogeneous Dirichlet boundary condition. For completeness, we present its modification including mixed Dirichlet-Neumann boundary conditions.

Proof

Starting from (20), adding the terms ±(Rhτ,v)K,m and ±(∇·σhτ,v)K,m, and using the integration by parts, we obtain43 ∑K,mbK,m(uhτ,v)=∑K,m(ϑ′(uhτ)-g,v)K,m-(gN,v)∂KN,m+(σ(uhτ,∇uhτ),∇v)K,m+({ϑ(uhτ)}m-1,v+m-1)K=∑K,m((Rhτ)′-∇·σhτ-g,v)K,m-∑K,m(σhτ-σ(uhτ,∇uhτ),∇v)K,m-∑K,m((Rhτ-ϑ(uhτ))′,v)K,m-({ϑ(uhτ)}m-1,v+m-1)K+∑K,m(σhτ·n-gN,v)∂KN,m=:ξ1+ξ2+ξ3+ξ4.

The terms ξi, i=1,⋯,4 are estimated separately.

Let vK,m∈P0(K×Im) be the piecewise constant projection of v∈Vτ given by the identity (vK,m,1)K,m=(v,1)K,m. Using the Cauchy–Schwarz inequality, assumption (33), the Poincarè inequality (35), and (22), we have44 |ξ1|≤∑K,m|((Rhτ)′-∇·σhτ-g,v)K,m|=∑K,m|((Rhτ)′-∇·σhτ-g,v-vK,m)K,m|≤∑K,mCP(Rhτ)′-∇·σhτ-gK,mhK2∇vK,m2+τm2v′K,m21/2=∑K,mCPdK,m(Rhτ)′-∇·σhτ-gK,mvVK,m=∑K,mCPηR,K,mvVK,m.

Furthermore, by the Cauchy–Schwarz inequality and (22), we obtain45 |ξ2|≤∑K,m|(σhτ-σ(uhτ,∇uhτ),∇v)K,m|≤∑K,mdK,mhKσhτ-σ(uhτ,∇uhτ)K,mhKdK,m∇vK,m=∑K,mηS,K,mhKdK,m∇vK,m.

The use of (29), and a similar manipulations as in (45), give46 |ξ3|≤∑K,m|((Rhτ-ϑ(uhτ))′,v)K,m-({ϑ(uhτ)}m-1,v+m-1)K|=∑K,m|(Rhτ-ϑ(uhτ),v′)K,m|≤∑K,mdK,mτmRhτ-ϑ(uhτ)K,mτmdK,mv′K,m=∑K,mηT,K,mτmdK,mv′K,m.

Hence, estimates (45)–(46), the Cauchy inequality and (22) imply47 |ξ2|+|ξ3|≤∑K,mηS,K,mhKdK,m∇vK,m+ηT,K,mτmdK,mv′K,m≤∑K,mηS,K,m2+ηT,K,m21/2vVK,m.

Furthermore, let vγ,m∈P0(γ×Im), γ⊂∂KN be the L2-orthogonal projection from Lemma 3. Then using assumption (34), the Cauchy inequality and the space-time trace inequality (36), we have48 |ξ4|=∑K,m∑γ⊂∂KN(σhτ·n-gN,v-vγ,m)γ,m≤∑K,m∑γ⊂∂KNσhτ·n-gNγ,mv-vγ,mγ,m≤CT∑K,m∑γ⊂∂KNmax(1,hγ-1/2)dK,mσhτ·n-gNγ,mvVK,m.

The particular estimates (44), (47), and (48), together with the discrete Cauchy–Schwarz inequality, imply (41). □

Remark 2

Obviously, if ∂K∩ΓN≠∅, then ηN,K,m=0.

Efficiency

The aim is to show that the local individual error estimators ηR,K,m, ηS,K,m and ηT,K,m from (41)–(42) are locally efficient, i.e., they provide local lower bounds to the error measure up to a generic constant C>0 which is independent of u, uhτ, h and τ, but may depend on data problems and the degrees of polynomial approximation p and q. A dependence of the estimate up to this generic constant we will denote by ≲.

In order to derive the local variants of the error measure, we denote by ωK the set of elements sharing at least a vertex with K∈Thm, i.e.,49 ωK=∪K′∩K≠0K′,K∈Thm,m=0,⋯,r.

Moreover, we define the functional sub-spaces VD,m={v∈Vτ:supp(v)⊂D×Im¯} for any set D⊂Ω (cf. (18)) and the corresponding error measures (cf. (23))50 RD,m(w)=sup{0≠v∈VD,m}1vVτ∑K,mbK,m(w,v).

Obviously, the definition of VD,m and RD,m(uhτ) together with the shape regularity implies51 ∑K,mRK,m(uhτ)≤∑K,mRωK,m(uhτ)≲R(uhτ).

Moreover, for each space-time element K×Im, K∈Thm, m=1,⋯,r, we introduce the L2(K×Im)-projection of the non-polynomial functions, namely52 ϑ′(uhτ)¯∈Pq(Im,PpK(K):(ϑ′(uhτ)¯,v)K,m=(ϑ′(uhτ),v)K,m∀v∈Pq(Im,PpK(K))g¯∈Pq(Im,PpK(K)):(g¯,v)K,m=(g,v)K,m∀v∈Pq(Im,PpK(K)).

Finally, for each vertex a of the mesh Thm, we denote by ωa a patch of elements K∈Thm that share this vertex. By pa=maxK∈ωapK we denote the maximal polynomial degree on ωa. Then, for each a of K∈Thm, we define a vector-valued function σ¯a=σ¯a(uhτ,∇uhτ)∈Pq(Im,RTNpa(K)) (cf. (32)) by53 (σ¯a·nK,v)γ,m=(ψaσ(uhτ,∇uhτ)·nK,v)γ,m∀v∈Pq(Im,Ppa(γ)),γ⊂K(σ¯a·v)K,m=(ψaσ(uhτ,∇uhτ),v)K,m∀v∈Pq(Im,Ppa-1(K)d),

where · denotes the mean value on γ⊂∂K and ψa is a continuous piecewise linear function such that ψa(a)=1 and it vanishes at the other vertices of K. Finally, we set σ¯|K×Im=∑a∈Kσ¯a.

The proof of the local efficiency of the error estimates presented is based on the choice of a suitable test function in (23). We set54 w(x,t)=dK,m2τmPh({ϑ(uhτ)}m-1)(x)χK(x)Φm(t).

where χK(x) is the standard bubble function on K, Φm(t) is the Legendre polynomial of degree q+1 on Im (and vanishing outside) and Ph({ϑ(uhτ)}m-1)∈PpK(K) is the L2(K)-projection weighted by χK(x) given by55 (Ph({ϑ(uhτ)}m-1),χKv)K=({ϑ(uhτ)}m-1,χKv)K∀v∈PpK(K).

We note that56 Ph({ϑ(uhτ)}m-1)≠{ϑ(uhτ)¯}m-1,

in general, compare with (52).

Using the inverse inequality, the polynomial function w given by (54) can be estimated as57 wVK,m2=dK,m-2hK2∇wK,m2+τm2w′K,m2≲dK,m-2wK,m2≤dK,m2τm2Ph({ϑ(uhτ)}m-1)K2∫ImΦm2(t)dt≲dK,m2τmPh({ϑ(uhτ)}m-1)K2.

Similarly as in [11] or [18], we introduce the oscillation terms58 ηG,K,m:=dK,mg¯-gK,m,ηϑ,K,m:=dK,mτm{ϑ(uhτ)}m-1-Ph({ϑ(uhτ)}m-1)K,ηϑ′,K,m:=dK,mϑ′(uhτ)¯-ϑ′(uhτ)K,m,ησ,K,m:=dK,mhKσ¯-σ(uhτ,∇uhτ)K,m+dK,m∇·σ¯-∇·σ(uhτ,∇uhτ)K,m.

The goal is to prove the lower bounds of the proposed error estimates, namely to estimate ηT,K,m, ηR,K,m and ηS,K,m by RK,m(uhτ) and the oscillation terms (58), K∈Th, m=1,⋯,r.

Theorem 2

Let ηT,K,m, K∈Thm, m=1,⋯,r be the error estimates given by (42), then59 ηT,K,m≲RK,m(uhτ)+ηG,K,m+ηϑ′,K,m+ηϑ,K,m+ηS,K,m.

where RK,m are the local error measures defined by (49)–(50) and the oscillation terms ηG,K,m, ηϑ,K,m and ηϑ′,K,m are given by (58).

Proof

We start the proof by the putting function w from (54) as the test function in (50), i.e.60 RK,m(uhτ)=sup0≠v∈VK,m∑K,mbK,m(uhτ,v)vVτ≥bK,m(uhτ,w)wVτ

since supp(w)=K×Im, cf. (54). Then, using (20) and the fact that w vanishes on ∂K, we have61 RK,m(uhτ)≥(ϑ′(uhτ)-g,w)K,m+(σ(uhτ,∇uhτ),∇w)K,m+({ϑ(uhτ)}m-1,w+m-1)KwVK,m=(ϑ′(uhτ)¯-g¯,w)K,m+(σhτ,∇w)K,mwVK,m+({ϑ(uhτ)}m-1,w+m-1)KwVK,m=:ξ1+ξ2+(g¯-g,w)K,m+(σ-σhτ,∇w)K,m+(ϑ′(uhτ)-ϑ′(uhτ)¯,w)K,mwVK,m=:ξ3+ξ4+ξ5.

The functions ϑ′(uhτ)¯, g¯ and σhτ are polynomials of degree q in time whereas w and ∇w are the (Legendre) polynomial of degree (q+1) in time, cf. (54). Due to the L2(Im)-orthogonality of the Legendre polynomials, we have ξ1=0, since62 (ϑ′(uhτ)¯-g¯,w)K,m+(σhτ,∇w)K,m=0

Moreover, using inequality (57), relations (54)-(55) and the equivalence of norms on finite dimensional spaces,

we obtain63 ξ2≳(Ph({ϑ(uhτ)}m-1),dK,m2τmPh({ϑ(uhτ)}m-1)χK)KdK,mτmPh({ϑ(uhτ)}m-1)K≳dK,mτm(Ph({ϑ(uhτ)}m-1),Ph({ϑ(uhτ)}m-1))KPh({ϑ(uhτ)}m-1)K=dK,mτmPh({ϑ(uhτ)}m-1)K.

Furthermore, let wK,m=1K×Im∫K×Imwdxdt be the mean value of w on the space-time element K×Im. Due to (52), the Cauchy–Schwarz inequality and (35), we have64 |ξ3|=|(g¯-g,w-wK,m)K,m|wVK,m≤g¯-gK,mw-wK,mK,mwVK,m≲dK,mg¯-gK,m=ηG,K,m,

and65 |ξ5|≲dK,mϑ′(uhτ)-ϑ′(uhτ)¯K,m=ηϑ′,K,m.

Similarly, the Cauchy–Schwarz inequality and (22) imply66 |ξ4|≤dK,mhKσ(uhτ,∇uhτ)-σhτK,mhK∇wK,mdK,mwVK,m≤dK,mhKσ(uhτ,∇uhτ)-σhτK,m=ηS,K,m.

Collecting (61)–(66), we have67 RK,m(uhτ)≳dK,mτmPh({ϑ(uhτ)}m-1)K-ηG,K,m-ηS,K,m-ηϑ′,K,m.

Moreover, using (42c), (26), integration by parts, the boundedness of the Radau polynomials, the triangle inequality and (58), we have68 ηT,K,m=dK,mτmRhτ-ϑ(uhτ)K,m=dK,mτm{ϑ(uhτ)}m-1rmK,m=dK,mτm{ϑ(uhτ)}m-1K∫Imrm2dt≲dK,mτm{ϑ(uhτ)}m-1K≤dK,mτmPh({ϑ(uhτ)}m-1)K+ηϑ,K,m.

Hence, (67) and (68)69 ηT,K,m≤RK,m(uhτ)+ηϑ,K,m+ηG,K,m+ηϑ′,K,m+ηS,K,m,

which proves the theorem. □

Theorem 3

Let ηS,K,m and ηR,K,m, K∈Thm, m=1,⋯,r be the error estimates given by (42), then70 ηR,K,m≲RωK,m(uhτ)+ηG,K,m+ησ,K,m+ηS,K,m,

71 ηS,K,m≲RωK,m(uhτ)+ηG,K,m+∑K⊂ωKησ,K,m,

where RωK,m is the local error measures defined by (49)–(50) and the oscillation terms ηG,K,m, ηϑ,K,m and ηϑ′,K,m are given by (58).

Proof

The proof is in principle identical with the proof [18, Lemmas 7-9], we present the main step for completeness. Let g¯ and σ¯ be the projection given by (52) and (53). Using the triangle inequality, the inverse inequality and (58), we obtain72 ηR,K,m=dK,m(Rhτ)′-∇·σhτ-gK,m≤dK,m(Rhτ)′-∇·σ¯-g¯K,m+dK,mg¯-gK,m+dK,m∇·σ¯-∇·σhτK,m≲dK,m(Rhτ)′-∇·σ¯-g¯K,m+ηG,K,m+dK,mhKσ¯-σhτK,m.

The first term on the right-hand side of (72) can be estimated as in [36, Theorem 4.10] by73 dK,m(Rhτ)′-∇·σ¯-g¯K,m≲ResωK,m(uhτ)+ηG,K,m+ησ,K,m,

where the resulting oscillation terms are estimated with the aid (58). Moreover, the last term on the right-hand side of (72) together with (42b) and assumption (58), reads74 dK,mhKσ¯-σhτK,m≤dK,mhKσ¯-σ(uhτ,∇uhτ)K,m+dK,mhKσ(uhτ,∇uhτ)-σhτK,m≤ησ,K,m+ηS,K,m,

which proves (70).

The proof of (71) is based on the decomposition75 σhτ-σ(uhτ,∇uhτ)K,m≤σhτ-σ¯K,m+σ¯-σ(uhτ,∇uhτ)K,m.

While the second term on the right-hand side of (75) can be estimated by assumption (58), the estimate of the first term is somewhat more technical. It depends on the flux reconstruction used. For the flux reconstruction in Sect. 5.2, the proof is identical to the proof of [18, Lemma 9], which mimics the stationary variant [24, Theorem 3.12]. On the other hand, using the flux reconstruction from Sect. 5.1, it is possible to apply the technique from [11, Lemma 7.5], where the final relation has to be integrated over Im. □

Spatial Flux Reconstructions and Stopping Criteria

We present two ways of reconstructing the spatial flux σhτ∈L2(0,T,H(div,Ω)) that satisfies the assumptions (33)–(34). The first one, proposed in [19] for the case of homogeneous Dirichlet boundary condition, is defined by the volume and edge momenta of the Raviart-Thomas-Nedelec (RTN) elements, cf. [7], and is easy to compute. The second approach is based on the solution of local Neumann problems on patches associated with each vertex of the mesh. This idea comes from, e.g., [24], its space-time variant was proposed in [18] for nonlinear convection-diffusion equations. Finally, in Sect. 5.3, we discuss the errors arising from the solution of algebraic systems and introduce a stopping criterion for the appropriate iterative solver.

Element-Wise Variant

We denote by pK,max the maximum polynomial degree over the element K and its neighbours that share the entire edge with K and pγ,max the maximum polynomial degree on neighbouring elements having a common edge γ. Let RTNpK,max(K) be the space of RTN finite elements of order pK,max for element K∈Thm, cf. (32), and uhτ∈Shpτq be the approximate solution. The spatial reconstruction σhτ is defined element-wise: for each K∈Thm, find σhτ|K×Im∈Pq(Im,RTNpK,max(K)) with σhτ·n|γ×Im∈Pq(Im,Ppγ,max(γ)) such that76 (σhτ·n,v)γ,m=(K(uhτ)∇uhτ·n-α[uhτ]·n,v)γ,m∀v∈Pq(Im,Ppγ,max(γ)),γ⊂∂K\ΓN(gN,v)γ∀v∈Pq(Im,Ppγ,max(γ)),γ⊂∂KN(σhτ,v)K,m=(K(uhτ)∇uhτ,∇v)K,m+(β-12)(K(uhτ)[uhτ],∇v)∂K\Γ,m+(2β-1)(K(uhτ)[uhτ],∇v)∂K∩ΓD,m∀v∈Pq(Im,PpK,max-1(K)d).

The edge momenta in (76) are uniquely defined and since pγ,max≤pK,max, σhτ in (76) is well defined as well. Here, the numerical flux K(uhτ)∇uhτ·n-α[uhτ]·n is conservative on interior edges, which implies that σhτ·n are the same on each interior edge γ and therefore the resulting reconstruction σhτ∈L2(0,T,H(div,Ω)) globally.

Obviously, the first relation in (76) with pK≤pγ,max directly implies assumption (34). Moreover, using the Green theorem, (76), (12), (31) and pK≤pγ,max≤pK,max, we obtain77 (∇·σhτ,v)K,m=-(σhτ,∇v)K,m+(σhτ·nK,v)∂K,m=-(K(uhτ)∇uhτ,∇v)K,m+(K(uhτ)∇uhτ·n-α[uhτ]·n,v)∂K\ΓN,m-(β-12)(K(uhτ)[uhτ],∇v)∂K\Γ,m-(2β-1)(K(uhτ)[uhτ],∇v)∂K∩ΓD,m+(gN,v)∂KN,m=-AK,m(uhτ,v)-(g,v)K,m=((Rhτ)′-g,v)K,m∀v∈Pq(Im,PpK(K)),K∈Thm,

which justifies the assumption (33).

Patch-Wise Variant

For each vertex a of the mesh Thm, we denote by ωa a patch of elements K∈Thm sharing this vertex. By pa=maxK∈ωapK we denote the maximal polynomial degree on ωa. Let Ppa∗(ωa) be the space of piecewise polynomial discontinuous functions of degree pa on ωa with mean value zero for a∉∂Ω. We define the space78 WRTN,paN(ωa)={v∈H(div,ωa);v|K∈RTNpa(K),v·n=0on∂ωa},a∉∂ΩWRTN,paN(ωa)={v∈H(div,ωa);v|K∈RTNpa(K),v·n=0on∂ωa\∂Ω,&(v·n,ϕ)γ,m=(gN,ϕ)γ,m∀ϕ∈Pq(Im,Ppa(γ))on∂ωa∩∂KN},a∈∂Ω.

We set the local problems on patches ωa for all vertices a: find σhτ∈Pq(Im,WRTN,paN(ωa)) and raτ∈Pq(Im,Ppa∗(ωa)) such that79 (σaτ,v)ωa,m-(raτ,∇·v)ωa,m=(ξa1,v)ωa,m∀v∈Pq(Im,WRTN,paN(ωa))(∇·σaτ,ϕ)ωa,m=(ξa2,ϕ)ωa,m∀ϕ∈Pq(Im,Ppa∗(ωa)),

where80 ξa1=ψaσ(uhτ,∇uhτ)ξa2=ψa(Rhτ)′-ψag+∇ψa·ξ(uhτ,∇uhτ),

with81 ξ(uhτ,∇uhτ)=σ(uhτ,∇uhτ)+(2β-1)∑γ⊄ΓNℓm,γ(uhτ),

and ℓm,γ:Shp,m→[Sh0,m]d is the lifting operator defined by82 ∫Ωℓm,γ(uhτ)·vdx=∫γ[uhτ]K(uhτ)vdx∀v∈[Sh0,m]d,γ⊄ΓN.

Then the final reconstructed flux is obtained by summing up σaτ on each element that contains vertex a, i.e.,83 σhτ|K,m=∑a∈Kσaτ|K.

The assumption (34) follows directly from (78) and pK≤pa. Inserting the hat function ψav for a∉∂Ω and v∈Pq(Im) in (17), using (5), (82) and omitting the zero terms, we have84 ∑K∈ThmAK,m(uhτ,ψav)=∑K∈Thm(K(uhτ)∇uhτ,∇ψav)K,m+(2β-1)∑γ⊄∂Ω([uhτ],K(uhτ)∇ψav)γ,m+(2β-1)∑γ⊂ΓD([uhτ],K(uhτ)∇ψav)γ,m-(g,ψav)Ω,m=(ξa2,v)ωa,m-(Rhτ,ψav)ωa,m

Applying (13) and (31), we gain for a∉∂Ω and v∈Pq(Im)85 (∇·σaτ,v)ωa,m=∑K⊂ωaAK,m(uhτ,ψav)+(Rhτ,ψav)K,m=(ξa2,v)ωa,m.

From this it follows that the second relation in (79) holds element-wise, i.e.86 (∇·σaτ,ϕ)K,m=(ξa2,ϕ)K,m,∀ϕ∈Pq(Im,Ppa(K)).

Then (33) follows from87 (∇·σhτ,ϕ)K,m=∑a⊂K(∇·σaτ,ϕ)K,m=∑a⊂K(ξa2,ϕ)K,m=((Rhτ)′-g,ϕ)K,m∀ϕ∈Pq(Im,Ppa(K))

and from pK≤pa.

Stopping Criteria for Iterative Solvers

The space-time discretization (13) leads to a system of nonlinear algebraic equations for each time level m=1,⋯,r. These systems have to be solved iteratively by a suitable solver, e.g., the Picard method, the Newton method or their variants. Therefore, it is necessary to set a suitable stopping criterion for the iterative solvers so that, on the one hand, the algebraic errors do not affect the quality of the approximate solution and, on the other hand, an excessive number of algebraic iterations is avoided.

However, the error estimates presented in Sect. 4 do not take into account errors arising from the inaccurate solution of these systems. Indeed, the aforementioned reconstructions fulfill assumption (33) only if the systems given by (13) are solved exactly. The full a posteriori error analysis including algebraic errors has been treated, e.g., in [8, 23, 29]. These error estimators are based on additional flux reconstructions that need to be evaluated at each iteration, and therefore, the overall computational time is increased.

To speed up the computations and control the algebraic errors, we adopt the technique of [17]. This approach offers (i) the measurement of algebraic errors by a quantity similar to the error measure (23), (ii) the setting of the stopping criterion for iterative solvers with one parameter corresponding to the relative error, and (iii) a fast evaluation of the required quantities.

For each m=1,⋯,r, we define the estimators (cf. (23))88 ηalgm(uhτ)=sup0≠v∈Shpτq∑K∈ThmbK,m(uhτ,v)vVτ,ηspam(uhτ)=sup0≠v∈Shp+1τq+1∑K∈ThmbK,m(uhτ,v)vVτ,

where the norm ·Vτ is given by (22),89 Shp+1τq+1={v∈L2(Ω×(0,T)):v|Im∈Pq+1(Im,Shp+1,m),m=1,⋯,r},andShp+1,m={v∈L2(Ω):v|K∈PpK+1(K),K∈Thm},m=0,⋯,r.

The space Shp+1τq+1 is an enrichment space of Shpτq by polynomials of the space degree pK+1 and the time degree q+1 for each K×Im, K∈Thm, m=1,⋯,r. Finally, we define the global in time quantities90 ηalg(uhτ)=∑m=1r(ηalgm(uhτ))21/2,ηspa(uhτ)=∑m=1r(ηspam(uhτ))21/2.

Obviously, if uhτ fulfills (13) exactly, then ηalgm(uhτ)=0 for all m=0,⋯,r. Moreover, if uhτ is the weak solution (6) then ηspam(uhτ)=0 for all m=0,⋯,r. Comparing (88) with (23), the quantity ηspa(uhτ) exhibits a variant of the error measure R(uhτ). Nevertheless, ηspa(uhτ) is neither lower nor upper bound of R(uhτ) since Shp+1τq+1⊄Vτ and Vτ⊄Shp+1τq+1.

The quantities (88) can be evaluated very fast since the suprema (maxima) are the sum of the suprema (maxima) for all space-time elements K×Im, K∈Thm, m=1,⋯,r, which are computable separately, cf. [17] for details. Hence, we prescribe the stopping criterion for the corresponding iterative solver as91 ηalgm(uhτ)≤cAηspam(uhτ),m=1,⋯,r,

where cA∈(0,1) is the user-dependent constant. The justification of this approach and the influence of algebraic errors on the error estimates are studied numerically in Sect. 6.1.1.

Numerical Experiments

We present numerical experiments that justify the a posteriori error estimates (41)–(42). Since the error measure (23) is the dual norm of the residual, it is not possible to evaluate the error even if the exact solution is known. Therefore, similarly to [18], we approximate the error by solving the dual problem given for each time interval Im,m=1,⋯,r: Find ψm∈Ymτ=L2(Im,H1(Ω)),92 (ψm,ϕ)Ymτ=∑K,mbK,m(uhτ,ϕ)∀ϕ∈Ymτ,

where (cf. (21a)–(22))93 (u,v)Ymτ=∑K∈ThmdK,m-2hK2(∇u,∇v)K,m+τm2(u′,v′)K,m,m=1,⋯,r.

Then we have R(uhτ)2=∑m=1rψYmτ2. We solve (92) for each m=1,⋯,r by linear conforming finite element on a global refinement of the space-time mesh Thm×Im which is proportional to the space and time polynomial approximation degrees. We denote this quantity by R~(uhτ). The second error contribution J given by (24) is computable, so the total error E (cf. (25)) is approximated by E~(uhτ):=R~(uhτ)2+J(uhτ)1/2.

Remark 3

Sometimes, this approximate evaluation of the (exact) error is not sufficiently accurate for fine grids and high polynomial approximation degrees. In this case, very fine global refinement is required and then the resulting algebraic systems are too large to be solved in a reasonable time.

All numerical experiments were carried out using the patch-wise reconstruction from Sect. 5.2 using the in-house code ADGFEM [10]. The arising nonlinear algebraic systems are solved iteratively by a Newton-like method, we refer to [14] for details. Each Newton-line iteration leads to a linear algebraic system that is solved by GMRES method with block ILU(0) preconditioner.

Barenblatt Problems

We consider two nonlinear variants of (3) following from the Barenblatt problem [4] where the analytical solution exists. The first variant reads94 ∂tϑ(u)-Δu=0,ϑ(u)=u1/m,m∈(0,1),

where the analytical solution is95 u(x1,x2,t)=11+t(⌊[1-m-14m2x12+x22(1+t)1/m⌋+)mm-1,⌊z⌋+:=max(z,0),z∈R

Using the substitution v:=u1/m, we have the second variant96 ∂tv-∇·(m|v|m-1∇v)=0,m>1,

having the solution97 v(x1,x2,t)={11+t(⌊1-m-14m2x12+x22(1+t)1/m⌋+)mm-1}1/m.

For both problems ((94)–(95) and (96)–(97)), the computational domain is Ω=(-6,6)2 and the Dirichlet boundary condition is prescribed on all boundaries by (95) or (97). The final time is T=1.

We carried out computation using a sequence of uniform triangular grids (having 288, 1152, 4608 and 18432 triangles) with several combinations of polynomial approximation degrees with respect to space (p) and time (q). The time step has been chosen constant τ=0.01. Besides the error quantities (R~(uhτ) and J(uhτ)) and its estimators η, ηR:=∑K,mηR,K,m, ηS:=∑K,mηS,K,m and ηT:=∑K,mηT,K,m, we evaluate the effectivity indices98 ieff=ηR~(uhτ),iefftot=η2+J(uhτ)1/2E~(uhτ).

In addition, we present the experimental order of convergence (EoC) of the errors and the estimators for each pair of successive meshes.

Tables 1–4 show the results for both Barenblatt problems ((94)–(95) with m=0.25 and (96)–(97) with m=2) with two variants of the scaling parameter dK,m, K∈Thm, m=1,⋯,r given by (21a) and (21b). The quantity #DoF represents the number of degrees of freedom in the space, that is, #DoF=dimShp,m, m=1,⋯,r. We observe a good correspondence between R~(uhτ) and η, the effectivity index ieff varies between 1 and 2.5 for all tested values of p and q and both variants of dK,m ((21a) and (21b)).Table 1 Barenblatt problem (94)–(95), m=0.25, scaling parameter dK,m given by (21a), approximation of the error and the error estimators, EOC in parenthesis

h	#DoF	R~(uhτ)	η	J(uhτ)	ηR	ηS	ηT	ieff	iefftot	
p=1 & q=1	
1.41	864	8.42×10-4	1.46×10-3	4.01×10-3	7.67×10-6	1.22×10-3	7.94×10-4	1.73	1.13	
0.71	3456	7.31×10-4	1.29×10-3	3.69×10-3	7.68×10-6	1.16×10-3	5.38×10-4	1.76	1.13	
(0.20)	(0.18)	(0.12)	(0.00)	(0.07)	(0.56)	
0.35	13824	4.98×10-4	1.04×10-3	2.95×10-3	8.56×10-6	1.02×10-3	1.55×10-4	2.09	1.16	
(0.55)	(0.30)	(0.33)	(-0.16)	(0.18)	(1.80)	
0.18	55296	4.40×10-4	1.01×10-3	2.81×10-3	9.00×10-6	1.01×10-3	3.22×10-5	2.31	1.18	
(0.18)	(0.04)	(0.07)	(-0.07)	(0.02)	(2.26)	
p=2 & q=2	
1.41	1728	2.06×10-4	3.49×10-4	1.32×10-3	5.60×10-7	3.17×10-4	1.46×10-4	1.70	1.09	
0.71	6912	1.16×10-4	1.86×10-4	7.91×10-4	4.59×10-8	1.74×10-4	6.64×10-5	1.60	1.08	
(0.82)	(0.91)	(0.74)	(3.61)	(0.86)	(1.14)	
0.35	27648	6.51×10-5	8.69×10-5	4.33×10-4	5.46×10-8	8.57×10-5	1.44×10-5	1.34	1.04	
(0.84)	(1.10)	(0.87)	(-0.25)	(1.02)	(2.20)	
0.18	110592	3.51×10-5	4.34×10-5	2.28×10-4	6.54×10-8	4.32×10-5	2.46×10-6	1.23	1.03	
(0.89)	(1.00)	(0.92)	(-0.26)	(0.99)	(2.56)	
p=3 & q=2	
1.41	2880	6.32×10-5	1.16×10-4	3.74×10-4	4.39×10-8	9.31×10-5	6.83×10-5	1.83	1.12	
0.71	11520	2.03×10-5	3.02×10-5	1.45×10-4	3.48×10-8	2.82×10-5	1.08×10-5	1.49	1.06	
(1.64)	(1.94)	(1.37)	(0.33)	(1.73)	(2.66)	
0.35	46080	5.77×10-6	8.16×10-6	4.11×10-5	5.31×10-8	8.03×10-6	1.20×10-6	1.41	1.05	
(1.81)	(1.89)	(1.81)	(-0.61)	(1.81)	(3.17)	
0.18	184320	1.35×10-6	2.09×10-6	9.47×10-6	6.53×10-8	2.04×10-6	8.52×10-8	1.55	1.07	
(2.10)	(1.97)	(2.12)	(-0.30)	(1.98)	(3.82)	

Table 2 Barenblatt problem (94)–(95), m=0.25, scaling parameter dK,m given by (21b), approximation of the error and the error estimators, EOC in parenthesis

h	#DoF	R~(uhτ)	η	J(uhτ)	ηR	ηS	ηT	ieff	iefftot	
p=1 & q=1	
1.41	864	2.34×10-1	4.06×10-1	1.12	2.14×10-3	3.39×10-1	2.21×10-1	1.73	1.13	
0.71	3456	1.02×10-1	1.79×10-1	5.14×10-1	1.07×10-3	1.62×10-1	7.50×10-2	1.76	1.13	
(1.20)	(1.18)	(1.12)	(1.00)	(1.07)	(1.56)	
0.35	13824	3.47×10-2	7.26×10-2	2.05×10-1	5.96×10-4	7.14×10-2	1.08×10-2	2.09	1.16	
(1.55)	(1.30)	(1.33)	(0.84)	(1.18)	(2.80)	
0.18	55296	1.53×10-2	3.54×10-2	9.80×10-2	3.14×10-4	3.51×10-2	1.12×10-3	2.31	1.18	
(1.18)	(1.04)	(1.07)	(0.93)	(1.02)	(3.26)	
0.09	221184	7.63×10-3	1.76×10-2	4.85×10-2	1.59×10-4	1.75×10-2	1.44×10-4	2.31	1.18	
(1.01)	(1.00)	(1.01)	(0.98)	(1.00)	(2.97)	
p=2 & q=2	
1.41	1728	5.73×10-2	9.73×10-2	3.68×10-1	1.56×10-4	8.83×10-2	4.08×10-2	1.70	1.09	
0.71	6912	1.62×10-2	2.60×10-2	1.10×10-1	6.40×10-6	2.43×10-2	9.24×10-3	1.60	1.08	
(1.82)	(1.91)	(1.74)	(4.61)	(1.86)	(2.14)	
0.35	27648	4.53×10-3	6.06×10-3	3.02×10-2	3.81×10-6	5.97×10-3	1.01×10-3	1.34	1.04	
(1.84)	(2.10)	(1.87)	(0.75)	(2.02)	(3.20)	
0.18	110292	1.22×10-3	1.51×10-3	7.96×10-3	2.28×10-6	1.51×10-3	8.56×10-5	1.23	1.03	
(1.89)	(2.00)	(1.92)	(0.74)	(1.99)	(3.55)	
p=3 & q=2	
1.41	2880	1.76×10-2	3.22×10-2	1.04×10-1	1.22×10-5	2.59×10-2	1.90×10-2	1.83	1.12	
0.71	11520	2.83×10-3	4.20×10-3	2.02×10-2	4.85×10-6	3.92×10-3	1.51×10-3	1.49	1.06	
(2.64)	(2.94)	(2.37)	(1.33)	(2.73)	(3.66)	
0.35	46080	4.42×10-4	5.68×10-4	2.87×10-3	3.70×10-6	5.59×10-4	8.38×10-5	1.29	1.04	
(2.68)	(2.89)	(2.81)	(0.39)	(2.81)	(4.17)	
0.18	184320	4.71×10-5	7.28×10-5	3.30×10-4	2.27×10-6	7.10×10-5	2.97×10-6	1.55	1.07	
(3.23)	(2.97)	(3.12)	(0.70)	(2.98)	(4.82)	

Table 3 Barenblatt problem (96)–(97), m=2, scaling parameter dK,m given by (21a), approximation of the error and the error estimators, EOC in parenthesis

h	#DoF	R~(uhτ)	η	J(uhτ)	ηR	ηS	ηT	ieff	iefftot	
p=1 & q=1	
1.41	864	3.35×10-3	5.76×10-3	1.59×10-2	6.44×10-7	4.58×10-3	3.50×10-3	1.72	1.13	
0.71	3456	1.71×10-3	3.19×10-3	7.67×10-3	5.90×10-7	2.74×10-3	1.64×10-3	1.87	1.16	
(0.97)	(0.85)	(1.05)	(0.13)	(0.74)	(1.09)	
0.35	13824	1.09×10-3	2.27×10-3	4.57×10-3	4.81×10-7	2.17×10-3	6.69×10-4	2.07	1.21	
(0.64)	(0.49)	(0.75)	(0.29)	(0.34)	(1.29)	
0.18	55296	9.17×10-4	2.04×10-3	3.60×10-3	3.62×10-7	2.02×10-3	2.41×10-4	2.22	1.25	
(0.26)	(0.15)	(0.34)	(0.41)	(0.10)	(1.48)	
p=2 & q=2	
1.41	1728	6.33×10-4	1.27×10-3	4.52×10-3	1.98×10-7	9.68×10-4	8.18×10-4	2.00	1.12	
0.71	6912	3.26×10-4	5.98×10-4	2.20×10-3	1.63×10-7	4.70×10-4	3.71×10-4	1.84	1.11	
(0.96)	(1.08)	(1.04)	(0.28)	(1.04)	(1.14)	
0.35	27648	1.91×10-4	2.99×10-4	1.17×10-3	1.24×10-7	2.60×10-4	1.47×10-4	1.57	1.08	
(0.77)	(1.00)	(0.91)	(0.40)	(0.85)	(1.33)	
0.18	110592	1.23×10-4	1.71×10-4	7.18×10-4	8.96×10-8	1.63×10-4	5.26×10-5	1.39	1.06	
(0.63)	(0.81)	(0.71)	(0.47)	(0.68)	(1.48)	
p=3 & q=2	
1.41	2880	2.57×10-4	5.65×10-4	2.00×10-3	1.26×10-7	3.69×10-4	4.28×10-4	2.20	1.14	
0.71	11520	1.51×10-4	2.70×10-4	1.06×10-3	8.66×10-8	2.07×10-4	1.74×10-4	1.79	1.10	
(0.77)	(1.06)	(0.91)	(0.54)	(0.84)	(1.30)	
0.35	46080	1.00×10-4	1.50×10-4	6.70×10-4	6.12×10-8	1.33×10-4	6.91×10-5	1.50	1.06	
(0.59)	(0.85)	(0.66)	(0.50)	(0.63)	(1.33)	
0.18	184320	6.92×10-5	9.37×10-5	4.48×10-4	4.30×10-8	9.04×10-5	2.46×10-5	1.35	1.05	
(0.53)	(0.68)	(0.58)	(0.51)	(0.56)	(1.49)	

Table 4 Barenblatt problem (96)–(97), m=2, scaling parameter dK,m given by (21b), approximation of the error and the error estimators, EOC in parenthesis

h	#DoF	R~(uhτ)	η	J(uhτ)	ηR	ηS	ηT	ieff	iefftot	
p=1 & q=1	
1.41	864	3.35×10-1	5.76×10-1	1.59	6.44×10-5	4.58×10-1	3.50×10-1	1.72	1.13	
0.71	3456	8.54×10-2	1.60×10-1	3.84×10-1	2.95×10-5	1.37×10-1	8.19×10-2	1.87	1.16	
(1.97)	(1.85)	(2.05)	(1.13)	(1.74)	(2.09)	
0.35	13824	2.74×10-2	5.68×10-2	1.14×10-1	1.20×10-5	5.43×10-2	1.68×10-2	2.07	1.21	
(1.64)	(1.49)	(1.75)	(1.29)	(1.34)	(2.29)	
0.18	55296	1.15×10-2	2.57×10-2	4.53×10-2	4.56×10-6	2.55×10-2	3.03×10-3	2.22	1.25	
(1.25)	(1.15)	(1.34)	(1.40)	(1.09)	(2.47)	
0.09	221184	5.54×10-3	1.27×10-2	2.13×10-2	1.64×10-6	1.26×10-2	6.05×10-4	2.29	1.27	
(1.06)	(1.02)	(1.09)	(1.48)	(1.01)	(2.32)	
p=2 & q=2	
1.41	1728	6.33×10-2	1.27×10-1	4.52×10-1	1.98×10-5	9.68×10-2	8.18×10-2	2.00	1.12	
0.71	6912	1.63×10-2	2.99×10-2	1.10×10-1	8.16×10-6	2.35×10-2	1.85×10-2	1.84	1.11	
(1.96)	(2.08)	(2.04)	(1.28)	(2.04)	(2.14)	
0.35	27648	4.78×10-3	7.49×10-3	2.94×10-2	3.10×10-6	6.52×10-3	3.69×10-3	1.57	1.08	
(1.77)	(2.00)	(1.90)	(1.40)	(1.85)	(2.33)	
0.18	110592	1.55×10-3	2.15×10-3	9.04×10-3	1.13×10-6	2.04×10-3	6.62×10-4	1.39	1.06	
(1.63)	(1.80)	(1.70)	(1.46)	(1.67)	(2.48)	
p=3 & q=2	
1.41	2880	2.57×10-2	5.65×10-2	2.00×10-1	1.26×10-5	3.69×10-2	4.28×10-2	2.20	1.14	
0.71	11520	7.54×10-3	1.35×10-2	5.31×10-2	4.33×10-6	1.03×10-2	8.70×10-3	1.79	1.10	
(1.77)	(2.06)	(1.91)	(1.54)	(1.84)	(2.30)	
0.35	46080	2.51×10-3	3.76×10-3	1.68×10-2	1.53×10-6	3.34×10-3	1.73×10-3	1.50	1.06	
(1.59)	(1.85)	(1.66)	(1.50)	(1.63)	(2.33)	
0.18	184320	8.71×10-4	1.18×10-3	5.64×10-3	5.42×10-7	1.14×10-3	3.09×10-4	1.36	1.05	
(1.53)	(1.67)	(1.57)	(1.50)	(1.55)	(2.48)	

Finally, we note that the experimental orders of convergence EoC in Tables 1–4) of the error R~(uhτ) and its estimate η are O(hp) for the choice (21b) of the scaling parameter dK,m but only O(hp-1) for the choice (21a). This follows from the fact that τm≪hK for the computations of the Barenblatt problem and then the dominant part of dK,m is τm-2TdϑduK,m,∞, cf. (21a), which implies that dK,m=O(h0) (the time step is the same for all computations). The dominant part of the error estimator is ηS,K,m, hence if σhτ-σ(uhτ,∇uhτ)K,m=O(hp) then ηS,K,m=O(hp-1), cf. (42b). Nevertheless, comparing the pairs of Tables 1–2 and Tables 3–4, we found that the effectivity indexes are practically independent of the choice of dK,m.

Justification of the Algebraic Stopping Criterion (91)

We present the numerical study of the stopping criterion (91) which is used in the iterative solution of algebraic systems given by (13). We consider again the Barenblatt problem (94)–(95) with m=0.25 and (96)–(97) with m=2. The user-dependent constant cA in (91) has been chosen as 10-1, 10-2, 10-3 and 10-4. Tables 5 and 6 show the estimators η, J(uhτ), ηalg and ηalg, cf. (90), for selected meshes and polynomial approximation degrees and the scaling parameter dK,m chosen by (21a).

Additionally, we present the total number of steps of the Newton-like solver Nnon, the total number of GMRES iterations Nlin and the computational time in seconds. The computational time has only an informative character.Table 5 Barenblatt problem (94)–(95), m=0.25, scaling parameter dK,m given by (21a), numerical study of the algebraic stopping criterion (91)

cA	η	J(uhτ)	ηalg	ηspa	Nnon	Nlin	time(s)	
h=0.35, p=1 & q=1, #DoF=13824	
1.0×10-1	1.2475×10-3	3.0760×10-3	8.1322×10-4	1.3804×10-2	202	14148	422.1	
1.0×10-2	1.0559×10-3	2.9565×10-3	7.6470×10-5	1.3483×10-2	362	21589	606.8	
1.0×10-3	1.0435×10-3	2.9468×10-3	7.9268×10-6	1.3458×10-2	529	26545	693.6	
1.0×10-4	1.0423×10-3	2.9457×10-3	7.3279×10-7	1.3456×10-2	579	27766	705.1	
h=0.35, p=2 & q=2, #DoF=27648	
1.0×10-1	1.0443×10-4	4.3586×10-4	5.7369×10-5	1.1019×10-3	406	10581	1968.3	
1.0×10-2	8.8249×10-5	4.3375×10-4	6.1984×10-6	1.0961×10-3	536	12059	2119.4	
1.0×10-3	8.7054×10-5	4.3350×10-4	6.0680×10-7	1.0956×10-3	576	12541	2030.1	
1.0×10-4	8.6948×10-5	4.3347×10-4	5.1172×10-8	1.0955×10-3	618	13580	2132.1	
h=0.35, p=3 & q=2, #DoF=46080	
1.0×10-1	9.9098×10-6	4.1201×10-5	5.3825×10-6	1.0693×10-4	534	10480	6610.2	
1.0×10-2	8.2946×10-6	4.1156×10-5	6.0342×10-7	1.0670×10-4	602	11479	6998.3	
1.0×10-3	8.1647×10-6	4.1150×10-5	4.5285×10-8	1.0669×10-4	636	12288	7181.7	
1.0×10-4	8.1577×10-6	4.1150×10-5	4.5439×10-9	1.0669×10-4	668	13178	7566.5	

Table 6 Barenblatt problem (96)–(97), m=2, scaling parameter dK,m given by (21a), numerical study of the algebraic stopping criterion (91)

cA	η	J(uhτ)	ηalg	ηspa	Nnon	Nlin	time(s)	
h=0.35, p=1 & q=1, #DoF=13824	
1.0×10-1	2.3055×10-3	4.6659×10-3	4.0757×10-4	9.3084×10-3	100	2199	224.6	
1.0×10-2	2.2689×10-3	4.5703×10-3	1.3712×10-5	9.2211×10-3	200	4957	284.5	
1.0×10-3	2.2688×10-3	4.5694×10-3	3.5878×10-6	9.2216×10-3	299	7715	352.6	
1.0×10-4	2.2688×10-3	4.5696×10-3	3.9773×10-7	9.2220×10-3	378	9790	413.7	
h=0.35, p=2 & q=2, #DoF=27648	
1.0×10-1	3.0332×10-4	1.1855×10-3	7.0279×10-5	1.7702×10-3	201	4065	1535.3	
1.0×10-2	2.9940×10-4	1.1753×10-3	7.2859×10-6	1.7634×10-3	286	5820	1675.9	
1.0×10-3	2.9916×10-4	1.1747×10-3	8.0717×10-7	1.7619×10-3	393	7803	1759.7	
1.0×10-4	2.9916×10-4	1.1747×10-3	8.9211×10-8	1.7620×10-3	529	10172	1984.6	
h=0.35, p=3 & q=2, #DoF=46080	
1.0×10-1	1.5523×10-4	6.8813×10-4	7.1710×10-5	1.1705×10-3	202	4222	5539.6	
1.0×10-2	1.5037×10-4	6.6961×10-4	5.7493×10-6	1.1586×10-3	316	6538	6068.2	
1.0×10-3	1.5026×10-4	6.6968×10-4	6.3387×10-7	1.1580×10-3	453	8880	6615.7	
1.0×10-4	1.5026×10-4	6.6967×10-4	5.9895×10-8	1.1580×10-3	591	11150	7101.3	

We observe that the error estimators η, J(uhτ) and also ηspa converge to the limit values for decreasing cA in (91) which mimic the case when the algebraic errors are negligible. Moreover, the relative differences between the actual values η and J(uhτ) and their limits correspond more or less to the value of cA. Obviously, smaller values of cA cause prolongation of the computational time, due to a higher number of iterations, with a negligible effect on accuracy. Thus, the choice cA=10-2 seems to be optimal in order to balance accuracy and efficiency.

The presented numerical experiments indicate that the estimator ηspa(uhτ) gives an upper bound of R(uhτ), however, this observation is not supported by the theory. The quantity ηspa(uhτ) is used only in the stopping criterion (91).

Tracy Problem

Tracy problem represents a standard benchmark, where the analytical solutions of the Richards equation are available [35]. We consider the Gardners constitutive relations [26]99 K(u)=Ksexp(-αψ)ifψ>0Ksifψ≤0,ϑ(u)=θr+(θs-θr)exp(-αψ)ifψ>0θsifψ≤0

where ψ=u-z is the pressure head, z is the vertical coordinate and the material parameters Ks=1.2I, θs=0.5, θr=0.0, and α=0.1 are the isotropic conductivity, saturated water content, residual water content, and the soil index parameter related to the pore-size distribution, respectively.

The computational domain is Ω=(0,1)2, the initial condition is set u=ur:=-10 in Ω where ur corresponds to the hydraulic head when the porous medium is dry. On the top part of the boundary Γ1:={(x,z),x∈(0,1),z=1}, we prescribe the boundary condition100 u(x)=1αlog(exp(αur)+(1-exp(αur)sin(πx)),x∈(0,1)

and on the rest of boundary Γ we set u=ur. We note that this benchmark poses an inconsistency between the initial and boundary conditions on Γ1. Hence, the most challenging part is the computation close to t=0. In order to avoid the singularity at t=0, we investigate the error only on the interval t∈[1.0×10-5,1.1×10-4] with the fixed time step τ is 1.0×10-6.

We perform a computation using a sequence of uniform triangular grids with several combinations of polynomial approximation degrees and the choice (21b), the results are shown in Table 7. We observe reasonable values of the effectivity indices except for the finest grids and the higher degrees of polynomial approximation, where the effectivity indices ieff are below 1. Based on the values of EoC, we suppose that ieff below 1 is not caused by the failure of the error estimator but due to an inaccurate approximation R~(uhτ) of the exact error; see Remark 3.Table 7 Tracy problem scaling parameter dK,m given by (21b), approximation of the error and the error estimators, EOC in parenthesis

h	#DoF	R~(uhτ)	η	J(uhτ)	ηR	ηS	ηT	ieff	iefftot	
p=1 & q=1	
0.18	384	2.43×10-1	2.95×10-1	6.62×102	2.25×10-3	2.90×10-1	4.83×10-2	1.22	1.00	
0.09	1536	8.77×10-2	1.19×10-1	2.47×102	9.92×10-4	1.10×10-1	4.36×10-2	1.35	1.00	
(1.47 )	(1.32)	(1.43)	(1.18)	(1.40)	(0.15)	
0.04	6144	1.50×10-2	2.51×10-2	4.68×101	1.33×10-4	2.39×10-2	7.31×10-3	1.67	1.00	
(2.55)	(2.24)	(2.40)	(2.90)	(2.20)	(2.58)	
0.02	24576	7.34×10-3	1.22×10-2	2.27×101	8.11×10-5	1.20×10-2	1.86×10-3	1.66	1.00	
(1.03)	(1.04)	(1.05)	(0.71)	(0.99)	(1.98)	
p=2 & q=2	
0.18	768	4.88×10-2	6.18×10-2	1.92×102	3.35×10-4	6.04×10-2	1.22×10-2	1.27	1.00	
0.09	3072	1.75×10-2	1.98×10-2	6.03×101	1.44×10-5	1.95×10-2	3.35×10-3	1.13	1.00	
(1.48)	(1.65)	(1.67)	(4.54)	(1.63)	(1.86)	
0.04	12288	5.59×10-3	6.28×10-3	1.94×101	4.68×10-6	6.08×10-3	1.56×10-3	1.12	1.00	
(1.64)	(1.65)	(1.64)	(1.62)	(1.68)	(1.11)	
0.02	49,152	1.90×10-3	1.33×10-3	4.37	2.41×10-6	1.31×10-3	1.87×10-4	0.70	1.00	
(1.56)	(2.24)	(2.15)	(0.96)	(2.21)	(3.06)	
p=3 & q=2	
0.18	1280	2.24×10-2	2.64×10-2	9.60×101	2.42×10-5	2.61×10-2	4.48×10-3	1.18	1.00	
0.09	5120	6.26×10-3	7.94×10-3	2.77×101	1.21×10-5	7.13×10-3	3.48×10-3	1.27	1.00	
(1.84)	(1.74)	(1.79)	(1.00)	(1.87)	(0.37)	
0.04	20480	1.40×10-3	4.60×10-4	1.63	2.87×10-6	4.49×10-4	8.90×10-5	0.33	1.00	
(2.16)	(4.11)	(4.08)	(2.08)	(3.99)	(5.29)	
0.02	81920	1.37×10-3	8.85×10-5	3.08×10-1	2.37×10-6	8.59×10-5	1.09×10-5	0.06	1.00	
(0.03)	(2.38)	(2.41)	(0.28)	(2.39)	(3.03)	

Mesh Adaptive Algorithm

We introduce the mesh adaptive algorithm which is based on the a posteriori error estimates η, cf. (41). Let δ>0 be the given tolerance, the goal of the algorithm is to define the sequence of time steps τm, meshes Thm and spaces Shp,m, m=1,⋯,r such that the corresponding approximate solution uhτ∈Shpτq given by (13) satisfies the condition101 η=η(uhτ)≤δ.

Another possibility is to require η2+J(uhτ)1/2≤δ, then the following considerations have to be modified appropriately.

The mesh adaptation strategy is built on the equi-distribution principle, namely the sequences {τm,Thm,Shp,m}m=1r should be generated such that 102a ηm≤δm:=δτm/T∀m=1,⋯,r,

102b ηK,m≤δK,m:=δm1/#Thm∀K∈Thm∀m=1,⋯,r,

where ηm:=(∑K∈ThmηK,m2)1/2 is the error estimate corresponding to the time interval Im, m=1,⋯,r and #Thm denotes the number of elements of Thm. Obviously, if all the conditions in (102) are valid, then the criterion (101) is achieved.

Algorithm 1 Space-time mesh adaptive algorithm.

Based on (101)–(102), we introduce the abstract Algorithm 1. The size of τm, m=1,⋯,r (step 8 of the algorithm) are chosen to equilibrate estimates of the spatial and temporal reconstruction, ηS,m:=(∑K∈Thm(ηS,K,m)2)1/2 and ηT,m:=(∑K∈Thm(ηT,K,m)2)1/2, cf. (42). Particularly, we set the new time step according to the formula103 τm+1=τmcFηS,mηT,m1/(q+1),m=1,⋯,r,

where cF∈(0,1) is the security factor and q≥0 is the polynomial degree with respect to time. Therefore, q+1 corresponds to the temporal order of convergence.

The construction of the new mesh (step 11 in Algorithm 1) is based on the modification of the anisotropic hp-mesh adaptation method from [15, 20]. Having the actual mesh Thm, for each K∈Thm we set the new volume of K according the formula104 νK=|K|Λ(δK,m/ηK,m),K∈Thm,

where δK,m is the local tolerance from (102b), |K| is the volume of |K| and Λ:R+→R+ is a suitable increasing function such that Λ(1)=1. For particular variants of Λ, we refer to [15, 20].

When the new volume of mesh elements is established by (104), the new shape of K and a new polynomial approximation degree pK are optimized by minimizing the interpolation error. This optimization is done locally for each mesh element. In one adaptation level, we admit the increase or decrease of pK by one. Setting the new area, shape, and polynomial approximation degree for each element of the current mesh, we define the continuous mesh model [16] and carry out a remeshing using the code ANGENER [9].

The generated meshes are completely non-nested and non-matching, hence the evaluation of the time-penalty term (cf. Remark 1) is delicate. We refer to [20] where this aspect is described in detail and numerically verified. The presented numerical analysis takes into account the errors arising from the re-meshing in the temporal reconstruction Rhτ, which contains term {ϑ(uhτ)}m-1, cf. (26). The following numerical experiments show that the error estimator is under the control also after each remeshing.

Barenblatt Problem

We apply Algorithm 1 to the Barenblatt problem (96) with m=2. Table 8 shows the error estimators obtained by adaptive computation for three different tolerances δ. Compared with the error estimators from Table 4, we observe that the adaptive computations achieve significantly smaller error estimates using a significantly smaller number of degrees of freedom. We note that we are not able to present the quantity R~ (cf. (92)–(93)) approximating the error since the finite element code used for the evaluation of R~ supports only uniform grids.Table 8 Barenblatt problem (96)–(97), scaling parameter dK,m given by (21b), the error estimators obtained by the adaptive computations using Algorithm 1

hp adaptation	
δ	#DoF	η	J(uhτ)1/2	ηR	ηS	ηT	
2.0E−03	4 543	7.82×10-4	1.69×10-3	2.18×10-4	5.40×10-4	3.23×10-4	
1.0E−03	6 244	4.57×10-4	1.17×10-3	1.48×10-4	3.13×10-4	1.43×10-4	
5.0E−04	9 071	2.10×10-4	7.02×10-4	6.75×10-5	1.38×10-4	7.79×10-5	
The quantity #DoF is the average number of space degrees of freedom per one time step

Figure 1 shows the performance of Algorithm 1, where each dot corresponds to one time step m=1,⋯,r. We plot the values of the accumulated estimators η¯m=∑i=1mηi for all m=1,⋯,r. The red nodes correspond to all computed time steps, including the rejected ones (steps 11–12 of Algorithm 1) whereas the blue nodes mark only the accepted time steps. The rejected time steps indicate the re-meshing. Moreover, we plot the “accumulated” tolerance δ(tm/T)1/2, cf. (101) and (102a). We observe that the resulting estimator η at t=T is below the tolerance δ by a factor of approximately 2.5 since conditions (102) are stronger than (101).Fig. 1 Barenblatt problem, (96)–(97), m=2, performance of Algorithm 1, accumulated error estimator η¯m and the “accumulated” tolerance δ(tm/T)1/2 for m=1,⋯,r

Figure 2, left, shows the hp-mesh obtained by Algorithm 1 at the final time T=1, each triangle is highlighted by a color corresponding to the polynomial degree used pK, K∈Thm. We observe a strong anisotropic refinement about the circular singularity of the solution when u→0+, see the analytical formula (97). Outside of this circle, large triangles with the smallest polynomial degree (p=1) are generated. On the other hand, due to the regularity of the solution in the interior of the circle, the polynomial degrees p=2 or p=3 are generated.Fig. 2 Barenblatt problem, hp-mesh obtained by Algorithm 1 (left) and the error estimators ηK,m, K∈Thm at T=1

Moreover, Fig. 2, right, shows the error estimator ηK,m, K∈Thm at T=1. The elements in the exterior of the circle have small values of ηK,m≈10-17 – 10-14 due to a constant solution and negligible errors. On the other hand, the values of ηK,m for the rest of elements K∈Thm are in the range 10-13–10-11 due to the equidistant principle used.

Single Ring Infiltration

We deal with the numerical solution of the single ring infiltration experiment, which is frequently used for the identification of saturated hydraulic conductivity, cf. [32, 39] for example. We consider the Richards equation (3) where the active pore volume ϑ is given by (2), the water content function θ is given by the van Genuchten’s law [27] and the conductivity K(u)=KsKr(u) is given by the Mualem function [31], namely105 θ(u)=θs-θr(1+-αψ)nm+θrforψ<0,θsforψ≥0,Kr(u)=1-(-αψ)mn(1+(-αψ)n)-m21+(-αψ)nm/2forψ<0,1forψ≥0,

where ψ=u-z is the pressure head, z is the vertical coordinate and the material parameters Ks=0.048Im·hours-1, θs=0.55, θr=0.0, α=0.8m-1, n=1.2, m=1/6 and Ss=10-3m-1 (cf. (2)).

The computational domain together with the boundary parts is sketched in Fig. 3a. On the boundary part ΓD we set the Dirichlet boundary condition u=1.05m, and on ΓN=Γ\ΓD we consider the homogeneous Neumann boundary condition. The smaller “magenta” vertical lines starting at ΓD belong to ΓN. At t=0, a dry medium with u=ψ+z=-2m is prescribed. We carried out the computation until the physical time T=2hours. The inconsistency of the initial and boundary condition on ΓD makes the computation quite difficult for t≈0.Fig. 3 Single ring infiltration problem

Figure 3b verifies the conservativity of the adaptive method. We plot the quantities106 F(t)=∫0t∫ΓK(u)∇u·ndSdt,ΔV(t)=V(t)-V(0),V(t)=∫Ωϑ(u(·,t))dx,t∈[0,T],

where F(t) is the total flux of the water through the boundary Γ till time t and ΔV(t) is the changes of the water content in the domain between times 0 and t. From equation (3) and the Stokes theorem, we have the conservation law F(t)=ΔV(t) for all t∈[0,T]. Therefore, we also show the relative difference between these quantities |F(t)-ΔV(t)|/ΔV(t) for t>0 in Fig. 3b the vertical label on the right. We observe that, except for the time close to zero, where the inconsistency between initial and boundary conditions is problematic, the relative difference is at the level of several percent.

Furthermore, Fig. 4 shows the accumulated estimators η¯m=∑i=1mηi for time levels tm, m=1,⋯,m. The red nodes correspond to all computed time steps, including the rejected steps whereas the blue line connects only the accepted time steps. The rejected time steps are followed by the remeshing which is carried out namely for small t. We observe that the elimination of the rejected time steps causes that the errors arising from the remeshing do not essentially affect the total error estimate η.Fig. 4 Single ring infiltration, performance of Algorithm 1, accumulated error estimator η¯m with respect to tm, m=1,⋯,r

Moreover, Fig. 5 shows the hp-meshes, the hydraulic head and the error estimator ηK,m, K∈Thm at selected time levels obtained from Algorithm 1 with δ=5.0×10-3. We observe the mesh adaptation namely at the (not sharp) interface between the saturated and non-saturated medium and also in the vicinity of the domain singularities. The error estimators ηK,m, K∈Thm indicate an equi-distribution of the error.Fig. 5 Single ring infiltration, hydraulic head (top ), the corresponding hp-meshes obtained by Algorithm 1 (center) and the error estimators ηK,m, K∈Thm (bottom) at t=0.4, t=0.8 and t=2 hours (from left to right)

Conclusion

We derived reliable and efficient a posteriori error estimates in the residual-based norm for the Richards equation discretized by the space-time discontinuous Galerkin method. The numerical verification indicates the effectivity indexes between 1 and 2.5 for the tested examples. Moreover, we introduced the hp-mesh adaptive method handling varying non-nested and non-matching meshes and demonstrated its efficiency for simple test benchmark and its applicability for the numerical solution of the single ring infiltration experiment.

It will be possible to generalize the presented approach to genuinely space-time hp-adaptive method, where the (local) polynomial order q in time is varied as well. However, the question is of potential benefit. Based on our experience, the setting q=1 gives sufficiently accurate approximation for the majority of tested problems.

On the other hand, the choice q=0 would be sufficient only in subdomains of Ω where the solution is almost constant in time. Therefore, we suppose that the benefit of local varying of polynomial order in time will be low.

Although the presented numerical examples are two-dimensional, it would be possible to apply the presented error estimates and mesh adaptation to three-dimensional problems as well. We refer, e.g., to [1] and the references therein, where the anisotropic mesh adaptation techniques are developed for time-dependent 3D problems.

Funding

Open access publishing supported by the National Technical Library in Prague. This work has been supported by the Czech Science Foundation Grant No. 20-01074 S (V.D.), the Charles University grant SVV-2023-260711, and the Grant Agency of Charles University Project No. 28122 (H.S.), European Development Fund-Project “Center for Advanced Aplied Science” No. CZ.02.1.01/0.0/0.0./16 019/0000778 (M.V.). V.D. acknowledges the membership in the Nečas Center for Mathematical Modeling ncmm.karlin.mff.cuni.cz.

Data Availability

No datasets were generated or analysed during the current study.

Declarations

Conflict of interest

The authors declare that they have no conflict of interest.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
==== Refs
References

1. Alauzet F Loseille A Olivier G Time-accurate multi-scale anisotropic mesh adaptation for unsteady flows in CFD J. Comput. Phys. 2018 373 28 63 10.1016/j.jcp.2018.06.043
Alauzet, F., Loseille, A., Olivier, G.: Time-accurate multi-scale anisotropic mesh adaptation for unsteady flows in CFD. J. Comput. Phys. 373, 28–63 (2018)
2. Alt HW Luckhaus S Quasilinear elliptic-parabolic differential equations Math. Z. 1983 183 3 311 341 10.1007/BF01176474
Alt, H.W., Luckhaus, S.: Quasilinear elliptic-parabolic differential equations. Math. Z. 183(3), 311–341 (1983)
3. Baradi M Difonzo FV Strong solutions for Richards’ equation with Cauchy conditions and constant pressure Environ. Fluid Mech. 2020 20 165 174 10.1007/s10652-019-09705-w
Baradi, M., Difonzo, F.V.: Strong solutions for Richards’ equation with Cauchy conditions and constant pressure. Environ. Fluid Mech. 20, 165–174 (2020)
4. Barenblatt GI On some unsteady motions of a liquid and gas in a porous medium Akad. Nauk SSSR. Prikl. Mat. Meh. 1952 16 67 78
Barenblatt, G.I.: On some unsteady motions of a liquid and gas in a porous medium. Akad. Nauk SSSR. Prikl. Mat. Meh. 16, 67–78 (1952)
5. Baron V Coudière Y Sochala P Adaptive multistep time discretization and linearization based on a posteriori error estimates for the Richards equation Appl. Numer. Math. 2017 112 104 125 10.1016/j.apnum.2016.10.005
Baron, V., Coudière, Y., Sochala, P.: Adaptive multistep time discretization and linearization based on a posteriori error estimates for the Richards equation. Appl. Numer. Math. 112, 104–125 (2017)
6. Bernardi C El Alaoui L Mghazli Z A posteriori analysis of a space and time discretization of a nonlinear model for the flow in partially saturated porous media IMA J. Numer. Anal. 2014 34 3 1002 1036 10.1093/imanum/drt014
Bernardi, C., El Alaoui, L., Mghazli, Z.: A posteriori analysis of a space and time discretization of a nonlinear model for the flow in partially saturated porous media. IMA J. Numer. Anal. 34(3), 1002–1036 (2014)
7. Brezzi F Fortin M Mixed and Hybrid Finite Element Methods, Springer Series in Computational Mathematics 1991 New York Springer
Brezzi, F., Fortin, M.: Mixed and Hybrid Finite Element Methods, Springer Series in Computational Mathematics, vol. 15. Springer, New York (1991)
8. Di Pietro DA Vohralík M Yousef S An a posteriori-based, fully adaptive algorithm with adaptive stopping criteria and mesh refinement for thermal multiphase compositional flows in porous media Compu. Math. Appl. 2014 68 12, Part B 2331 2347 10.1016/j.camwa.2014.08.008
Di Pietro, D.A., Vohralík, M., Yousef, S.: An a posteriori-based, fully adaptive algorithm with adaptive stopping criteria and mesh refinement for thermal multiphase compositional flows in porous media. Compu. Math. Appl. 68(12, Part B), 2331–2347 (2014)
9. Dolejší, V.: ANGENER – Anisotropic mesh generator, in-house code. Charles University, Prague, Faculty of Mathematics and Physics (2000). https://msekce.karlin.mff.cuni.cz/~dolejsi/angen/
10. Dolejší, V.: ADGFEM – Adaptive discontinuous Galerkin finite element method, in-house code. Charles University, Prague, Faculty of Mathematics and Physics (2020). https://msekce.karlin.mff.cuni.cz/~dolejsi/adgfem/
11. Dolejší V Ern A Vohralík M A framework for robust a posteriori error control in unsteady nonlinear advection-diffusion problems SIAM J. Numer. Anal. 2013 51 2 773 793 10.1137/110859282
Dolejší, V., Ern, A., Vohralík, M.: A framework for robust a posteriori error control in unsteady nonlinear advection-diffusion problems. SIAM J. Numer. Anal. 51(2), 773–793 (2013)
12. Dolejší V Ern A Vohralík M hp-adaptation driven by polynomial-degree-robust a posteriori error estimates for elliptic problems SIAM J. Sci. Comput. 2016 38 5 A3220 A3246 10.1137/15M1026687
Dolejší, V., Ern, A., Vohralík, M.: -adaptation driven by polynomial-degree-robust a posteriori error estimates for elliptic problems. SIAM J. Sci. Comput. 38(5), A3220–A3246 (2016)
13. Dolejší V Feistauer M Discontinuous Galerkin Method: Analysis and Applications to Compressible Flow. Springer Series in Computational Mathematics 2015 Cham Springer
Dolejší, V., Feistauer, M.: Discontinuous Galerkin Method: Analysis and Applications to Compressible Flow. Springer Series in Computational Mathematics, vol. 48. Springer, Cham (2015)
14. Dolejší V Kuráž M Solin P Adaptive higher-order space-time discontinuous Galerkin method for the computer simulation of variably-saturated porous media flows Appl. Math. Model. 2019 72 276 305 10.1016/j.apm.2019.02.037
Dolejší, V., Kuráž, M., Solin, P.: Adaptive higher-order space-time discontinuous Galerkin method for the computer simulation of variably-saturated porous media flows. Appl. Math. Model. 72, 276–305 (2019)
15. Dolejší V May G Anisotropic hp-Mesh Adaptation Methods 2022 Cham Birkhäuser
Dolejší, V., May, G.: Anisotropic -Mesh Adaptation Methods. Birkhäuser, Cham (2022)
16. Dolejší V May G Rangarajan A A continuous hp-mesh model for adaptive discontinuous Galerkin schemes Appl. Numer. Math. 2018 124 1 21 10.1016/j.apnum.2017.09.015
Dolejší, V., May, G., Rangarajan, A.: A continuous -mesh model for adaptive discontinuous Galerkin schemes. Appl. Numer. Math. 124, 1–21 (2018)
17. Dolejší V Roskovec F Vlasák M Residual based error estimates for the space-time discontinuous Galerkin method applied to the compressible flows Comput. Fluids 2015 117 304 324 10.1016/j.compfluid.2015.05.027
Dolejší, V., Roskovec, F., Vlasák, M.: Residual based error estimates for the space-time discontinuous Galerkin method applied to the compressible flows. Comput. Fluids 117, 304–324 (2015)
18. Dolejší V Roskovec F Vlasák M A posteriori error estimates for higher order space-time Galerkin discretizations of nonlinear parabolic problems SIAM J. Numer. Anal. 2021 59 3 1486 1509 10.1137/18M117594X
Dolejší, V., Roskovec, F., Vlasák, M.: A posteriori error estimates for higher order space-time Galerkin discretizations of nonlinear parabolic problems. SIAM J. Numer. Anal. 59(3), 1486–1509 (2021)
19. Dolejší, V., Shin, H.G.: A posteriori error estimate and mesh adaptation for the numerical solution of the richards equation. In: J.M. Melenk, et al. (eds.) Spectral and High Order Methods for Partial Differential Equations ICOSAHOM 2020+1, Lecture Notes in Computational Science and Engineering 137 (2023)
20. Dolejší, V., May, G.: An anisotropic -mesh adaptation method for time-dependent problems based on interpolation error control. J. Sci. Comput. 95(2) (2023)
21. Ern A Guermond JL Theory and Practice of Finite Elements 2004 Berlin Springer
Ern, A., Guermond, J.L.: Theory and Practice of Finite Elements. Springer, Berlin (2004)
22. Ern A Smears I Vohralík M Guaranteed, locally space-time efficient, and polynomial-degree robust a posteriori error estimates for high-order discretizations of parabolic problems SIAM J. Numer. Anal. 2017 55 6 2811 2834 10.1137/16M1097626
Ern, A., Smears, I., Vohralík, M.: Guaranteed, locally space-time efficient, and polynomial-degree robust a posteriori error estimates for high-order discretizations of parabolic problems. SIAM J. Numer. Anal. 55(6), 2811–2834 (2017)
23. Ern A Vohralík M Adaptive inexact Newton methods with a posteriori stopping criteria for nonlinear diffusion PDEs SIAM J. Sci. Comput. 2013 35 4 A1761 A1791 10.1137/120896918
Ern, A., Vohralík, M.: Adaptive inexact Newton methods with a posteriori stopping criteria for nonlinear diffusion PDEs. SIAM J. Sci. Comput. 35(4), A1761–A1791 (2013)
24. Ern A Vohralík M Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations SIAM J. Numer. Anal. 2015 53 2 1058 1081 10.1137/130950100
Ern, A., Vohralík, M.: Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations. SIAM J. Numer. Anal. 53(2), 1058–1081 (2015)
25. Farthing M Ogden F Numerical solution of Richards’ equation: a review of advances and challenges Soil Sci. Soc. Am. J. 2017 81 6 1257 1269 10.2136/sssaj2017.02.0058
Farthing, M., Ogden, F.: Numerical solution of Richards’ equation: a review of advances and challenges. Soil Sci. Soc. Am. J. 81(6), 1257–1269 (2017)
26. Gardner WR Some steady state solutions of the unsaturated moisture flow equation with application to evaporation from a water table Soil Sci. 1958 85 228 232 10.1097/00010694-195804000-00006
Gardner, W.R.: Some steady state solutions of the unsaturated moisture flow equation with application to evaporation from a water table. Soil Sci. 85, 228–232 (1958)
27. van Genuchten MT Closed-form equation for predicting the hydraulic conductivity of unsaturated soils Soil Sci. Soc. Am. J. 1980 44 5 892 898 10.2136/sssaj1980.03615995004400050002x
van Genuchten, M.T.: Closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44(5), 892–898 (1980)
28. Kim IC Požár N Nonlinear elliptic-parabolic problems Arch. Ration. Mech. Anal. 2013 210 975 1020 10.1007/s00205-013-0663-3
Kim, I.C., Požár, N.: Nonlinear elliptic-parabolic problems. Arch. Ration. Mech. Anal. 210, 975–1020 (2013)
29. Mallik, G., Vohralík, M., Yousef, S.: Goal-oriented a posteriori error estimation for conforming and nonconforming approximations with inexact solvers. J. Comput. Appl. Math. 366 (2020)
30. Mitra, K., Vohralík, M.: A posteriori error estimates for the Richards equation. Tech. Rep. hal-03328944v2, INRIA (2022)
31. Mualem Y A new model for predicting the hydraulic conductivity of unsaturated porous media Water Resour. Res. 1976 12 3 513 522 10.1029/WR012i003p00513
Mualem, Y.: A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12(3), 513–522 (1976)
32. Nakhaei, M., Šimnek, J.: Parameter estimation of soil hydraulic and thermal property functions for unsaturated porous media using the hydrus-2d code. J. Hydrol. Hydromech. 62(1) (2014)
33. Richards LA Capillary conduction of liquids through porous mediums J. Appl. Phys. 1931 1 5 318 333
Richards, L.A.: Capillary conduction of liquids through porous mediums. J. Appl. Phys. 1(5), 318–333 (1931)
34. Rudin W Real and Complex Analysis 1987 3 New York McGraw-Hill
Rudin, W.: Real and Complex Analysis, 3rd edn. McGraw-Hill, New York (1987)
35. Tracy, F.T.: Clean two- and three-dimensional analytical solutions of Richards equation for testing numerical solvers. Water Resour Res 42(8) (2006)
36. Verfürth R A Posteriori Error Estimation Techniques for Finite Element Methods. Numerical Mathematics and Scientific Computation 2013 Oxford Oxford University Press
Verfürth, R.: A Posteriori Error Estimation Techniques for Finite Element Methods. Numerical Mathematics and Scientific Computation. Oxford University Press, Oxford (2013)
37. Vohralík, M.: A posteriori error estimates for efficiency and error control in numerical simulations. Universite Pierre et Marie Curie – Paris 6 (2018). Lecture notes Course NM497
38. Vohralík M Yousef S A simple a posteriori estimate on general polytopal meshes with applications to complex porous media flows Comput. Methods Appl. Mech. Eng. 2018 331 728 760 10.1016/j.cma.2017.11.027
Vohralík, M., Yousef, S.: A simple a posteriori estimate on general polytopal meshes with applications to complex porous media flows. Comput. Methods Appl. Mech. Eng. 331, 728–760 (2018)
39. Xu X Lewis C Liu W Albertson J Kiely G Analysis of single-ring infiltrometer data for soil hydraulic properties estimation: comparison of BEST and Wu methods Agric. Water Manag. 2012 107 34 41 10.1016/j.agwat.2012.01.004
Xu, X., Lewis, C., Liu, W., Albertson, J., Kiely, G.: Analysis of single-ring infiltrometer data for soil hydraulic properties estimation: comparison of BEST and Wu methods. Agric. Water Manag. 107, 34–41 (2012)
