An adjoint-based a posteriori analysis of numerical approximation of Richards equation

1.   Introduction

The subject of investigation in this paper is numerical solutions of the Richards Equation [23]. This equation is a governing mathematical principle for modeling the water flow in an unsaturated porous medium that is driven by the gravity and capillarity that disregards air flow. Since ability to construct closed form solutions to this equation is very limited (see for example [25,18,22,24] for some related effort on the subject), a reliance on numerical approximations is a necessity. However, even with the emergence of many advances of computing technology, this equation remains one of the most challenging problems in porous media flow and transport. Recent review on its numerical solutions can be found in [17].

There are several outstanding issues attributed to the challenge. Richards Equation is strongly nonlinear, which appears as the dependence of the soil unsaturated hydraulic conductivity ($\kappa$) and the water content ($\vartheta$) to the pressure head ($u$). Note that the presence of the water content in the equation is in terms of its temporal rate of change. Inclusion of gravity in the Darcy's velocity $q$, written as $q = -\kappa(u) \partial_{z} (u-z)$, can potentially create instability in the numerical solutions, in particular, when simulating dry soil conditions. The variable $z$ in the expression of $q$ denotes the vertical spatial coordinate, which is positive in the downward direction. It represents the influence of gravity to the flow. Furthermore, some of the more realistic scenario requires taking into account the soil heterogeneity in the simulations.

The mixed-form (or coupled-form) of Richards Equation in one dimensional soil column is written as follows:

Two typical boundary conditions are

Here we assume that $\vartheta : (-\infty, 0] \to (0, 1)$ and $\kappa : (-\infty, 0] \to (\kappa_{{{\rm{min}}}}, \kappa_s)$ with $\kappa_s>\kappa_{{{\rm{min}}}} >0$. The choice $(-\infty, 0]$ as the domain of $\vartheta$ and $\kappa$ is done to reflect the physical relevance that the pressure head $u$ is always nonpositive and Richards Equation typically models unsaturated flow.

The major theme in this paper is two-fold. One aspect centers on the development of a numerical approximation of Richards Equation in space-time finite element methods obtained from an appropriate variational formulation. Space-time finite element methods have been previously used for parabolic equations (see for example [20]) and for reaction-diffusion system (see for example [15]). A recent work on application of control volume finite element in combination with method of lines for solving Richards Equation is recorded in [10]. To the best of the author's knowledge, there has not been any attempt to apply space-time finite element methods technique to Richards Equation. In particular, a finite volume element spatial discretization (with linear finite element) is employed due to its inherent local mass conservation property. This is an important trait commonly desired and in some cases imperative in order to produce reliable numerical simulations of flow and transport in porous media (see for example [5,8,16]). The space-time variational formulation in combination with a certain variational crime in the form of numerical integration techniques would in turn yield implementable time marching schemes for approximate solution of Richards Equation.

The other aspect is concerned with an a posteriori error estimation of the resulting numerical approximation. In this regard, some investigations on a posteriori error analysis of numerical methods for Richards Equation are already available. Bause and Knabner [3] use adaptive mixed hybrid finite element discretizations to solve Richards Equation, where the adaptivity is performed under an a posteriori error indicator that is based on either superconvergence or residual of the approximation. Baron et al. [2] employ Discrete Duality Finite Volume (DDFV) scheme along with second-order backward differentiation formula to solve the equation. They derive an a posteriori error bound of the approximation using the equilibrated fluxes method. Bernardi et al. [4] perform a semi discretization of Richards Equation by finite element method and apply Backward Euler scheme to get the full discretization. Then a posteriori error bounds are derived that aim at distinguishing components of contribution of spatial discretization from temporal discretization.

In many practical situations, it may not be necessary to measure global property of the approximate solution. More often, an accuracy is desired only for some specified quantities of interest associated with the numerical approximation, which is usually expressed as a functional of the approximate solution. For this purpose, a suitable a posteriori analysis is based on duality, adjoint operators governing the generalized Green's function, and a variational formulation. This approach is adopted in this paper. It is also suitable because, as mentioned before, the numerical approximation of Richards Equation is based on certain variational equations. Utilizations of adjoint equations are not new (see for example [21] for an extensive exposition on their applications). On various roles of adjoint methodologies in performing a posteriori error estimations of numerical methods for differential equations, one can consult [19,15,11,13,12,14,1] and references therein.

