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

Abstract

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.

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:

$$\frac{\partial u}{\partial t}=\frac{\partial }{\partial x}\left(K(x)a(u)\left(\frac{\partial u}{\partial x}+g\right)\right)\text{,}$$

where $u=u(t\text{,}x)\text{,}(t\text{,}x)\in [0\text{,}\infty )\times [0\text{,}2\alpha ]$ for some $\alpha >0$, $K(x)$ is a piecewise constant positive function,

$$K(x)=\left\{\begin{array}{c}{K}_1\text{,}x\in [0\text{,}\alpha ]\text{,}{K}_2\text{,}x\in (\alpha \text{,}2\alpha ]\text{,}\hfill \end{array}\right.$$

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 $\mathbb{R}$ satisfying the ellipticity condition

$$\frac{1}{1+|u{|}^{\sigma }}⩽a(u)⩽1\text{,}u\in \mathbb{R}\text{,}\sigma >0\text{.}$$

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, 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

$$\frac{\partial u}{\partial t}=\frac{\partial }{\partial x}\left(u^s\frac{\partial u}{\partial x}\right)\text{,}s>0\text{,}x⩾0\text{,}$$

subject to the homogeneous initial condition $u(0\text{,}x)=0$ and the boundary conditions $u(t\text{,}0)={\mathit{ct}}^{\frac{1}{s}}$ and $u(t\text{,}+\infty )=0$. This equation has a unique solution in the form of a traveling wave

$$\left\{\begin{array}{c}u=(-\mathit{vs}(x-\mathit{vt}){)}^{\frac{1}{s}}\text{,}x>\mathit{vt}\text{,}v=\sqrt{\frac{c^s}{s}}\text{;}\hfill \\ u=0\text{,}x⩽\mathit{vt}\text{.}\hfill \end{array}\right.$$

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 Fig. 7, Fig. 8, Fig. 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 Fig. 7, Fig. 8, Fig. 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), (37)). These conditions are given in terms of the nonnegative parameter $0⩽\gamma ⩽1$, which measures deviation from an upwind scheme. We select different flux formulas, based on the range of values of $\gamma$ that provide a monotone scheme.

Next, we generalize the necessary and sufficient conditions for the multi-dimensional Richards’ equation

$$\frac{\partial u}{\partial t}=\mathrm{div}\left(K(x)a(u)\left(\nabla u+g\overrightarrow{e_1}\right)\right)\text{,}$$

where $\mathrm{\Omega }={\mathrm{\Omega }}_1\cup {\mathrm{\Omega }}_2\subset {\mathbb{R}}^d\text{,}d⩾1$, with ${\mathrm{\Omega }}_1$ and ${\mathrm{\Omega }}_2$ being connected domains with Lipschitz boundaries, $\overrightarrow{e_1}≔(1\text{,}0\text{,}\dots \text{,}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\in {\mathrm{\Omega }}_i$.

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).