Characterizing the impact of model error in hydrogeologic time series recovery inverse problems

TL;DR

This paper investigates how systematic model errors affect hydrogeologic time series reconstruction, deriving bounds and diagnostics to improve understanding and handling of such errors in inverse problems.

Contribution

It introduces a systematic analysis of model error impact in hydrogeologic inverse problems, including new bounds and diagnostic criteria for transfer function uncertainties.

Findings

Derived upper and lower error bounds for model error effects.

Monte Carlo simulations validate the bounds and reveal error behavior.

Identified a diagnostic criterion for transfer function characterization.

Abstract

Hydrogeologic models are commonly over-smoothed relative to reality, owing to the difficulty of obtaining accurate high-resolution information about the subsurface. When used in an inversion context, such models may introduce systematic biases which cannot be encapsulated by an unbiased "observation noise" term of the type assumed by standard regularization theory and typical Bayesian formulations. Despite its importance, model error is difficult to encapsulate systematically and is often neglected. Here, model error is considered for a hydrogeologically important class of inverse problems that includes interpretation of hydraulic transients and contaminant source history inference: reconstruction of a time series that has been convolved against a transfer function (i.e., impulse response) that is only approximately known. Using established harmonic theory along with two results established here regarding triangular Toeplitz matrices, upper and lower error bounds are derived for the effect of systematic model error on time series recovery for both well-determined and over-determined inverse problems. A Monte Carlo study of a realistic hydraulic reconstruction problem is presented, and the lower error bound is seen informative about expected behavior. A possible diagnostic criterion for blind transfer function characterization is also uncovered.

Figures

Equations

$${\varphi }_n(t)=\frac{e^{\frac{t}{2}}}{n!}\frac{d^n}{d{t}^n}(e^{-t}{t}^n),$$

$$h(0,t)$$

$$b(\mathbf{x},t)$$

$$h(\mathbf{x},t)$$

$$\mathbf{c}=\mathbf{Ba},$$

$$\mathbf{B}=T\left[\begin{array}{ccccc}{b}_0& 0& 0& \dots & 0\\ b_1-b_0& b_0& 0& \dots & 0\\ b_2-b_1& b_1-b_0& b_0& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ b_{N-1}-b_{N-2}& b_{N-2}-b_{N-3}& b_{N-3}-b_{N-4}& \dots & b_0\end{array}\right].$$

$${\mathbf{B}}^{-1}\mathbf{c}=\mathbf{a}.$$

$${\left[\begin{array}{cc}A& 0\\ C& D\end{array}\right]}^{-1}=\left[\begin{array}{cc}{A}^{-1}& 0\\ -D^{-1}C{A}^{-1}& D^{-1}\end{array}\right],$$

$${\mathbf{F}}_{ik}$$

$${\mathbf{G}}_{kj}$$

$$(\mathbf{FG}{)}_{ij}$$

$$=$$

$$\mathbf{c}=\tilde{\mathbf{B}}\tilde{\mathbf{a}}.$$

$${\int }_0^{\infty }(h(0,t)-\tilde{h}(0,t){)}^{2}dt\approx {\Vert \mathbf{a}-\tilde{\mathbf{a}}\Vert }_2^2,$$

$${\Vert \mathbf{a}-\tilde{\mathbf{a}}\Vert }_2={\Vert (\mathbf{I}-{\tilde{\mathbf{B}}}^{-1}\mathbf{B})\mathbf{a}\Vert }_2.$$

$${\mathbf{e}}_0$$

$${\mathbf{e}}_1$$

$$|a_0|\left|1-\frac{b_0}{{\tilde{b}}_0}\right|\le {\Vert \mathbf{a}-\tilde{\mathbf{a}}\Vert }_2,$$

$$\left|1-\frac{b_0}{{\tilde{b}}_0}\right|\le \frac{{\Vert \mathbf{a}-\tilde{\mathbf{a}}\Vert }_2}{{\Vert \mathbf{a}\Vert }_2}.$$

