Second-order accurate monotone finite volume scheme for Richards’ equation

Journal of Computational Physics 239 (2013) 123–137

Contents lists available at SciVerse ScienceDirect

Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp

Second-order accurate monotone finite volume scheme for Richards’ equation Oleksandr Misiats a,⇑, Konstantin Lipnikov b a b

Department of Mathematics, The Pennsylvania State University, United States MS B284, Los Alamos National Laboratory, Los Alamos, NM 87544, United States

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 9 February 2012 Received in revised form 27 August 2012 Accepted 3 September 2012 Available online 16 October 2012

In this work we perform a theoretical and numerical analysis of Richards’ equation. For certain types of nonlinearities we provide explicit analytical solutions. These solutions are used to show that conventional unconditionally monotone finite volume schemes have only first-order accuracy. We derive necessary and sufficient conditions for the monotonicity of finite volume discretizations and use these conditions to construct a monotone finite volume discretization accurate to second-order. Ó 2012 Elsevier Inc. All rights reserved.

Keywords: Richards’ equation Finite volume Monotonicity Accuracy

1. Introduction Richards’ equation plays a crucial role in the analysis of saturated–unsaturated fluid flows through a porous medium. The underlying physical model makes the fundamental assumption that the movement of the gas phase in a two-phase flow can be neglected. The focus of this article is on the numerical and theoretical analysis of steady-state solutions to Richards’ equation in heterogeneous media. Therefore, without loss of generality, we consider a slightly modified equation:

   @u @ @u ¼ KðxÞaðuÞ þg ; @t @x @x

ð1Þ

where u ¼ uðt; xÞ; ðt; xÞ 2 ½0; 1Þ  ½0; 2a for some a > 0, KðxÞ is a piecewise constant positive function,

KðxÞ ¼ f K 1 ;

x 2 ½0; a; K 2 ;

x 2 ða; 2a;

and g > 0 is a constant. Some of the results obtained will be extended to a two-dimensional Richards’ equation. We assume that the nonlinearity aðuÞ is an absolutely continuous Lipschitz function on R satisfying the ellipticity condition

1 6 aðuÞ 6 1; 1 þ jujr

u 2 R;

r > 0:

ð2Þ

We shall refer to u and KðxÞ as the potential and the diffusion coefficients, respectively. Originally suggested by Lorenzo Richards in 1931 [20], this model has been extensively studied by a variety of authors both in one and multiple dimensions. Several models for the nonlinear function aðuÞ have been suggested by Corey–Brooks [7] as well as Mualem [27] and van Genuchten [23]. The quasilinear evolution Eq. (1) presents a challenge to the analysis, ⇑ Corresponding author. E-mail address: [email protected] (O. Misiats). 0021-9991/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2012.09.004

124

O. Misiats, K. Lipnikov / Journal of Computational Physics 239 (2013) 123–137

since it has a closed-form analytical solution only in a limited number of cases. Some particular cases of (1) which admit analytical traveling front solutions have been studied in [11]. The idea can be illustrated with the one-dimensional Zeldovich equation

  @u @ @u ; ¼ us @t @x @x

s > 0; x P 0; 1

subject to the homogeneous initial condition uð0; xÞ ¼ 0 and the boundary conditions uðt; 0Þ ¼ ct s and uðt; þ1Þ ¼ 0. This equation has a unique solution in the form of a traveling wave