The rest of this paper is organized as follows. A space-time variational equation governing the numerical approximation of Richards Equation is derived in Section 2. The description includes examples of application of numerical integration to produce the time marching schemes. Section 3 carries out the formulation of an adjoint-based a posteriori error analysis of the quantities of interest calculated using the approximate solution. Since the variational equation of the solution is nonlinear, an appropriate linearization is conducted that would make construction of the corresponding adjoint equation amenable. Some numerical examples to demonstrate performance of the numerical methods and the error estimators are shown in Section 4. Here much of the effort is devoted to illustrate global accuracy of numerical methods and reliability of the error estimators in terms of their capability to decompose the total error into relevant components. A comparison to the actual error in some specified quantities of interest is conducted. Finally, the conclusion and future work is discussed in Section 5.

2.   Finite volume element in space and finite element in time

In what follows, we assume that

are supplied to (1.1). Denoting

the solution of (1.1) supplied with (2.1) satisfies

Here $\langle \cdot, \cdot \rangle$ is the usual scalar product in $L^2(a, b)$, and $A: C[a, b] \times H^1_{\rm{D}} \times H^1_{\rm{D}} \to \mathbb{R}$ is defined as

The spatial domain $(a, b)$ is partitioned into a collection of $M$ subintervals ${\mathcal T}_h$, such that $\tau_j = (z_{j-1}, z_j) \in \mathcal{T}_h$, with length $h_j = z_j - z_{j-1}$, for $j = 1, \cdots, M$ and $(a, b) = \cup_{j = 1}^M \tau_j$, where $h = \max \{ h_j : 1\le j \le M \}$. On this $\mathcal{T}_h$, let

where $\phi_j(z)$ is the usual 'hat' function such that $\phi_j(z_i) = \delta_{ij}$.

2.1.   A brief excursion to finite volume element

The finite volume approximations rely on a local conservation property associated with the governing equation, in particular with respect to the second order differential operator in (1.1). To fix the idea, consider $\tau^* = (z_l, z_r) \subset (a, b)$ and apply fundamental theorem of calculus to get

The above $\tau^*$ is called a control volume. We choose $M$ control volumes. Specifically, given $z_j$, for $j = 1, \cdots, M-1$, we set $\tau^*_j = (z_{j-1/2}, z_{j+1/2})$, where $z_{j-1/2}$ is the mid-point of $\tau_{j}$ and $z_{j+1/2}$ is the mid-point of $\tau_{j+1}$. For $z_0$, we set $\tau^*_0 = (z_0, z_{1/2})$. Collection of these control volumes is denoted by $\mathcal{T}_h^*$. On this $\mathcal{T}_h^*$, let

where $\phi^*_j(z)$ is a piecewise constant function such that it is equal to 1 in $\tau^*_j$ and zero elsewhere. Set $J_h : H^1_{\rm{D}} \to \mathcal{Y}_h$ such that $[J_h v](z_j) = v(z_j)$, i.e., given $v \in H^1_{\rm{D}}$, $J_h v$ is a piecewise constant function over ${\mathcal{T}}_h^*$. A standard interpolation estimate suggests

where $\| \cdot \| = \sqrt{\langle \cdot, \cdot \rangle}$. By the Cauchy-Schwarz inequality, this implies

To express (2.3) in a variational setting, define $A_h : C[a, b] \times H^1_{\rm{D}} \times H^1_{\rm{D}} \to \mathbb{R}$ as

Exclusion of $a$ in the above equation is due to the Neumann boundary condition at that point. Next, we quantify the discrepancy of $A_h(\cdot ; \cdot, \cdot)$ from $A(\cdot ; \cdot, \cdot)$.

Proposition 2.1. Let $w \in C[a, b]$, $u, v \in H^1_\mathit{\rm D}$. Then

where

Furthermore, when $u, w \in \mathcal{X}_h$ and $\kappa^\prime(r)$ is bounded for every $r \in (-\infty, 0]$, then

where $C_\kappa = \sup_{u \in (-\infty, 0]} |\kappa^\prime(u)|$.