$$b_0={\int }_0^{\infty }{e}^{-\frac{t}{2}}b(0,t)dt,$$

$${\left|a_0\left(1-\frac{b_0}{{\tilde{b}}_0}\right)\right|}^2+{\left|a_0\left(\frac{1}{{\tilde{b}}_0}\left[b_1-{\tilde{b}}_1\left(\frac{b_0}{{\tilde{b}}_0}\right)\right]\right)+a_1\left(1-\frac{b_0}{{\tilde{b}}_0}\right)\right|}^2\le {\Vert \mathbf{a}-\tilde{\mathbf{a}}\Vert }_2^2.$$

$$\sum _{l=1}^M\Vert {\mathbf{c}}_{\mathbf{l}}-\widetilde{{\mathbf{c}}_{\mathbf{l}}}{\Vert }_2^2=\underset{\tilde{\mathbf{a}}}{\min }\sum _{l=1}^M\Vert {\mathbf{c}}_{\mathbf{l}}-\widetilde{{\mathbf{c}}_{\mathbf{l}}}{\Vert }_2^2.$$

$${\mathbf{B}}_{\otimes }$$

$${\tilde{\mathbf{B}}}_{\otimes }$$

$$\mathbf{D}=\left[\begin{array}{c}{\mathbf{I}}_{\mathbf{N}}\\ {\mathbf{I}}_{\mathbf{N}}\\ ⋮\\ {\mathbf{I}}_{\mathbf{N}}\end{array}\right].$$

$$\sum _{l=1}^M\Vert {\mathbf{c}}_{\mathbf{l}}-\widetilde{{\mathbf{c}}_{\mathbf{l}}}{\Vert }_2^2={\mathbf{a}}^{\mathbf{T}}{\mathbf{D}}^{\mathbf{T}}{\mathbf{B}}_{\otimes }^{\mathbf{T}}{\mathbf{B}}_{\otimes }\mathbf{Da}-2{\mathbf{a}}^{\mathbf{T}}{\mathbf{D}}^{\mathbf{T}}{\mathbf{B}}_{\otimes }^{\mathbf{T}}{\tilde{\mathbf{B}}}_{\otimes }\mathbf{D}\tilde{\mathbf{a}}+{\tilde{\mathbf{a}}}^{\mathbf{T}}{\mathbf{D}}^{\mathbf{T}}{\tilde{\mathbf{B}}}_{\otimes }^{\mathbf{T}}{\tilde{\mathbf{B}}}_{\otimes }\mathbf{D}\tilde{\mathbf{a}}.$$

$$\frac{d}{d\tilde{\mathbf{a}}}\sum _{l=1}^M\Vert {\mathbf{c}}_{\mathbf{l}}-\widetilde{{\mathbf{c}}_{\mathbf{l}}}{\Vert }_2^2$$

$$=$$

$$\mathbf{D}\tilde{\mathbf{a}}={\tilde{\mathbf{B}}}_{\otimes }^{-1}{\mathbf{B}}_{\otimes }\mathbf{Da},$$

$$\tilde{\mathbf{a}}=\frac{1}{M}{\mathbf{D}}^{\mathbf{T}}{\tilde{\mathbf{B}}}_{\otimes }^{-1}{\mathbf{B}}_{\otimes }\mathbf{Da}.$$

$${\Vert \mathbf{a}-\tilde{\mathbf{a}}\Vert }_2$$

$$=$$

$${\Vert \mathbf{a}-\tilde{\mathbf{a}}\Vert }_2=\frac{1}{M}\sum _{l=1}^M{\Vert ({\mathbf{I}}_{\mathbf{N}}-{\tilde{\mathbf{B}}}_{\mathbf{l}}^{-1}{\mathbf{B}}_{\mathbf{l}})\mathbf{a}\Vert }_2.$$

$$\frac{{\Vert \mathbf{a}-\tilde{\mathbf{a}}\Vert }_2}{{\Vert \mathbf{a}\Vert }_2}\le {\Vert {\mathbf{I}}_{\mathbf{MN}}-{\tilde{\mathbf{B}}}_{\otimes }^{-1}{\mathbf{B}}_{\otimes }\Vert }_2.$$

