HIGH-ORDER MIMETIC DIFFERENCE SIMULATION OF UNSATURATED FLOW USING RICHARDS EQUATION

. The vadose zone is the portion of the subsurface above the water table and its pore space usually contains air and water. Due to the presence of inﬁltration, erosion, plant growth, microbiota, contaminant transport, aquifer recharge, and discharge to surface water, it is crucial to predict the transport rate of water and other substances within this zone. However, ﬂow in the vadose zone has many complications as the parameters that control it are extremely sensitive to the saturation of the media, leading to a nonlinear problem. This ﬂow is referred as unsaturated ﬂow and is governed by Richards equation. Analytical solutions for this equation exists only for simpliﬁed cases, so most practical situations require a numerical solution. Nevertheless, the nonlinear nature of Richards equation introduces challenges that causes numerical solutions for this problem to be computationally expensive and, in some cases, unreliable. High-order mimetic ﬁnite diﬀerence operators are discrete analogs of the continuous diﬀerential operators and have been extensively used in the ﬁelds of ﬂuid and solid mechanics. In this work, we present a numerical approach involving high-order mimetic operators along with a Newton root-ﬁnding algorithm for the treatment of the nonlinear component. Fully-implicit time discretization scheme is used to deal with the problem’s stiﬀness.


Introduction
In numerical analysis, classical numerical differentiation approaches begin by discretizing the corresponding system of partial differential equations of the specific problem to solve. This process leads to an approximation, at some order of accuracy, that may not result in a stable numerical scheme. In contrast, mimetic or compatible methods are formulated in such a way that allow them to discretely mimic and preserve the physical properties of the vector calculus operators used to describe continuous problems. These discrete operators are then employed to discretize the given problem. Typically, the mimetic differentiation scheme will achieve stability if the continuous problem is also stable.
Castillo-Grone's (CG) mimetic differential operators [4] are discrete approximations of their continuous counterparts and have been broadly used with success in multiple fields of physics and engineering. Some of the applications include: wave propagation, fluid dynamics, seismic studies, electromagnetism and image processing [2,4,9,10,11]. Compared to common numerical analysis techniques such as Finite Elements and Discontinuous Galerkin (DG), CG mimetic methods are computationally less expensive while achieving the same order of accuracy at the domain interior as well as at the boundary, which is something the Staggered Summation by Parts methods, for instance, cannot achieve. Recent developments have pushed the operators into solving for more realistic problems with complex coefficients [3] and non-trivial geometries [1].
Derived from CG operators, Castillo-Corbino's (CC) high-order mimetic difference operators provide the same features of the original operators with a more intuitive approach as they do not rely on free parameters which need to be selected for their construction. This results in discrete operators that have an optimal bandwidth and generally performed better, in terms of accuracy, when compared to CG counterparts [8]. The new operators are capable of achieving a uniform high-order of accuracy (up to eight order) in one-, two-, and three-dimensional space.
Unsaturated flow in porous media is a common physical phenomenon. One instance of this phenomenon is the flow of water in dry soils, which is of high importance in the fields of hydrology and agriculture. The Richards equation [6] is a nonlinear partial differential equation (PDE) that describes the flow of water in unsaturated soils. Most analytical solutions to Richards equation are impractical since they rely on simplifications and are limited to lower dimensions. There are multiple numerical approaches to tackle this problem, yet they all turn out to be computationally expensive and in certain circumstances, unreliable.
The purpose of this paper is to explain the use of a relatively new numerical approach (mimetic finite differences) to a well known and highly nonlinear problem (Richards equation). The work here presented describes and verifies the employment and accuracy of high-order mimetic approximations in conjunction with Newton's method to simulate flow in unsaturated porous media. The advantage of using our approach is that it allows for higher order approximations that are not only computationally cheaper but also intuitive to implement. This paper is organized as follows: In Section 2, we talk about the problem's governing equation and its parameters. In Section 3, we give a brief description of mimetic operators and staggered grids. Section 4 is about time discretization, and combining the operators presented in Section 3 with the equation shown in Section 2. In Section 5, we solve multiple test cases, and we compare against results obtained by [7]. Finally, in Section 6 we present our conclusions and proposed future work.

