A numerical method for efficient 3D inversions using Richards equation

TL;DR

This paper introduces an efficient numerical method for 3D inversion of hydraulic properties using Richards equation, enabling large-scale, deterministic parameter estimation with reduced computational resources.

Contribution

The paper presents a Jacobian-free, vector product-based inversion algorithm for Richards equation, scalable to large 3D problems and compatible with modest computational resources.

Findings

Successfully applied to 3D saturated hydraulic conductivity inversion

Allows larger problems than traditional Jacobian-based methods

Operates efficiently on modest computational hardware

Abstract

Fluid flow in the vadose zone is governed by Richards equation; it is parameterized by hydraulic conductivity, which is a nonlinear function of pressure head. Investigations in the vadose zone typically require characterizing distributed hydraulic properties. Saturation or pressure head data may include direct measurements made from boreholes. Increasingly, proxy measurements from hydrogeophysics are being used to supply more spatially and temporally dense data sets. Inferring hydraulic parameters from such datasets requires the ability to efficiently solve and deterministically optimize the nonlinear time domain Richards equation. This is particularly important as the number of parameters to be estimated in a vadose zone inversion continues to grow. In this paper, we describe an efficient technique to invert for distributed hydraulic properties in 1D, 2D, and 3D. Our algorithm does not store the Jacobian, but rather computes the product with a vector, which allows the size of the inversion problem to become much larger than methods such as finite difference or automatic differentiation; which are constrained by computation and memory, respectively. We show our algorithm in practice for a 3D inversion of saturated hydraulic conductivity using saturation data through time. The code to run our examples is open source and the algorithm presented allows this inversion process to run on modest computational resources.

Tables

Table 1. Table 1: Fictitious source test for Richards equation in 1D using the mixed-form Newton iteration.

Mesh Size (n)	${\Vert \Psi (x,0.5)-{\Psi }_{\text{true}}(x,0.5)\Vert }_{\infty }$	Order Decrease, $𝒪\left(\delta \right)$
64	5.485569e+00
128	2.952912e+00	0.894
256	1.556827e+00	0.924
512	8.035072e-01	0.954
1024	4.086729e-01	0.975
2048	2.060448e-01	0.988
4096	1.034566e-01	0.994
8192	5.184507e-02	0.997

Table 2. Table 2: Comparison of total solve time and iteration of numerical schemes ( Δ ​ t = 10 Δ 𝑡 10 \Delta t=10 s)

Method Name	Solve Time (s)	Total Iterations
Mixed-Picard	1.574	479
Mixed-Newton	1.381	112

Table 3. Table 3: Comparison of the memory necessary for storing the dense explicit sensitivity matrix compared to the peak memory used for the implicit sensitivity calculation excluding the matrix solve. The calculations are completed on a variety of mesh sizes for a single distributed parameter ( K s subscript 𝐾 𝑠 K_{s} ) as well as for five distributed van Genuchten model parameters ( K s , α , n , θ r subscript 𝐾 𝑠 𝛼 𝑛 subscript 𝜃 𝑟 K_{s},\alpha,n,\theta_{r} , and θ s subscript 𝜃 𝑠 \theta_{s} ). Values are reported in gigabytes (GB).

	Explicit Sensitivity		Implicit Sensitivity
Mesh Size	1 parameter	5 parameters	1 parameter	5 parameters
$32\times 32\times 32$	1.31	6.55	0.136	0.171
$64\times 64\times 64$	10.5	52.4	0.522	0.772
$128\times 128\times 128$	83.9	419	3.54	4.09

Equations

$$\frac{\partial \theta (\psi )}{\partial t}-\nabla \cdot k(\psi )\nabla \psi -\frac{\partial k(\psi )}{\partial z}=0\psi \in \Omega$$