Full text

1 Introduction

The effect of systematic model error is being increasingly studied in the hydrological literature, specifically in the context of forward modeling (Gupta et al, 2008; Vrugt et al, 2008; Lin and Beck, 2012; Gong et al, 2013; Vrugt and Sadegh, 2013; White et al, 2014) and data assimilation (Del Giudice et al, 2015b), with systematic treatments presented by Refsgaard et al (2006) and Gupta et al (2012). Some attention has also been directed explicitly at the effect of model error in inverse modeling (Hansen and Vesselinov, 2016).

Inverse analyses now form a major part of hydrogeologic practice, being of relevance to all forms of model parametrization, including hydraulic tomography and contaminant source identification. In the inverse problems literature, it is common to assume a perfect model, with all divergence between model prediction and observed data vector attributable to “noise” drawn from a symmetric, zero-mean probability distribution function. This theoretical approach underlies classical regularization methods such as Tikhonov and TSVD techniques (Hansen, 1992), and is also used for specifying the likelihood function in the Bayesian inversion paradigm (BuiThanh, 2012). That the approach of encoding all errors as unbiased parametric uncertainties may not be appropriate in hydrogeologic inverse modeling has been recognized. However in the absence of a paradigm that captures model error in a systematic fashion, the perfect model assumption remains common in practice (Lin and Beck, 2012; Del Giudice et al, 2015a). Thus, it is timely to consider formal analyses of the systematic model errors on inverse-model estimates.

Current approaches to quantifying the effect of model error are typically probabilistic, treating the impact of the model uncertainty on output with Bayesian (Krzysztofowicz, 1999) or information theoretic (Gong et al, 2013) formalisms. The uncertainty about model structure is modeled by parameterizing the model itself as a probability distribution function (pdf) linking inputs and outputs, or as a deterministic numerical model with pdfs defined on its state variables (see summary in Renard et al, 2010). To quantify uncertainty using either approach, some type of Monte Carlo computational analysis is indicated. Since computational exploration cannot proceed in infinite dimensions, it is of course unavoidable that explorations will impose some sort of coarse-grained parametrization which is at best approximately valid.

This paper follows a somewhat different path. The focus is on a specific class of inverse problems commonly faced by hydrogeologists: time series recovery problems with a temporal convolution structure. More concretely, this means the recovery of an input signal from one or more remote output signal measurements, where each output signal is generated by temporal convolution of the (shared) input signal with a (unique) transfer function, which is only approximately known. Hydrogeologic examples of such problems include the inference of hydraulic head history at some location of interest from available time series obtained at remote monitoring wells, and the inference of contaminant source histories from breakthrough curves.

For these problems, it is shown in Sect. 2 how it is possible to formally decompose the transfer function(s) as well as the input and output signals into generalized Fourier series. Some apparently new results concerning triangular Toeplitz matrices are established. Techniques of matrix analysis are then employed to derive concrete error bounds on the L2 signal reconstruction as a function of the error in dominant components of the transfer function(s). In Sect. 3, a Monte Carlo study of hydraulic inversion is presented which contextualizes the theoretical results shown in Sect. 2 and some empirical observations that go beyond the theoretical work are noted. Section 4 summarizes what has been learned and suggests interesting future research directions.

2 Derivation of error bounds

2.1 Laguerre expansion method

Consider a linear system in which it is possible to express a transient output signal, $h(\bm{x},t)$, at a location $\bm{x}$, resulting from a transient input signal at location $\bm{0}$, $h(\bm{0},t)$, by means of the convolution $h(\bm{x},t)=b(\bm{x},t)*h(\bm{0},t)$, where $b$ is a Green’s function (i.e., transfer function) representing the response to an instantaneous Dirac input signal at the origin (i.e. $h(\bm{0},t)=\delta(t)$). In general, for some fixed $\bm{x}$, the input signal, transfer function, and output signal can be expanded as a generalized Fourier series in a basis of Laguerre functions, ${\phi_{n}(\cdot)}$ ($n\geq 0$). Each Laguerre function is defined according to the formula:

[TABLE]

and together they form an orthonormal basis on $[0,\infty)$ (Abate et al, 1996). The series expansions are written as follows:

[TABLE]

where $T$ is a characteristic time of the problem, chosen to accelerate convergence. Let $\mathbf{c}$ be a vector of $N$ Laguerre coefficients, such that its $n$-th entry, $\mathbf{c}_{n}=c_{n}$. Similarly, define $\mathbf{a}$ to be a vector of $N$ Laguerre coefficients, such that $\mathbf{a}_{n}=a_{n}$. It has been shown (Hansen and Kueper, 2009) that, in general, these vectors of coefficients can be related by the matrix operation

[TABLE]

where $\mathbf{B}$ is the following lower triangular Toeplitz (LTT) matrix:

[TABLE]

Because $\mathbf{B}$ is a full-rank square matrix, if it is known perfectly then inversion is well defined (although not necessarily numerically stable):

[TABLE]

2.2 Results on triangular Toeplitz matrix manipulation

To continue the analysis of the last section, we need to establish some properties of LTT matrices. First, it is proven that the inverse of a (lower or upper) triangular Toeplitz matrix is itself a (lower or upper) triangular Toeplitz matrix. Second, it is proven that the product of two (both lower or both upper) triangular Toeplitz matrices is similarly a (lower or upper) triangular Toeplitz matrix. Without loss of generality, it is assumed the matrices are LTT in both proofs.

2.2.1 Lemma 1: Inversion of triangular Toeplitz matrices

Let $\mathbf{M_{N}}$ be an $N\times N$ LTT matrix, $N$ arbitrary, and let $\mathbf{M_{N}^{-1}}$ be its inverse. It may be shown by induction that $\mathbf{M_{N}^{-1}}$ is LTT. This argument makes repeated use of the following identity for block-triangular matrices (Bernstein, 2005, p. 71):

[TABLE]

where $A$, [math], $C$, and $D$ are compatibly-shaped sub-matrices.

The base case is trivial: note that that for any $2\times 2$ LTT matrix, $\mathbf{M_{2}}$, Eq. (8) implies directly that $\mathbf{M_{2}^{-1}}$ is LTT.

For the inductive step, assume that it has been established for $(N-1)\times(N-1)$ LTT matrices that their inverses are LTT. Define $\mathbf{M_{N-1}}$ to be the sub-matrix consisting of the first $N-1$ rows and first $N-1$ columns of an arbitrary LTT matrix, $\mathbf{M_{N}}$. Note that $\mathbf{M_{N-1}}$ is also LTT, and by our inductive assumption so is $\mathbf{M_{N-1}^{-1}}$. It is valid to apply Eq. (8) in two different ways. First, make the assignment $A\equiv\mathbf{M_{N-1}}$ and apply Eq. (8). This implies that $A^{-1}$ is LTT, and also that $\mathbf{M_{N}^{-1}}$ is lower triangular. This analysis has accounted for all but the $N$-th row of $\mathbf{M_{N}^{-1}}$. To see that the constant descending diagonals continue into the bottom row, note that the sub-matrix consisting of the last $N-1$ rows and last $N-1$ columns of $\mathbf{M_{N}}$ is also $\mathbf{M_{N-1}}$. Make the assignment $D\equiv\mathbf{M_{N-1}}$ and apply Eq. (8) again, implying that $D^{-1}$ is LTT, and also that $A^{-1}=D^{-1}$ (Note that these are both $(N-1)\times(N-1)$ matrices which are largely overlapping, and do not participate in the same block partitioning of $\mathbf{M_{N}^{-1}}$). It is thus shown that all descending diagonals of $\mathbf{M_{N}^{-1}}$ are constant (the single element $(\mathbf{M_{N}^{-1}})_{N1}$ can have any value without affecting this). It has thus been shown that, subject to our inductive assumption, $\mathbf{M_{N}^{-1}}$ is LTT, for arbitrary LTT $\mathbf{M_{N}}$.