Proof. For any $\tau_j \in \mathcal{T}_h$, integration by parts gives

which when applied to $A(\cdot ; \cdot, \cdot)$ gives

For $j = 0$, set $K_j = \tau_j^*$. For $j = 1, \cdots, M-1$, fix a $\tau^*_j \in \mathcal{T}_h^*$ and suppose $\tau_{j}, \tau_{j+1} \in \mathcal{T}_h$ are such that $K_j = \tau_j \cap \tau^*_j = (x_{j-1/2}, x_j)$, $K_{j+1} = \tau_{j+1} \cap \tau^*_j = (x_j, x_{j+1/2})$. By fundamental theorem of calculus,

for $e = j, j+1$. Recognizing that

we may apply (2.10) in $A_h(\cdot ; \cdot, \cdot)$ to get

Subtraction of (2.9) from (2.11) and recalling that $J_h v (z_j) = v(z_j)$ gives (2.6).

Furthermore, when $w, u \in \mathcal{X}_h$, product rule of differentiation for $z \in \tau_j$ gives

so using this identity and the Cauchy-Schwarz inequality,

This completes the proof.

Remark 2.1. The foregoing exposition gives an indication that $\langle w, J_h v - v \rangle \to 0$ and $A_h(w ; u, v) \to A(w ; u, v)$ as $h \to 0$. This will play a role later on in the a posteriori error analysis. Various estimates such as described in (2.4), (2.5), and (2.8) have been established in several literatures on finite volume element methods (see for example [6,7,9] and references therein).

2.2.   A variational equation for the approximation

In a similar fashion to the spatial variable, we partition $[0, T]$ into a collection of subintervals $\mathcal{I}_k$, such that $I_n = [t_{n-1}, t_n] \in \mathcal{I}_k$ with time step $k_n = t_n - t_{n-1}$ and $[0, T] = \cup_{n = 1}^N I_n$, where $k = \max \{ k_n: 1\le n \le N \}$. We denote the jump of a function $w(\cdot, t)$ across $t_n$ by $[w]_n = w_n^+ - w_n^-$, where $w_n^+ = \lim_{s\rightarrow t_n^+} w(\cdot, s)$ and $w_n^- = \lim_{s\rightarrow t_n^-} w(\cdot, s)$. On every space-time slab $[a, b] \times I_n$, the approximate solution is sought in a functional space that contains functions that are piecewise linear polynomial in spatial variable and polynomial of degree $q$ in temporal variable. In particular, we define

We denote by $\mathcal{W}^q_h$ the space of functions defined on $[a, b] \times [0, T]$ such that restriction of $w \in \mathcal{W}^q_h$ to $[a, b] \times I_n$ belongs to $\mathcal{W}_h^q(I_n)$. The approximation amounts to finding $\widetilde{u} \in \mathcal{W}_h^q$ that is governed by

where

Notice that the first equation in (2.13) is a global formulation in that the integration is over $(0, T)$. While an implementation based on this formulation is possible, it is perhaps more amenable to construct an implementation that is local over $I_n$ for $n = 1, \cdots, N$. The corresponding equation for this formulation can be derived from (2.13) by choosing $w \in \mathcal{W}_h^q$ such that $w |_{(a, b)\times I_n} = v \in \mathcal{W}^q_h(I_n)$ and it is zero everywhere else, which yields

2.3.   Some examples

In what follows, we describe two specific examples that transform (2.15) into computable algebraic schemes.

2.3.1.   FVEM in space - dG0 in time

Here $\tilde{u} \in \mathcal{W}^0_h$, i.e.,

which for every $w_0 \in \mathcal{X}_h$ is governed by

for $n = 1, \cdots, N$. Notice that (2.17) is mimicking the Backward Euler difference scheme with $v_{0, n} \in \mathcal{X}_h$ being the unknown function to be solved. In particular, setting $(U_{0, n}, U_{1, n}, \cdots, U_{M-1, n}) = \boldsymbol{U}_n \in \mathbb{R}^{M}$, such that

then $\boldsymbol{U}_n$ is governed by