Richards equation
Richards equation describes the flow of water through unsaturated or partially saturated porous media. It is derived by combining Darcy's law [6] with the equation of mass conservation. There are two different forms of Richards equation that will be considered in this paper, they differ on how to deal with the non-linearity in the transient term.
• Head-based form: where C(ψ) = ∂θ ∂ψ is the specific moisture capacity function, K is the hydraulic conductivity, ψ is the pressure head, θ is the water content, and z is the depth. There are many empirical formulations for the hydraulic conductivity (K) and water content (θ) functions. The most popular model is the van Genuchten [6], which describes these functions without discontinuities. It is worth noting that both, the hydraulic conductivity and the water content depend on the pressure head (ψ). Indeed, both functions are highly nonlinear since they can dramatically change over a small range of ψ [6]. A version of the van Genuchten model that can be found at Celias paper [6] is: where K s corresponds to the saturated hydraulic conductivity, and θ s and θ r represent the saturated and residual moisture content, respectively. Richards equation typical use is to simulate infiltration experiments (in both lab and field scale). These experiments begin with a dry soil and then water is poured on top of the ground surface, making the connection with Darcy's law pretty obvious.

Mimetic operators
The gradient (∇) and divergence (∇·) are vector calculus operators that perform linear transformations (differentiation), and as such, they have matrix representation. These matrices can be constructed using standard finite differences [6] to approximate the derivative of the unknown function (in this case ψ). A special subset of these matrices that is closer to the continuum counterpart is known as mimetic. Mimetic finite-differences (MFD) have been widely and effectively used in many fields of applied mathematics and physics [5] and [8]. MFD provide high-order uniform accuracy without compromising physical coherence. Represented by sparse matrices, the main idea behind the construction of these operators is to find high-order approximations that satisfy the Extended Gauss Divergence theorem in the discrete sense [5]: where G ≡ ∇ and D ≡ ∇· are the mimetic gradient and divergence operators, respectively, B is a mimetic boundary operator, and the weight matrices P , Q and I are self-adjoint. In particular, the Q inner product accounts for the scalar inner product in cell centers, the P inner product accounts for a vector-field inner product at the cell faces, and the I inner product is at the boundaries. Notice that weighted inner products are defined in the standard form, Finally, as discrete counterparts, mimetic operators "mimic" the following vector calculus identities of their continuous analogs making them more faithful to the physics: with L ≡ ∇ 2 and C ≡ ∇× as the mimetic Laplacian and curl operators, respectively.
3.1. Staggered Grids. Mimetic operators are defined over staggered grids. Scalar variables such as density, pressure, and temperature are stored at the cells' centers, while vector variables (velocity components, electric conductivity, hydraulic conductivity, etc.) are stored at the edges (or faces in 3D domains).

3.2.
One-dimensional Mimetic Operators. Considering improvements with respect to the original Castillo-Grone mimetic Operators in terms of accuracy and optimal bandwidth, we follow the Castillo and Corbino [8] approach for the construction of our operators. We illustrate the second-order onedimensional mimetic divergence (D) and gradient (G) operators, which are the foundations of mimetic operators in higher dimensions and higher order. In our one-dimensional staggered grid discretization (depicted in Figure 1), the mimetic divergence operator acts on vector components (v-values) defined at m + 1 nodes, with x i = i∆x, i = 0, 1, ..., m.
These v-values are regarded as an (m + 1)-tuple. Conversely, the mimetic gradient operator acts on u-values defined at both boundary nodes (x 0 on the left, and x m on the right), as well as the m cellcenters, x i+ 1 /2 = (i + 1 /2)∆x, i = 0, 1, ..., m − 1. Therefore, u-values are regarded as (m + 2)-tuple. D is then an (m + 2) × (m + 1) sparse matrix with first and last row as zero vectors (required since the divergence is calculated at cell-centers). The gradient operator, G is a (m + 1) × (m + 2) matrix. The one-dimensional mimetic divergence operator is given by: and the mimetic gradient: Note there is a minimum number of cells needed to construct these mimetic operators. The gradient requires at least 2k cells, whereas the divergence requires at least 2k + 1, where k is the desired order of accuracy. Finally, these operators can also be extended to higher dimensions as seen in [8].

Compact
Operators. An important feature of mimetic operators is that they provide uniform order of accuracy all the way to the boundary. High-order (m ≥ 4) mimetic operators can be represented in a "compact way" by factoring the original matrices (i.e. the 2 nd -order matrix). By doing this, higher orders of accuracy can be attained using only the minimum number of points from the mesh. For instance, a k th -order mimetic gradient operator can be constructed as: where L kth represents the left factor matrix of k th -order. Conversely, for a k th -order mimetic divergence operator: Here R kth is the right factor matrix of k th -order. Finally, we write the k th -order mimetic Laplacian operator in compact form as: Our implementation uses compact representations of the fourth-order mimetic gradient and divergence for Equation (2.1).