(

1

u ¼ ðv sðx  v tÞÞ s ; u ¼ 0;

x > v t;

v¼

qffiffiffi cs ; s

x 6 v t:

Various analytical and quasi-analytical solutions to the Richards’ equations can be found in the literature. A large class of solutions was obtained using the exponential function, aðuÞ, which allows for the linearization of the flow equation (see [21,18,26] and references therein). Many solutions were limited to homogeneous hydraulic conductivity (see [3,6] and references therein). A steady-state solution to (1) with a piecewise-constant KðxÞ can always be expressed via quadratures [13,25]; however, for the typical form of the hydraulic conductivity, the quadrature integrals cannot be expressed via elementary functions. One of the goals of this work is to find closed-form analytical steady-state solutions for a two-layered soil. These analytical solutions will be used to analyze the convergence of the finite volume (FV) schemes. Due to the essential nonlinear nature of Eq. (1), finding a monotone discretization of second-order accuracy is a non-trivial issue. In the last few decades, several finite difference, finite element, and finite volume methods have been developed (see [9] and references therein). The existence and uniqueness of FV solutions to Richards’ equation and their convergence to the continuum solutions have been established in [19]. The FV method proved to be the most successful in capturing the steep gradients that typically arise in the solutions of (1). Nevertheless, on coarse grids that are often used in sensitivity and uncertainty quantification analysis, the standard FV schemes may produce non-monotone solutions, as illustrated in Figs. 7–9. Such overshoots and undershoots have been observed by several authors, e.g. [2,14,12]. The latter works focused on finding a weighted estimate of the hydraulic conductivity between adjacent control volumes in order to get the most accurate discretizations. However, to the best of our knowledge, rigorous analysis of the monotonicity on under-resolved meshes has not been performed. It is natural to expect that the steeper the solution gradients are, the finer the mesh required to resolve them must be. As shown in Figs. 7–9, any discrete solution becomes monotone on sufficiently fine meshes. However, in practice, the mesh may be given, while the solution gradients are a priori unknown and may be arbitrarily large. Therefore, it is important to have a discretization that is monotone for an arbitrarily coarse mesh. It is known that a monotone discretization may be obtained via an upwind scheme [24,5,8]. The simplest donor-type upwind scheme is inexpensive but only first-order. On the other hand, there exists upwind schemes accurate to second-order (e.g. see [22] for conservation laws, [4,15,17] for scalar advection equations, and [10] for singularly perturbed problems), which are monotone but expensive, especially in multi-dimensions. The present work focuses on studying a monotone FV scheme that is second-order accurate and has complexity comparable to that of the donor-type upwind scheme. After proving a local monotonicity result, we use it to derive necessary and sufficient conditions for monotonicity (see formulas (36) and (37)). These conditions are given in terms of the nonnegative parameter 0 6 c 6 1, which measures deviation from an upwind scheme. We select different flux formulas, based on the range of values of c that provide a monotone scheme. Next, we generalize the necessary and sufficient conditions for the multi-dimensional Richards’ equation

   @u ! ¼ div KðxÞaðuÞ ru þ ge1 ; @t

ð3Þ

! where X ¼ X1 [ X2  Rd ; d P 1, with X1 and X2 being connected domains with Lipschitz boundaries, e1 :¼ ð1; 0; . . . ; 0Þ is the direction of the gravity (which is essentially one-dimensional), and K is the piecewise constant function, KðxÞ ¼ K i when x 2 Xi . The outline of the paper is as follows. In Section 2, we show the well-posedness of the quasilinear Eq. (1). In Section 3, we derive an explicit solution for the steady state problem and use it to test the accuracy of various FV schemes in Section 4. In Section 5, we derive necessary and sufficient monotonicity conditions and propose a monotone scheme with second-order accuracy. Finally, in Section 6 we provide necessary and sufficient monotonicity conditions for a FV solution of Eq. (3). 2. Existence and uniqueness of a weak solution We start with the analysis of the solutions to the evolution Eq. (1) subject to the time-independent Dirichlet boundary conditions,

uðt; 0Þ ¼ 0;

uðt; 2aÞ ¼ ur ;

and the initial condition

t P 0;

ð4Þ

125

O. Misiats, K. Lipnikov / Journal of Computational Physics 239 (2013) 123–137

uð0; xÞ ¼ u0 ðxÞ 2 L2 ð0; 2aÞ:

ð5Þ 1;2

Clearly, if K 1 – K 2 , the boundary value problem (1), (4), (5) does not have a classical solution u 2 C ð½0; 1Þ  ½0; 2aÞ. Thus, the natural question is its well-posedness, i.e. the existence of a weak solution. In order to simplify the presentation, we perform the change of variables uðt; xÞ :¼ uðt; xÞ þ gx and rewrite the problem as follows:

8 @u   @u @ > < @t ¼ @x KðxÞ aðu  gxÞ @x ; uðt; 0Þ ¼ 0; uðt; 2aÞ ¼ ur þ 2ag :¼ ur ; > : uð0; xÞ ¼ u0 ðxÞ þ gx:

ð6Þ

r Let u0 ðxÞ ¼ u x and X ¼ ð0; 2aÞ. 2a

Definition 1. The problem (6) is said to have a weak solutionu 2 L2 ð½0; T; H1 ðXÞÞ if

Z

T

Z 0

0

for all

2a

ðuðx; tÞ  u0 ðxÞÞ

@v  @t

Z

T

Z

0

2a

KðxÞ aðu  gxÞ

0

@u @v dx ¼ 0 @x @x

v 2 L2 ð½0; T; H10 ðXÞÞ such that v t 2 L1 ð½0; T  XÞ and v ðx; TÞ ¼ 0.

Remark 1. In one-dimension any weak solution of (6) is an absolutely continuous function of x. Recall that absolute continuity of functions is a smoothness property which is stricter than continuity and uniform continuity. At this moment, we need the following theorem that follows immediately from the main result of [16]. Therefore, we omit its proof. Theorem 1. Let us assume that 9c > 0 such that

aðu  gxÞ P c;

8x; u 2 R:

ð7Þ

Then Eq. (6) has a unique weak solution on any time interval ½0; T; T > 0. Proposition 1. Let us assume that function aðuÞ satisfies the condition (2). Then problem (6) has a unique weak solution u 2 L2 ð½0; T; H1 ðXÞÞ on any time interval ½0; T; T > 0. Proof. We cannot apply Theorem 1 directly since the assumption (2) is weaker then the condition (7). To this end, we fix a constant M > 0 and define

8 aðv Þ; if jv j 6 M; > < aM ðv Þ :¼ aðMÞ; if v P M; > : aðMÞ; if v 6 M: Consider the problem

  8 @u @ uM @ M > ¼ KðxÞ ½ a M ðuM  gxÞ @x ; > @x @t > <

ð8Þ

uM ðt; 0Þ ¼ 0; uM ðt; 2aÞ ¼ ur ; > > > : uM ð0; xÞ ¼ u0 ðxÞ þ gx:

Theorem 1 guarantees the existence and uniqueness of the solution uM to (8). Moreover, uM satisfies (8) with a uniformly elliptic right-hand side. Thus, by the maximum principle,

sup ðt;xÞ2½0;T½0;2a

juM ðt; xÞj 6 supju0 ðxÞ þ gxj :¼ c0 ; x

with c0 independent of M. Choosing M :¼ c0 þ 2ajgj þ 1, we see that

aM ðuM ðt; xÞ  gxÞ  aðuM ðt; xÞ  gxÞ;

8ðt; xÞ 2 ½0; T  ½0; 2a:

Thus, uM ðt; xÞ  uðt; xÞ solves also Eq. (6). The uniqueness of uðt; xÞ is proved using the same idea. Let us assume that u1 ðt; xÞ and u2 ðt; xÞ are two distinct solutions to (6). Then u1 ðt; xÞ and u2 ðt; xÞ are two distinct solutions of (8) with

M :¼

sup ðt;xÞ2½0;T½0;2a

ju1 ðt; xÞj þ

which contradicts Theorem 1. h

sup

ju2 ðt; xÞj þ 1

ðt;xÞ2½0;T½0;2a

126

O. Misiats, K. Lipnikov / Journal of Computational Physics 239 (2013) 123–137

3. Analysis of the steady state problem We proceed to look for the steady state solutions of (1):

(

   KðxÞaðuÞ @u þ g ¼ 0; @x

@ @x

ð9Þ

uðt; 0Þ ¼ 0; uðt; 2aÞ ¼ ur : Integrating this equation, we obtain

  @u þ g  C; KðxÞaðuÞ @x

ð10Þ

where C is a constant flux, unknown at this moment and uniquely determined by the boundary condition ur . The flux is positive for evaporation and negative for infiltration. First, we note the critical case C ¼ 0 (later it is referred to as the first critical case) in which u ¼ gx and ur ¼ 2g a. It also follows from (10) that C > 0 iff ur > 2g a. The following example shows that the value of C is crucial in determining the behavior of the solutions to (9). It also illustrates the variety of solutions that can be used for verifying numerical codes. Let

aðuÞ ¼

1 : 1 þ u2

ð11Þ

This choice of nonlinearity is relatively simple but preserves the main features of the problem and shows the complexity of the possible solutions. In order to slightly reduce the number of cases to be considered, we assume that the parameters K 1 ; K 2 ; a and g satisfy the technical condition

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi K1

!

p2 1 < K2: 1þ 2 2 þ 4 4g a

ð12Þ

The analysis in the case when (12) does not hold may be performed analogously. Depending on the value of C, we will analyze all possible scenarios in full generality and illustrate them with the following choice of problem parameters:

K 1 ¼ 0:5;

K 2 ¼ 2;

g ¼ 2 and a ¼ 5:

ð13Þ

As mentioned above, there is a one-to-one correspondence between the flux C and the value ur at the right endpoint of the computational interval. In what follows, it is convenient to treat the value ur as a function of C. If C – 0, we can separate variables in (10):

dx ¼

C KðxÞ



du 1  gKðxÞ þ u2 C

:

ð14Þ

Let us denote u1 ðxÞ :¼ uðxÞ when x 2 ½0; a and u2 ðxÞ :¼ uðxÞ when x 2 ½a; 2a. To simplify formulas, we introduce two numbers

U1 ¼ 1 

g K1 C

and U2 ¼ 1 

g K2 : C

Case 1. Let C < 0. It follows from (14) and the boundary condition uð0Þ ¼ 0 that

u1 ðxÞ ¼ u2 ðxÞ ¼

pffiffiffiffiffiffi

U1 tan

pffiffiffiffiffiffi

U2 tan



 Cx pffiffiffiffiffiffi U1 ; K1

ð15Þ

Cðx  aÞ pffiffiffiffiffiffi U2 þ arctan K2

sffiffiffiffiffiffi

 pffiffiffiffiffiffi!! U1 Ca U1 tan : U2 K1

ð16Þ

The value ur is expressed from (16) as

ur ¼

pffiffiffiffiffiffi

C a pffiffiffiffiffiffi U2 tan U2 þ arctan K2

sffiffiffiffiffiffi



U1 C a pffiffiffiffiffiffi tan U1 U2 K1

!! :

ð17Þ

Formulas (15) and (16) impose natural restrictions on C : u1 and u2 have no vertical asymptotes on ½0; a and ½a; 2a, respectively, iff 0 > C > C 1 ¼ C 1 ðK 1 ; K 2 ; g; aÞ. The expression (17) provides a one-to-one correspondence between C 2 ðC 1 ; 0Þ and ur 2 ð1; 2g aÞ and thus C can be uniquely solved for given ur . For the sample choice of parameters (13), we find that C 1  0:014939. The solutions for C  0:00399, corresponding to ur ¼ 25, are illustrated in Fig. 2.

O. Misiats, K. Lipnikov / Journal of Computational Physics 239 (2013) 123–137

127

Case 2. Let 0 < C < gK 1 . Since uð0Þ ¼ 0, we obtain

 pffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffi Cx U1 ; u1 ðxÞ ¼  U1 tanh K1 pffiffiffiffiffiffiffiffiffiffi Cðx  aÞ pffiffiffiffiffiffiffiffiffiffi U2 þ arctanh u2 ðxÞ ¼  U2 tanh K2

sffiffiffiffiffiffi



U1 C a pffiffiffiffiffiffiffiffiffiffi U1 tanh U2 K1

!! :

ð18Þ

Similarly to Case 1,

pffiffiffiffiffiffiffiffiffiffi C a pffiffiffiffiffiffiffiffiffiffi ur ¼ ur ðCÞ ¼  U2 tanh U2 þ arctanh K2

sffiffiffiffiffiffi



U1 C a pffiffiffiffiffiffiffiffiffiffi tanh U1 U2 K1

!! ;

providing a one-to-one correspondence between C 2 ð0; gK 1 Þ and ur 2 ð2aK 1 ; ur ðgK 1 ÞÞ. For instance, the choice of parameters (13) gives ur ðgK 1 Þ  1:73. Case 3. Let C ¼ gK 1 which is the second critical case. We obtain

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi! K2 gK 1 K2 u2 ¼   1 tanh ðx  aÞ 1 ; K1 K2 K1

u1 ¼ 0; which gives

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi! K2 gK 1 a K 2 ur ¼   1 tanh 1 : K1 K2 K1 Parameters (13) give ur  1:73. Case 4. Let C > K 1 g and

pffiffiffiffiffiffi

U1 tan

 pffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffi Ca U1 < U2 : K1

ð19Þ

Then

u1 ðxÞ ¼

pffiffiffiffiffiffi

U1 tan



 Cx pffiffiffiffiffiffi U1 ; K1

ð20Þ

pffiffiffiffiffiffiffiffiffiffi Cðx  aÞ pffiffiffiffiffiffiffiffiffiffi U2  arctanh u2 ðxÞ ¼  U2 tanh K2

sffiffiffiffiffiffiffiffiffiffi U1

U2

tan

 pffiffiffiffiffiffi!! Ca K1

U1

;

and

pffiffiffiffiffiffiffiffiffiffi C a pffiffiffiffiffiffiffiffiffiffi U2  arctanh ur ¼  U2 tanh K2

sffiffiffiffiffiffiffiffiffiffi U1

 pffiffiffiffiffiffi!! Ca : U1 tan U2 K1

ð21Þ

Note that solution (20) does not have asymptotes on ½0; a if and only if

C<

K1 2

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

p2 2 gK 1 þg þ : a2 2

ð22Þ

The choice of coefficients (13) satisfies the inequality (12) and the condition (19) becomes C < 1:0216096. Now, formula (21) is a bijection between C 2 ð1; 1:021696Þ and ur 2 ð1:73; 1:71Þ. Case 5. Let C > K 1 g and

pffiffiffiffiffiffi

U1 tan

  pffiffiffiffiffiffiffiffiffiffi Cf a pffiffiffiffiffiffi U1 ¼ U2 ; K1

ð23Þ

which is the third critical case. Then u1 is given by (20) and

u2 ¼

pffiffiffiffiffiffiffiffiffiffi U2 :

For the choice of parameters (13), we obtain u2  ur  1:71.

ð24Þ

128

O. Misiats, K. Lipnikov / Journal of Computational Physics 239 (2013) 123–137

Remark 2. Assuming inequality (12), the function u given by (20) and (24) with C given by (23) is the unique solution of the boundary value problem with a homogeneous Neumann boundary condition at the right endpoint:

(

@ @x

   KðxÞaðuÞ @u þ g ¼ 0; @x

½1exuðt; 0Þ ¼ 0;

@uðt;2aÞ @x

¼ 0:

Case 6. Let C > K 1 g and

pffiffiffiffiffiffi

U1 tan

 pffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffi Ca U1 > U2 : K1

Then u1 is given by formula (20),

pffiffiffiffiffiffiffiffiffiffi Cðx  aÞ pffiffiffiffiffiffiffiffiffiffi U2  arccoth u2 ðxÞ ¼  U2 coth K2

sffiffiffiffiffiffiffiffiffiffi U1



C a pffiffiffiffiffiffi U1 tan U2 K1

!! ;

ð25Þ

and

pffiffiffiffiffiffiffiffiffiffi C a pffiffiffiffiffiffiffiffiffiffi ur ¼  U2 coth U2  arccoth K2

sffiffiffiffiffiffiffiffiffiffi U1

 pffiffiffiffiffiffi!! Ca tan U1 : U2 K1

In order for both solutions (20) and (25) not to have asymptotes on the intervals ð0; aÞ and ða; 2aÞ, respectively, in addition to (22), we require

C a pffiffiffiffiffiffiffiffiffiffi U2  arccoth K2

sffiffiffiffiffiffiffiffiffiffi U1

U2

 tan

C a pffiffiffiffiffiffi U1 K1

! < 0:

This implies that C < C 2 ðK 1 ; K 2 ; g; aÞ. For the choice of parameters (13), we obtain C 2 ¼ 1:0216104. Thus, the bijection is between C 2 ð1:0216096; 1:0216104Þ and ur 2 ð1:71; þ1Þ. The analysis of this example is complete. To summarize, we illustrate the variety of solutions u1 and u2 in Figs. 1–4.

Fig. 1. Structure of possible solutions for various values of ur and parameters (13).

O. Misiats, K. Lipnikov / Journal of Computational Physics 239 (2013) 123–137

129

Fig. 2. Left: case 1 with ur ¼ 25. Right: case 2 with ur ¼ 10.

Fig. 3. Left: case 3 with ur  1:73. Right: case 4 with ur ¼ 0.

Fig. 4. Left: case 5 with ur ¼ 1:71. Right: case 6 with ur ¼ 10.

4. Finite volume discretizations In this section, we use of the explicit steady state solutions from Section 3 to analyze various numerical discretizations for accuracy and monotonicity. Recall that the nonlinearity aðuÞ is given by (11). To this end, we will use the following two choices of parameters:

K 1 ¼ 0:5;

K 2 ¼ 2;

a ¼ 5; g ¼ 2; ur ¼ 0;

ð26Þ

K 1 ¼ 0:5;

K 2 ¼ 2;

a ¼ 5; g ¼ 10; ur ¼ 0:

ð27Þ

and

In both cases, the solutions are tangents on ½0; 5 and hyperbolic tangents on ½5; 10. The parameters (26) will be used to test the accuracy of various finite volume discretizations, while the parameters (27) will be used to determine the monotonicity of these schemes. The reason for working with two different reference solutions

130

O. Misiats, K. Lipnikov / Journal of Computational Physics 239 (2013) 123–137

is that the one induced by (27) has larger derivatives around the interface point x ¼ a and near the boundary x ¼ 2a and requires finer meshes for the convergence analysis. In turn, larger derivatives enable us to test the schemes for monotonicity rather accurately, since the lack of monotonicity makes a more profound impact on a numerical solution. We now proceed with discretizing Eq. (1) using a finite volume scheme. Denoting um n :¼ uðt m ; xn Þ, we get

umþ1  um n n

s

¼

1 h



 m   m  unþ1  um un  um n n1 þ g  jm þ g : 1 n2 2 h h

jmnþ1

ð28Þ

The numerical integration in time is performed up to T ¼ 500 using the time step s ¼ 0:0001. The initial condition is uð0; xÞ ¼ 0. The spatial discretization is done on a uniform mesh fxi ¼ hi; 0 6 i 6 Ng occupying the interval ½0; 10 with N ¼ 24; 36; 48; 96; 192 and 384. The values N ¼ 24; 36 and 48 correspond to coarse grids and are used to test the monotonicity of finite volume schemes. The larger values (N ¼ 48 and up) are used to determine their order of accuracy. Remark 3. The scheme (28) is an explicit finite volume scheme for the Richards’ equation. In terms of stability, this scheme is more restrictive than the implicit scheme, which is commonly used to solve the steady state Richards’ equation. However, for the purpose of this paper, we can use either of the methods, explicit or implicit. In order to simplify the theoretical 10 analysis of monotonicity, we use the explicit scheme, provided that we guarantee its stability. Since the mesh size is h ¼ N1 , the Courant number m for the finest mesh, N ¼ 384, satisfies the sufficient conditions for stability:

m :¼ max jKðxÞ aðuðxÞÞj x

s 2

h

 0:295 <

1 : 2

Thus, the explicit numerical scheme (28) is stable for all N 6 384. Also, an important observation is that the difference between the numerical solution of the stationary problem (9) and the numerical solution of the evolution Eq. (1) at T ¼ 500 is negligible. To verify this, the scheme (28) has been tested on time intervals T ¼ 1000 and T ¼ 2000. The difference between the obtained solutions was too small to be detected in the first six significant digits. The finite volume scheme requires the definition of transmissibility coefficients jm and jm , which are essentially the nþ12 n12 approximations of the nonlinear term KðxÞ aðuðt m ; xÞÞ at the points xnþ1 :¼ xn þ 2h and xn1 :¼ xn  2h, respectively. The choice of 2 2 this approximation is crucial, since it determines the properties of the scheme. We tested the following five choices of jm nþ12 summarized in Table 1. 4.1. Accuracy of FV schemes Tables 2 and 3 illustrate the l2 and l1 norms of errors for g ¼ 2 and g ¼ 10, respectively:

kykl1 :¼ maxjyi j and kykl2 :¼ 16i6n

n 1X y2 n i¼1 i

!1=2 ;

where y 2 Rn is an error vector. The larger the value of g is, the steeper the gradients of the solutions are; thus, the finer mesh is required to solve the problem accurately. Note that the schemes (HU) and (A) are only accurate to first-order, while the schemes (H), (HA) and (HMP) provide second-order accuracy. The scheme (HMP) is the most accurate in both l2 and l1 norms. The experimental order of convergence is calculated using the linear regression algorithm. As expected it is close to the theoretical one (in the l2 -norm) in Table 2. In Table 3, it deviates from the theoretical order of convergence due to the fact that the meshes N ¼ 192 and N ¼ 384 are not fine enough to capture the asymptotic (as N ! 1) convergence rates for g ¼ 10. However, these meshes rather accurately describe the asymptotic convergence rates for a smaller value of g : g ¼ 2, see Table 2. It is worth mentioning that the solution is non-smooth around the finitely many interface points,

Table 1 Five choices of transmissibility coefficients. Name

Symbol

Harmonic average

H

Formula Kðx ÞKðx

n jm ¼ 2 Kðxnþ1 Þaðunþ1 m nþ1 2

Arithmetic average

A

jm ¼ nþ1 2

Harmonic in K, arithmetic in a

HA

Harmonic in K, upwind in a

HU

Harmonic in K, midpoint in a

HMP

Þaðum Þ aðum n Þ nþ1 ÞþKðxn Þ aðum n Þ

nþ1

Kðxnþ1 Þ aðum ÞþKðxn Þ aðum n Þ nþ1 2

Kðxn ÞKðxnþ1 Þ m jm ¼ Kðx ðaðum nþ1 Þ þ aðun ÞÞ nþ1 n ÞþKðxnþ1 Þ 2

2Kðxn ÞKðxnþ1 Þ jm ¼ Kðx aðum n Þ nþ1 n ÞþKðxnþ1 Þ 2

um þum 

2Kðxn ÞKðxnþ1 Þ jm ¼ Kðx a nþ1 n ÞþKðxnþ1 Þ 2

n

nþ1

2

131

O. Misiats, K. Lipnikov / Journal of Computational Physics 239 (2013) 123–137 Table 2 l2 and l1 norms of error for the FV schemes (28) with g ¼ 2 and their order of accuracy. Scheme

Norm

H

l2

1.2 10

l1

8.11 102

3.15 102

2

2

A

HA

HU

HMP

N ¼ 48

N ¼ 96 2

N ¼ 192

3.37 10

3

N ¼ 384

Rate 4

1.93

1.03 102

3.03 103

1.58

3

3

1.02

8.72 10

4

2.19 10

l2

3.87 10

1.95 10

l1

2.39 101

1.47 101

8.28 102

4.39 102

0.82

3

3

4

4

1.84

9.57 10

4.69 10

l2

5.97 10

1.73 10

l1

3.59 102

1.15 102

3.22 103

9.01 104

1.78

2

2

2

3

0.92

4.85 10

1.29 10

l2

4.81 10

2.60 10

l1

1.49 101

8.76 102

4.87 102

2.60 102

0.84

l2

3.71 103

1.08 103

2.77 104

6.93 105

1.92

l1

2.06 102

8.21 103

2.63 103

7.60 104

1.59

1.38 10

7.16 10

Table 3 l2 and l1 norms of error for the FV schemes (28) with g ¼ 10.

H

A

HA

HU

HMP

Norm

N ¼ 48

N ¼ 96

N ¼ 192

N ¼ 384

Rate

l2

5.87 102

2.21 102

7.68 103

2.25 103

1.56

l1

2.81 101

1.92 101

1.02 101

4.28 102

0.91

2

2

2

2

0.83

l2

6.12 10

3.73 10

l1

4.02 101

3.58 101

2.67 101

1.73 101

0.41

2

2

3

3

2.12 10

1.09 10

l2

3.22 10

1.19 10

1.14 10

1.61

l1

1.78 101

1.08 101

4.80 102

1.65 102

1.15

2

2

2

2

3.81 10

l2

8.90 10

4.93 10

1.45 10

0.87

l1

4.44 101

2.86 101

1.74 101

1.05 101

0.70

2

3

3

4

7.31 10

1.52

1.10 102

0.99

l2

1.78 10

6.37 10

l1

9.42 102

4.06 102

2.68 10

2.40 10

2.52 102

and thus, the finer mesh is needed to see the asymptotically optimal (theoretical) convergence rate in l1 norm than in l2 norm. 4.2. Monotonicity of FV schemes It worth noting that the steady-state solution is always monotone on ½0; a and ½a; 2a for any choice of nonlinearity aðuÞ. In each subinterval we have @u ¼ Fðu; C; K i Þ where C is an integration constant, which implies that @u must have the constant @x @x sign on ½0; a and ½a; 2a. Indeed, if we assume that F changes sign, e.g. on ð0; aÞ, there exists a nonconstant solution u1 ðxÞ and a point x1 2 ð0; aÞ such that Fðu1 ðx1 Þ; C; K 1 Þ ¼ 0. This fact leads to a contradiction with Picard’s theorem, which guarantees the uniqueness of solutions, since there is another, trivial constant solution u2 ðxÞ  u1 ðx1 Þ passing through the point ðx1 ; u1 ðx1 ÞÞ. The formal proof of this fact is presented in the next section. Thus, the time-dependent solution of (1) is also monotone in x for sufficiently large t. The question of monotonicity (in x) for arbitrary t P 0 (provided that the initial data is a monotone function of x) has been studied in [1]. The simplest way to measure loss of monotonicity in each FV scheme is to identify the largest undershoot and overshoot. In our case, the undershoot is defined as the negative value of the discrete solution, while the overshoot is the difference between the maximal value of the discrete solution and the value at x ¼ 5. For a monotone discretization, both the overshoot and the undershoot must be zero. The results are presented in Table 4. The loss of monotonicity is also illustrated in Figs. 5–7. Additionally, we see that the boundary layer (the steep gradients) is resolved relatively accurately by the schemes (A), (HA) and (HMP) on a mesh with 36 grid points. Resolution with the scheme (H) requires 48 grid points, and for the scheme (HU) the boundary layer is resolved even on a coarse mesh with 24 grid points. The interior layer is resolved relatively accurately on a 36-point grid for schemes (A) and (HA), on a 48-point grid for the scheme (H), and on a 24-point grid for the schemes (HU) and (HMP). As expected, the only monotone scheme is the upwind (HU). However, it has only the first-order of accuracy (OðhÞ error), 2 while the schemes (H), (HA), and (HMP) provide second-order accuracy (Oðh Þ error). In the subsequent sections we propose an inexpensive FV scheme that remains monotone when solution gradients are under-resolved and has asymptotic secondorder accuracy.

132

O. Misiats, K. Lipnikov / Journal of Computational Physics 239 (2013) 123–137

Table 4 Overshoots and undershoots for the FV schemes (28) with g ¼ 2. N ¼ 24 H

A

HA

HU HMP

Overshoot Undershoot

N ¼ 36

N ¼ 48

1:51  10

1:70  10

2

0

2.76 101

1.09 101

0

1

3

2

Overshoot

7:29  10

1:59  10

0

Undershoot

2.39 101

1.47 101

0

1

3

0

Overshoot

7.15 10

Undershoot

5.86 102

1.58 10 0

Overshoot Undershoot

0 0

0 0

0 0

Overshoot

3.50 101 0

7.84 103 0

0

Undershoot

0

0

Fig. 5. Finite volume schemes H (left) and A (right).

Fig. 6. Finite volume schemes HA (left) and HU (right).

5. Necessary and sufficient conditions for monotonicity Since the mesh required for resolving the boundary and interior layers is a priori unknown, a mathematical criterion is developed in the next two subsections that allows us to formulate an adaptive monotone FV scheme. 5.1. One-dimensional problem In this section, we derive necessary and sufficient conditions for obtaining a monotonic discrete solution of the stationary problem (1) on the interval ½0; a. The other interval is analyzed analogously. In what follows, for definiteness we assume that the nonlinearity aðuÞ 2 C 1 ðRÞ and is a non-increasing function of u, i.e.

O. Misiats, K. Lipnikov / Journal of Computational Physics 239 (2013) 123–137

133

Fig. 7. Finite volume scheme HMP.

da 6 0: du A similar but technically more involved condition may be obtained in the general case. It is worth mentioning that in the majority of practically relevant cases the function aðuÞ is non-increasing, e.g.

( aðuÞ ¼

1 1þu2

1;

;

u > 0;

u 6 0:

The steady-state solution of the finite volume scheme (28) satisfies

jnþ12 ðunþ1  un þ ghÞ  jn12 ðun  un1 þ ghÞ ¼ 0: Without loss of generality, we look for

ð29Þ

jnþ12 and jn12 on ½0; a as follows:

jnþ12 ¼ K 1 ðcaðun Þ þ ð1  cÞaðunþ1 ÞÞ; jn12 ¼ K 1 ðcaðun1 Þ þ ð1  cÞaðun ÞÞ;

ð30Þ

where c ¼ cn 2 ½0; 1. The value c ¼ 1 gives the FV scheme (HU). The value c ¼ 0:5 gives the FV scheme (HA). Definition 2. The discretization (29) is called monotone on interval ½0; a if any discrete solution that satisfies (29) has the following property:

minfun1 ; unþ1 g 6 un 6 maxfun1 ; unþ1 g

ð31Þ

for all n P 1 such that xnþ1 2 ½0; a. Proposition 2. The discretization (29) with

jnþ12 and jn12 given by (30) is monotone for any mesh size h if and only if c ¼ 1.

Proof. Without loss of generality, we assume that K 1 ¼ 1. The condition (31) is clearly equivalent to two conditions:

if un1 6 un ; then unþ1 P un ;

ð32Þ

if un1 P un ; then unþ1 6 un :

ð33Þ

and

Assume un1 6 un . Then aðun1 Þ P aðun Þ and thus

jn12 ðun  un1 þ ghÞ ¼ ðaðun Þ þ cðaðun1 Þ  aðun ÞÞÞðun  un1 þ ghÞ P aðun Þgh: It follows from (29) that

ðaðunþ1 Þ þ cðaðun Þ  aðunþ1 ÞÞÞðunþ1  un þ ghÞ P aðun Þgh: If c ¼ 1, (34) directly implies (32). Assume now that c < 1. Expanding

aðunþ1 Þ ¼ aðun Þ þ

@aðhn Þ ðunþ1  un Þ; @u

hn 2 ðun ; unþ1 Þ;

ð34Þ

134

O. Misiats, K. Lipnikov / Journal of Computational Physics 239 (2013) 123–137

the inequality (34) reads as

  @aðhn Þ @aðhn Þ gh þ ð1  cÞ ðunþ1  un Þ P 0; ðunþ1  un Þ aðun Þ þ ð1  cÞ @u @u

ð35Þ

nÞ which is a quadratic inequality in unþ1  un . Since ð1  cÞ @aðh 6 0, the roots of the quadratic inequality (35) are positive, i.e. @u (32) holds, if and only if

aðun Þ þ ð1  cÞ

@aðhn Þ gh P 0: @u

ð36Þ

Assume now that un1 P un . Repeating the above calculations, we may conclude as before that if c ¼ 1, (33) holds for any h. If c < 1, we arrive at a similar quadratic inequality

  @aðhn Þ @aðhn Þ gh þ ð1  cÞ ðunþ1  un Þ 6 0: ðunþ1  un Þ aðun Þ þ ð1  cÞ @u @u If

@aðhn Þ @u

¼ 0, the inequality (33) clearly holds. Otherwise, we may guarantee that (33) holds if

8 nÞ < aðun Þ þ ð1  cÞ @aðh gh P 0; @u : junþ1  un j 6 

aðun Þ @aðhn Þ @u

ð1cÞ

ð37Þ

 gh:

Note that the second condition in (37) imposes the restriction only on the steepness of the discrete gradient, without any a priori knowledge of the gradient direction. The proof of the proposition is complete. h Corollary 1. The condition (36) is the necessary condition for the discretization to be monotone. The conditions (37) are the sufficient conditions for the discretization to be monotone. The conditions (37) enable us to develop an adaptive monotone FV scheme for (1). This scheme is constructed as follows. In practice, for every n ¼ 1; . . . ; N  1, we check the conditions (37) with c ¼ 12 (arithmetic average) and hn ¼ un at n  1; n and n þ 1. This leads to considering the following four cases. 1. If the conditions (37) hold for n  1; n, and n þ 1 choose

jnþ12 ¼

Kðxn ÞKðxnþ1 Þ ðaðunþ1 Þ þ aðun ÞÞ; Kðxn Þ þ Kðxnþ1 Þ

jn12 ¼

Kðxn ÞKðxn1 Þ ðaðun1 Þ þ aðun ÞÞ: Kðxn Þ þ Kðxn1 Þ

2. If the conditions (37) fail for n  1 but hold for n and n þ 1 choose

jnþ12 ¼

Kðxn ÞKðxnþ1 Þ ðaðunþ1 Þ þ aðun ÞÞ; Kðxn Þ þ Kðxnþ1 Þ

jn12 ¼

2Kðxn ÞKðxn1 Þ aðun1 Þ: Kðxn Þ þ Kðxn1 Þ

3. If the conditions (37) hold for n  1 and n but fail for n þ 1 choose

jnþ12 ¼

2Kðxn ÞKðxnþ1 Þ aðun Þ; Kðxn Þ þ Kðxnþ1 Þ

jn12 ¼

Kðxn ÞKðxn1 Þ ðaðun1 Þ þ aðun ÞÞ: Kðxn Þ þ Kðxn1 Þ

jn12 ¼

Kðxn ÞKðxn1 Þ aðun1 Þ: Kðxn Þ þ Kðxn1 Þ

4. In all other possibilities choose

jnþ12 ¼

2Kðxn ÞKðxnþ1 Þ aðun Þ; Kðxn Þ þ Kðxnþ1 Þ

Convergence of the adaptive FV scheme is illustrated in Table 5. Note that while the discretization remains monotone on an extremely coarse grid, N ¼ 8, for N ¼ 96 and higher the errors in the above discretization coincide with the errors produced by the scheme (HA). This is natural to expect since for any nonlinearity a there exists h0 > 0 such that 8h < h0 the conditions (37) hold true. The experimental convergence rate indicates that the scheme is of the second order of accuracy;

Table 5 Accuracy and monotonicity of the adaptive FV scheme for (28) with g ¼ 2 and g ¼ 10. g

Monotone for N ¼ 8

Norm

2

Yes

l2

3.14 10

l1

1.53 101

1.15 102

2

2

10

Yes

N ¼ 48

N ¼ 96 2

l2

9.88 10

l1

3.79 101

1.73 10

1.19 10

N ¼ 192 3

1.08 101

N ¼ 384

Rate

4.85 10

4

1.29 10

2.56

3.22 103

9.01 104

2.41

3

3

4

3.81 10

1.14 10

2

4.80 102

1.65 102

1.47

O. Misiats, K. Lipnikov / Journal of Computational Physics 239 (2013) 123–137

135

in fact, in some instances the experimental convergence rate is greater than 2, which is a artificial effect due to the switch in schemes for the lower values of N. Let us now illustrate how the above idea can be used to obtain a monotone FV scheme for the van Genutchen–Mualem model. Consider the van Genutchen–Mualem model with parameters n ¼ 2 and m ¼ 1=n, which give the following form of nonlinearity:

( aðuÞ ¼

5

ð1þu2 Þ4 ð

1;

p1 ffiffiffiffiffiffiffiffi

1þu2 þuÞ2

;

u > 0; ð38Þ

u 6 0:

Since the function aðuÞ is concave when u is positive, we may verify the sufficient conditions only at the points n; n  1 and n þ 1. Fig. 8 compares the FV scheme (HA) to the adaptive monotone FV scheme described above on a mesh with 24 points. The scheme (HA) produces a solution that has both overshoots and undershoots. The adaptive FV scheme fixes this problem. Note that both solutions are linear when x P 8:4 where Richards’ equation reduces to a linear elliptic equation. 5.2. Multi-dimensional problems Now, we generalize the results in Proposition 2 to Eq. (3) with d > 1. In order to simplify the presentation and illustrate the numerical results, we assume that d ¼ 2; however, similar reasoning will hold true for any d P 2). Along with the evolution Eq. (3), we consider the stationary Richards’ equation

  ! div KðxÞaðuÞðru þ ge1 Þ ¼ 0:

ð39Þ

Proofs of the existence and uniqueness results for Eqs. (3) and (39) use the same arguments as that in the one-dimensional case (see Section 2). Lemma 1. Let KðxÞ be constant in subdomains K 1 and K 2 . Let u be the solution of (39). Then u does not have local maxima or local minima inside X1 or X2 .

@a Proof. Recall that aðuÞ is absolutely continuous, so that @u exists almost everywhere. Then u solves the following equation in

X1 and X2 : div ðaðu Þru Þ þ g

@a  @u ðu Þ ¼ 0: @u @x1

For u being a fixed function, consider the linear equation

div ðaðu Þrv Þ þ g

@a  @ v ¼ 0: ðu Þ @u @x1

ð40Þ

Eq. (40) has a unique solution v  in X1 and X2 satisfying the maximum principle in these domains. Consequently, v  cannot have local maxima or minima in X1 and X2 . On the other hand, the uniqueness of the solution of the linear Eq. (40) implies that v   u and thus u cannot have local maxima or minima. h

Fig. 8. Harmonic and combined schemes for N ¼ 24 and aðuÞ given by (38).

136

O. Misiats, K. Lipnikov / Journal of Computational Physics 239 (2013) 123–137

Fig. 9. Left: discrete local minimum for Scheme I (harmonic average). Right: monotone discretization with Scheme II (mixed harmonic-upwind).

Consider a square mesh and denote ukn :¼ uðxn1 ; xk2 Þ. A FV scheme for (39) results in a five-point stencil in X:

 kþ1   k   k   k  k1 unþ1  ukn u  ukn1  ukn kþ1 u k1 u  un  jn 2 n ¼0 þ g  jkn1 n þ g þ jn 2 n 2 2 h h h h

jknþ1

ð41Þ

with

jknþ1 ¼ 2

2Kðxkn ÞKðxknþ1 Þ Kðxkn Þ þ Kðxknþ1 Þ

ðcaðukn Þ þ ð1  cÞaðuknþ1 ÞÞ

2Kðxkn ÞKðxkþ1 n Þ Kðxkn Þ þ Kðxkþ1 n Þ

ðgaðukn Þ þ ð1  gÞaðukþ1 n ÞÞ

and kþ1

jn 2 ¼

for 0 6 k; n 6 N. Our objective is to study monotonicity of the discretization (41) in the domains, where the analytic solution of (39) is monotone, namely, in X1 and X2 . Definition 3. The discretization (41) is called monotone if for all 1 6 k; n 6 N  1 k1 k k k kþ1 k1 minfuknþ1 ; ukn1 ; ukþ1 n ; un g 6 un 6 maxfunþ1 ; un1 ; un ; un g:

The following proposition provides the necessary and sufficient conditions for the monotonicity of the discretization (41). Proposition 3. The condition

aðukn Þ þ ð1  cÞ

@aðhkn Þ gh P 0; @u

hkn 2 ðukn ; uknþ1 Þ

ð42Þ

is the necessary condition for monotonicity of (41). The conditions

8 @aðhk Þ > < aðukn Þ þ ð1  cÞ @un gh P 0; k k > : junþ1  un j 6 

aðukn Þ @aðhk Þ ð1cÞ @un

 gh;

hkn 2 ðukn ; uknþ1 Þ; @aðhkn Þ @u

– 0;

ð43Þ

are the sufficient conditions for monotonicity of (41). Proof. Assume that K 1 ¼ 1. We first derive the condition that guarantees the absence of discrete local minima in (41). Taking the symmetry into account, the question of the absence of local minima reduces to the following two auxiliary problems. Case 1. Assuming that ukn1 P ukn ; unk1 P ukn and ukþ1 P ukn ; 1 6 k; n 6 N, find the condition under which uknþ1 6 ukn . n Case 2. Assuming that ukn1 P ukn ; unk1 P ukn and uknþ1 P ukn ; 1 6 k; n 6 N, find the condition under which unkþ1 6 ukn . Arguing as in the proof of Proposition 2, Case 1 leads to the conditions (43). Case 2 does not impose any additional conditions on c and g. Indeed, using the monotonicity of aðuÞ, it is straightforward to verify that under the assumptions of Case 2, ukþ1 6 ukn for all 0 6 c; g 6 1 and g; h > 0. The derivation of the condition that guarantees the absence of discrete local n maxima follows the same lines and leads to the condition (42). h

O. Misiats, K. Lipnikov / Journal of Computational Physics 239 (2013) 123–137

137

The following numerical simulations support the above arguments. Let X ¼ ½0; 102 with X1 ¼ ½0; 5  ½0; 10 and X2 ¼ ½5; 10  ½0; 10. As before, we set K 1 ¼ 0:5 and K 2 ¼ 2 in X1 and X2 respectively, and g ¼ 10. We proceed with implementing two schemes on a 24  24 point mesh: Scheme I: harmonic averaging in both x1 and x2 :

jknþ1 ¼ 2

2Kðxkn ÞKðxknþ1 Þaðukn Þaðuknþ1 Þ ; Kðxkn Þaðukn Þ þ Kðxknþ1 Þaðuknþ1 Þ

1

2 jkþ ¼ n

k kþ1 2Kðxkn ÞKðxkþ1 n Þaðun Þaðun Þ : k k kþ1 Kðxn Þaðun Þ þ Kðxn Þaðukþ1 n Þ

Scheme II: harmonic averaging in x2 and harmonic - upwind scheme in the direction x1 :

jknþ1 ¼ 2

2Kðxkn ÞKðxknþ1 Þ aðukn Þ; Kðxkn Þ þ Kðxknþ1 Þ

1

2 jkþ ¼ n

k kþ1 2Kðxkn ÞKðxkþ1 n Þaðun Þaðun Þ : kþ1 Þ Kðxkn Þaðukn Þ þ Kðxkþ1 Þaðu n n

The first scheme clearly has a discrete local minimum in X1 (undershoot), while the second scheme is monotone in both

X1 and X2 . 6. Conclusion We performed the theoretical and numerical analysis of the one-dimensional and two-dimensional Richards’ equation. This analysis allowed us to develop a computationally efficient monotone FV scheme of second-order accuracy. The scheme calculates transmissibility coefficients based on derived sufficient conditions for monotonicity. We also presented a few analytical solutions (in a closed-form) for a two-layer media that can be used to verify numerical codes. 7. Acknowledgments This work was carried out under the auspices of the National Nuclear Security Administration of the U.S. Department of Energy at Los Alamos National Laboratory under Contract No. DE-AC52-06NA25396 and the DOE Office of Science Advanced Scientific Computing Research (ASCR) Program in Applied Mathematics. References [1] A. Shcheglov, On the monotonicity of the solution of a mixed problem for a quasilinear heat equation with a discontinuous coefficient, J. Comput. Math. Math. Phys. 36 (6) (1996) 86–94. [2] A. Warrick, Numerical approximation of Darcian flow through unsaturated soil, Water Resour. Res. 27 (1991) 215–222. [3] D.A. Barry, J.-Y. Parlange, G.C. Sander, M. Sivaplan, A class of exact solutions for Richards’s equation, J. Hydrol. 142 (1993) 29–46. [4] T.J. Barth, D.C. Jespersen, The Design and Application of Upwind Schemes on Unstructured Meshes, AIAA Paper 89G0366, 1989. [5] H. Berninger, Domain Decomposition Methods for Ellipitc Problems with Jumping Nonlinearities and Application to the Richards Equation, Ph.D. Thesis, 2007. http://page.mi.fu-berlin.de/haiko/dissBerninger.pdf. [6] P. Broadbridge, I. White, Constant rate rainfall infiltration – A versatile nonlinear model. 1. Analytical solution, Water Resour. Res. 24 (1) (1988) 145– 154. [7] R. Brooks, A. Corey, Hydraulic Properties of Porous Media, Hydrol. Pap. 3, Colorado State Univ., 1964. [8] Raimund Bürger, Anı´bal Coronel, Mauricio Sepúlveda, On an upwind difference scheme for strongly degenerate parabolic equations modelling the settling of suspensions in centrifuges and non-cylindrical vessels, Appl. Numer. Math. 56 (10–11) (2006) 1397–1417. [9] G. Manzini, S. Ferraris, Mass-conservative finite volume methods on 2-D unstructured grids for the Richards equation, Adv. Water Res. 27 (2004) 1199– 1215. [10] H.-G. Roos, A second order monotone upwind scheme, Computing 36 (1–2) (1986) 57–67. [11] J. Caputo, Y. Stepanyants, Front solutions of Richards equation, Transp. Porous Med. 74 (1) (2008) 1–20. [12] J. Gasto, J.M. Gripoll, Y. Cohen, Estimation of internodal permeabilities for numerical simulation of unsaturated flows, Water Resour. Res. 38 (2002) 1– 10. [13] J. Vanderborght, R. Kasteel, M. Herbert, M. Javaux, D. thiery, M. Vanclooster, C. Mouvet, H. Vereecken, A set of analytical benchmarks to test numerical models of flow and transport in soils, Vadoze Zone J. 4 (2005) 206–221. [14] J. Zaidel, D. Russo, Estimation of finite difference interblock conductivities for simulation of infiltration into variably dry soils, Water Resour. Res. 28 (1992) 2285–2295. [15] D. Kuzmin, A vertex-based hierarchica slope limiter for p-adaptive discontinuous Galerkin methods, J. Comput. Phys. 233 (2010) 3077–3085. [16] F. Otto, L1 -contraction and uniqueness for quasilinear elliptic–parabolic equations, J. Diff. Equ. 131 (1) (1996) 20–38. [17] P. Smolarkiewicz, Multidimensional positive definite advection transport algorithm: an overview, Int. J. Numer. Methods Fluids 50 (10) (2006). [18] A.J. Pullan, The quasilinear approximation for unsaturated porous media flow, Water Resour. Res. 26 (6) (1990) 1219–1234. [19] R. Eymard, M. Gutnic, D. Hilhorst, The finite volume method for Richards equation, Comput. Geosci. 3 (1999) 259–294. [20] L. Richards, Capillary conduction of liquids through porous mediums, Physics 1 (5) (1931) 318–333. [21] R. Shrivastava, T.-C. Jim Yeh, Analytical solutions for one-dimensional, transient infiltration towards the water table in homogeneous and layered soils, Water Resour. Res. 27 (1991) 753–762. [22] S. Spekreijse, Multigrid solution of monotone second-order discretization of hyperbolic conservation laws, Math. Comput. 49 (179) (1987) 135–155. [23] M. van Genuchten, A closed-form equation for predicting the hydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. J. 44 (1980) 892–898. [24] B. van Leer, Upwind and high-resolution methods for compressible flow: from donor cell to residual-distribution schemes, Commun. Comput. Phys. 1 (2) (2006) 192–206. [25] A.W. Warrick, T.-C. Jim Yeh, One-dimensional, steady vertical flow in a layered solid profile, Adv. Water Res. 13 (4) (1990) 207–210. [26] T.-C. Jim Yeh, One-dimensional steady state infiltration in heterogeneous soils, Water Resour. Res. 25 (10) (1989) 2149–2158. [27] Y. Mualem, A new model for predicting the hydraulic conductivity of unsaturated porous media, Water Resour. Res. 12 (1976) 513–522.