where $G : \mathbb{R}^M \to \mathbb{R}^M$ such that $G_i : \mathbb{R}^M \to \mathbb{R}$, for $i = 0, 1, \cdots, M-1$, is constructed from the left hand side of (2.17) by replacing $w_0$ by $\phi_i$. Here, the dependence on $\boldsymbol{U}_n$ is realized through (2.18). Clearly (2.19) is a nonlinear algebraic system of equations governing $\boldsymbol{U}_n$.

2.3.2.   FVEM in space - dG1 in time

Here $\tilde{u} \in \mathcal{W}^1_h$, i.e.,

which implies that $\widetilde{u}_{n-1}^+ = v_{0, n} + t_{n-1} v_{1, n}$ and $\widetilde{u}_{n}^- = v_{0, n} + t_n v_{1, n}$. We may equivalently write

Choosing $w = \frac{t_n - t}{k_n} \, \psi_{n-1}^+$ with $\psi_{n-1}^+ \in \mathcal{X}_h$, and using integration by parts along with acknowledging some cancellations,

In a similar fashion, using $w = \frac{t - t_{n-1}}{k_n} \, \psi_{n}^-$ with $\psi_{n}^- \in \mathcal{X}_h$ yields

Thus $\widetilde{u} \in \mathcal{W}^1_h$ satisfies (2.21) and for every $n = 1, \cdots, N$, it is governed by

where

Setting

and

and $\boldsymbol{U}_n = (\boldsymbol{U}^+_{n-1}, \boldsymbol{U}^-_n) \in \mathbb{R}^{2M}$, then (2.9) yields

where $G:\mathbb{R}^{2M} \to \mathbb{R}^{2M}$, with $G = (G^+, G^-)$, and $G^+ : \mathbb{R}^{2M} \to \mathbb{R}^M$ and $G^- : \mathbb{R}^{2M} \to \mathbb{R}^M$ such that $G^+_i : \mathbb{R}^{2M} \to \mathbb{R}$ and $G^-_i : \mathbb{R}^{2M} \to \mathbb{R}$ are respectively constructed from the left hand side of (2.9) by setting $\psi^+_{n-1} = \psi^-_n = \phi_i$, for $i = 0, 1, \cdots, M-1$.

2.4.   A variational crime by numerical integrations

The preceding description is a derivation of algebraic equations governing the approximation that is faithful to the variational equation (2.15) and the choice of polynomial degree of the temporal variable. Still, for a completely implementable scheme, one must rely on further approximation of the integrations appeared in (2.17) and (2.9). In the current setting, there are two integrations that need to be approximated: the spatial integration $\langle \cdot , \cdot \rangle$ and the temporal integration $\int_{I_n} \cdot \, {\rm d} t$. Utilization of various numerical integration techniques are pretty common. In particular, in the standard finite element methods for typical steady state problems, forms/functionals in the variational equations, which are expressed as integrations of spatial variables, are approximated by various Gaussian quadratures. This is clearly applicable for $\langle \cdot, \cdot \rangle$. To minimize the associated pollution to the global accuracy of the approximation, the numerical integrations must be chosen such that the degree of their errors is of similar order to the errors corresponding to $W^q_h$.

Furthermore, what is more crucial in this case is how $\int_{I_n} \cdot \, {\rm d} t$ is to be approximated. We note that the only temporal integration in (2.17) is the one associated with the Neumann condition $g_a(t)$, and for this a right hand point rule resulting in

is adequate.

Derivation of numerical integrations for (2.9) is a bit more involved. A viable option is the following two point Gaussian quadrature

With this, set

where $\gamma = \frac{1+\sqrt{3}}{2\sqrt{3}}$. The resulting approximations $G^+_{n, i} \approx G^+_i$ and $G^-_{n, i} \approx G^-_i$ are expressed as

where $c^+_l = c^-_r = \gamma$, $c^+_r = c^-_l = 1-\gamma$, and $\delta_{ij}$ is the usual Kronecker delta.

The approach proceeds with the construction of algebraic equations for $(\widetilde{u}_{1, n}, \widetilde{u}_{2, n})$ where $\widetilde{u}^-_n$ appearing on the second equation above is represented as

