Acta Polytechnica https://doi.org/10.14311/AP.2021.61.0059 Acta Polytechnica 61(SI):59–67, 2021 © 2021 The Author(s). Licensed under a CC-BY 4.0 licence Published by the Czech Technical University in Prague NUMERICAL SIMULATION OF A SINGLE RING INFILTRATION EXPERIMENT WITH hp-ADAPTIVE SPACE-TIME DISCONTINUOUS GALERKIN METHOD Vít Dolejšíb, Michal Kurážc, Pavel Solina,∗ a University of Nevada, College of Science, Department of Mathematics and Statistics, 1664 N Virginia St, NV 89557 Reno, USA b Charles University, Faculty of Mathematics and Physics, Sokolovská 83, 186 75 Prague, Czech Republic c Czech University of Life Sciences Prague, Faculty of Environmental Sciences, Department of Water Resources and Environmental Modeling, Kamýcká 129, 165 00 Prague, Czech Republic ∗ corresponding author: solin@unr.edu Abstract. We present a novel hp-adaptive space-time discontinuous Galerkin (hp-STDG) method for the numerical solution of the nonstationary Richards equation equipped with Dirichlet, Neumann and seepage face boundary conditions. The hp-STDG method presented in this paper is a generalization of a hp-STDG method which was developed for time dependent non-linear convective-diffusive problems. We describe the method and the single ring experiment, and then we present a numerical experiment which clearly demonstrates the superiority of the hp-STDG method over a discontinuous Galerkin method based on a static fine mesh. Keywords: Richards equation, porous media flow, space-time discontinuous Galerkin method. 1. Introduction The single ring (hereafter SR) infiltration experiment is a standard and robust dynamic field experiment. The steady part of this experiment is traditionally used for the identification of saturated hydraulic con- ductivity. More recently, the unsteady part of this experiment was used to analyze various hydraulic properties of soils such as van Genuchten’s α (see, for example, [1] and [2]). The SR experiment can be mathematically de- scribed by the Richards equation [3], which is a govern- ing equation for many different applications ranging from hydrology, soil and snow physics [4–6], contami- nant transport simulations [7], etc. The Richards equation belongs to the class of quasi- linear degenerate parabolic equations since the space and time derivative terms can vanish or blow up. The analytical results concerning the existence and unique- ness were published, e.g., in [8–10]. First attempts of solving the Richards equation numerically date back to [11] in the 1970s. Various methods for the nu- merical solution of the Richards equation have been studied and improved in the last decades [12–20]. There are two main reasons why hp-adaptive dis- continuous Galerkin methods are naturally suitable for the numerical simulation of the SR experiment: (1.) The flow pattern is typically formed out of a mov- ing wetting front. This wetting front is the source of considerable computational difficulties due to steep gradients in the solution and abrupt changes in nonlinear functions. (2.) The flow domain is effectively discontinuous, since the thickness of the steel ring (≈ 1 mm) is negligible compared to the scale of the flow pattern, as shown in Fig. 1. In this paper we use the discontinuous Galerkin (DG) method to solve the nonstationary Richards equation describing the SR experiment. The Richards equation is equipped with both Dirichlet and Neu- mann boundary conditions (BCs), as well as with a seepage face boundary condition (BC). We combine the DG method with an anisotropic mesh adapta- tion algorithm and present a simulation of the SR experiment which demonstrates the full potential of the proposed numerical treatment for this class of problems. 2. Governing equations We will model the single ring infiltration experiment using the nonstationary planar Richards equation in the form C(ψ)∂Ψ ∂t −∇ · (K(ψ)∇Ψ) = 0, (2.1) where ψ denotes the pressure head [L] and Ψ is the hydraulic head [L]. They are connected by the simple relation Ψ = ψ + z, (2.2) where z is the vertical distance from the reference level (the geodetic head [L]). Moreover, the ten- sor function K(ψ) is the unsaturated hydraulic con- ductivity [L.T−1], and it is defined as the product 59 https://doi.org/10.14311/AP.2021.61.0059 https://creativecommons.org/licenses/by/4.0/ https://www.cvut.cz/en V. Dolejší, M. Kuráž, P. Solin Acta Polytechnica Figure 1. Scheme of the single ring infiltration experiment. of the the relative hydraulic conductivity Kr(ψ) [-] and the (tensor value) saturated hydraulic conductiv- ity Ks [L.T−1], K(ψ) = Kr(ψ)Ks. (2.3) Furthermore, the scalar function C(ψ) denotes the water retention capacity [L−1], expressed in the form C(ψ) = dθ(ψ) dψ + θ(ψ) θS Ss, (2.4) where the symbol θ(ψ) denotes the water content function [-], and θS is the limited saturated water content [-] and Ss is the specific storativity (aquifer storage) [L−1]. The function θ(ψ) is given by the van Genuchten’s law [21], which can be written in the form θ(ψ) =  θs − θr (1 + (−αψ)n)m + θr, for ψ < 0, θs for ψ ≥ 0, (2.5) where θr[-] is the residual water content, θs[-] is the saturated water content, n [-] and m [-] are pore size distribution parameters, and α [L−1] is given as the inverse of the air entry value. The relative conductivity Kr(ψ) appearing in (2.3) is given by the Mualem function [22] which in combi- nation with (2.5) reads Kr(ψ) = { (1−(−αψ)m n(1+(−αψ)n)−m)2 (1+(−αψ)n)m/2 for ψ < 0, 1 for ψ ≥ 0, (2.6) In order to close the Richard’s equation (2.1) we add the usual Dirichlet and Neumann BCs. Moreover, in order to simulate the seepage of the fluid from the medium through its boundary, we consider also the seepage face boundary condition which reads ψ := Ψ− z ≤ 0, −K(Ψ− z)∇Ψ · n ≥ 0, (2.7) ψ(∇Ψ · n) = 0. The first condition in (2.7) means that the pressure head ψ must be non-positive and the second one means that the fluid cannot flow into the medium because there is no fluid outside the medium. Finally, the last relation in (2.7) gives that the fluid leaves the medium if and only if the medium becomes saturated (ψ = 0). For the purpose of numerical implementation we have to reformulate (2.7). According to [23], we treat the seepage BC as a Signorini condition. This means that if the medium on the boundary is unsaturated (ψ < 0) then this boundary behaves as a no-flow (homogeneous Neumann) boundary. On the other hand, if the medium on the boundary is saturated then the boundary changes into a Dirichlet one with a zero pressure head. These conditions can be then written as ∇Ψ · n = 0, if ψ < 0, (2.8a) ψ = 0, if ∇Ψ · n < 0. (2.8b) 2.1. Initial boundary value problem In order to formulate (2.1) in a form more suitable for numerical discretization, and taking into account (2.2), we write the first term in (2.1) as C(ψ)∂Ψ ∂t = C(Ψ− z)∂Ψ ∂t = C(Ψ− z)∂(Ψ− z) ∂t = ∂ϑ(Ψ− z) ∂t = ∂ϑ(ψ) ∂t . (2.9) Here, in virtue of (2.4), we define the quantity ϑ(ψ) := θ(ψ) + Ss θS ∫ ψ −∞ θ(s)ds (2.10) which represents the active pore volume affected by the changing pressure (hydrostatic and capillary). For vanishing specific storage Ss, the active pore volume ϑ(ψ) coincides with the water content θ(ψ). Now we are ready to formulate the initial boundary value problem. The porous medium occupies a domain Ω ⊂ R2. If x ∈ Ω then we write either x = (x1, x2) or x = (x1, z). Furthermore, the symbol Γ denotes the boundary of Ω, T > 0 is the final physical time of our interest and we put I := (0, T ). Moreover, symbols 60 vol. 61 Special Issue/2021 Numerical simulation of a single ring infiltration experiment ΓD, ΓN and ΓE , respectively, denotes the correspond- ing parts of Γ, where the Dirichlet, Neumann and seepage face BCs are prescribed, respectively. Obvi- ously, these parts have to be mutually disjoint and ΓD ∪ ΓN ∪ ΓE = Γ. In virtue of (2.1) and (2.8) – (2.9), we define the following problem. Find a function Ψ = Ψ(x, t) : Ω× I → R such that ∂ϑ(ψ) ∂t −∇ · (K(ψ)∇Ψ) = 0 in Ω× I, (2.11a) Ψ = ΨD on ΓD × I, (2.11b) K(ψ)∇Ψ · n = qN on ΓN × I, (2.11c) ∇Ψ · n = 0 if Ψ ≤ 0 Ψ = 0 if ∇Ψ · n ≤ 0 } on ΓE × I, (2.11d) Ψ(x, t = 0) = Ψ0 in Ω. (2.11e) 3. Space-time discontinuous Galerkin discretization We solve (2.11) numerically using the space-time dis- continuous Galerkin (STDG) method. We refer to [24] where it was developed for the numerical solution of a time-dependent partial differential equation. Here, we present its adaptation for the Richards equation. 3.1. Space-time partition First, using a finite sequence 0 = t0 < t1 < . . . < tr = T we split the time interval I into sub-intervals intervals Im = (tm−1, tm], m = 1, . . . , r of lengths |Im| = τm. We define τ = maxm=1,...,rτm. For each interval tm, m = 1, . . . , r, we consider a triangulation Th,m of Ω consisting of a finite number of closed mutually disjoint triangles K covering Ω, i.e., Ω = ∪K∈Th,m K, m = 0, . . . , r. The meshes Th,m can differ for each time level. By ∂K and hK we denote the boundary and the diameter, respectively, of K ∈ Th,m. Furthermore, by Γh,m we denote the set of edges of Th,m. We define several subsets of Γh,m, namely ΓIh,m = {γ ∈ Γh,m, γ ⊂ Ω} (interior edges), (3.1) ΓDh,m = {γ ∈ Γh,m, γ ⊂ ΓD} (Dirichlet edges), ΓNh,m = {γ ∈ Γh,m, γ ⊂ ΓN} (Neumann edges), ΓEh,m = {γ ∈ Γh,m, γ ⊂ ΓE} (seepage face edges), ΓBh,m = ΓDh,m ∪ ΓNh,m ∪ ΓEh,m, ΓIDh,m = ΓIh,m ∪ ΓDh,m. Finally, a unit normal to γ ∈ Γh,m is denoted by nγ . We assume that vectors nγ , γ ∈ ΓBh,m are outer normals to ∂Ω. The orientation of nγ of interior edges γ ∈ ΓIh,m is arbitrary but fixed. 3.2. Finite element space We will seek the approximate solution in a space of piecewise polynomial but discontinuous functions. We use a fixed polynomial degree q ∈ N wrt. t ∈ I but a possibly varying polynomial degree pK wrt. x ∈ Ω for each mesh element K ∈ Th,m, m = 0, . . . , r. We define the spaces of space-time discontinuous piecewise polynomial functions on Ω× I as Sτ,qh,p := {ϕ : Ω× I → R; ϕ|K×Im ∈ P pK ,q(K × Im), K ∈ Th,m, m = 1, . . . , r}, (3.2) where P pK ,q(K × Im) is the space of space-time poly- nomial functions on K × Im having degree ≤ pK wrt. x ∈ K and degree ≤ q wrt. t ∈ Im for K ∈ Th,m and m = 1, . . . , r. Moreover, a restriction of Sτ,qh,p on the time layer Im, m = 1, . . . , r is given by Sτ,qm,h,p := {ϕ : Ω× Im → R; (3.3) ϕ|K×Im ∈ P pK ,q(K × Im), K ∈ Th,m}. Finally, we introduce the jumps and mean val- ues of functions from Sτ,qh,p. For γ ∈ ΓIh,m and v ∈ Sτ,qm,h,p, m = 0, . . . , r we define [v]γ = v|(+) γ − v|(−) γ , (3.4) 〈v〉γ = 1 2 ( v|(+) γ + v|(−) γ ) , (3.5) where v|(±) γ = limδ→0± v(x + δnγ), x ∈ γ. For γ ∈ ΓBh,m, we set [v]γ = 〈v〉γ = v|γ . In the case that nγ , [·]γ and 〈·〉γ are arguments of ∫ γ . . . dS, γ ∈ Γh,m we omit the subscript γ. Additionally, we set {ϕ}m := ϕ|+m − ϕ|−m, ϕ|±m := lim δ→0± ϕ(tm + δ), (3.6) defining a jump of ϕ ∈ Sτ,qh,p wrt. time at tm, m = 0, . . . , r. We observe that the finite element space Sτ,qh,p is independent on the type of boundary conditions since they are incorporated directly in the numerical scheme. This is different from conforming finite element meth- ods. 3.3. Space discretization We discretize problem (2.11) in space using the incom- plete interior penalty Galerkin (IIPG) method. The main advantage of the IIPG technique is the absence of additional stabilization terms, cf. [24, Chapter 2]. In order to simplify an exposition of IIPG method, we assume in this section that ΓE = ∅. The general case will be discussed in Section 3.4. Let Ψ be a sufficiently smooth exact solution of (2.11). Let m = 1, . . . , r be arbitrary but fixed. We multiply (2.11a) with ϕ ∈ Sτ,qh,p, integrate over K ∈ Th,m, sum over all elements K ∈ Th,m, and carry out some manipulation. The Green’s theorem applied to 61 V. Dolejší, M. Kuráž, P. Solin Acta Polytechnica the conductivity term and (2.11c) yields − ∑ K∈Th,m ∫ K ∇ · (K(Ψ− z)∇Ψ)ϕdx (3.7) = ∑ K∈Th,m (K(Ψ− z)∇Ψ,∇ϕ)K − ∑ γ∈ΓID h,m (〈K(Ψ− z)∇Ψ〉 · n, [ϕ])γ − ∑ γ∈ΓN h,m (qN , ϕ)γ . In order to ensure the coercivity of the resulting nu- merical scheme, we include additional terms to (3.7) which are vanishing if Ψ is smooth and satisfies (2.11b). Namely,∑ γ∈ΓI h,m (σ[Ψ], [ϕ])γ + ∑ γ∈ΓD h,m (σ(Ψ−ΨD), ϕ)γ , (3.8) where σ is the penalty parameter. The terms in (3.8) are called the interior and boundary penalties and they replace the interelement continuity of approximate solution from Sτ,qh,p. Taking into account (3.7)–(3.8), we obtain the rela- tion (∂tϑ(Ψ), ϕ)Ω + ah,m(Ψ, ϕ) = 0 ∀ϕ ∈ Sτ,qh,p, (3.9) where ah,m(Ψ, ϕ) := ∑ K∈Th,m (K(Ψ− z)∇Ψ,∇ϕ)K (3.10) − ∑ γ∈ΓID h,m (〈K(Ψ− z)∇Ψ〉 · n, [ϕ])γ + ∑ γ∈ΓI h,m (σ[ϕ], [Ψ])γ + ∑ γ∈ΓD h,m (σϕ,Ψ−ΨD)γ − (qN , ϕ)ΓN . Let us note that Ψ and ϕ are functions of x = (x1, x2) = (x1, z) ∈ Ω and t ∈ I. 3.4. Numerical realization of seepage face boundary (2.8) In this section, we extend the IIPG discretization of (2.1) from Section 3.3 to a general case, i.e., ΓE 6= ∅. Similarly as in [23], we implement condition (2.8) using the scheme if ψ < 0 and ∇Ψ · n > 0 then ∇Ψ · n = 0, if ψ < 0 and ∇Ψ · n < 0 then ∇Ψ · n = 0, if ψ > 0 and ∇Ψ · n > 0 then ∇Ψ · n = 0, if ψ > 0 and ∇Ψ · n < 0 then ψ = 0, (3.11) The quantities before “then” are taken from the inte- rior of Ω and the equalities after “then” correspond to the prescribed BCs. Let us note that practical implementation of (3.11) in the framework of conforming finite element methods is quite difficult since the finite element space, where the approximate solution is sought, depends on the Dirichlet part of boundary ΓD and then the realization of (3.11) requires frequent changes of finite element space. The realization of (3.11) using the DG method is much easier. On the seepage face part of boundary ΓE , we define a scalar function σ∗ : ΓE → {0, 1} whose role is to switch between (2.8a) (Neumann BC) for σ∗ = 0 and (2.8b) (Dirichlet BC) for σ∗ = 1. Then we supplement the form (3.10) by additional integrals over ΓE in the form ah,m(Ψ, ϕ) := ∑ K∈Th,m (K(Ψ− z)∇Ψ,∇ϕ)K (3.12) − ∑ γ∈ΓID h,m (〈K(Ψ− z)∇Ψ〉 · n, [ϕ])γ − ∑ γ∈ΓE h,m (σ∗K(Ψ− z)∇Ψ · n, ϕ)γ + ∑ γ∈ΓI h,m (σ[ϕ], [Ψ])γ + ∑ γ∈ΓD h,m (σϕ,Ψ−ΨD)γ + ∑ γ∈ΓE h,m (σ∗σϕ,Ψ− z)γ − (qN , ϕ)ΓN . Form ah,m from (3.12) behaves in the following way: Let gE ⊂ ΓE be an arbitrary but fixed part of the seepage face boundary. If σ∗|gE = 1 then Dirichlet BC Ψ = z is enforced on gE . On the other hand, if σ∗|gE = 0 then the homogeneous Neumann BC is used on gE . Now, let us discuss the choice of the function σ∗ in (3.12). In virtue of (3.11), we set σ∗ according to the signs of ψ and −K(Ψ− z)∇Ψ · n on ΓE . Therefore, we define σ∗(Ψ(x, t)) = { 1 if ψ(x, t) ≥ 0 and ∇Ψ · n ≤ 0, 0 otherwise . (3.13) However, this choice is not suitable from the practical point of view because the function σ∗ only admits two distinct values 0 and 1. As a result, iterative solvers for the resulting algebraic system will fail to converge. This unpleasant effect can be suppressed via smoothing of σ∗ near the transition from σ∗ = 0 and σ∗ = 1. Hence, we introduce transition regions of small widths sph > 0 and sfl > 0, and define a function σ∗h(Ψ(x, t)) which is continuously differentiable with respect to Ψ. It is equal to σ∗ if ψ(x, t) 6∈ [0, sph] and K(Ψ− z)∇Ψ · n 6∈ [0, sfl]. In order to ensure the consistency of the resulting numerical scheme, the transition regions must van- ish for a vanishing mesh size. Therefore we define 62 vol. 61 Special Issue/2021 Numerical simulation of a single ring infiltration experiment sph|K := c1 hK and sfl|K := c2 hK , where c1 > 0 and c2 > 0 are empirically chosen constants. Typically we use the values c1 = c2 = 0.1. Numerical experiments show that the numerical scheme only depends weakly on sph|K and sfl|K . Finally, we replace σ∗ in (3.12) with σ∗h and obtain ah,m(Ψ, ϕ) := ∑ K∈Th,m (K(Ψ− z)∇Ψ,∇ϕ)K (3.14) − ∑ γ∈ΓID h,m (〈K(Ψ− z)∇Ψ〉 · n, [ϕ])γ − ∑ γ∈ΓE h,m (σ∗h(Ψ)K(Ψ− z)∇Ψ · n, ϕ)γ + ∑ γ∈ΓI h,m (σ[ϕ], [Ψ])γ + ∑ γ∈ΓD h,m (σϕ,Ψ−ΨD)γ + ∑ γ∈ΓE h,m (σ∗h(Ψ)σϕ,Ψ− z)γ − (qN , ϕ)ΓN . 3.5. Time discretization The time discretization is carried out using the space- time discontinuous Galerkin method (STDGM) (see [24]), which reads Ah,m(Ψ, ϕ) = 0 ∀ϕ ∈ Sτ,qh,p, m = 1, . . . , r, (3.15) where Ah,m(Ψ, ϕ) = ∫ Im ( (∂tϑ(Ψ− z), ϕ)Ω + ah,m(Ψ, ϕ) ) dt + ({ϑ(Ψ− z)}m−1, ϕ| + m−1)Ω. (3.16) Now, we define the final hp-space-time discontinu- ous Galerkin (hp-STDG) solution of problem (2.11): Definition 3.1. Let Ah,m be the form (3.16). The function Ψhτ ∈ Sτ,qh,p is called the approximate solution of (2.11) if Ah,m(Ψhτ , ϕ) = 0 ∀ϕ ∈ Sτ,qh,p, m = 1, . . . , r (3.17) where Ψhτ |−0 := Ψ0, cf. (2.11e). The term ({ϑ(Ψ− z)}m−1, ϕ| + m−1)Ω in (3.16) re- places in some sense the interelement continuity in time. The approximate solutions on two subsequent time layers are connected together only by this term. A major advantage of this treatment is that differ- ent grids (and polynomial degrees pK , K ∈ Th,m) can be used in each time interval separately. This is advantageous from the point of view of adaptive methods. 3.6. Solution of the discrete problem The approximate solution (3.17) represents a sequence of r nonlinear algebraic systems of equations, one system for each m = 1, . . . , r. The size of the system (= number of degrees of freedom) is equal to the dimension of the space Sτ,qm,h,p and it changes for m = 1, . . . , r. The system is solved using the Anderson acceleration method [25]. The stopping criteria and the choice of the time step are controlled in such a way that the errors arising due to the space discretization, time discretization, and iterative solution of nonlinear algebraic systems are balanced. The meshes and polynomial degrees are chosen adaptively using a technique introduced in [26, 27]. Using a higher order reconstruction of the computed approximate solution, this technique proposes new sizes and shapes of mesh elements, as well as the corresponding polynomial degrees. The basic idea of the mesh adaptation algorithm is to minimize the number of degrees of freedom on the prescribed toler- ance level for the interpolation error measured in the H1(Ω)-seminorm. 4. SR experiment setup The proposed SR experiment is depicted in Fig. 2. The ring is located in the vicinity of a hillside, which is a typical application for the seepage face boundary discussed in this paper. The porous environment is anisotropic, with layers inclined into the hillside. Such experimental setup was chosen to combine the effects of a moving wetting front, seepage face boundary, and discontinuous flow domain. The domain is illustrated in the left part of Fig. 2. The initial and boundary conditions are given as: Γ1 (=: ΓN ) no flow boundary ∂Ψ ∂~n = 0, ∀(x, t) ∈ Γ1 × t ∈ [0, T ) Γ2 (=: ΓD) ponding depth 5 cm Ψ(x, t) = 0.05 + z, ∀(x, t) ∈ Γ2 × t ∈ [0, T ) Γ3 (=: ΓE) seepage face, cf. (2.7) and Section 3.4. This is shown in the right part of Fig. 2. The origin of the coordinate system is located at the left-bottom corner of the domain. The initial condition is a hydrostatic state with zero capillary pressure at the impermeable bottom, Ψ(x, t0) = −2 m, ∀x ∈ Ω. (4.1) The soil hydraulic properties were assumed as fol- lows: α = 0.8 m−1, n = 1.2, m = 0.167, θs = 0.55, θr = 0, Ss = 10−3 m−1. Soil properties were assumed here to be anisotropic with rotated main axes of anisotropy as depicted on Fig. 3. The rotation angle αxx′ was -30◦, and so the layers were inclined into the hillslope. Kx′ was 63 V. Dolejší, M. Kuráž, P. Solin Acta Polytechnica Figure 2. Scheme of the single ring infiltration experiment. Figure 3. Example of rotated main axes of anisotropy of hydraulic conductivity, obtained from [28]. assumed to be ten times greater than Kz′ . The re- sulting tensor of anisotropy of saturated hydraulic conductivity for our global system of coordinates was Ks = ( 1.031× 10−4 −5.187× 10−5 −5.187× 10−5 4.323× 10−5 ) m.s−1. Please note that the example presented here is purely synthetic. 4.1. Numerical results Fig. 4 presents the approximate solution (hydraulic head) for different time levels as well as the dynam- ically adapted meshes. Fig. 5 shows the dynamic development of the mesh and the number of degrees of freedom. It should be noted that our case study is sensitive to discretization qualities as explained in [29]. The numerical solution exhibits physically justifiable characteristics while maitaining a low number of de- grees of freedom. The maximum number of degrees of freedom was roughly 1000. Taking into account the domain size and dimension, such automatically adapted mesh is highly optimal. Finally we also present a quantitative comparison of an adaptive computation and a computation on a fixed fine mesh. Namely, we compare the actual water flux through the input boundary F (t) = − ∫ Γ2 K(ψ)∇Ψ(x, t) · n dx (4.2) and the total (accumulated flux) F (t) = − ∫ t 0 ∫ Γ2 K(ψ)∇Ψ(x, t′) · n dxdt′. (4.3) The fine mesh was chosen in such a way that both computations have roughly the same level of accuracy. Fig. 6 shows the values of the actual and total fluxes in dependence on time. We observe that the actual flux is the highest at the beginning of the computation and starting from t ≈ 0.2 hours, it is approximately constant. Finally, when the seepage face is active, the actual flux decreases (the total flux stagnates). Moreover, we find that both computations (adaptive and fixed meshes) give very similar behaviours. The main difference is at the end of the computation where the actual flux is already negligible. The difference is caused by a large final time step chosen by the adap- tive method. Furthermore, the adaptive computation uses approximately 200 elements and 10,000 degrees of freedom. It takes 74,646 seconds. The fixed mesh computation uses 2279 elements and 45,580 degrees of freedom, and it takes 401,769 seconds. Hence, clearly, the adaptive computation is much more efficient. References [1] X. Xu, C. Lewis, W. Liu, et al. Analysis of single-ring infiltrometer data for soil hydraulic properties estimation: Comparison of best and wu methods. Agricultural Water Management 107:34 – 41, 2012. doi:10.1016/j.agwat.2012.01.004. [2] M. Nakhaei, J. Šim�unek. Parameter estimation of soil hydraulic and thermal property functions for unsaturated porous media using the hydrus-2d code. Journal of Hydrology and Hydromechanics 62(1):7–15, 2014. doi:10.2478/johh-2014-0008. [3] L. A. Richards. Capillary conduction of liquids through porous mediums. Journal of Applied Physics 1(5):318–333, 1931. doi:10.1063/1.1745010. [4] S. C. Iden, J. R. Blöcher, E. Diamantopoulos, et al. Numerical test of the laboratory evaporation method using coupled water, vapor and heat flow modelling. Journal of Hydrology 570:574 – 583, 2019. doi:10.1016/j.jhydrol.2018.12.045. 64 http://dx.doi.org/10.1016/j.agwat.2012.01.004 http://dx.doi.org/10.2478/johh-2014-0008 http://dx.doi.org/10.1063/1.1745010 http://dx.doi.org/10.1016/j.jhydrol.2018.12.045 vol. 61 Special Issue/2021 Numerical simulation of a single ring infiltration experiment Figure 4. Hydraulic head and a detail of the dynamically adapted mesh. 1e+02 1e+03 1e+04 1e+05 0 0.5 1 1.5 2 2.5 3 3.5 4 σ m time [hours] number of mesh elements number of DoF Figure 5. Adaptively updated mesh: number of mesh elements and number of degrees of freedom w.r.t. physical time. 65 V. Dolejší, M. Kuráž, P. Solin Acta Polytechnica 1.0e-03 1.0e-02 1.0e-01 1.0e+00 0 0.5 1 1.5 2 2.5 3 3.5 4 ac tu al flu x F (t ) = − ∫ Γ 2 K ∇ Ψ (x ,t )· nd x time [hours] flux, no adaptation flux, hp adaptation 1.0e-03 1.0e-02 1.0e-01 1.0e+00 0 0.5 1 1.5 2 2.5 3 3.5 4 -1.0e-03 -5.0e-04 0.0e+00 5.0e-04 1.0e-03 1.5e-03 2.0e-03 to ta lfl ux F( t) = − ∫ t 0 ∫ Γ 2 K ∇ Ψ (x ,t ′ )· nd xd t′ di ff er en ce hp ad ap t- no ad ap t[ m · da y− 1 ] time [hours] total flux, no adaptation total flux, hp adaptation difference: hp adapt - no adapt Figure 6. Comparison of the actual flux F (t) (up) and the total (accumulated) flux (down) for the two computations, cf. (4.2)–(4.3). [5] A. Binley, K. Beven. Vadose zone flow model uncertainty as conditioned on geophysical data. Ground Water 41(2):119–127, 2003. doi:10.1111/j.1745-6584.2003.tb02576.x. [6] S. Würzer, N. Wever, R. Juras, et al. Modelling liquid water transport in snow under rain-on-snow conditions – considering preferential flow. Hydrology and Earth System Sciences 21(3):1741–1756, 2017. doi:10.5194/hess-21-1741-2017. [7] M. Kuraz, P. Mayer, V. Havlicek, et al. Dual permeability variably saturated flow and contaminant transport modeling of a nuclear waste repository with capillary barrier protection. Applied Mathematics and Computation 219(13):7127 – 7138, 2013. ESCO 2010 Conference in Pilsen, June 21- 25, 2010, doi:10.1016/j.amc.2011.08.109. [8] H. Alt, S. Luckhaus. Quasilinear elliptic-parabolic differential equations. Mathematische Zeitschrift 183(3):311–341, 1983. doi:10.1007/BF01176474. [9] F. Otto. L1-contraction and uniqueness for quasilinear elliptic-parabolic equations. Journal of Differential Equations 131(1):20–38, 1996. [10] F. Otto. L1–contraction and uniqueness for unstationary saturated-unsaturated porous media flow. Adv Math Sci Appl 7(2):537–553, 1997. [11] L. Lam, D. Fredlund. Saturated-unsaturated transient finite element seepage model for geotechnical engineering. Advances in Water Resources 7(3):132 – 136, 1984. doi:10.1016/0309-1708(84)90042-3. [12] C. Kees, M. Farthing, C. Dawson. Locally conservative, stabilized finite element methods for variably saturated flow. Computer Methods in Applied Mechanics and Engineering 197(51):4610 – 4625, 2008. doi:10.1016/j.cma.2008.06.005. [13] J. Šembera, M. Beneš. Nonlinear Galerkin method for reaction-diffusion systems admitting invariant regions. Journal of Computational and Applied Mathematics 136(1):163 – 176, 2001. doi:10.1016/S0377-0427(00)00582-3. 66 http://dx.doi.org/10.1111/j.1745-6584.2003.tb02576.x http://dx.doi.org/10.5194/hess-21-1741-2017 http://dx.doi.org/10.1016/j.amc.2011.08.109 http://dx.doi.org/10.1007/BF01176474 http://dx.doi.org/10.1016/0309-1708(84)90042-3 http://dx.doi.org/10.1016/j.cma.2008.06.005 http://dx.doi.org/10.1016/S0377-0427(00)00582-3 vol. 61 Special Issue/2021 Numerical simulation of a single ring infiltration experiment [14] P. Solin, M. Kuraz. Solving the nonstationary Richards equation with adaptive hp-FEM. Advances in Water Resources 34(9):1062 – 1081, 2011. New Computational Methods and Software Tools, doi:10.1016/j.advwatres.2011.04.020. [15] M. Tocci, C. Kelley, C. Miller. Accurate and economical solution of the pressure-head form of Richards’ equation by the method of lines. Advances in Water Resources 20(1):1 – 14, 1997. doi:10.1016/S0309-1708(96)00008-5. [16] C. Miller, C. Abhishek, M. Farthing. A spatially and temporally adaptive solution of Richards’ equation. Advances in Water Resources 29(4):525 – 545, 2006. doi:10.1016/j.advwatres.2005.06.008. [17] M. Kuraz, P. Mayer, V. Havlicek, P. Pech. Domain decomposition adaptivity for the Richards equation model. Computing 95(1):501–519, 2013. doi:10.1007/s00607-012-0279-8. [18] M. Kuraz, P. Mayer, P. Pech. Solving the nonlinear Richards equation model with adaptive domain decomposition. Journal of Computational and Applied Mathematics 270:2 – 11, 2014. Fourth International Conference on Finite Element Methods in Engineering and Sciences (FEMTEC 2013), doi:10.1016/j.cam.2014.03.010. [19] M. Kuraz, P. Mayer, P. Pech. Solving the nonlinear and nonstationary Richards equation with two-level adaptive domain decomposition (dd-adaptivity). Applied Mathematics and Computation 267:207 – 222, 2015. The Fourth European Seminar on Computing (ESCO 2014), doi:10.1016/j.amc.2015.03.130. [20] V. Dolejší, M. Kuráž, P. Solin. Adaptive higher-order space-time discontinuous galerkin method for the computer simulation of variably-saturated porous media flows. Applied Mathematical Modelling 72:276 – 305, 2019. doi:10.1016/j.apm.2019.02.037. [21] M. T. van Genuchten. Closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal 44(5):892–898, 1980. doi:10.2136/sssaj1980.03615995004400050002x. [22] Y. Mualem. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resources Research 12(3):513–522, 1976. https://agupubs.onlinelibrary.wiley.com/doi/ pdf/10.1029/WR012i003p00513 doi:10.1029/WR012i003p00513. [23] M. Kuraz, P. Mayer, M. Leps, D. Trpkosova. An adaptive time discretization of the classical and the dual porosity model of Richards’ equation. Journal of Computational and Applied Mathematics 233(12):3167 – 3177, 2010. Finite Element Methods in Engineering and Science (FEMTEC 2009), doi:10.1016/j.cam.2009.11.056. [24] V. Dolejší, M. Feistauer. Discontinuous Galerkin Method – Analysis and Applications to Compressible Flow. Springer Series in Computational Mathematics 48. Springer, Cham, 2015. [25] H. F. Walker, P. Ni. Anderson acceleration for fixed-point iterations. SIAM J Numer Anal 49(4):1715–1735, 2011. doi:10.1137/10078356X. [26] V. Dolejší. Anisotropic hp-adaptive method based on interpolation error estimates in the Lq-norm. Appl Numer Math 82:80–114, 2014. doi:10.1016/j.apnum.2014.03.003. [27] V. Dolejší. Anisotropic hp-adaptive method based on interpolation error estimates in the H1-seminorm. Appl Math 60(6):597–616, 2015. doi:10.1007/s10492-015-0113-7. [28] M. Kuraz, J. R. Bloecher. Hydrodynamic of porous media. CULS in Prague, 2017. Http://drutes.org/documents/notes.pdf. [29] J. Dusek, M. Dohnal, T. Vogel. Numerical analysis of ponded infiltration experiment under different experimental conditions. Soil and Water Research 4(SPECIAL ISSUE 2):22–27, 2009. doi:10.17221/1368-SWR. 67 http://dx.doi.org/10.1016/j.advwatres.2011.04.020 http://dx.doi.org/10.1016/S0309-1708(96)00008-5 http://dx.doi.org/10.1016/j.advwatres.2005.06.008 http://dx.doi.org/10.1007/s00607-012-0279-8 http://dx.doi.org/10.1016/j.cam.2014.03.010 http://dx.doi.org/10.1016/j.amc.2015.03.130 http://dx.doi.org/10.1016/j.apm.2019.02.037 http://dx.doi.org/10.2136/sssaj1980.03615995004400050002x https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/WR012i003p00513 https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/WR012i003p00513 http://dx.doi.org/10.1029/WR012i003p00513 http://dx.doi.org/10.1016/j.cam.2009.11.056 http://dx.doi.org/10.1137/10078356X http://dx.doi.org/10.1016/j.apnum.2014.03.003 http://dx.doi.org/10.1007/s10492-015-0113-7 http://dx.doi.org/10.17221/1368-SWR Acta Polytechnica 61(SI):59–67, 2021 1 Introduction 2 Governing equations 2.1 Initial boundary value problem 3 Space-time discontinuous Galerkin discretization 3.1 Space-time partition 3.2 Finite element space 3.3 Space discretization 3.4 Numerical realization of seepage face boundary (2.8) 3.5 Time discretization 3.6 Solution of the discrete problem 4 SR experiment setup 4.1 Numerical results References