By combination of base case and inductive step it follows that if $\mathbf{M_{N}}$ is an LTT $N\times N$ matrix then so is $\mathbf{M_{N}^{-1}}$, $\forall\ N\geq 2$. $\square$

2.2.2 Lemma 2: Multiplication of triangular Toeplitz matrices

For any two $N\times N$ matrices $\mathbf{F}$ and $\mathbf{G}$, it is true that the element $(\mathbf{FG})_{ij}=\sum_{k=1}^{N}\mathbf{F}_{ik}\mathbf{G}_{kj}$. If the matrices are also LTT, it follows that

[TABLE]

where $\mathbf{f}_{n}$ and $\mathbf{g}_{n}$ are the elements on the $n$-th diagonal of $\mathbf{F}$ and $\mathbf{G}$, respectively (with the main diagonal having index 0, the sub-diagonal having index 1, and so on). Then it follows that

[TABLE]

$(\mathbf{FG})_{ij}$ is thus a function only of $i-j$ and is zero for $i<j$. Thus, $\mathbf{FG}$ is LTT. $\square$

2.3 Effect of imperfect model: single observation location

Assume perfect knowledge of $h(\bm{x},t)$, but imperfect knowledge of $b(\bm{x},t)$, and a need to infer $h(\bm{0},t)$. The imperfect knowledge of $b$ will lead to an approximate solution $\tilde{h}(\bm{0},t)$, whose Laguerre coefficients lie in vector $\mathbf{\tilde{a}}$.

In matrix form, this can be written by distinguishing the (unknown) true matrix, $\mathbf{B}$, from the approximate matrix, $\mathbf{\tilde{B}}$, resulting from our imperfect knowledge of the Green’s function, $b$. The matrix inverse problem that is being solved is thus

[TABLE]

It is also true (by Parseval’s theorem) that the squared error of our source history estimate can be expressed in vector form via

[TABLE]

with equality in the limit $N\rightarrow\infty$. Although the error introduced by spectral leakage (i.e., series truncation) has been recognized as important in some geophysical inversion (Sneider and Trampert, 1999), many transfer functions and input signals are smooth in hydrology, and the approximate equality will be taken to be exact in subsequent analysis. The error analysis can thus be performed in the matrix domain.

Because $\mathbf{\tilde{B}}$ and $\mathbf{B}$ are invertible, there is a unique solution to Eq. (13) and thus:

[TABLE]

By application of Lemma 1 and Lemma 2, it follows that $\mathbf{I-\tilde{B}^{-1}B}$ is an LTT matrix. Let $\mathbf{\tilde{b}^{-1}}$, $\mathbf{b}$, and $\mathbf{e}$ be the vectors of coefficients on the diagonals of $\mathbf{\tilde{B}^{-1}}$, $\mathbf{B}$, and $\mathbf{I-\tilde{B}^{-1}B}$, respectively, indexed as in the proof of Lemma 2. By applying Eq. (8) with $A$ defined as the upper right $2\times 2$ sub-matrix of $\mathbf{\tilde{B}}$, it follows immediately that $\mathbf{\tilde{b}^{-1}}_{0}=\frac{1}{\tilde{b}_{0}}$ and $\mathbf{\tilde{b}^{-1}}_{1}=-\frac{\tilde{b}_{1}-\tilde{b}_{0}}{\tilde{b}_{0}^{2}}$. Also, from inspection of Eq. (6), $\mathbf{b}_{0}=b_{0}$ and $\mathbf{b}_{1}=b_{1}-b_{0}$. Then, by applying Eq. (8) and Eq. (12), it follows that

[TABLE]

Note that these elements (like all of $\mathbf{e}$) are zero when $\mathbf{\tilde{B}=B}$.

The following lower error bound follows from consideration of the first element of $\mathbf{(I-\tilde{B}^{-1}B)a}$:

[TABLE]

In the useful special case in which the input signal is an arbitrary decaying exponential (with rate constant by $1/2T$, noting that $T$ is a free parameter), the only nonzero term of its Laguerre series is $a_{0}$ and a lower bound on the relative error follows immediately:

[TABLE]

The coefficient $b_{0}$ is computed

[TABLE]

and similarly for $\tilde{b}_{0}$. It should be clear that if $b(x,t)$ and $\tilde{b}(x,t)$ have different shapes, particularly if $\tilde{b}(x,t)$ represents transmission through a homogeneous medium, and the true Green’s function, $b(x,t)$, is characterizes a medium that is heterogeneous or homogeneous with substantially different properties, then it is possible to have $b_{0}/\tilde{b}_{0}\gg 1$. In such cases, the error due to fitting of the inaccurate model overwhelms the signal and the errors of signal measurement (detection).

It is generally possible to use the approach developed here to generate lower error bounds relating the first $k$ terms of the sequences $\{a_{n}\}$, $\{b_{n}\}$, and $\{\tilde{b}_{n}\}$, for arbitrary $k$, depending on the amount of information available. For instance, for $k=2$:

[TABLE]

It is also possible to derive an upper bound, which does not depend on $\{a_{n}\}$, but which requires $k=N$ terms of the other sequences. This is seen in the next section.

2.4 Effect of imperfect model: multiple observation locations

In the case of $M$ monitoring locations, the problem is generally over-determined, and instead of directly computing the inverse, one may define the optimal solution, $\mathbf{\tilde{a}}$, as the one which minimizes the sum of squared residuals at each location, i.e., satisfies the following condition:

[TABLE]

This problem may be placed in a matrix form by defining the block diagonal matrices

[TABLE]

It is also useful to define the following block-columnar matrix of $M$, $N\times N$ identity matrices:

[TABLE]

Using this notation, Eq. (22) can be represented as:

[TABLE]

Differentiating,

[TABLE]

The optimal solution will be when this quantity equals zero, which is satisfied when

[TABLE]

or

[TABLE]

It is clear that

[TABLE]

By inspection of the last equality, it is apparent that

[TABLE]

This is simply the average of the individual model errors if only a single monitoring location were to be used (see Eq. (15)), implying that the error bound theory developed above can be carried over straightforwardly. If $\mathbf{a}$ is viewed as a random variable, the expected error is also not reduced by incorporating additional measurements, unlike the scenario of uncorrelated random noise.

One can also derive an upper bound on the relative error of $\mathbf{\tilde{a}}$ from Eq. (32) by noting that the Hölder norm $\lVert\cdot\rVert_{2}$ is submultiplicative, and that $\frac{1}{M}\lVert\mathbf{Da}\rVert_{2}=\lVert\mathbf{a}\rVert_{2}$. From this, it follows immediately that

[TABLE]

Naturally, this applies to single measurement location reconstruction as a special case. Practically, this depends on the full sets of coefficients $\{b_{n}\}$ and $\{\tilde{b}_{n}\}$, which is a greater information demand than for the lower bounds, which only involved relationships of the dominant components.

3 Monte Carlo study: reconstruction of hydraulic transients

In this section, the inference of a hydraulic head transient history along an aquifer boundary (which might be interpreted as a river stage transient, where the river is in a hydraulic connection with the aquifer) based on a time series of measurements made at a single nearby groundwater monitoring well is considered. This represents both an application of the above ideas, and also a study of independent interest.

3.1 Procedure

It is assumed here that the specific storage is known and spatially uniform, and that the log hydraulic conductivity is defined by a multi-Gaussian spatially random field whose mean is known, but which is otherwise unknown. Assuming flow is described by the groundwater flow equation on this heterogeneous conductivity field, $b(\bm{x},t)$ is determined as the head history at a fixed location, $\bm{x}$. A natural interpretive model, $\tilde{b}(\bm{x},t)$, is selected: the same groundwater flow equation, but solved on a homogeneous conductivity field that is everywhere equal to the mean of the true log hydraulic conductivity field.