which is obtained from (2.26). Thus, with $(G^+_{n, i}, G^-_{n, i})$ replacing $(G^+_i, G^-_i)$, (2.24) is solved to get $(\widetilde{u}_{1, n}, \widetilde{u}_{2, n})$, after which $\widetilde{u}^-_n$ is recovered from (2.27).

3.   An adjoint-based a posteriori error analysis

In many realistic situations, it is often desirable to achieve an acceptable level of accuracy of a numerical approximation in some quantities of interest. Relevant examples include average water content over a certain region and at some time instances or the water content at some locations. Along this line of argument, it may be computationally infeasible as well as very inefficient to attempt to control the error in a global fashion when all that is required is accuracy on those aforementioned quantities. A practical alternative is to estimate the error of the numerical approximation in the specified quantity of interest, whose representation is expressed as a functional of $u$:

for given data $\psi_T:(a, b) \to \mathbb{R}$, $\psi:(a, b) \times (0, T) \to \mathbb{R}$, and $\psi_a : (0, T) \to \mathbb{R}$. If one wants to quantify the (averaged) water content at time $t = T$, then (3.1) uses $\psi = 0$, $\psi_a = 0$, and $\psi_T$ is set to be a piecewise constant function in the spatial variable that reflects the desired nature of the average quantity. On the other hand, if an accumulated water content is the quantity of interest to be approximated, then $\psi_T = 0$, $\psi_a = 0$, and $\psi = 1$ in (3.1).

To derive the error in approximating $Q(u)$, we use a generalized Green's function that solves the adjoint problem corresponding to a special choice of (adjoint) data $\psi_T$, $\psi$, and $\psi_a$, as illustrated by the description in the previous paragraph. As alluded to earlier, the formulation of an adjoint problem broadens the applications of Green's functions (see for example [19,15,1] and references therein). However, an adjoint operator formally corresponds to a linear operator. Since Richards Equation and the associated numerical approximation are nonlinear problem, we must perform a linearization, after which the adjoint problem is built to correspond to that linearized representation.

The nonlinearity in Richards Equation stems from $\theta(u)$ and $\kappa(u)$, and that is where the linearization effort is concentrated on. To this end, letting

the Mean Value Theorem for integral gives

Furthermore, setting $F : H^1(a, b) \to \mathbb{R}$ by

its Fréchet derivative is $F^\prime(w) : H^1(a, b) \to \mathbb{R}$ such that

Utilizing again the Mean Value Theorem for integral, one gets

where

At this stage, we are in a position to formulate the adjoint problem. For $t \in [T, 0]$, let $\varphi(\cdot, t) \in H^1_{\rm{D}}$ satisfy

for every $w(\cdot, t) \in H^1_{\rm{D}}$. Here $\varphi$ is solution to the adjoint problem, which is governed by a linear problem as stated in (3.8). The two theorems below state the quantification of error in the approximation of $Q(u)$, which is written in terms of residuals of $\widetilde{u}$ weighted against $\varphi$. In what follows, we use

Theorem 3.1. For $\widetilde{u} \in \mathcal{W}^q_h$ in (2.13) and $u$ in (2.2),

where

with $R_{h, n}(\widetilde{u} ; \widetilde{u}, \varphi)$ as expressed in (2.14), and

Proof. Substitute $w = e = u - \widetilde{u}$ in (3.8) so that

where we have used the first equation in (2.2). Since

integration of (3.8) over $I_n$ yields

where $R_n$ is as expressed in (3.9). Next we sum up (3.15) over $n = 1, \cdots, N$ and take advantage of the continuity of $\varphi$ in $t$ to get

The residual $R_n(\widetilde{u} ; \widetilde{u}, \varphi)$ can be decomposed into three components, which is obtained by adding and subtracting $R_{h, n}(\widetilde{u} ; \widetilde{u}, \varphi)$:

where

Putting(3.18) to (3.17) and in turn to (3.16) completes the proof.

Theorem 3.2. For $\widetilde{u} \in \mathcal{W}^q_h$ in (2.15) and $u$ in (2.2),

where

and $\pi_h^q \varphi \in \mathcal{W}_h^q$ is the usual projection of $\varphi$ onto $\mathcal{W}_h^q$.

Proof. Most derivation steps follow the proof of Theorem 3.1 up to (3.16):