$$\theta (\psi )=\{\begin{array}{rlr}{\theta }_r& +\frac{{\theta }_s-{\theta }_r}{(1+|\alpha \psi {|}^n{)}^m}& \psi <0\\ {\theta }_s& & \psi \ge 0\end{array}$$

$$k(\psi )=\{\begin{array}{rlr}{K}_s& {\theta }_e(\psi {)}^l(1-(1-{\theta }_e(\psi {)}^{-m}{)}^m{)}^2& \psi <0\\ K_s& & \psi \ge 0\end{array}$$

$${\theta }_e(\psi )=\frac{\theta (\psi )-{\theta }_r}{{\theta }_s-{\theta }_r},m=1-\frac{1}{n},n>1$$

$$\frac{\partial \theta (\psi )}{\partial t}-\nabla \cdot \vec{f}-\frac{\partial k(\psi )}{\partial z}=0$$

$$k(\psi {)}^{-1}\vec{f}=\nabla \psi$$

$$\frac{d\mathbf{\theta }(\mathbf{\psi })}{dt}-\mathbf{D}\text{ diag}\left({\mathbf{k}}_{Av}({\mathbf{\psi }}^{n+1})\right)\mathbf{G}\mathbf{\psi }-{\mathbf{G}}_z\left({\mathbf{k}}_{Av}({\mathbf{\psi }}^{n+1})\right)=0$$

$${\mathbf{k}}_{Av}(\mathbf{\psi })=\frac{1}{{\mathbf{A}}_v\frac{1}{\mathbf{k}(\mathbf{\psi })}}$$

$$F({\mathbf{\psi }}^{n+1},{\mathbf{\psi }}^n)=\frac{\mathbf{\theta }({\mathbf{\psi }}^{n+1})-\mathbf{\theta }({\mathbf{\psi }}^n)}{\Delta t}-\mathbf{D}\text{ diag}\left({\mathbf{k}}_{Av}({\mathbf{\psi }}^{n+1})\right)\mathbf{G}{\mathbf{\psi }}^{n+1}-{\mathbf{G}}_z\left({\mathbf{k}}_{Av}({\mathbf{\psi }}^{n+1})\right)=0$$

$$F({\mathbf{\psi }}^{n+1},{\mathbf{\psi }}^n)\approx F({\mathbf{\psi }}^{n,m},{\mathbf{\psi }}^n)+{\mathbf{J}}_{{\mathbf{\psi }}^{n+1,m}}\delta \mathbf{\psi }$$

$${\mathbf{J}}_{{\psi }^{n+1,m}}={\frac{\partial F(\mathbf{\psi },{\mathbf{\psi }}^n)}{\partial \mathbf{\psi }}|}_{{\mathbf{\psi }}^{n+1,m}}$$

$${\mathbf{J}}_{{\mathbf{\psi }}^{n+1,m}}=\frac{1}{\Delta t}\frac{d\mathbf{\theta }({\mathbf{\psi }}^{n+1,m})}{d{\mathbf{\psi }}^{n+1,m}}-\frac{d}{d{\mathbf{\psi }}^{n+1,m}}\left(\mathbf{D}\mathsf{diag}\left({\mathbf{k}}_{Av}({\mathbf{\psi }}^{n+1,m})\right)\mathbf{G}{\mathbf{\psi }}^{n+1,m}\right)-{\mathbf{G}}_z\frac{d{\mathbf{k}}_{Av}({\mathbf{\psi }}^{n+1,m})}{d{\mathbf{\psi }}^{n+1,m}}$$

$$\frac{d{\mathbf{k}}_{Av}(\mathbf{\psi })}{d\mathbf{\psi }}=\mathrm{diag}\left(({\mathbf{A}}_v{\mathbf{k}}^{-1}(\mathbf{\psi }){)}^{-2}\right){\mathbf{A}}_v\mathsf{diag}\left({\mathbf{k}}^{-2}(\mathbf{\psi })\right)\frac{d\mathbf{k}(\mathbf{\psi })}{d\mathbf{\psi }}$$

$$\frac{\partial }{\partial \mathbf{\psi }}\left(\mathbf{D}\text{ diag}\left({\mathbf{k}}_{Av}(\mathbf{\psi })\right)\mathbf{G}\mathbf{\psi }\right)=\mathbf{D}\text{ diag}\left({\mathbf{k}}_{Av}(\mathbf{\psi })\right)\mathbf{G}+\mathbf{D}\text{ diag}\left(\mathbf{G}\mathbf{\psi }\right)\frac{\partial {\mathbf{k}}_{Av}(\mathbf{\psi })}{\partial \mathbf{\psi }}$$

$${\mathbf{J}}_{{\psi }^{n+1,m}}\delta \mathbf{\psi }=-F({\mathbf{\psi }}^{n+1,m},{\mathbf{\psi }}^n)$$

$${\mathbf{\psi }}^{n+1,m+1}={\mathbf{\psi }}^{n+1,m}+\alpha \delta \psi$$

$$\Vert F({\mathbf{\psi }}^{n,m+1},{\mathbf{\psi }}^n)\Vert <\Vert F({\mathbf{\psi }}^{n,m},{\mathbf{\psi }}^n)\Vert$$

$$\hat{\mathbf{m}}=\arg \underset{\mathbf{m}}{\min }\Phi (\mathbf{m})=\frac{1}{2}{\Vert {\mathbf{W}}_d({\mathbf{d}}_{\mathrm{pred}}(\mathbf{m})-{\mathbf{d}}_{\mathrm{obs}})\Vert }_2^2+\frac{\beta }{2}{\Vert {\mathbf{W}}_m(\mathbf{m}-{\mathbf{m}}_{\mathrm{ref}})\Vert }_2^2.$$

$${\nabla }_{\mathbf{m}}\Phi (\mathbf{m})={\mathbf{J}}^{\top }{\mathbf{W}}_d^{\top }{\mathbf{W}}_d({\mathbf{d}}_{\mathrm{pred}}(\mathbf{m})-{\mathbf{d}}_{\mathrm{obs}})+\beta {\mathbf{W}}_m^{\top }{\mathbf{W}}_m(\mathbf{m}-{\mathbf{m}}_{\mathrm{ref}})=0$$

$${\nabla }_{\mathbf{m}}^2\Phi (\mathbf{m})\approx \mathbf{H}={\mathbf{J}}^{\top }{\mathbf{W}}_d^{\top }{\mathbf{W}}_d\mathbf{J}+\beta {\mathbf{W}}_m^{\top }{\mathbf{W}}_m=0$$

$$\mathbf{H}\delta \mathbf{m}=-{\nabla }_{\mathbf{m}}\Phi (\mathbf{m})$$

$${\mathbf{m}}_{k+1}={\mathbf{m}}_k+{\alpha }_k\delta \mathbf{m}$$

$$F({\mathbf{\psi }}^{n+1}(\mathbf{m}),{\mathbf{\psi }}^n(\mathbf{m}),\mathbf{m})=\frac{{\mathbf{\theta }}^{n+1}({\mathbf{\psi }}^{n+1})-{\mathbf{\theta }}^n({\mathbf{\psi }}^n)}{\Delta t}$$

$$-\mathbf{D}\text{ diag}\left({\mathbf{k}}_{Av}({\mathbf{\psi }}^{n+1},\mathbf{m})\right)\mathbf{G}{\mathbf{\psi }}^{n+1}-{\mathbf{G}}_z{\mathbf{k}}_{Av}({\mathbf{\psi }}^{n+1},\mathbf{m})=0$$

$${\nabla }_{\mathbf{m}}F({\mathbf{\psi }}^n,{\mathbf{\psi }}^{n+1},\mathbf{m})=\frac{\partial F}{\partial {\mathbf{k}}_{Av}}\frac{\partial {\mathbf{k}}_{Av}}{\partial \mathbf{m}}+\frac{\partial F}{\partial {\mathbf{\psi }}^n}\frac{\partial {\mathbf{\psi }}^n}{\partial \mathbf{m}}+\frac{\partial F}{\partial {\mathbf{\psi }}^{n+1}}\frac{\partial {\mathbf{\psi }}^{n+1}}{\partial \mathbf{m}}=0$$

$$\frac{1}{\Delta t}\left(\frac{\partial {\mathbf{\theta }}^{n+1}}{\partial {\mathbf{\psi }}^{n+1}}\frac{\partial {\mathbf{\psi }}^{n+1}}{\partial \mathbf{m}}-\frac{\partial {\mathbf{\theta }}^n}{\partial {\mathbf{\psi }}^n}\frac{\partial {\mathbf{\psi }}^n}{\partial \mathbf{m}}\right)-\mathbf{D}\text{ diag}\left(\mathbf{G}{\mathbf{\psi }}^{n+1}\right)\left(\frac{\partial {\mathbf{k}}_{Av}}{\partial \mathbf{m}}+\frac{\partial {\mathbf{k}}_{Av}}{\partial {\mathbf{\psi }}^{n+1}}\frac{\partial {\mathbf{\psi }}^{n+1}}{\partial \mathbf{m}}\right)$$

$$-\mathbf{D}\text{ diag}\left({\mathbf{k}}_{Av}({\mathbf{\psi }}^{n+1})\right)\mathbf{G}\frac{\partial {\mathbf{\psi }}^{n+1}}{\partial \mathbf{m}}-{\mathbf{G}}_z\left(\frac{\partial {\mathbf{k}}_{Av}}{\partial \mathbf{m}}+\frac{\partial {\mathbf{k}}_{Av}}{\partial {\mathbf{\psi }}^{n+1}}\frac{\partial {\mathbf{\psi }}^{n+1}}{\partial \mathbf{m}}\right)$$

$$\stackrel{{\mathbf{A}}_0({\mathbf{\psi }}^{n+1})}{\overbrace{\left[\frac{1}{\Delta t}\frac{\partial {\mathbf{\theta }}^{n+1}}{\partial {\mathbf{\psi }}^{n+1}}-\mathbf{D}\text{ diag}\left(\mathbf{G}{\mathbf{\psi }}^{n+1}\right)\frac{\partial {\mathbf{k}}_{Av}}{\partial {\mathbf{\psi }}^{n+1}}-\mathbf{D}\text{ diag}\left({\mathbf{k}}_{Av}({\mathbf{\psi }}^{n+1},\mathbf{m})\right)\mathbf{G}-{\mathbf{G}}_z\frac{\partial {\mathbf{k}}_{Av}}{\partial {\mathbf{\psi }}^{n+1}}\right]}}\frac{\partial {\mathbf{\psi }}^{n+1}}{\partial \mathbf{m}}$$

$$+\underset{{\mathbf{A}}_{-1}({\mathbf{\psi }}^n)}{\underbrace{\left[-\frac{1}{\Delta t}\frac{\partial {\mathbf{\theta }}^n}{\partial {\mathbf{\psi }}^n}\right]}}\frac{\partial {\mathbf{\psi }}^n}{\partial \mathbf{m}}=\underset{\mathbf{B}({\psi }^{n+1})}{\underbrace{\left[-\mathbf{D}\text{ diag}\left(\mathbf{G}{\mathbf{\psi }}^{n+1}\right)\frac{\partial {\mathbf{k}}_{Av}}{\partial \mathbf{m}}-{\mathbf{G}}_z\frac{\partial {\mathbf{k}}_{Av}}{\partial \mathbf{m}}\right]}}$$

$$\overbrace{\left[\begin{array}[]{cccc}\mathbf{A}_{0}(\boldsymbol{\psi}_{1})\\ \mathbf{A}_{-1}(\boldsymbol{\psi}_{1})&\mathbf{A}_{0}(\boldsymbol{\psi}_{2})\\ &\mathbf{A}_{-1}(\boldsymbol{\psi}_{2})&\mathbf{A}_{0}(\boldsymbol{\psi}_{3})\\ &\qquad\ddots&\qquad\ddots\\ &&\mathbf{A}_{-1}(\boldsymbol{\psi}_{n_{t}-1})&\mathbf{A}_{0}(\boldsymbol{\psi}_{n_{t}})\\ \end{array}\right]}^{\mathbf{A}(\boldsymbol{\Psi},\mathbf{m})}\overbrace{\left[\begin{array}[]{c}\frac{\partial\boldsymbol{\psi}_{1}}{\partial\mathbf{m}}\\ \frac{\partial\boldsymbol{\psi}_{2}}{\partial\mathbf{m}}\\ \vdots\\ \frac{\partial\boldsymbol{\psi}_{n_{t}-1}}{\partial\mathbf{m}}\\ \frac{\partial\boldsymbol{\psi}_{n_{t}}}{\partial\mathbf{m}}\\ \end{array}\right]}^{\frac{\partial\boldsymbol{\Psi}}{\partial\mathbf{m}}}=\overbrace{\left[\begin{array}[]{c}\mathbf{B_{1}}(\boldsymbol{\psi}_{1})\\ \mathbf{B_{2}}(\boldsymbol{\psi}_{2})\\ \vdots\\ \mathbf{B_{n-1}}(\boldsymbol{\psi}_{n_{t}-1})\\ \mathbf{B_{n}}(\boldsymbol{\psi}_{n_{t}})\\ \end{array}\right]}^{\mathbf{B}(\boldsymbol{\Psi},\mathbf{m})}$$

$$\mathbf{J}=\mathbf{PA}(\mathbf{\Psi },\mathbf{m}{)}^{-1}\mathbf{B}(\mathbf{\Psi },\mathbf{m})$$

$${\mathbf{J}}^{\top }=\mathbf{B}(\mathbf{\Psi },\mathbf{m}{)}^{\top }\mathbf{A}(\mathbf{\Psi },\mathbf{m}{)}^{-\top }{\mathbf{P}}^{\top }$$

$${\Psi }_{\text{true}}(x,t)=-20\arctan (20((x-0.25)-t))-40$$

Full text

††: Inverse Problems

Keywords: Richards equation, parameter estimation, hydrogeophysics, inversion, finite volume

1 Introduction

Studying the processes that occur in the vadose zone, the region between the earth’s surface and the fully saturated zone, is of critical importance for understanding our groundwater resources. Fluid flow in the vadose zone is described by Richards equation and parameterized by hydraulic conductivity, which is a nonlinear function of pressure head [1, 2]. Typically, hydraulic conductivity is heterogeneous and can have a large dynamic range. In any site characterization, the spatial estimation of the hydraulic conductivity function is an important step. Achieving this, however, requires the ability to efficiently solve and optimize the nonlinear, time-domain Richards equation. Rather than working with a full, implicit, 3D time-domain system of equations, simplifications are consistently used to avert the conceptual, practical, and computational difficulties inherent in the parameterization and inversion of Richards equation. These simplifications typically parameterize the conductivity and assume that it is a simple function in space, often adopting a homogeneous or one dimensional layered soil profile (cf. [3, 4, 5, 6]). Due to the lack of constraining hydrologic data, such assumptions are often satisfactory for fitting observed measurements, especially in two and three dimensions as well as in time. However, as more data become available, through spatially extensive surveys and time-lapse proxy measurements (e.g. direct current resistivity surveys and distributed temperature sensing), extracting more information about subsurface hydrogeologic parameters becomes a possibility. The proxy data can be directly incorporated through an empirical relation (e.g. [7]) or time-lapse estimations can be structurally incorporated through some sort of regularization technique [8, 9, 5]. Recent advances have been made for the forward simulation of Richards equation in a computationally-scalable manner [10]. However, the inverse problem is non-trivial, especially in three-dimensions [11], and must be considered using modern numerical techniques that allow for spatial estimation of hydraulic parameters.

Inverse problems in space and time are often referred to as history matching problems (see [12, 13, 14, 15, 16] and reference within). Inversions use a flow simulation model, combined with some a-priori information, in order to estimate a spatially variable hydraulic conductivity function that approximately yields the observed data. The literature shows a variety of approaches for this inverse problem, including trial-and-error, stochastic methods, and various gradient based methods [17, 3, 18, 19, 20, 21, 22]. The way in which the computational complexity of the inverse method scales becomes important as problem size increases [11]. Computational memory and time often become a bottleneck for solving the inverse problem, both when the problem is solved in 2D and, particularly, when it is solved in 3D [23]. To solve the inverse problem, stochastic methods are often employed, which have an advantage in that they can examine the full parameter space and give insights into non-uniqueness [24]. However, as the number of parameters we seek to recover in an inversion increases, these stochastic methods require that the forward problem be solved many times, which often makes these methods impractical. This scalability, especially in the context of hydrogeophysics has been explicitly noted in the literature (cf. [3, 4, 11, 25]).

Derivative-based optimization techniques become a practical alternative when the forward problem is computationally expensive or when there are many parameters to estimate (i.e. thousands to millions). Inverse problems are ill-posed and thus to pose a solvable optimization problem, an appropriate regularization is combined with a measure of the data misfit to state a deterministic optimization problem [26]. Alternatively, if prior information can be formulated using a statistical framework, we can use Bayesian techniques to obtain an estimator through the Maximum A Posteriori model (MAP) [27]. In the context of Bayesian estimation, gradient based methods are also important, as they can be used to efficiently sample the posterior [28].

A number of authors have sought solutions for the inverse problem, where the forward problem is Richards equation (cf. [29, 30, 16] and references within). Since the problem is parabolic (therefore, stiff), most work discretizes the forward problem by an implicit method in time and a finite volume or finite element in space. Most work uses a Newton-like method for the resulting nonlinear system, which arises from the discretization of the forward problem. For the deterministic inverse problem using Richards equation, previous work uses some version of a Gauss-Newton method (e.g. Levenberg-Marquardt), with a direct calculation of the sensitivity matrix [24, 22, 29]. However, while these approaches allow for inversions of moderate scale, they have one major drawback: the sensitivity matrix is large and dense; its computation requires dense linear algebra and a non-trivial amount of memory (cf. [11]). Previous work used either external numerical differentiation (e.g. PEST) or automatic differentiation in order to directly compute the sensitivity matrix [20, 29, 31, 11]. These techniques can generate inaccuracies in the sensitivity matrix and, consequently, tarry the convergence of the optimization algorithm. Furthermore, external numerical differentiation is computationally intensive and limits the number of model parameters that can be estimated.

The goal of this paper is to suggest a modern numerical formulation that allows the inverse problem to be solved without explicit computation of the sensitivity matrix by using exact derivatives of the discrete formulation [23]. Our technique is based on the discretize-then-optimize approach, which discretizes the forward problem first and then uses a deterministic optimization algorithm to solve the inverse problem [32]. To this end, we require the discretization of the forward problem. Similar to the work of [2], we use an implicit Euler method in time and finite volume in space. Given the discrete form, we show that we can analytically compute the derivatives of the forward problem with respect to distributed hydraulic parameters and, as a result, obtain an implicit formula for the sensitivity. The formula involves the solution of a linear time-dependent problem; we avoid computing and storing the sensitivity matrix directly and, rather, suggest a method to efficiently compute the product of the sensitivity matrix and its adjoint times a vector. Equipped with this formulation, we can use a standard inexact Gauss-Newton method to solve the inverse problem for distributed hydraulic parameters in 3D. This large-scale distributed parameter estimation becomes computationally tractable with the technique presented in this paper and can be employed with any iterative Gauss-Newton-like optimization technique.

This paper is structured as follows: in Section 2, we discuss the discretization of the forward problem on a staggered mesh in space and backward Euler in time; in Section 3, we formulate the inverse problem and construct the implicit functions used for computations of the Jacobian-vector product. In Section 4, we demonstrate the validity of the implementation of the forward problem and sensitivity calculation. Section 5 compares two iterative approaches, Picard and Newton, for handling the nonlinearity in Richards equation. Finally, in Section 6, we show an example of a 3D inversion for hydraulic conductivity and discuss extensions for inverting for multiple distributed hydraulic parameters from Richards equation and contrast the scalability of our methodology to standard techniques.

To accelerate both the development and dissemination of this approach, we have built these tools on top of an open source framework for organizing simulation and inverse problems in geophysics [33]. We have released our numerical implementation under the permissive MIT license. Our implementation of the implicit sensitivity calculation for Richards equation and associated inversion implementation is provided and tested to support 1D, 2D, and 3D forward and inverse simulations with respect to custom empirical relations and sensitivity to parameters within these functions. The source code can be found at https://github.com/simpeg/simpeg and may be a helpful resource for researchers looking to use or extend our implementation.

2 Forward problem

In this section, we describe Richards equations and its discretization [1]. Richards equation is a nonlinear parabolic partial differential equation (PDE) and we follow the so-called mixed formulation presented in [2] with some modifications. In the derivation of the discretization, we give special attention to the details used to efficiently calculate the effect of the sensitivity on a vector, which is needed in any derivative based optimization algorithm.

2.1 Richards equation

The parameters that control groundwater flow depend on the effective saturation of the media, which leads to a nonlinear problem. The groundwater flow equation has a diffusion term and an advection term which is related to gravity and only acts in the $z$-direction. There are two different forms of Richards equation; they differ in how they deal with the nonlinearity in the time-stepping term. Here, we use the most fundamental form, referred to as the ‘mixed’-form of Richards equation [2]:

[TABLE]

where $\psi$ is pressure head, $\theta(\psi)$ is volumetric water content, and $k(\psi)$ is hydraulic conductivity. This formulation of Richards equation is called the ‘mixed’-form because the equation is parameterized in $\psi$ but the time-stepping is in terms of $\theta$. The hydraulic conductivity, $k(\psi)$, is a heterogeneous and potentially anisotropic function that is assumed to be known when solving the forward problem. In this paper, we assume that $k$ is isotropic, but the extension to anisotropy is straightforward [33, 34]. The equation is solved in a domain, $\Omega$, equipped with boundary conditions on $\partial\Omega$ and initial conditions, which are problem-dependent.

An important aspect of unsaturated flow is noticing that both water content, $\theta$, and hydraulic conductivity, $k$, are functions of pressure head, $\psi$. There are many empirical relations used to relate these parameters, including the Brooks-Corey model [35] and the van Genuchten-Mualem model [21, 36]. The van Genuchten model is written as:

[TABLE]

where

[TABLE]

Here, $\theta_{r}$ and $\theta_{s}$ are the residual and saturated water contents, $K_{s}$ is the saturated hydraulic conductivity, $\alpha$ and $n$ are fitting parameters, and, $\theta_{e}(\psi)\in[0,1]$ is the effective saturation. The pore connectivity parameter, $l$, is often taken to be $\frac{1}{2}$, as determined by [21]. Figure 1 shows the functions over a range of negative pressure head values for four soil types (sand, loam, sandy clay, and clay). The pressure head varies over the domain $\psi\in(-\infty,0)$. When the value is close to zero (the left hand side), the soil behaves most like a saturated soil where $\theta=\theta_{s}$ and $k=K_{s}$. As the pressure head becomes more negative, the soil begins to dry, which the water retention curve shows as the function moving towards the residual water content ($\theta_{r}$). Small changes in pressure head can change the hydraulic conductivity by several orders of magnitude; as such, $k(\psi)$ is a highly nonlinear function, making Richards equation a nonlinear PDE.

2.2 Discretization

Richards equation is parameterized in terms of pressure head, $\psi$. Here, we describe simulating Richards equation in one, two, and three dimensions. We start by discretizing in space and then we discretize in time. This process yields a discrete, nonlinear system of equations; for its solution, we discuss a variation of Newton’s method.

2.2.1 Spatial Discretization

In order to conservatively discretize Richards equation, we introduce the flux ${\vec{f}}$ and rewrite the equation as a first order system of the form:

[TABLE]

We then discretize the system using a standard staggered finite volume discretization (cf. [37, 38, 34]). This discretization is a natural extension of mass-conservation in a volume where the balance of fluxes into and out of a volume are conserved [39]. Here, it is natural to assign the entire cell one hydraulic conductivity value, $k$, which is located at the cell center. Such assigning leads to a piecewise constant approximation for the hydraulic conductivity and allows for discontinuities between adjacent cells. From a geologic perspective, discontinuities are prevalent, as it is possible to have large differences in hydraulic properties between geologic layers in the ground. The pressure head, $\psi$, is also located at the cell centers and the fluxes are located on cell faces, which lead to the usual staggered mesh or Marker and Cell (MAC) discretization in space [40]. We demonstrate the discretization in 1D, 2D and 3D on the tensor mesh in Figure 2. We discretize the function, $\psi$, on a cell-centered grid, which results in a grid function, $\boldsymbol{\psi}$. We use bold letters to indicate other grid functions.

The discretization of a diffusion-like equation on an orthogonal mesh is well-known (see [41, 40, 42, 43] and reference within). We discretize the differential operators by using the usual mass balance consideration and the elimination of the flux, $\mathbf{f}$ 111Here we assume an isotropic conductivity that leads to a diagonal mass matrix and this yields easy elimination of the fluxes.. This spatial discretization leads to the following discrete nonlinear system of ordinary differential equations (assuming homogeneous Dirichlet boundary conditions):

[TABLE]

Here, $\mathbf{D}$ is the discrete divergence operator and $\mathbf{G}$ is the discrete gradient operator. The discrete derivative in the $z$-direction is written as $\mathbf{G}_{z}$. The values of $\psi$ and $k(\psi)$ are known on the cell-centers and must be averaged to the cell-faces, which we complete through harmonic averaging [41].

[TABLE]

where $\mathbf{A}_{v}$ is a matrix that averages from cell-centers to faces and the division of the vector is done pointwise; that is, we use the vector notation, $(1/{\mathbf{v}})_{i}=1/\mathbf{v}_{i}$. We incorporate boundary conditions using a ghost-point outside of the mesh [44].

2.2.2 Time discretization and stepping

Richards equation is often used to simulate water infiltrating an initially dry soil. At early times in an infiltration experiment, the pressure head, $\psi$, can be close to discontinuous. These large changes in $\psi$ are also reflected in the nonlinear terms $k(\psi)$ and $\theta(\psi)$; as such, the initial conditions imposed require that an appropriate time discretization be chosen. Hydrogeologists are often interested in the complete evolutionary process, until steady-state is achieved, which may take many time-steps. Here, we describe the implementation of a fully-implicit backward Euler numerical scheme. Higher-order implicit methods are not considered here because the uncertainty associated with boundary conditions and the fitting parameters in the Van Genuchten models (eq. 2) have much more effect than the order of the numerical method used.

The discretized approximation to the mixed-form of Richards equation, using fully-implicit backward Euler, reads:

[TABLE]

This is a nonlinear system of equations for $\boldsymbol{\psi}^{n+1}$ that needs to be solved numerically by some iterative process. Either a Picard iteration (as in [2]) or a Newton root-finding iteration with a step length control can be used to solve the system. Note that to deal with dependence of $\theta$ with respect to $\psi$ in Newton’s method, we require the computation of $\frac{d\boldsymbol{\theta}}{d\boldsymbol{\psi}}$. We can complete this computation by using the analytic form of the hydraulic conductivity and water content functions (e.g. derivatives of eq. 2). We note that a similar approach can be used for any smooth curve, even when the connection between $\theta$ and $\psi$ are determined empirically (for example, when $\theta(\psi)$ is given by a spline interpolation of field data).

2.3 Solving the nonlinear equations

Regardless of the empirical relation chosen, we must solve 7 using an iterative root-finding technique. Newton’s method iterates over $m=1,2,\dots$ until a satisfactory estimation of $\boldsymbol{\psi}^{n+1}$ is obtained. Given $\boldsymbol{\psi}^{n+1,m}$, we approximate $F(\boldsymbol{\psi}^{n+1},\boldsymbol{\psi}^{n})$ as:

[TABLE]

where the Jacobian for iteration, $m$, is:

[TABLE]

The Jacobian is a large dense matrix, and its computation necessitates the computation of the derivatives of $F(\boldsymbol{\psi}^{n+1,m},\boldsymbol{\psi}^{n})$. We can use numerical differentiation in order to evaluate the Jacobian (or its product with a vector). However, in the context of the inverse problem, an exact expression is preferred. Given the discrete forward problem, we obtain that:

[TABLE]

Here, recall that $\mathbf{k}_{Av}$ is harmonically averaged and its derivative can be obtained by the chain rule:

[TABLE]

Similarly, for the second term in (10) we obtain:

[TABLE]

Here the notation n+1,m has been dropped for brevity. For the computations above, we need the derivatives of functions $\mathbf{k}(\boldsymbol{\psi})$ and $\boldsymbol{\theta}(\boldsymbol{\psi})$; note that, since the relations are assumed local (point wise in space) given the vector, $\boldsymbol{\psi}$, these derivatives are diagonal matrices. For Newton’s method, we solve the linear system:

[TABLE]

For small-scale problems, we can solve the linear system using direct methods; however, for large-scale problems, iterative methods are more commonly used. The existence of an advection term in the PDE results in a non-symmetric linear system. Thus, when using iterative techniques to solve this system, an appropriate iterative method, such as bicgstab or gmres [45], must be used.

At this point, it is interesting to note the difference between the Newton iteration and the Picard iteration suggested in [2]. We can verify that the Picard iteration uses an approximation to the Jacobian $\mathbf{J}_{\psi^{n+1,m}}\,\delta\boldsymbol{\psi}$ given by dropping the second term from (12). This term can have negative eigenvalues and dropping it is typically done when considering the lagged diffusivity method [46]. However, as discussed in [46], ignoring this term can slow convergence.

Finally, a new iterate is computed by adding the Newton update to the last iterate:

[TABLE]

where $\alpha$ is a parameter that guarantees that

[TABLE]

To obtain $\alpha$, we perform an Armijo line search [47]. In our numerical experiments, we have found that this method can fail when the hydraulic conductivity is strongly discontinuous and changes rapidly. In such cases, Newton’s method yields a poor descent direction. Therefore, if the Newton iteration fails to converge to a solution, the update is performed with the mixed-form Picard iteration. Note that Picard iteration can be used, even when Newton’s method fails, because Picard iteration always yields a descent direction [46].

At this point, we have discretized Richards equation in both time and space while devoting special attention to the derivatives necessary in Newton’s method and the Picard iteration as described in [2]. The exact derivatives of the discrete problem will be used in the following two sections, which outline the implicit formula for the sensitivity and its incorporation into a general inversion algorithm.

3 Inverse Problem

The location and spatial variability of, for example, an infiltration front over time is inherently dependent on the hydraulic properties of the soil column. As such, direct or proxy measurements of the water content or pressure head at various depths along a soil profile contain information about the soil properties. We pose the inverse problem, which is the estimation of distributed hydraulic parameters, given either water content or pressure data. We frame this problem under the assumption that we wish to estimate hundreds of thousands to millions of distributed model parameters. Due to the large number of model parameters that we aim to estimate in this inverse problem, Bayesian techniques or external numerical differentiation, such as the popular PEST toolbox [31], are not computationally feasible. Instead, we will employ a direct method by calculating the exact derivatives of the discrete Richards equation and solving the sensitivity implicitly. For brevity, we show the derivation of the sensitivity for an inversion model of only saturated hydraulic conductivity, $K_{s}$, from pressure head data, $\mathbf{d}_{\rm obs}$. This derivation can be readily extended to include the use of water content data and inverted for other distributed parameters in the heterogeneous hydraulic conductivity function. We will demonstrate the sensitivity calculation for multiple distributed parameters in the numerical examples (Section 6).

The Richards equation simulation produces a pressure head field at all points in space as well as through time. Data can be predicted, $\mathbf{d}_{\rm pred}$, from these fields and compared to observed data, $\mathbf{d}_{\rm obs}$. To be more specific, we let $\boldsymbol{\Psi}=[(\boldsymbol{\psi}^{1})^{\top},\ldots,(\boldsymbol{\psi}^{n_{t}})^{\top}]^{\top}$ be the (discrete) pressure field for all space and $n_{t}$ time steps. When measuring pressure head recorded only in specific locations and times, we define the predicted data, $\mathbf{d}_{\rm pred}$, as $\mathbf{d}_{\rm pred}=\mathbf{P}\boldsymbol{\Psi}(\mathbf{m})$. Here, the vector $\mathbf{m}$ is the vector containing all of the parameters which we are inverting for (e.g. $K_{s},\alpha,n,\theta_{r}$, or $\theta_{s}$ when using the van Genuchten empirical relation). The matrix, $\mathbf{P}$, interpolates the pressure head field, $\boldsymbol{\Psi}$, to the locations and times of the measurements. Since we are using a simple finite volume approach and backward Euler in time, we use linear interpolation in both space and time to compute $\mathbf{d}_{\rm pred}$ from $\boldsymbol{\Psi}$. Thus, the entries of the matrix $\mathbf{P}$ contain the interpolation weights. For linear interpolation in 3D, $\mathbf{P}$ is a sparse matrix that contains up to eight non-zero entries in each row. Note that the time and location of the data measurement is independent and decoupled from the numerical discretization used in the forward problem. A water retention curve, such as the van Genuchten model, can be used for computation of predicted water content data, which requires another nonlinear transformation, $\mathbf{d}_{\rm pred}^{\theta}=\mathbf{P}\boldsymbol{\theta}(\boldsymbol{\Psi}(\mathbf{m}),\mathbf{m})$. Note here that the transformation to water content data, in general, depends on the model to be estimated in the inversion, which will be addressed in the numerical examples. For brevity in the derivation that follows, we will make two simplifying assumptions: (1) that the data are pressure head measurements, which requires a linear interpolation that is not dependent on the model; and, (2) that the model vector, $\mathbf{m}$, describes only distributed saturated hydraulic conductivity. Our software implementation does not make these assumptions; our numerical examples will use water content data, a variety of empirical relations, and calculate the sensitivity to multiple heterogeneous empirical parameters.

3.1 Solution through optimization

We can now formulate the discrete inverse problem to estimate saturated hydraulic conductivity, $\mathbf{m}$, from the observed pressure head data, $\mathbf{d}_{\rm obs}$. We frame the inversion as an optimization problem, which minimizes a data misfit and a regularization term.

[TABLE]

The first term in the objective function is the data misfit; it contains a weighted difference between the observed and predicted data. Assuming that the observed data is noisy, with independent distributed Gaussian noise having standard deviation $\boldsymbol{\sigma}$, the weighting matrix, $\mathbf{W}_{d}$, is a diagonal matrix that contains the entries $\boldsymbol{\sigma}_{i}^{-1}$ on its diagonal. The matrix, $\mathbf{W}_{m}$, is a regularization matrix that is chosen based on a-priori information and assumptions about the geologic setting [48]. In a regularized inversion framework, if we estimate only the log of the static conductivity, $K_{s}$, in (2), then we typically choose $\mathbf{W}_{m}$ as a combination of the first order gradient matrix and an identity matrix. In a Bayesian framework, where one assumes that $\mathbf{m}$ is Gaussian, the matrix $\mathbf{W}_{m}^{\top}\mathbf{W}_{m}$ is the inverse covariance matrix. If $\mathbf{m}$ contains more than a single distributed hydrogeologic parameter, then a correlation term is typically added [8]. The reference model, $\mathbf{m}_{\rm ref}$, is chosen based on a-priori information about the background, $\mathbf{m}$. Finally, we add the misfit and regularization terms with a regularization parameter, $\beta$, which balances the a-priori information, or smoothness, with the fit to noisy data. This is the standard approach in geophysical inversions [26, 48, 49, 38] where hundreds of thousands to millions of distributed parameters are commonly estimated in a deterministic inversion. The hydrogeologic literature also shows the use of these techniques; however, there is also a large community advancing stochastic inversion techniques and geologic realism (cf. [50]). Regardless of the inversion algorithm used, an efficient method to calculate the sensitivity is crucial; this method is the focus of our work.

To minimize the objective function, we use an inexact Gauss-Newton method that is well-suited for large-scale problems [47]. For the Gauss-Newton method, the gradient of Equation 14 is computed and iteratively driven towards [math] by using an approximate Hessian, $\mathbf{H}$. For the problem at hand, the gradient can be expressed as:

[TABLE]

where the sensitivity matrix, $\mathbf{J}=\nabla_{\mathbf{m}}\mathbf{d}_{\rm pred}(\mathbf{m})$, is the derivative of the predicted data with respect to the model parameters. The Hessian of the objective function, $\mathbf{H}$, is approximated with the first order terms and is guaranteed to be positive definite.

[TABLE]

Finally, the (inexact) Gauss-Newton update, $\delta\mathbf{m}$, is computed by (approximately) solving the system

[TABLE]

using some iterative technique. In this work, we use the preconditioned conjugate gradient (PCG) method, allowing us to work with large-scale problems. To precondition the conjugate gradient algorithm, we use a standard limited-memory BFGS algorithm initiated with the inverse of $\mathbf{W}_{m}^{\top}\mathbf{W}_{m}$ instead of the identity. For more details on the preconditioner, see [51].

It is important to note that the sensitivity matrix, $\mathbf{J}$, as well as the approximate Hessian, $\mathbf{H}$, are large, dense matrices and their explicit computations result in an algorithm that is constrained by computational memory. However, as we show in the following section, one can avoid this computation and efficiently compute the products of $\mathbf{J}$ and $\mathbf{J}^{\top}$ with a vector. This computation, coupled with an iterative optimization method, results in an efficient algorithm.

Once the model update, $\delta\bf m$, is found, a Gauss-Newton step is taken:

[TABLE]

where $\alpha_{k}$ is the current line search parameter. We use a backtracking Armijo line search [52] to enforce sufficient decrease in our objective function [47]. We repeat the minimization procedure, which updates the linearization until the norm of the gradient falls below a certain tolerance, the model update is sufficiently small, or until we have exhausted the maximum number of iterations. Values for these stopping criteria are often problem-dependent. Running synthetic models with elements similar to observed data can often give insight into appropriate stopping criteria [53].

3.2 Implicit sensitivity calculation

The optimization problem requires the derivative of the pressure head with respect to the model parameters, $\frac{\partial\boldsymbol{\Psi}}{\partial\mathbf{m}}$. We can obtain an approximation of the sensitivity matrix through a finite difference method on the forward problem [22, 24, 20]. One forward problem, or two, when using central differences, must be completed for each column in the Jacobian at every iteration of the optimization algorithm. This style of differentiation proves advantageous in that it can be applied to any forward problem; however, it is highly inefficient and introduces errors into the inversion that may slow the convergence of the scheme [31]. Automatic differentiation (AD) can also be used [47]. However, AD does not take the structure of the problem into consideration and requires that the dense Jacobian be explicitly formed. [29] presents three algorithms (finite difference, adjoint, and direct) to directly compute the elements of the dense sensitivity matrix for Richards equation. As problem size increases, the memory required to store this dense matrix often becomes a practical computational limitation [54, 11]. As we show next, it is possible to explicitly write the derivatives of the Jacobian and evaluate their products with vectors using only sparse matrix operations. This algorithm is much more efficient than finite differencing or AD, especially for large scale simulations, since it does not require explicitly forming and storing a large dense matrix. Rather, the algorithm efficiently computes matrix-vector and adjoint matrix-vector products with sensitivity. We can use these products for the solution of the Gauss-Newton system when using the conjugate gradient method, which bypasses the need for the direct calculation of the sensitivity matrix and makes solving large-scale inverse problems possible. Other geophysical inverse problems have used this idea extensively, especially in large-scale electromagnetics (cf. [23]). The challenge in both the derivation and implementation for Richards equation lies in differentiating the nonlinear time-dependent forward simulation with respect to multiple distributed hydraulic parameters.

The approach to implicitly constructing the derivative of Richards equation in time involves writing the whole time-stepping process as a block bi-diagonal matrix system. The discrete Richards equation can be written as a function of the model. For a single time-step, the equation is written:

[TABLE]

In this case, $\mathbf{m}$ is a vector that contains all the parameters of interest. Note that $\boldsymbol{\psi}^{n+1}$ and $\boldsymbol{\psi}^{n}$ are also functions of $\mathbf{m}$. In general, $\boldsymbol{\theta}^{n+1}$ and $\boldsymbol{\theta}^{n}$ are also dependent on the model; however, for brevity, we will omit these derivatives. The derivatives of $F$ to the change in the parameters $\mathbf{m}$ can be written as:

[TABLE]

or, in more detail:

[TABLE]

The above equation is a linear system of equations and, to solve for $\frac{d\boldsymbol{\Psi}}{d\mathbf{m}}$, we rearrange the block-matrix equation:

[TABLE]

Here, we use the subscript notation of $\mathbf{A}_{0}(\boldsymbol{\psi}^{n+1})$ and $\mathbf{A}_{-1}(\boldsymbol{\psi}^{n})$ to represent two block-diagonals of the large sparse matrix $\mathbf{A}({\boldsymbol{\Psi}},\mathbf{m})$. Note that all of the terms in these matrices are already evaluated when computing the Jacobian of Richards equations in Section 2 and that they contain only basic sparse linear algebra manipulations without the inversion of any matrix. If $\boldsymbol{\psi}_{0}$ does not depend on the model, meaning the initial conditions are independent, then we can formulate the block system as:

[TABLE]

This is a block matrix equation; both the storage and solve will be expensive if it is explicitly computed. Indeed, its direct computation is equivalent to the adjoint method [29, 12].

Since only matrix vector products are needed for the inexact Gauss-Newton optimization method, the matrix $\mathbf{J}$ is never needed explicitly and only the products of the form $\mathbf{J}\mathbf{v}$ and $\mathbf{J}^{\top}\mathbf{z}$ are needed for arbitrary vectors $\mathbf{v}$ and $\mathbf{z}$. Projecting the full sensitivity matrix onto the data-space using $\mathbf{P}$ results in the following equations for the Jacobian:

[TABLE]

In these equations, we are careful to not write $\frac{d\boldsymbol{\Psi}}{d\mathbf{m}}$, as it is a large dense matrix which we do not want to explicitly compute or store. Additionally, the matrices $\bf A(\boldsymbol{\Psi},\mathbf{m})$ and $\bf B(\boldsymbol{\Psi},\mathbf{m})$ do not even need to be explicitly formed because the matrix $\bf A(\boldsymbol{\Psi},\mathbf{m})$ is a triangular block-system, which we can solve using forward or backward substitution with only one block-row being solved at a time (this is equivalent to a single time step). To compute the matrix vector product, $\mathbf{J}\mathbf{v}$, we use a simple algorithm:

Given the vector $\mathbf{v}$ calculate $\mathbf{y}=\mathbf{Bv}$ 2. 2.

Solve the linear system $\mathbf{Aw}=\mathbf{y}$ for the vector $\mathbf{w}$ 3. 3.

Set $\mathbf{J}\mathbf{v}=\mathbf{P}\mathbf{w}$

Here, we note that we complete steps (1) and (2) using a for-loop with only one block-row being computed and stored at a time. As such, only the full solution, $\boldsymbol{\Psi}$, is stored and all other block-entries may be computed as needed. There is a complication here if data is in terms of water content or effective saturation, as the data projection is no longer linear and may have model dependence. These complications can be dealt with using the chain rule during step (3). Similarly, to compute the adjoint $\mathbf{J}^{\top}\mathbf{z}$ involves the intermediate solve for $\mathbf{y}$ in $\mathbf{A}^{\top}\mathbf{y}=\mathbf{P}^{\top}\mathbf{z}$ and then computation of $\mathbf{B}^{\top}\mathbf{y}$. Again, we solve the block-matrix via backward substitution with all block matrix entries being computed as needed. Note that the backward substitution algorithm can be viewed as time stepping, which means that it moves from the final time back to the initial time. This time stepping is equivalent to the adjoint method that is discussed in [12] and references within. The main difference between our approach and the classical adjoint-based method is that our approach yields the exact gradient of the discrete system; no such guarantee is given for adjoint-based methods.

The above algorithm and the computations of all of the different derivatives summarizes the technical details of the computations of the sensitivities. Equipped with this “machinery”, we now demonstrate that validity of our implementation and our ability to solve large-scale problems using modest computational resources.

4 Validation

4.1 Forward problem

Richards equation has no analytic solution, which makes testing the code more involved. Code-to-code comparisons have been completed for comparison to [2], which can be found in [55]. Here we have chosen to use a fictitious source experiment to rigorously test the code. In this experiment, we approximate an infiltration front by an arctangent function in one dimension, which is centered over the highly nonlinear part of the van Genuchten curves, with $\psi\in[-60,-20]$ centimeters. The arctangent curve advects into the soil column with time. The advantage of using an analytic function is that the derivative is known explicitly and can be calculated at all times. However, it should be noted that Richards equation does not satisfy the analytic solution exactly, but differs by a function, $S(x,t)$. We refer to this function as the fictitious source term. The analytic function that we used has similar boundary conditions and shape to an example in [2] and is considered over the domain $x\in[0,1]$.

[TABLE]

This analytic function is shown at times 0 and 0.5 in Figure 3 and has a pressure head range of $\psi\in[-60,-20]$. We can compare these values to the van Genuchten curves in Figure 1. We can then put the known pressure head into Richards equation (1) and calculate the analytic derivatives and equate them to a source term, $S(x,t)$. Knowing this source term and the analytic boundary conditions, we can solve discretized form of Richards equation, which should reproduce the analytic function in Equation 23. Table 1 shows the results of the fictitious source test when the number of mesh-cells is doubled and the time-discretization is both fixed and equivalent to the mesh size (i.e. $k=h$). In this case, we expect that the backward-Euler time discretization, which is $\mathcal{O}(\delta)$, will limit the order of accuracy. The final column of Table 1 indeed shows that the order of accuracy is $\mathcal{O}(\delta)$. The higher errors in the coarse discretization are due to high discontinuities and changes in the source term, which the coarse discretization does not resolve.

4.2 Inverse problem

In order to test the implicit sensitivity calculation, we employ derivative and adjoint tests as described in [38]. Given that the Taylor expansion of a function $f(\mathbf{m}+h\Delta\mathbf{m})$ is

[TABLE]

for any of the model parameters considered, we see that our approximation of $f(\mathbf{m}+h\Delta\mathbf{m})$ by $f(\mathbf{m})+h\mathbf{J}\Delta\mathbf{m}$ should converge as $\mathcal{O}(h^{2})$ as $h$ is reduced. This allows us to verify our calculation of $\mathbf{J}\mathbf{v}$. To verify the adjoint, $\mathbf{J}^{\top}\mathbf{v}$, we check that

[TABLE]

for any two random vectors, $\mathbf{w}$ and $\mathbf{v}$. These tests are run for all of the parameters considered in an inversion of Richards equation. Within our implementation, both the derivative and adjoint tests are included as unit tests which are run on any updates to the implementation (https://travis-ci.org/simpeg/simpeg).

5 Comparison of iterative methods

With the goal of extending these simple one-dimensional examples into three dimensions, the efficiency of the iterative method used is of critical importance for computational solve time. Figure 4 shows a comparison of the iterative methods used for the numerical example described above (where $\Delta t=10$s) . Here, both of the Picard schemes take between 10 to 25 iterations per time-step, while the Newton scheme generally uses less than five iterations (Figure 4a). The difference between iteration counts occurs due to the use of the gradient information in the Newton method, which speeds convergence when the initial guess is close to the solution. When the initial guess is not close to the solution, for example at $t=0$s, the Newton method is slower to converge than the Picard method. The first iteration of Newton has both a higher number of iterations and solve time (Figure 4b). In both Figure 4a & 4b, there is an overall decrease in the number of inner iterations and the time taken per iteration as the system reaches steady-state. The iteration count has previously been used to inform an adaptive time-stepping scheme [10]. Table 2 shows the overall results of solve time and total iterations. The Newton method is marginally faster than the Picard method because it takes the fewest number of iterations.

6 Three dimensional inversion

In this section, we turn our attention to recovering a three-dimensional soil structure given water content data. The example, motivated by a field experiment introduced in [56], shows a time-lapse electrical resistivity tomography survey completed in the base of a managed aquifer recharge pond. The goal of this management practice is to infiltrate water into the subsurface for storage and subsequent recovery. Such projects require input from geology, hydrology, and geophysics to map the hydrostratigraphy, to collect and interpret time-lapse geophysical measurements, and to integrate all results to make predictions and decisions about fluid movement at the site. As such, the hydraulic properties of the aquifer are important to characterize, and information from hydrogeophysical investigations has been demonstrated to inform management practices [56]. We use this context to motivate both the model domain setup of the following synthetic experiment and the subsequent inversion.

For the inverse problem solved here, we assume that time-lapse water-content information is available at many locations in the subsurface. In reality, water content information may only be available through proxy techniques, such as electrical resistivity methods. These proxy data can be related to hydrogeologic parameters using inversion techniques based solely on the geophysical inputs (cf. [57]). For the following numerical experiments, we do not address complications in empirical transformations, such as Archie’s equation [7]. The synthetic numerical model has a domain with dimensions 2.0 m $\times$ 2.0 m $\times$ 1.7 m for the $x$, $y$, and $z$ dimensions, respectively. The finest discretization used is 4 cm in each direction. We use padding cells to extend the domain of the model (to reduce the effect of boundary conditions in the modelling results). These padding cells extend at a factor of 1.1 in the negative $z$ direction. We use an exponentially expanding time discretization with 40 time steps and a total time of 12.3 hours. This choice in discretization leads to a mesh with 1.125$\times 10^{5}$ cells in space ($50\times 50\times 45$). To create a three-dimensionally varying soil structure, we construct a model for this domain using a three-dimensional, uniformly random field, $\in[0,2]$, that is convolved with an anisotropic smoothing kernel for a number of iterations. We create a binary distribution from this random field by splitting the values above and below unity. Figure 5 shows the resulting model, which reveals potential flow paths. We then map van Genuchten parameters to this synthetic model as either a sand or a loamy-sand. The van Genuchten parameters for sand are: $K_{s}$: 5.83e-05m/s, $\alpha$: 13.8, $\theta_{r}$: 0.02, $\theta_{s}$: 0.417, and $n$: 1.592; and for loamy-sand are: $K_{s}$: 1.69e-05m/s, $\alpha$: 11.5, $\theta_{r}$: 0.035, $\theta_{s}$: 0.401, and $n$: 1.474. For this inversion, we are interested in characterizing the soil in three dimensions.

For calculation of synthetic data, the initial conditions are a dry soil with a homogeneous pressure head ($\psi=-30$cm). The boundary conditions applied simulate an infiltration front applied at the top of the model, $\psi=-10$cm $\in\delta\Omega^{\text{top}}$. Neumann (no-flux) boundary conditions are used on the sides of the model domain. Figure 8 shows the pressure head and water content fields from the forward simulation. Figure 6a and 6b show two cross-sections at time 5.2 hours and 10.3 hours of the pressure head field. These figures show true soil type model as an outline, where the inclusions are the less hydraulically-conductive loamy sand. The pressure head field is continuous across soil type boundaries and shows the infiltration moving vertically down in the soil column. We can compute the water content field from the pressure head field using the nonlinear van Genuchten model chosen; Figure 6c and 6d show this computation at the same times. The loamy sand has a higher relative water content for the same pressure head and the water content field is discontinuous across soil type boundaries, as expected.

The observed data, which will be used for the inversion, is collected from the water content field at the points indicated in both Figure 6c and 6d. The sampling location and density of this three-dimensional grid within the model domain is similar to the resolution of a 3D electrical resistivity survey. Our implementation supports data as either water content or pressure head; however, proxy water content data is more realistic in this context. Similar to the field example in [56], we collect water content data every 18 minutes over the entire simulation, leading to a total of 5000 spatially and temporally extensive measurements. The observed water content data for a single infiltration curve is plotted through depth in Figure 7. The green circles in Figure 6 show the locations of these water content measurements. The depth of the observation is colour-coded by depth, with the shallow measurements being first to increase in water content over the course of the infiltration experiment. To create the observed dataset, $\mathbf{d}_{\rm obs}$, from the synthetic water content field, 1% Gaussian noise is added to the true water content field. This noise is below what can currently be expected from a proxy geophysical measurement of the water content. However, with the addition of more noise, we must reduce our expectations of our ability to recover the true parameter distributions from the data. In this experiment, we are interested in examining what is possible to recover under the best of circumstances, and therefore have selected a low noise level.

For the inverse problem, we are interested in the distribution of soil types that fits the measured data. We parameterize these soil types using the van Genuchten empirical model (2) with as least five spatially distributed properties. Inverting for all 5.625$\times 10^{5}$ parameters in this simulation with only 5000 data points is a highly underdetermined problem, and thus there are many possible models that may fit those data. In this 3D example, we will invert solely for saturated hydraulic conductivity and assume that all other van Genuchten parameters are equivalent to the sand; that is, they are parameterized to the incorrect values in the loamy sand. Note that this assumption, while reasonable in practice, will handicap the results, as the van Genuchten curves between these two soils differ. Better results can, of course, be obtained if we assume that the van Genuchten parameters are known; this assumption is unrealistic in practice, which means that we will not be able to recreate the data exactly. However, the distribution of saturated hydraulic conductivity may lead to insights about soil distributions in the subsurface. Figure 8 shows the results of the inversion for saturated hydraulic conductivity as a map view slice at 66 cm depth and two vertical sections through the center of the model domain. The recovered model shows good correlation to the true distribution of hydraulic properties, which is superimposed as a dashed outline. Figure 7 shows the predicted data overlaid on the true data for five water content measurement points through time; these data are from the center of the model domain and are colour-coded by depth. As seen in Figure 7, we do a good job of fitting the majority of the data. However, there is a tendency for the predicted infiltration front to arrive before the observed data, which is especially noticeable at deeper sampling locations. The assumptions put on all other van Genuchten parameters to act as sand, rather than loamy sand, lead to this result.

6.1 Scalability of the implicit sensitivity

For the forward simulation presented, the Newton root-finding algorithm took 4-12 iterations to converge to a tolerance of $1\times 10^{-4}$m on the pressure head. The inverse problem took 20 iterations of inexact Gauss-Newton with five internal CG iterations used at each iteration. The inversion fell below the target misfit of 5000 at iteration 20 with $\phi_{d}=4.893\times 10^{3}$; this led to a total of 222 calls to functions to solve the products $\mathbf{J}\mathbf{v}$ and $\mathbf{J}^{\top}\mathbf{z}$. Here, we again note that the Jacobian is neither computed nor stored directly, which makes it possible to run this code on modest computational resources; this is not possible if numerical differentiation or direct computation of the Jacobian is used. For these experiments, we used a single Linux Debian Node on Google Compute Engine (Intel Sandy Bridge, 16 vCPU, 14.4 GB memory) to run the simulations and inversion. The forward problem takes approximately 40 minutes to solve. In this simulation, the dense Jacobian matrix would have 562.5 million elements. If we used a finite difference algorithm to explicitly calculate each of the 1.125e5 columns of the Jacobian with a simple forward difference, we would require a calculation for each model parameter – or approximately 8.5 years of computational time. Furthermore, we would need to recompute the Jacobian at each iteration of the optimization algorithm. In contrast, if we use the implicit sensitivity algorithm presented in this paper, we can solve the entire inverse problem in 34.5 hours.

Table 3 shows the memory required to store the explicit sensitivity matrix for a number of mesh sizes and contrasts them to the memory required to multiply the implicit sensitivity by a vector. These calculations are modifications on the example presented above and use 5000 data points. The memory requirements are calculated for a single distributed parameter ($K_{s}$) as well as five spatially distributed parameters ($K_{s},\alpha,n,\theta_{r}$, or $\theta_{s}$). Neither calculation includes the memory required to solve the matrix system, as such, the reported numbers underestimate the actual memory requirements for solving the inverse problem. The aim of this comparison is to demonstrate how the memory requirements scale, an appropriate solver must also be chosen for either method to solve the forward problem. When using an explicitly calculated sensitivity matrix to invert for additional physical properties, the memory footprint increases proportionally to the number of distributed physical properties; this is not the case for the implicit sensitivity calculation. For example, on a $128\times 128\times 128$ mesh, the explicit formation of the sensitivity requires 419 GB for five spatially distributed model parameters, which is five times the requirement for a single distributed model parameter (83.9 GB). For the implicit sensitivity on the same mesh, only 4.09 GB of memory is required, which is 1.2 times the requirement for a single distributed model parameter (3.54 GB). For this mesh, inverting for five spatially distributed parameters requires over 100 times less memory when using the implicit sensitivity algorithm, allowing these calculations to be run on modest computational resources.

7 Conclusions

The number of parameters that are estimated in Richards equation inversions has grown and will continue to grow as time-lapse data and geophysical data integration become standard in site characterizations. The increase in data quantity and quality provides the opportunity to estimate spatially distributed hydraulic parameters from Richards equation; doing so requires efficient simulation and inversion strategies. In this paper, we have shown a computationally efficient derivative-based optimization algorithm that does not store the Jacobian, but rather computes its effect on a vector (i.e. $\bf Jv$ or $\bf J^{\top}z$). By not storing the Jacobian, the size of the problem that we can invert becomes much larger. We have presented efficient methods to compute the Jacobian that can be used for all empirical hydraulic parameters, even if the functional relationship between parameters is obtained from laboratory experiments.

Our technique allows a deterministic inversion, which includes regularization, to be formulated and solved for any of the empirical parameters within Richards equation. For a full 3D simulation, as many as ten spatially distributed parameters may be needed, resulting in a highly non-unique problem; as such, we may not be able to reasonably estimate all hydraulic parameters. Depending on the setting, amount of a-priori knowledge, quality and quantity of data, the selection of which parameters to invert for may vary. Our methodology enables practitioners to experiment in 1D, 2D and 3D with full simulations and inversions, in order to explore the parameters that are important to a particular dataset. Our numerical implementation is provided in an open source repository (http://simpeg.xyz).

Acknowledgments

The funding for this work is provided through the Vanier Canada Graduate Scholarships Program and grants from The University of British Columbia (NSERC 22R47082). Thank you also to contributors of the SimPEG project (http://simpeg.xyz).