The study is then performed according to the following basic procedure: First, a true, exponentially decaying, transient in the river stage on the aquifer domain boundary is specified, along with the location of a monitoring well at which a time series of measurements is to be made. The accuracy of reconstruction of the river stage transient from the transient at the monitoring well is studied, given an overly smooth model of the subsurface. For simplicity, the free parameter, $T$, is selected so that the Laguerre decomposition of the true transient is the vector $\mathbf{a}=<1,0,0,\dots,0>$. Next, 500 two-dimensional subsurface realizations are generated with different heterogeneous log-hydraulic conductivity fields, all of which have the same multi-Gaussian statistical correlation structure and geometric mean conductivity. Subsequently, using finite element analysis, head time series are computed at the monitoring well for each of the 500 subsurface realizations resulting from a Dirac head impulse at the river stage. Each impulse response (Green’s function) is decomposed as a vector of Laguerre coefficients, $\mathbf{b}$. Again, using finite element analysis, the impulse response at the well is computed, but assuming a uniform hydraulic conductivity field with the same geometric mean hydraulic conductivity as used in each of the heterogeneous realizations. This impulse response is decomposed as a vector of Laguerre coefficients, $\mathbf{\tilde{b}}$. Finally, for each realization, the reconstruction error, $\left\lVert\mathbf{(I-\tilde{B}^{-1}B)a}\right\rVert_{2}$, is computed and compared with the analytical lower bound in Eq. (19). Statistics about this quantity are tabulated so that its relationship to qualitative features of the inverse model discrepancy may be studied.

Both the true solution and the interpretive model are described by the following system of equations:

[TABLE]

solved subject to the initial and boundary conditions

[TABLE]

where the only difference enters due to different choices for $K(\bm{x})$. In the above equations, $\bm{x}\>[\mathrm{L}]$ represents the location, $t\>[\mathrm{T}]$ represents time, $S_{s}\>[\mathrm{L^{-1}}]$ represents specific storage, $h\>[\mathrm{L}]$ represents hydraulic head, $\bm{q}\>[\mathrm{LT^{-1}}]$ represents groundwater flux, $K\>[\mathrm{LT^{-1}}]$ represents hydraulic conductivity, and $\bm{n}\>[1]$ is the outward-facing unit normal vector. For vector quantities, the units reported are for each of their components.

More concretely, a two-dimensional model of saturated flow in a heterogeneous porous medium is defined over the rectangular domain $\bm{D}=(0,L_{1})\times(0,L_{2})$, where $L_{1}=10\ [L]$ and $L_{2}=4\ [L]$ ($L$ is any consistent length unit), and specific storage, $S_{s}=1$ (Fig. 1). $\Gamma_{L}$ represents the left boundary of $\bm{D}$ (at $x=0$), $\Gamma_{R}$ represents the right boundary of $\bm{D}$ (at $x=10$), and $\Gamma_{S}$ represents the union of the other two sides of $\bm{D}$ (at $y=0$ and $y=4$, respectively). Let $Y(\bm{x},\omega)=\ln[K(\bm{x},\omega)]$ be a random field, where $\omega$ belongs to the space of random events $\bm{\Omega}$. Assuming $Y(\bm{x},\omega)$ is Gaussian with zero mean and a separable exponential covariance function,

[TABLE]

where $\sigma^{2}_{Y}=2$, $\eta_{1}=4$ and $\eta_{2}=2$ are the variance and the correlation lengths of the random field.