At this stage, we intend to insert (2.15), which is valid when the test function is $\pi_h^q \varphi \in \mathcal{W}^q_h$. To do so, add and subtract $R_n(\widetilde{u} ; \widetilde{u}, \pi_h^q \varphi)$ so that

where similar equation to (3.18) has been used, and $\varepsilon_{A, n}$ and $\varepsilon_{h, n}$ are as in (3.12). Putting all these results back to (3.21) completes the proof.

4.   Numerical examples

Several numerical examples are presented in this section to achieve two goals: 1) to investigate the global/norm-based accuracy of the proposed approximation, and 2) to validate the robustness of error estimators that are derived from Theorem 3.1 and Theorem 3.2. While the former cannot satisfactorily substitute for a rigorous a priori error analysis, at least it should give an illustrative indicator on the global convergence property of the approximation. With respect to the latter, we only concentrate on the estimators' accuracy in predicting error and their capability to decompose it into relevant components. Various pertinent applications of the proposed error estimators to other aspects in numerical simulation of Richards Equation, such as its role in adaptivity, will be a subject of future work.

A uniform set of discretization parameters $h_i = h = (b-a)/M$ and $k_n = k = T/N$ is used to construct the algebraic equations (2.19). The time marching is executed by solving this system using the standard Newton's method of iteration.

4.1.   A problem with closed form solution

While the proposed procedures enjoy a flexibility in their implementation, a closed form solution of Richards Equation is needed for the purpose of assessing their performance. As alluded to in the introduction, it is only on a very rare occasion that a closed form solution of Richards Equation is available. One such instance is when the constitutive relations are expressed as

where $\kappa_s$ is the saturated hydraulic conductivity, $\vartheta_r$ and $\vartheta_s$ are respectively the residual and saturated water content, and $\alpha$ is the reciprocal of vertical height associated with the capillary fringe. When $g_\star(t) = g_\star = {\rm{constant}}$, $\star = a, b$, then the closed form solution can be expressed as a series representation:

where

The pair $\{\phi_n, \lambda_n\}_{n = 1}^\infty$ constitutes an eigenfunction and an eigenvalue that satisfies

where $\tilde{\mathcal{B}}_\star$ are boundary conditions for $w$, which are appropriately derived from $\mathcal{B}_\star$ via the relation $w(z, t) = \kappa_s e^{\alpha u(z, t)}$. The constants $\{C_1, C_2\}$ are obtained from imposing the boundary conditions $\mathcal{B}_\star u = g_\star$.

Two examples are considered in the numerical experiments whose data are listed in Table 4.1. Solution profiles of these examples are shown in Figure 4.1 and Figure 4.2. The axes on these figures are flipped to follow the plotting style for profiles associated with Richards Equation (see for example [25,24]). The initial condition is

where for Ex. 1,

and for Ex. 2,

4.2.   An accuracy assessment of the approximation

In this subsection, a set of numerical experiments to investigate accuracy of the approximation is presented. We solve the two examples whose data are listed in Table 4.1.

Table 4.2 and Table 4.3 list the errors of approximation $\vartheta(\widetilde{u}(T))$ in $L^2(a, b)$-norm for Ex. 1 and Ex. 2, respectively. Four different number of elements ($M = 12, 24, 48, 96$) and four different number of time steps ($N = 1, 2, 4, 8$) are used to collect the error data in Table 4.2, while for error data in Table 4.3, $M = 5, 10, 20, 40$ and $N = 4, 8, 16, 32$ are used.

First, the accuracy of FVEM-dG1 generally outperforms that of FVEM-dG0, which is especially evident when solving Ex. 1 (see Table 4.2). The approximation error for Ex. 1 seems to be dominated by component of the temporal discretization. For a fixed $N$, refining $M$ does not quite improve the accuracy. However, for a fixed $M$, refining $N$ by two roughly reduces the error by two for FVEM-dG0 and by seven to ten for FVEM-dG1, especially for larger $M$. A strikingly different result is observed for Ex. 2 (see Table 4.3), for which the spatial discretization error component is more dominant. For a fixed $N$, the error of FVEM-dG1 shows a quadratic convergence with respect to $M$. On the other hand, when $N$ is still small, the error of FVEM-dG0 resembles a first order convergence with respect to $M$. As $N$ is increased, a better convergence rate is obtained.