Temporal Discretization
It is important to mention that mimetic operators are only for spatial discretization. Considering that we want to solve a time-dependent problem, a time discretization scheme has to be chosen. Both, Celia [6] and [7] use a first order forward scheme for time discretization: Since we want to replicate [6] and [7] results, we opted for Equation (4.1) as our time discretization scheme, and the mixed form of Richards Equation (2.2). Putting all together: where D and G are the mimetic divergence and gradient, respectively. The derivative in the z-direction is written D z , which in one-dimension is the same as the divergence matrix. Equation (4.2) is the equation we want to solve, and since it involves multiple ψ n+1 we must solve it using a root-finding method, we opted for Newton's method: where the superscript m indicates the Newton's iteration, α is the step length of the descent, and δψ (update) is obtained by solving the system of linear equations, where J ψ n+1,m is the Jacobian of the system, and F (ψ n+1,m , ψ n ) = R n+1,m is the Newton's residual (Equation (4.2)). Recall that:

Test case
In order to verify the accuracy of our scheme, we make a direct comparison with results obtained by Celia et al. [6] (they only implemented a Picard's iteration model of the head-based formulation), and Cockett in [7] (with both, Picard and Newton's iteration implementation of the mixed formulation). For this, we use a set of second-and fourth-order mimetic operators based on Castillo-Corbino's formulation, while employing a mixed formulation of the Richards equation (using a Newton's iteration model). The van Genuchten model (Equation (2.3)) for this experiment has the following parameters: α = 1.611 × 10 6 , θ s = 0.287, θ r = 0.075, β = 3.96, A = 1.175 × 10 6 , γ = 4.74, K s = 9.44 × 10 −5 m /s. For this experiment, the initial condition is a pressure head of ψ 0 (x, 0) = −61.5cm for a 40cm high one dimensional soil column (initially dry). Inhomogeneous Dirichlet boundary conditions are considered. Here, the bottom of the soil column is consistent with the initial condition: ψ(0cm, t) = −61.5cm. On the other hand, the top of the soil column is given by ψ(40cm, t) = −20.7cm. Naturally, this irregularity causes a boundary layer and a steep gradient in the pressure head at early times [7]. Although it can be expected that the Newton's iteration method will fail to converge at early times, experimentally, we found success with the method by using the initial condition as the starting guess at every time step. As in [6,7], the spatial grid resolution is a fix value of 1cm long (note that ours is a staggered grid), and we chose the same time steps to make a fair comparison. Figure 2 and Figure 3 depict the numerical solution of the Richards equation obtained by [6,7]. Furthermore, Figure 4 shows our numerical solution with high-order (k = 4) mimetic differences for the mentioned equation using a mixed formulation. Visually, our numerical approximation is congruent with the ones produced by [6] and [7]. We estimated the order of convergence of our implementation by comparing successively finer grids around the point x = 0.7cm, and using ∆t = 0.25s. Table 1 displays the numerical results for this experiment, based on this we can confirm a second-order convergence for our set of second-order mimetic operators while attaining better than O(h 3.5 ) for our set of fourth-order operators. These results corroborate better convergence for our operators when compared to the first order convergence obtained by [7], while correctly capturing the physics of the nonlinear problem.     Table 1. Numerical convergence for the mixed formulation of the Richards equation using our second-, and fourth-order mimetic operators. Where U k ≈ ψ with a k-order of accuracy.

Conclusions and future work
In this work we solved Richards equation in its mixed formulation by using second-and fourth-order mimetic versions of the gradient and divergence operators. The main reason why we study this equation is due to the highly nonlinear nature of the problem as well as the functions involved (hydraulic conductivity and water content) which are of great importance in mathematical biology.
The discretization scheme (first order in time, second-and fourth-order in space) was fully-implicit and involved the use of Newtons iteration to deal with the non-linearity. Experimentally, we have found that the best results are produced by having the initial condition as starting guess at each time step (it proved a faster convergence with respect to a random guess). Numerical results are coherent with those obtained in [6,7] via standard finite difference approaches. Furthermore, the implementation presented here not only allows closer approximations, but also a computationally cheaper and intuitive way to solve the same problem. Future work should include: • Extension of the model to higher dimensions.
• Use of a high-order scheme for time discretization.
• Include scenarios where the physical domain is irregular or non trivial (can be solved with overlapping grids).