For the Monte Carlo study, a set of 500 realizations of $b(t)$ is generated by setting $c(t)=\delta(t)$, generating 500 $\ln[K(\bm{x})]$ fields (Fig. 1, bottom), using a 100-term truncated Karhunen-Loève expansions (KLE) to represent the field as weighted sums of predefined spatially variable orthonormal functions (Zhang and Lu, 2004), and solving Eqs. (35-40) on each. The numerical solution is evaluated using the FEniCS (Logg et al, 2012) package to discretize Eqs. (35-36), using finite element methods in space and an implicit Euler method in time. Simulated hydraulic head data is recorded at the point (4, 2) for each realization of $K(\bm{x})$ (Fig. 2), and the 50-term LEM is used along with $T=100\ [T]$ to reconstruct the time series of hydraulic head values on the left boundary (Fig. 3). The same procedure is followed to generate the interpretive model, $\tilde{b}(t)$ except $K(\bm{x})=1\>\forall\bm{x}$ (i.e., $Y(\bm{x})=0\>\forall\bm{x}$) is employed.

3.2 Discussion of results

Given that the random hydraulic conductivity fields chosen for the forward modeling were only moderately heterogeneous, and choice of a spatially uniform interpretive model is a natural response to unresolved heterogeneity, the wide array of possible reconstructions is notable. In particular, a bifurcation of the response classes was noted based on whether the peaks of the forward model impulse response, $b$, preceded or lagged that of the interpretive model, $\tilde{b}$. In case in which the interpretive model predicted a faster response than existed in reality, the reconstruction of the decaying exponential boundary condition was typically smooth, with its peak at a time significantly greater than zero: the delayed reconstruction of the boundary condition compensated for the over-rapid model. However, because signal causality is enforced, the model cannot respond to an earlier-than-anticipated arrival with a non-zero signal at negative time. Instead, the optimal reconstruction features a large peak at time zero, followed by decaying corrective oscillations. This bifurcation of behavior is potentially useful as a model diagnostic tool that does not require any a priori knowledge of the true model (other than that it possesses a unimodal structure): multiple candidate interpretive models could be tested with peaks at different locations, and the true peak location pinpointed by the disappearance of the spurious oscillations.

In Fig. 4 (top), the empirical pdf for the L2 error, normalized by the L2 norm of the signal being reconstructed is shown. It is apparent that even for moderate heterogeneity, reconstruction error on the same magnitude as the signal itself is to be expected. In Fig. 4 (bottom), the L2 error of approximation in the reconstruction of the boundary condition for each of the 500 realizations is plotted against the zero-order lower error bound in Eq. (19).

4 Summary and conclusion

Systematic model error was considered in the context of inverse problems in systems whose output signal is determined by convolution of an input signal with a transfer function, or impulse response, which describes system behavior. Using a generalized Fourier series expansion in Laguerre basis functions, it was possible to translate the signal reconstruction inverse problem into a matrix inverse problem whose structure may be analyzed using some classic, and some apparently new, results in matrix algebra. It is thus seen possible to place upper and lower bounds on the L2 signal reconstruction error as in terms of the transfer function infidelity.

The inverse problem of recovering river level history from remote measurements at a well, which has convolution structure, was chosen for Monte Carlo study. Forward predictions were generated by solving the groundwater flow equation on a mildly heterogeneous hydraulic conductivity field, and these were interpreted using the groundwater flow equation, assuming a homogeneous hydraulic conductivity field. The L2 reconstruction error for the river level was established for all realizations, and this error was compared with the error bound developed above. The simple lower bound derived here (Sect. 2) was found to be informative regarding the reconstruction error in the specific realizations. A qualitative bifurcation in the reconstructed signal was discovered, depending on the location of the peak of the interpretation model transfer function relative to that of the true model. Looking forward, this may prove to be a useful tool for transfer function identification.

The Laguerre expansion approach, because of its high degree of structure, relative simplicity and computational efficiency, may also prove to be a profitable foundation for further analysis of systematic model error. The systematic model errors are commonly ignored in the theoretical and practical inverse-model analyses. The matrix transformation of the inverse problem developed here can be also applied for other problems of interest such as groundwater contaminant transport, propagation of low-frequency seismic waves, heat flow, infectious disease transmission, population dynamics, spreading chemical/biochemical substances in atmosphere, and many others.

Acknowledgements

The authors acknowledge the support of the LANL Environmental Programs. All data is synthetic; the authors maintain an archive of codes and simulation output employed in the paper.