4.3.   Performance of the a posteriori error estimators

As mentioned, the error equation for a quantity of interest $Q$ stated in Theorems 3.1 and 3.2 can be used to derive fully computable error estimators for the approximate solution $\widetilde{u}$. Notice that adjoint equation (3.8) is formulated based on linearization that utilizes Mean Value Theorem on the path $\widetilde{u}_\sigma = \widetilde{u} +\sigma(u-\widetilde{u})$, for $\sigma \in [0, 1]$, cf. (3.2). Since in reality $u$ is not available, calculation of solution to the adjoint equation (3.8) must be done using the only available information, namely, the approximate solution $\widetilde{u}$, so in practice, $\widetilde{u}_\sigma \approx \widetilde{u}$. The adjoint $\varphi$ is approximated preferably using higher order approximation than the one used to produce $\widetilde{u}$.

To test the proposed error estimators, we consider two quantities of interest:

● The spatial average of water content at time $T$, which is represented as

To calculate the adjoint solution associated with this quantity, the corresponding adjoint data is $\psi_T = 1/(b-a)$, and the rest are zero.

● The total average of water content over $(a, b)\times(0, T)$, which is expressed as

The corresponding adjoint data is $\psi = 1/(b-a)/T$ and the rest are zero.

The true values of these quantities of interest for the two examples are listed in Table 4.4.

The approximate solution is $\widetilde{u} \in \mathcal{W}_h^0$ (FVEM-dG0). The adjoint solution is solved by continuous and piecewise quadratic finite element in spatial variable and continuous piecewise linear in temporal variable. Profiles of $\varphi$ corresponding to each of these quantities of interest for each of the problems are shown in Figure 4.3 and Figure 4.4, respectively.

Table 4.5 and Table 4.6 demonstrate performance of the error estimator for the above quantities of interest when solving Ex. 1. In these tables, the error has been decomposed according to the components of error listed in Theorem 3.1. As in the accuracy assessment results, four different number of elements ($M = 12, 24, 48, 96$) and four different number of time steps ($N = 1, 2, 4, 8$) are used. The last column, which is labeled by Eff. denotes the ratio of error estimator (Err. Est.) to the actual error, so the closer Eff. is to 1 indicates a more accurate error estimator.

Notice that values in the tables give an indication that the error in $Q$ is dominated by the contribution from temporal discretization, as refinement of the spatial mesh for a fixed time step only causes negligible reduction in Err. Est. Reducing the time step by two seems to reduce the error by two, which demonstrates an asymptotic first order convergence with respect to the time step, i.e., ${\rm{Err.\;\;Est.}} = \mathcal{O}(k)$. Prominence of temporal discretization effect to the total error makes sense due to the longer time simulation.

The component $E_0$ quantifies the quality of representation of the true initial condition $u_0$ in the simulation. Representation of $u_0$ is realized through the projection of $u_0$ onto $\mathcal{X}_h$. Thus, $E_0$ measures the discrepancy attributed to this choice, which shows a second order convergence with respect to $h$ (i.e., $E_0 = \mathcal{O}(h^2)$) for a fixed time step. Comparison of values in the tables indicates that relative contribution of this component to the overall error is less significant.

The component $E_1$, which measures the residual of the finite volume element discretization weighted by the adjoint solution and integrated over time, shows a decrease with respect to time step as well. In fact, $E_1$ is clearly the main contributor to the total error with asymptotic behavior $E_1 = \mathcal{O}(k)$. The components $E_2$ and $E_3$ measures the discrepancy between the variational setting associated with finite volume element and that of standard continuous Galerkin finite element. In the realm of a priori error analysis, these two components are bounded in terms of the spatial mesh size $h$ (see Proposition 2.1 and (2.5)). For every fixed $N$, there is a reduction of $E_2$ and $E_3$ as $h$ is refined, with mostly $E_2 = \mathcal{O}(h)$ in Table 4.5, $E_2 = \mathcal{O}(h^2)$ in Table 4.6, and $E_3 = \mathcal{O}(h^2)$ in both tables. However, these two components are less dominant in relative comparison to $E_1$. Notice also that $E_0$ has a different sign than the rest of components. A capability to distinguish components of error and to recognize potential cancellation is arguably one of the strong advantages of the adjoint-based error estimation techniques. Finally, as illustrated by the Eff., as refinement is performed, the error estimator gives a more accurate prediction.

Next we utilize the proposed error estimator to the numerical solution of Ex. 2. Table 4.7 and Table 4.8 show the breakdown of error for each of the quantities of interest. Again, following the accuracy assessment results, four different number of elements ($M = 5, 10, 20, 40$) and four different number of time steps ($N = 4, 8, 16, 32$) are used. As in the results for Ex. 1, the proposed estimator performs really well in predicting the error. However, upon a closer observation, the detailed situation is quite different from what happened in Ex. 1. In Table 4.7 (error associated with $Q$ in (4.9)), we notice that the error caused by discretization of the spatial variable is more dominant, so refining the spatial mesh reduces the error. In particular, $E_3$ is seen to be the main contributor to the total error with asymptotic behavior $E_3 = \mathcal{O}(h^2)$, which results in ${\rm{Err.\;Est.}} = \mathcal{O}(h^2)$. Notice also that time step refinement does not seem to improve the accuracy.

The result in Table 4.8 shows that $E_1$ and $E_3$ are two competing error components with $E_1 = \mathcal{O}(k)$ and $E_3 = \mathcal{O}(h^2)$. Due to their different sign, they tend to cancel each other, especially when $M$ and $N$ are still smaller. This in turn lowers the magnitudes of total error. However, as $N$ is increased, $E_3$ tends to be more dominant than $E_1$, especially for small $M$. Since there is an intertwinement of the error components stemming from temporal and spatial variables, Err. Est. in this table does not indicate a clear pattern of convergence with respect to any discretization parameter. However, as the two discretizations are simultaneously refined, there is an observable reduction.

Tables 4.9 to 4.12 present the decomposition of error for Ex. 1 and Ex. 2 into various components as dictated by Theorem 3.2. The columns for $\mathcal{E}_0 = E_0$, Err. Est., and Eff. are not included since they are the same as in the corresponding columns in Tables 4.5 to 4.8, respectively. The component $\mathcal{E}_1$ seems to be the main contributor to the total error. For results associated with Ex. 1 (see Tables 4.9 and 4.10), the asymptotic behavior is roughly $\mathcal{E}_1 = \mathcal{O}(k)$. For Ex. 2 with $Q$ as stated in (4.9) (see Table 4.11), the asymptotic behavior is $\mathcal{E}_1 = \mathcal{O}(h^2)$. However, it is not quite the case for $Q$ in (4.10) (see Table 4.12).

5.   Conclusion and future work

This paper investigates the application of adjoint-based a posteriori error analysis for numerical approximation of Richards Equation. Construction of the approximate solution is cast into space-time variational formulation, specifically using the finite volume element in spatial variable and discontinuous Galerkin finite element in temporal variable. The resulting error estimators have the capability to predict components of error in certain quantities of interest that are expressed as a functional of the solution. The two examples give a strong indication that the error estimators are robust and capable to predict the error satisfactorily.

As for future work, we are interested in exploring further applications of the error estimators, in particular as to how they are applied to the setting of multidimensional problems. Owing to the various challenges persistent in the approximations of Richards Equation, a utilization of adaptivity is perhaps the only judicious route. Here the adaptivity is multi-faceted, not only as it pertains to local spatial refinement and dynamic time stepping, but also as it relates to determining optimal number of iterations when solving the nonlinear algebraic system. In this regard, the prospect of adjoint-based approach to estimate the components of error seems to be very promising.

Another interesting subject, which is not pursued in the present work, is a rigorous mathematical analysis of the proposed approximation. It must begin with establishing the existence of an approximate solution of (2.15). Here a potentially useful tool is either the Banach Fixed Point Theorem or the Brouwer Fixed Point Theorem. It should then be followed by a careful convergence analysis with the ultimate goal of showing the existence of a limit of the sequence of approximate solutions as $(h, k) \to (0, 0)$, and confirming that the limit satisfies a weak formulation of the Richards Equation. This can then be supplied with a study convergence rate of the approximate solution with respect to $h$ and $k$.