Data assimilation with multiple types of observation boreholes via the ensemble Kalman filter embedded within stochastic moment equations

We employ an approach based on the ensemble Kalman filter coupled with stochastic moment equations (MEs-EnKF) of groundwater flow to explore the dependence of conductivity estimates on the type of available information about hydraulic heads in a three-dimensional randomly heterogeneous field where convergent flow driven by a pumping well takes place. To this end, we consider three types of observation devices corresponding to (i) multi-node monitoring wells equipped with packers (Type A) and (ii) partially (Type B) and (iii) fully (Type C) screened wells. We ground our analysis on a variety of synthetic test cases associated with various configurations of these observation wells. Moment equations are approximated at second order (in terms of the standard deviation of the natural logarithm, Y, of conductivity) and are solved by an efficient transient numerical scheme proposed in this study. The use of an inflation factor imposed to the observation error covariance matrix is also analyzed to assess the extent at which this can strengthen the ability of the MEs-EnKF to yield appropriate conductivity estimates in the presence of a simplified modeling strategy where flux exchanges between monitoring wells and aquifer are neglected. Our results show that (i) the configuration associated with Type A monitoring wells leads to conductivity estimates with the (overall) best quality, (ii) conductivity estimates anchored on information from Type B and C wells are of similar quality, (iii) inflation of the measurement-error covariance matrix can improve conductivity estimates when a simplified flow model is adopted, and (iv) when compared with the standard Monte Carlo-based EnKF method, the MEs-EnKF can efficiently and accurately estimate conductivity and head fields.

1 Introduction

Parameter estimation for groundwater system modeling is a key and important challenge due to our incomplete knowledge of the spatial distributions of hydrogeological attributes, such as hydraulic conductivity. The ensemble Kalman filter (EnKF; Evensen, 1994) is a powerful approach to parameter estimation in subsurface flow (Hendricks Franssen and Kinzelbach, 2008; Zheng et al., 2019) and solute transport (Liu et al., 2008; Li et al., 2012; Chen et al., 2018; Xu and Gomez-Hernandez, 2018) scenarios. Estimated system parameters can include conductivity (Botto et al., 2018), permeability (Zovi et al., 2017), porosity (Li et al., 2012), specific storage (Hendricks Franssen et al., 2011), dispersivity (Liu et al., 2008), riverbed conductivity (Kurtz et al., 2014), or unsaturated flow characteristic quantities (Zha et al., 2019; Li et al., 2020).

The EnKF can assimilate data sequentially through a real-time updating process. Alternatively, all collected measurements can be assimilated simultaneously, for example, within a typical model calibration framework. With reference to the latter aspect, the EnKF becomes an ensemble smoother (ES, Van Leeuwen and Evensen, 1996), as it is associated with a smoothing probability density function (PDF) rather than a filtering PDF (Jazwinski, 1967). With reference to the ES, observations in the past and current stages are assimilated only once, thus yielding increased efficiency with respect to the EnKF (Skjervheim et al., 2011). Iterative forms of the EnKF and ES, usually denoted by IEnKF (Gu and Oliver, 2007; Sakov et al., 2012) and IES (Chen and Oliver, 2013; Emerick and Reynolds, 2013; Luo et al., 2015; Chang et al., 2017; Li et al., 2018), have been developed to improve assimilation performance in scenarios characterized by strongly nonlinear behaviors. A variety of studies investigate challenges linked to such (ensemble) data assimilation algorithms, including, e.g., the possibility of coping with non-Gaussian model parameter distributions (Zhou et al., 2011; Li et al., 2018), physical unphysical results stemming from the estimation workflow (Wen and Chen, 2006; Song et al., 2014), or spurious correlations (Panzeri et al., 2013; Bauser et al., 2018; Luo et al., 2019; Soares et al., 2019). All of these works contribute to improve the robustness of these algorithms for parameter estimation in complex environmental systems.

Recent studies include the work of Xia et al. (2018), who tackle conductivity estimation in a two-dimensional variable-density flow setting using a localized IEnKF to balance central processing unit (CPU) time and estimation accuracy. Bauser et al. (2018) develop an adaptive covariance inflation method for the EnKF to reduce the negative effect of spurious correlations and illustrate an application of the method in a soil hydrology field context. Mo et al. (2019) use a deep-learning-based model as a surrogate of a solute transport model to reduce the CPU time associated with ensemble-based data assimilation through an iterative local update ensemble smoother in a contaminant identification problem considering a synthetic two-dimensional heterogeneous conductivity field. Li et al. (2020) compare benefits and drawbacks of embedding machine-learning-based (artificial neural network, ANN) and physics-based models into an IES for a set of synthetic unsaturated flow scenarios and find that (a) the performance of an IES relying on the Richards' equation is significantly impacted by soil heterogeneity, initial, and boundary conditions, and (b) an IES based on either ANN or Richards' equation can be notably affected by the quality of the measurements.

In this broad framework, it is noted that the accuracy of parameter estimation for a given environmental system is jointly determined by the ability of the mathematical model to describe the system of interest (Sakov and Bocquet, 2018; Alfonzo and Oliver, 2020; Luo, 2019; Evensen, 2019), the ability of the assimilation algorithm used (Emerick and Reynolds, 2013; Bocquet and Sakov, 2014), as well as by the quantity and quality of available observations (Zha et al., 2019; Xia et al., 2018, and references therein).

With reference to a groundwater system, data that are commonly collected in a borehole and then employed for parameter estimation include head (water level or pressure), solute concentration, and/or in some cases fluxes. A well screen opened at multiple depths can provide information associated with preferential pathways of flow and/or solute transport. Hydraulic heads observed in such a setting can be considered to constitute an integrated type of information and to be representative of an average system state (Elci et al., 2001, 2003; Konikow et al., 2009; Zhang et al., 2019). Elci et al. (2001, 2003) conclude that the use of long-screen wells to collect measurements should be approached with caution, as these can yield misleading and ambiguous information concerning, e.g., hydraulic head, solute concentration, location of contaminant source, and plume geometry. These types of monitoring wells can be found in a variety of field settings where head and/or solute concentration data are collected (see, e.g., Elci et al. (2001, 2003), Post et al. (2007), Konikow et al. (2009), Zhang et al. (2019), and references therein). As an alternative, somehow localized information could be provided through the use of packers. Installing the latter can be costly and in some cases impractical.

Here, we aim to explore the effect that assimilating hydraulic head information collected over time within wells equipped with screens of differing lengths can have on our ability to characterize the spatial distribution of conductivity of a three-dimensional fully saturated heterogeneous aquifer. We consider multi-node wells (Konikow et al., 2009) to represent observation boreholes that can be (a) equipped with packers to mimic pointlike measurements, (b) fully screened, or (c) partially penetrating. To this end, we focus on a convergent flow scenario driven by a partially penetrating pumping well operating in a three-dimensional heterogeneous conductivity field. Hydraulic head information is collected at a network of multi-node wells to represent data associated with screened intervals of differing lengths along the vertical. We consider synthetic scenarios to provide transparent comparative analyses of the extent at which the quality of the estimated conductivity fields is influenced by the type of multi-node wells considered.

Data assimilation is performed by relying on an EnKF coupled with stochastic moment equations (MEs) of transient groundwater flow (e.g., Tartakovsky and Neuman, 1998a, b; Zhang, 2002; Ye et al., 2004). The latter are approximated at second order (in terms of the standard deviation of the natural logarithm of hydraulic conductivity) and are solved by an efficient numerical scheme proposed in this study.

While we refer to Zhang (2002) and Winter et al. (2003) for reviews of MEs in heterogeneous conductivity fields, we recall that MEs of groundwater flow have been previously incorporated into geostatistical inverse modeling approaches (e.g., Hernandez et al., 2003) or stochastic pumping test interpretation (Neuman et al., 2004, 2007) and have been considered in field settings (Riva et al., 2009; Bianchi Janetti et al., 2010; Panzeri et al., 2015). More recent developments have allowed embedding stochastic MEs of steady-state groundwater flow in model reduction strategies (Xia et al., 2020). MEs of transient groundwater flow have also been framed in the context of data assimilation or parameter estimation approaches based on the EnKF approach (Li and Tchelepi, 2006; Panzeri et al., 2013, 2014).

Panzeri et al. (2013, 2014, 2015) present an approach for data assimilation (hereafter termed the MEs-EnKF) that relies on embedding MEs of groundwater flow within an EnKF framework. They (a) demonstrate that the method does not suffer from spurious correlation, thus avoiding resorting to any localization or inflation techniques; (b) document the computational feasibility and accuracy of the approach in two-dimensional synthetic log-conductivity domains; and then (c) explore a first field application to estimate log-transmissivity through assimilation of drawdown data collected during a series of cross-hole pumping tests.

An aspect that still somehow limits the advantages of the MEs-EnKF is related to the formulation of MEs in terms of a Green's function approach (see also Ye et al., 2004). One is then required to solve the equation satisfied by a (zero-order mean) Green's function for each node of the numerical grid employed to discretize the computational domain. While one can take advantage of symmetries related to the evaluation of the Green's function, Panzeri et al. (2014) show in their illustrative examples that the CPU time required by the MEs-EnKF is equivalent to performing a classical EnKF relying on a collection of 35 000 Monte Carlo (MC) realizations. The negative impact of this computational scheme could be aggravated in three-dimensional scenarios. Here, we circumvent this issue by solving MEs for three-dimensional transient groundwater flow by relying on the (second-order accurate) approximations of MEs presented by Zhang (2002).

The remainder of the work is structured as follows. Section 2 details the main elements associated with the mathematical background of MEs and multi-node wells. Section 3 introduces the coupling between MEs and the EnKF approach. Section 4 illustrates the synthetic settings we analyze together with the criteria according to which the performance of the MEs-EnKF and the standard Monte Carlo-based EnKF (MC-EnKF) is assessed. Section 5 is devoted to the presentation and analysis of the key results. The main conclusions of this work are presented in Sect. 6.

2 Theoretical background

2.1 Stochastic moment equations for groundwater flow

We consider transient groundwater flow in a three-dimensional domain Ω described by

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& S_{\mathrm{s}}\left(\mathit{x}\right){\displaystyle \frac{\partial h(\mathit{x},t)}{\partial t}}+{\mathrm{\nabla }}_x\cdot \mathit{q}(\mathit{x},t)=f(\mathit{x},t)\mathrm{with}\\ & \mathit{q}(\mathit{x},t)=-K\left(\mathit{x}\right)\mathrm{\nabla }h(\mathit{x},t),\end{array}\end{array}$$

subject to initial and boundary conditions

$$\begin{array}{}\text{(2)}& h\left(\mathit{x},t_{\mathrm{0}}\right)=H_{\mathrm{0}}\left(\mathit{x}\right)\mathit{x}\in \mathrm{\Omega }\end{array}$$

$$\begin{array}{}\text{(3)}& {\displaystyle }& {\displaystyle }h(\mathit{x},t)=H(\mathit{x},t)\mathit{x}\in {\mathrm{\Gamma }}_D\text{(4)}& {\displaystyle }& {\displaystyle }[-\mathit{q}(\mathit{x},t\left)\right]\cdot \mathit{n}\left(\mathit{x}\right)=\mathbf{Q}(\mathit{x},t)\mathit{x}\in {\mathrm{\Gamma }}_N,\end{array}$$

where x denotes the vector of Euclidian coordinates, t is time, K is hydraulic conductivity, S_s is specific storage, h is hydraulic head, q is Darcy flux, f is a forcing term, H₀(x) is initial hydraulic head, H(x,t) is head along the Dirichlet boundary, and Q is a prescribed flux along the Neuman boundary. In the following, we consider S_s(x), H(x,t), and Q(x,t) as deterministic, while H₀(x), f(x,t) and K(x) are taken to be random quantities.

The natural logarithm of hydraulic conductivity, Y(x)=ln K(x), is assumed to be a second-order stationary process correlated in space with mean 〈Y(x)〉 and variance ${\mathit{\sigma }}_Y^{\mathrm{2}}$. Tartakovsky and Neuman (1998a, b) derive integro-differential MEs to compute space–time dynamics of (ensemble) means and covariances of hydraulic heads and fluxes. They then resort to a perturbation approach to derive recursive approximations of these otherwise exact integro-differential MEs. Ye et al. (2004) solve second-order (in the standard deviation of Y, σ_Y) approximations of these MEs by finite elements for superimposed mean-uniform and convergent flows for two-dimensional settings. Since numerical solutions of moment equations are heavy in terms of computational resources (Zhang, 2002; Ye et al., 2004), in the following subsections we illustrate a workflow that enables us to evaluate all quantities of interest (up to second order in σ_Y) with reduced computational efforts.

2.1.1 Mean head and flux

We start by expressing a given random quantity, ℐ, as the sum of its (ensemble) mean, 〈ℐ〉, and a zero-mean random fluctuation, ℐ^′. Here, and in the following, mean head and flux are approximated up to second order in σ_Y as

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& \langle h(\mathit{x},t)\rangle \approx \langle h^{\left(\mathrm{0}\right)}(\mathit{x},t)\rangle +\langle h^{\left(\mathrm{2}\right)}(\mathit{x},t)\rangle ;\\ & \langle \mathit{q}(\mathit{x},t)\rangle \approx \langle {\mathit{q}}^{\left(\mathrm{0}\right)}(\mathit{x},t)\rangle +\langle {\mathit{q}}^{\left(\mathrm{2}\right)}(\mathit{x},t)\rangle ,\end{array}\end{array}$$

and the following equations hold

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }\left\{\begin{array}{ll}{S}_{\mathrm{s}}\left(\mathit{x}\right)\frac{\partial \langle h^{\left(\mathrm{0}\right)}(\mathit{x},t)\rangle }{\partial t}+{\mathrm{\nabla }}_x\cdot \langle {\mathit{q}}^{\left(\mathrm{0}\right)}(\mathit{x},t)\rangle =\langle f^{\left(\mathrm{0}\right)}(\mathit{x},t)\rangle & \mathit{x}\in \mathrm{\Omega }\\ \langle {\mathit{q}}^{\left(\mathrm{0}\right)}(\mathit{x},t)\rangle =-K_{\mathrm{G}}\left(\mathit{x}\right)\mathrm{\nabla }\langle h^{\left(\mathrm{0}\right)}(\mathit{x},t)\rangle & \mathit{x}\in \mathrm{\Omega }\\ \langle h^{\left(\mathrm{0}\right)}(\mathit{x},t)\rangle =\langle H_{\mathrm{0}}^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)\rangle & \mathit{x}\in \mathrm{\Omega }\\ \langle h^{\left(\mathrm{0}\right)}(\mathit{x},t)\rangle =H(\mathit{x},t)& \mathit{x}\in {\mathrm{\Gamma }}_D\\ -\langle {\mathit{q}}^{\left(\mathrm{0}\right)}(\mathit{x},t)\rangle \cdot \mathit{n}\left(\mathit{x}\right)=\mathbf{Q}(\mathit{x},t)& \mathit{x}\in {\mathrm{\Gamma }}_N\end{array}\right.,\text{(7)}& {\displaystyle }& {\displaystyle }\left\{\begin{array}{ll}{S}_{\mathrm{s}}\left(\mathit{x}\right)\frac{\partial \langle h^{\left(\mathrm{2}\right)}(\mathit{x},t)\rangle }{\partial t}+{\mathrm{\nabla }}_x\cdot \langle {\mathit{q}}^{\left(\mathrm{2}\right)}(\mathit{x},t)\rangle =\langle f^{\left(\mathrm{2}\right)}(\mathit{x},t)\rangle & \mathit{x}\in \mathrm{\Omega }\\ \langle {\mathit{q}}^{\left(\mathrm{2}\right)}(\mathit{x},t)\rangle =-K_{\mathrm{G}}\left(\mathit{x}\right)& \\ \left({\mathrm{\nabla }}_x\langle h^{\left(\mathrm{2}\right)}(\mathit{x},t)\rangle +\frac{{\mathit{\sigma }}_Y^{\mathrm{2}}}{\mathrm{2}}{\mathrm{\nabla }}_x\langle h^{\left(\mathrm{0}\right)}(\mathit{x},t)\rangle \right)+\mathit{r}(\mathit{x},t)& \mathit{x}\in \mathrm{\Omega }\\ \langle h^{\left(\mathrm{2}\right)}(\mathit{x},t)\rangle =\langle H_{\mathrm{0}}^{\left(\mathrm{2}\right)}\left(\mathit{x}\right)\rangle & \mathit{x}\in \mathrm{\Omega }\\ \langle h^{\left(\mathrm{2}\right)}(\mathit{x},t)\rangle =\mathrm{0}& \mathit{x}\in {\mathrm{\Gamma }}_D\\ -\langle {\mathit{q}}^{\left(\mathrm{2}\right)}(\mathit{x},t)\rangle \cdot \mathit{n}\left(\mathit{x}\right)=\mathrm{0}& \mathit{x}\in {\mathrm{\Gamma }}_N\end{array}\right..\end{array}$$

Here, superscript (i) indicates terms that are strictly of order i (in terms of powers of σ_Y), K_G(x)=e^〈Y(x)〉 is the geometric mean of K(x), and r(x,t) is the second-order residual flux evaluated as $\mathit{r}(\mathit{x},t)=\underset{y\to x}{lim}[-{\mathrm{\nabla }}_{x}u(\mathit{y},\mathit{x},t\left)\right]$ (e.g., Xia et al. (2019) and references therein), where $u(\mathit{y},\mathit{x},t)=\langle K^{\prime }\left(\mathit{y}\right)h^{\prime }(\mathit{x},t){\rangle }^{\left(\mathrm{2}\right)}$ is the second-order approximation of the cross-covariance between hydraulic head and conductivity, computed as detailed in Sect. 2.1.2.

2.1.2 Cross-covariance between hydraulic head and conductivity

Multiplying Eqs. (1)–(4) by K(y) and taking expectation yield the following equation governing the evolution of $u(\mathit{y},\mathit{x},t)$

$$\begin{array}{}\text{(8)}& \left\{\begin{array}{ll}{S}_{\mathrm{s}}\left(\mathit{x}\right)\frac{\partial u(\mathit{y},\mathit{x},t)}{\partial t}=& \\ ={\mathrm{\nabla }}_x\cdot \left[K_{\mathrm{G}}\left(\mathit{x}\right){\mathrm{\nabla }}_{x}u(\mathit{y},\mathit{x},t)-K_{\mathrm{G}}\left(\mathit{y}\right)C_Y(\mathit{x},\mathit{y})\right.& \\ \left.\langle {\mathit{q}}^{\left(\mathrm{0}\right)}(\mathit{x},t)\rangle \right]+C_{\mathrm{f}K}(\mathit{y},\mathit{x},t)& \\ & \mathit{x}\in \mathrm{\Omega }\\ u\left(\mathit{y},\mathit{x},t_{\mathrm{0}}\right)=U_{\mathrm{0}}(\mathit{y},\mathit{x})& \mathit{x}\in \mathrm{\Omega }\\ u(\mathit{y},\mathit{x},t)=\mathrm{0}& \mathit{x}\in {\mathrm{\Gamma }}_D\\ \left[K_{\mathrm{G}}\left(\mathit{x}\right){\mathrm{\nabla }}_{x}u(\mathit{y},\mathit{x},t)-K_{\mathrm{G}}\left(\mathit{y}\right)C_Y(\mathit{x},\mathit{y})\langle {\mathit{q}}^{\left(\mathrm{0}\right)}(\mathit{x},t)\rangle \right]& \\ \cdot \mathit{n}\left(\mathit{x}\right)=\mathrm{0}& \mathit{x}\in {\mathrm{\Gamma }}_N\end{array}\right..\end{array}$$

Here, $C_Y(\mathit{x},\mathit{y})=\langle Y^{\prime }\left(\mathit{x}\right)Y^{\prime }\left(\mathit{y}\right)\rangle$ is the covariance of the log-conductivity field, $C_{\mathrm{f}K}(\mathit{y},\mathit{x},t)=\langle K^{\prime }\left(\mathit{y}\right)f^{\prime }(\mathit{x},t){\rangle }^{\left(\mathrm{2}\right)}$ is the second-order cross-covariance between conductivity K(y) and forcing term f(x,t), and $U_{\mathrm{0}}(\mathit{y},\mathit{x})=\langle K^{\prime }\left(\mathit{y}\right)H_{\mathrm{0}}^{\prime }\left(\mathit{x}\right){\rangle }^{\left(\mathrm{2}\right)}$ is the second-order approximation of the cross-covariance between K and initial hydraulic head. Note that U₀(y,x) vanishes when H₀ is deterministic.

2.1.3 Head covariance

The equation governing the evolution of the (second-order) head covariance between space–time locations (y, τ) and (x, t), $C_{\mathrm{h}}(\mathit{y},\mathit{x},\mathit{\tau },t)=\langle h^{\prime }(\mathit{y},\mathit{\tau })h^{\prime }(\mathit{x},t){\rangle }^{\left(\mathrm{2}\right)}$, is given by (Zhang, 2002)

$$\begin{array}{}\text{(9)}& \left\{\begin{array}{ll}{S}_{\mathrm{s}}\left(\mathit{x}\right)\frac{\partial C_{\mathrm{h}}(\mathit{y},\mathit{x},\mathit{\tau },t)}{\partial t}=& \\ ={\mathrm{\nabla }}_x\cdot \left[K_{\mathrm{G}}\left(\mathit{x}\right){\mathrm{\nabla }}_x{C}_{\mathrm{h}}(\mathit{y},\mathit{x},\mathit{\tau },t)+u(\mathit{x},\mathit{y},\mathit{\tau })\right.& \\ \left.{\mathrm{\nabla }}_x\langle h^{\left(\mathrm{0}\right)}(\mathit{x},t)\rangle \right]+C_{\mathrm{fh}}(\mathit{y},\mathit{x},\mathit{\tau },t)& \\ & \mathit{x}\in \mathrm{\Omega }\\ C_{\mathrm{h}}\left(\mathit{y},\mathit{x},\mathit{\tau },t_{\mathrm{0}}\right)=\langle h^{\prime }(\mathit{y},\mathit{\tau })h^{\prime }\left(\mathit{x},t_{\mathrm{0}}\right){\rangle }^{\left(\mathrm{2}\right)}& \mathit{x}\in \mathrm{\Omega }\\ C_{\mathrm{h}}(\mathit{y},\mathit{x},\mathit{\tau },t)=\mathrm{0}& \mathit{x}\in {\mathrm{\Gamma }}_D\\ \left[K_{\mathrm{G}}\left(\mathit{x}\right){\mathrm{\nabla }}_x{C}_{\mathrm{h}}(\mathit{y},\mathit{x},\mathit{\tau },t)+u(\mathit{x},\mathit{y},\mathit{\tau }){\mathrm{\nabla }}_x\langle h^{\left(\mathrm{0}\right)}(\mathit{x},t)\rangle \right]& \\ \cdot \mathit{n}\left(\mathit{x}\right)=\mathrm{0}& \mathit{x}\in {\mathrm{\Gamma }}_N\end{array}\right.,\end{array}$$

where $u(\mathit{x},\mathit{y},\mathit{\tau })=\langle K^{\prime }\left(\mathit{x}\right)h^{\prime }(\mathit{y},\mathit{\tau }){\rangle }^{\left(\mathrm{2}\right)}$ is given by Eq. (8) and $C_{\mathrm{fh}}(\mathit{y},\mathit{x},\mathit{\tau },t)=\langle h^{\prime }(\mathit{y},\mathit{\tau })f(\mathit{x},t){\rangle }^{\left(\mathrm{2}\right)}$ is the second-order cross-covariance between forcing term f(x,t) and hydraulic head h(y,τ). To minimize redundancy, hereinafter we omit stating that all cross-/auto-covariances of quantities of interest appearing in our formulations are to be considered as second-order approximations. It is worth noting that spatially heterogenous conductivities of aquifer systems are often modeled through a single, in some cases multimodal, distribution (Winter et al., 2003). This approach corresponds to a homogenization of conductivity values that might be associated with diverse geomaterials within a unique system. Otherwise, the domain can be conceptualized as composed by zones, each associated with a given geomaterial and hydrogeological attributes. This leads to modeling the system under investigation as composed by a collection of disjoint blocks, whose location might be uncertain and within which a quantity such as conductivity can be spatially heterogeneous (see, e.g., Winter and Tartakovsky (2000, 2002), Winter et al. (2002, 2003), Guadagnini et al. (2004), Short et al. (2010), Perulero Serrano et al. (2014), Bianchi Janetti et al. (2019), and references therein). In this framework one can represent conductivity within each block upon relying on a distribution associated with low to mild variance, which is compatible with the order of approximation associated with the groundwater flow moment equations we consider (Winter and Tartakovsky (2002), Winter et al. (2002, 2003), and references therein). The scenario we investigate can then be seen as corresponding to the type of internal variability associated with a given geologic unit.

2.2 Monitoring wells

We consider three kinds of observation wells, leading to three diverse types of hydraulic head information (see Fig. 1). Type A wells are characterized by packers located at three depths, where pointwise hydraulic head observations are collected. Otherwise, Type B and/or C wells represent partially and fully penetrating wells, respectively, and provide hydraulic head values that are averaged along the corresponding screened intervals. Note that, even though there is no pumping from B- and C-wells, there are flux exchanges between these wells and the surrounding aquifer system, as opposed to the setting associated with packers (A-wells). Such flow is related to the difference between the water level within the well and hydraulic head values along the borehole.

Figure 1 Type of monitoring wells: pointwise (Type A), partially (Type B) and fully penetrating (Type C) observation boreholes.

Following Konikow et al. (2009), neglecting linear (due to skin effects) and nonlinear (due to turbulent flow) well loss terms, the water level at well I, $h_I^{\mathrm{w}}$, (Type B and/or C) at a given time t (omitted in the following equations for brevity) can be evaluated through

$$\begin{array}{}\text{(10)}& h_I^{\mathrm{w}}={\displaystyle \frac{\sum _{i=\mathrm{1}}^n{b}_i{K}_i{h}_i}{\sum _{i=\mathrm{1}}^n{b}_i{K}_i}},\end{array}$$

where n is the number of nodes in the multi-node observation well I, i.e., the number of cells according to which the well screen is discretized; h_i, b_i, and K_i are the hydraulic head, thickness, and conductivity of the cell of the numerical grid whose centroid corresponds to the ith node in the multi-node well, respectively. Note that Eq. (10) has been derived assuming that the flux exchange, i.e., the flow into (or out of) the monitoring well at the ith node, Q_i, depends linearly on the product b_iK_i (see also Sect. 2.2.2). Numerical evaluation of h_i at a given time t requires evaluating Q_i, as shown in Sect. 2.2.2.

2.2.1 Moments for hydraulic head at observation wells

Mean head at well I is approximated (at second order in σ_Y) as $\langle h_I^{\mathrm{w}}\rangle \approx \langle h_I^{\mathrm{w}\left(\mathrm{0}\right)}\rangle +\langle h_I^{\mathrm{w}\left(\mathrm{2}\right)}\rangle$ where, starting from Eq. (10), one can obtain the zero-, $\langle h_I^{\mathrm{w}\left(\mathrm{0}\right)}\rangle$, and second-, $\langle h_I^{\mathrm{w}\left(\mathrm{2}\right)}\rangle$, order components as

$$\begin{array}{}\text{(11)}& {\displaystyle }\langle h_I^{\mathrm{w}\left(\mathrm{0}\right)}\rangle ={\displaystyle \frac{\sum _{i=\mathrm{1}}^n{T}_{\mathrm{G},i}\langle h_i^{\left(\mathrm{0}\right)}\rangle }{\sum _{i=\mathrm{1}}^n{T}_{\mathrm{G},i}}}\text{(12)}& {\displaystyle }\begin{array}{rl}\langle h_I^{\mathrm{w}\left(\mathrm{2}\right)}\rangle \sum _{i=\mathrm{1}}^n{T}_{\mathrm{G},i}& =\sum _{i=\mathrm{1}}^n{T}_{\mathrm{G},i}\left\{{\displaystyle \frac{u_{ii}}{K_{\mathrm{G},i}}}-{\displaystyle \frac{\langle K_i^{\prime }{h}_I^{{}^{\prime }\mathrm{w}}{\rangle }^{\left(\mathrm{2}\right)}}{K_{\mathrm{G},i}}}\right.\\ & \left.+\langle h_i^{\left(\mathrm{2}\right)}\rangle +{\displaystyle \frac{{\mathit{\sigma }}_{Y,i}^{\mathrm{2}}}{\mathrm{2}}}\left(\langle h_i^{\left(\mathrm{0}\right)}\rangle -\langle h_I^{\mathrm{w}\left(\mathrm{0}\right)}\rangle \right)\right\}.\end{array}\end{array}$$

Here, $T_{\mathrm{G},i}=b_i{K}_{\mathrm{G},i}$, $\langle h_i^{\left(\mathrm{0}\right)}\rangle$, $\langle h_i^{\left(\mathrm{2}\right)}\rangle$, K_G,i, and ${\mathit{\sigma }}_{Y,i}^{\mathrm{2}}$ correspond to the zero- (evaluated by Eq. 6) and the second- (evaluated by Eq. 7) order mean head, geometric mean of conductivity, and variance of log-conductivity at the ith cell of a multi-node monitoring well, respectively; $u_{ii}=\langle K_i^{\prime }{h}_i^{\prime }{\rangle }^{\left(\mathrm{2}\right)}$ is the cross-covariance between conductivity and head at the ith cell; $\langle K_i^{\prime }{h}_I^{{}^{\prime }\mathrm{w}}{\rangle }^{\left(\mathrm{2}\right)}$ is the cross-covariance between well head and conductivity at the ith cell (evaluated as detailed below; see Eq. 15).

The covariance between water levels at wells I (i.e., $h_I^{\mathrm{w}}$) and J (i.e., $h_J^{\mathrm{w}}$), $C_{{\mathrm{h}}^{\mathrm{w}}}=\langle h_I^{{}^{\prime }\mathrm{w}}{h}_J^{{}^{\prime }\mathrm{w}}{\rangle }^{\left(\mathrm{2}\right)}$, can be evaluated as (see also Appendix A, Eqs. A1–A3)

$$\begin{array}{}\text{(13)}& \begin{array}{rl}& C_{{\mathrm{h}}^{\mathrm{w}}}\sum _{i=\mathrm{1}}^n{T}_{\mathrm{G},i}\sum _{j=i}^m{T}_{\mathrm{G},j}=\sum _{j=\mathrm{1}}^m\sum _{i=\mathrm{1}}^n{T}_{\mathrm{G},i}{T}_{\mathrm{G},j}\\ & \left\{\begin{array}{l}\left(\langle h_i^{\left(\mathrm{0}\right)}\rangle -\langle h_I^{\mathrm{w}\left(\mathrm{0}\right)}\rangle \right)\left(\langle h_j^{\left(\mathrm{2}\right)}\rangle -\langle h_J^{\mathrm{w}\left(\mathrm{2}\right)}\rangle +\frac{u_{jj}}{K_{\mathrm{G},j}}+\frac{u_{ij}}{K_{\mathrm{G},i}}\right)\\ +\left(\langle h_j^{\left(\mathrm{0}\right)}\rangle -\langle h_J^{\mathrm{w}\left(\mathrm{0}\right)}\rangle \right)\left(\langle h_i^{\left(\mathrm{2}\right)}\rangle -\langle h_I^{\mathrm{w}\left(\mathrm{2}\right)}\rangle +\frac{u_{ji}}{K_{\mathrm{G},j}}+\frac{u_{ii}}{K_{\mathrm{G},i}}\right)\\ +\left(\frac{{\mathit{\sigma }}_{Y,i}^{\mathrm{2}}}{\mathrm{2}}+\frac{{\mathit{\sigma }}_{Y,j}^{\mathrm{2}}}{\mathrm{2}}+C_{Y,ij}\right)\left(\langle h_i^{\left(\mathrm{0}\right)}\rangle -\langle h_I^{\mathrm{w}\left(\mathrm{0}\right)}\rangle \right)\left(\langle h_j^{\left(\mathrm{0}\right)}\rangle -\langle h_J^{\mathrm{w}\left(\mathrm{0}\right)}\rangle \right)\\ +C_{\mathrm{h},ij}\end{array}\right\},\end{array}\end{array}$$

where $T_{\mathrm{G},i}=b_i{K}_{\mathrm{G},i}$, $T_{\mathrm{G},j}=b_j{K}_{\mathrm{G},j}$, I, and J are indices ranging from 1 to N_w (i.e., the total number of monitoring wells of Type B and C); n and m correspond to the total number of cells according to which the screens of boreholes I and J are discretized, respectively; C_Y,ij is the log-conductivity covariance between the ith and jth cells of boreholes I and J, respectively; u_ji (or u_ij) is the cross-covariance between the conductivity of the jth cell of well J (or the ith cell of well I) and head of the ith cell of well I (or the jth cell of well J); and C_h,ij is the head covariance between the ith cell of well I and the jth cell of well J. Terms u_ji (or u_ij) and C_h,ij can be readily obtained by solving Eqs. (8) and (9).

The cross-covariance between $h_I^{\mathrm{w}}$ at a given time t and aquifer head h at time τ (at a given location; omitted for brevity), i.e., $C_{{\mathrm{h}}_I^{\mathrm{w}}{\mathrm{h}}^{\mathit{\tau }}}=\langle h_I^{{}^{\prime }\mathrm{w}}\left(t\right)h^{\prime }\left(\mathit{\tau }\right){\rangle }^{\left(\mathrm{2}\right)}$ is given by

$$\begin{array}{}\text{(14)}& C_{{\mathrm{h}}_I^{\mathrm{w}}{\mathrm{h}}^{\mathit{\tau }}}\sum _{i=\mathrm{1}}^n{T}_{\mathrm{G},i}=\sum _{i=\mathrm{1}}^n{T}_{\mathrm{G},i}\left[C_{\mathrm{h},i,\mathit{\tau }}+{\displaystyle \frac{u_{i,\mathit{\tau }}}{K_{\mathrm{G},i}}}\left(\langle h_i^{\left(\mathrm{0}\right)}\rangle -\langle h_I^{\mathrm{w}\left(\mathrm{0}\right)}\rangle \right)\right].\end{array}$$

Here, $u_{i,\mathit{\tau }}=\langle K_i^{\prime }{h}^{\prime }\left(\mathit{\tau }\right){\rangle }^{\left(\mathrm{2}\right)}$ represents the cross-covariance between conductivity at the ith cell of well I and aquifer head (at a given location) at time τ. Likewise, $C_{\mathrm{h},i,\mathit{\tau }}=\langle h_i^{\prime }\left(t\right)h^{\prime }\left(\mathit{\tau }\right){\rangle }^{\left(\mathrm{2}\right)}$ represents head covariance between head at time t at the ith cell of the monitoring well I and head at time τ at a given location in the aquifer.

The cross-covariance between $h_I^{\mathrm{w}}$ at time t and conductivity at a given location in the aquifer, $C_{{\mathrm{h}}_I^{\mathrm{w}}K}=\langle h_I^{{}^{\prime }\mathrm{w}}\left(t\right)K^{\prime }{\rangle }^{\left(\mathrm{2}\right)}$, can be expressed as

$$\begin{array}{}\text{(15)}& C_{{\mathrm{h}}_I^{\mathrm{w}}K}\sum _{i=\mathrm{1}}^n{T}_{\mathrm{G},i}=\sum _{i=\mathrm{1}}^n{T}_{\mathrm{G},i}\left[u_i+K_{\mathrm{G}}{C}_{Y,i}\left(\langle h_i^{\left(\mathrm{0}\right)}\rangle -\langle h_I^{\mathrm{w}\left(\mathrm{0}\right)}\rangle \right)\right],\end{array}$$

where $C_{Y,i}=\langle Y_i^{\prime }{Y}^{\prime \prime }\rangle$ is the covariance between log-conductivity at the ith cell along the monitoring borehole I and log-conductivity at a given location in the domain, and $u_i=\langle K^{\prime }{h}_i^{\prime }\left(t\right){\rangle }^{\left(\mathrm{2}\right)}$ is cross-covariance between conductivity at a given point in the domain and aquifer head at time t at the ith cell along the monitoring borehole I.

It is worthwhile to note that covariances and cross-covariances evaluated in Eqs. (13)–(15) depend explicitly on the difference between the mean water level at the monitoring well and the mean hydraulic head along the well screen.

2.2.2 Moments for flux between a monitoring well and the aquifer

Assuming that the evolution of head at the observation borehole I can be conceptualized as a sequence of temporal events (each associated with the attainment of instantaneous equilibrium conditions), following Konikow et al. (2009), the link between $h_I^{\mathrm{w}}$, h_i, and Q_i can then be obtained by relying on the Thiem (1906) formulation as

$$\begin{array}{}\text{(16)}& h_I^{\mathrm{w}}=h_i+{\displaystyle \frac{Q_i}{a{b}_i{K}_i}},\mathrm{with}a={\displaystyle \frac{\mathrm{2}\mathit{\pi }}{\ln \left(r_{\mathrm{0}}/r_{\mathrm{w}}\right)}},\end{array}$$

where r₀ and r_w are the effective (i.e., the radius of a well that would be associated with the same head as that calculated at the node of the cell that contains the well) and the actual well radius, respectively. The mean flux exchange is approximated as $\langle Q_i\rangle \approx \langle Q_i^{\left(\mathrm{0}\right)}\rangle +\langle Q_i^{\left(\mathrm{2}\right)}\rangle$ and from Eq. (16) one can write

$$\begin{array}{}\text{(17)}& {\displaystyle }\langle Q_i^{\left(\mathrm{0}\right)}\rangle =a{T}_{\mathrm{G},i}\left(\langle h_I^{\mathrm{w}\left(\mathrm{0}\right)}\rangle -\langle h_i^{\left(\mathrm{0}\right)}\rangle \right);\text{(18)}& {\displaystyle }\begin{array}{rl}\langle Q_i^{\left(\mathrm{2}\right)}\rangle & =a{T}_{\mathrm{G},i}\left\{\langle h_I^{\mathrm{w}\left(\mathrm{2}\right)}\rangle -\langle h_i^{\left(\mathrm{2}\right)}\rangle +{\displaystyle \frac{{\mathit{\sigma }}_{Y,i}^{\mathrm{2}}}{\mathrm{2}}}\left(\langle h_I^{\mathrm{w}\left(\mathrm{0}\right)}\rangle \right.\right.\\ & \left.\left.-\langle h_i^{\left(\mathrm{0}\right)}\rangle \right)+{\displaystyle \frac{\langle h_I^{{}^{\prime }\mathrm{w}}{K}_i^{\prime }{\rangle }^{\left(\mathrm{2}\right)}}{K_{\mathrm{G},i}}}-{\displaystyle \frac{\langle h_i^{\prime }{K}_i^{\prime }{\rangle }^{\left(\mathrm{2}\right)}}{K_{\mathrm{G},i}}}\right\}.\end{array}\end{array}$$

The cross-covariance between Q_i at time t and aquifer head at time τ, $C_{Q_i{\mathrm{h}}^{\mathit{\tau }}}=\langle Q_i^{\prime }\left(t\right)h^{\prime }\left(\mathit{\tau }\right){\rangle }^{\left(\mathrm{2}\right)}$, can be expressed as

$$\begin{array}{}\text{(19)}& C_{Q_i{\mathrm{h}}^{\mathit{\tau }}}=a{T}_{\mathrm{G},i}\left\{C_{{\mathrm{h}}_I^{\mathrm{w}}{\mathrm{h}}^{\mathit{\tau }}}-C_{{\mathrm{h}}_i{\mathrm{h}}^{\mathit{\tau }}}+{\displaystyle \frac{u_{i,\mathit{\tau }}}{K_{\mathrm{G},i}}}\left(\langle h_I^{\mathrm{w}\left(\mathrm{0}\right)}\rangle -\langle h_i^{\left(\mathrm{0}\right)}\rangle \right)\right\},\end{array}$$

where $C_{{\mathrm{h}}_I^{\mathrm{w}}{\mathrm{h}}^{\mathit{\tau }}}=\langle h_I^{{}^{\prime }\mathrm{w}}\left(t\right)h^{\prime }\left(\mathit{\tau }\right){\rangle }^{\left(\mathrm{2}\right)}$, $C_{{\mathrm{h}}_i{\mathrm{h}}^{\mathit{\tau }}}=\langle h_i^{\prime }\left(t\right)h^{\prime }\left(\mathit{\tau }\right){\rangle }^{\left(\mathrm{2}\right)}$. Finally, one can obtain the following expression for the cross-covariance between Q_i and K, $C_{Q_{i}K}=\langle Q_{i^{\prime }}^{\prime }{K}^{\prime }{\rangle }^{\left(\mathrm{2}\right)}$,

$$\begin{array}{}\text{(20)}& C_{Q_{i}K}=a{T}_{\mathrm{G},i}\left\{C_{{\mathrm{h}}_I^{\mathrm{w}}K}-u_i+K_{\mathrm{G}}{C}_{Y,i}\left(\langle h_I^{\mathrm{w}\left(\mathrm{0}\right)}\rangle -\langle h_i^{\left(\mathrm{0}\right)}\rangle \right)\right\}.\end{array}$$

2.3 Numerical solution strategy

We solve Eqs. (6)–(9) numerically by approximating the spatial derivatives through a finite element approach and the temporal derivatives through an implicit method. As in Xia et al. (2019), moments 〈h⁽⁰⁾〉, u, and 〈h⁽²⁾〉 are sequentially obtained by solving Eqs. (6)–(8), respectively. Details associated with the evaluation of C_h, which requires 〈h⁽⁰⁾〉 and u as inputs, are illustrated in the following.

For the purpose of our data assimilation workflow, we start by noting that we are interested in computing C_h associated with two identical time coordinates, i.e., $C_{\mathrm{h}}(\mathit{y},\mathit{x},\mathit{\tau }=t,t)=\langle h^{\prime }(\mathit{y},\mathit{\tau }=t)h^{\prime }(\mathit{x},t){\rangle }^{\left(\mathrm{2}\right)}$. We then recall that Zhang (2002) computes $C_{\mathrm{h}}(\mathit{y},\mathit{x},\mathit{\tau }=t,t)$ for each time t (while $C_{\mathrm{h}}(\mathit{y},\mathit{x},\mathit{\tau }=t,t-\mathrm{\Delta }t)$ is also unknown, Δt being a constant temporal step size) by solving for $C_{\mathrm{h}}(\mathit{y},\mathit{x},\mathit{\tau }=t,t^{\prime })$ from $t^{\prime }=\mathrm{0}$ to $t^{\prime }=t$. While this procedure can be computationally heavy for long times, Zhang (2002) points out that when flow changes only mildly, $C_{\mathrm{h}}(\mathit{x},\mathit{y},\mathit{\tau }=t,t-\mathrm{\Delta }t)\approx C_{\mathrm{h}}(\mathit{x},\mathit{y},\mathit{\tau }=t-\mathrm{\Delta }t,t-\mathrm{\Delta }t)$, an approximation whose general validity is still not completely explored.

Here, we circumvent this issue and obtain high computational efficiency by directly evaluating $C_{\mathrm{h}}(\mathit{y},\mathit{x},\mathit{\tau }=t,t)$ from $C_{\mathrm{h}}(\mathit{y},\mathit{x},\mathit{\tau }=t-\mathrm{\Delta }t,t-\mathrm{\Delta }t)$ via (i) computing $C_{\mathrm{h}}(\mathit{y},\mathit{x},\mathit{\tau }=t,t-\mathrm{\Delta }t)$ through the solution of the equation obtained by considering Eq. (9) where the space and time derivatives operate on τ and y (instead of t and x) from time t−Δt to t using $C_{\mathrm{h}}(\mathit{y},\mathit{x},\mathit{\tau }=t-\mathrm{\Delta }t,t-\mathrm{\Delta }t)$ as initial condition and then (ii) assessing $C_{\mathrm{h}}(\mathit{y},\mathit{x},\mathit{\tau }=t,t)$ by solving Eq. (9) using $C_{\mathrm{h}}(\mathit{y},\mathit{x},\mathit{\tau }=t,t-\mathrm{\Delta }t)$ as the initial condition.

It is further noted that Eqs. (6)–(9) are characterized by the same format, their discretization leading to a system of equations where the coefficients of the unknown quantities are identical, the corresponding right-hand-side terms (i.e., the forcing terms) being a function of the (ensemble) moment to be solved. In this context, one can resort to a direct solver for each time step. Thus, factorization of the matrix containing the coefficients of the system of equations is performed only once, resulting in a high computational efficiency because only the right-hand-side term needs to be updated, depending on the moment of interest.

With reference to the forcing terms 〈f⁽⁰⁾〉, 〈f⁽²⁾〉, C_fK, and C_fh in Eqs. (6)–(9), we note that these vanish for Type A wells and when one disregards flux exchanges between Type B (or C) wells and the aquifer. In these instances, mean head values and the associated covariance are simply obtained upon evaluating Eqs. (6)—(9)numerically. Thus, when considering a time interval [$t-\mathrm{\Delta }t,t$], the main computational cost stems from the evaluation of $u(\mathit{y},\mathit{x},t)$, $C_{\mathrm{h}}(\mathit{y},\mathit{x},\mathit{\tau }=t,t-\mathrm{\Delta }t)$, and $C_{\mathrm{h}}(\mathit{y},\mathit{x},\mathit{\tau }=t,t)$, each of these requiring N times the computational cost (hereafter denoted as $C_{\mathrm{c}}^{\mathrm{MEs}}$) associated with the solution of the system of N equations resulting after discretization. Therefore, the total computational effort required for solving Eqs. (6)—(9) at each time step is 3 $N{C}_{\mathrm{c}}^{\mathrm{MEs}}$. Note that the computational effort is reduced to 2 $N{C}_{\mathrm{c}}^{\mathrm{MEs}}$ for the first time interval, when the initial head is deterministic, or for a steady-state flow scenario (see Xia et al., 2019).

Otherwise, considering flux exchange processes when representing Type B (or C) wells entails evaluation of the source terms in Eqs. (6)–(9) as $\langle f^{\left(\mathrm{0}\right)}\rangle =\langle Q_i^{\left(\mathrm{0}\right)}\rangle$, $\langle f^{\left(\mathrm{2}\right)}\rangle =\langle Q_i^{\left(\mathrm{2}\right)}\rangle$, $C_{\mathrm{f}K}=C_{Q_{i}K}$, and $C_{\mathrm{fh}}=C_{Q_i{\mathrm{h}}^{\mathit{\tau }}}$. The evaluation of the (ensemble) moments of interest across time interval [$t-\mathrm{\Delta }t,t$] is then performed through the workflow depicted in Fig. 2. In this case, we note that convergence of the iterative procedure is attained when the absolute difference between mean well heads at iteration iter + 1, $\langle h_I^{\mathrm{w}}{\rangle }_{\mathrm{iter}+\mathrm{1}}$, and iter, $\langle h_I^{\mathrm{w}}{\rangle }_{\mathrm{iter}}$, is lower than a preset value ε. The main computational effort required for these evaluations corresponds to 3 (iter + 1)$N{C}_{\mathrm{c}}^{\mathrm{MEs}}$ for each time step.

Figure 2 Workflow for the numerical solution of MEs within time interval [t–Δt, t] when flux between monitoring wells and the aquifer is considered.

3 Ensemble Kalman filter coupled with moment equations

We start by introducing the mean system state vector 〈φ〉 as

$$\begin{array}{}\text{(21)}& \langle \mathit{\phi }\rangle ={\left[\begin{array}{lll}\langle \mathit{h}\rangle & \langle {\mathit{h}}^{\mathrm{w}}\rangle & \langle \mathit{Y}\rangle \end{array}\right]}^{\mathrm{T}},\end{array}$$

where 〈h〉, 〈Y〉 correspond, respectively, to N-dimensional vectors of mean head and mean log-conductivity, and 〈h^w〉 is a N_w-dimensional mean well head vector, subscript “T” representing transpose.

Each data assimilation cycle, corresponding to time interval [$t-\mathrm{\Delta }t,t$] comprises a forecast (or forward propagation) step and an update (or analysis) step. The forecast step is implemented by solving the moment equations described in Sect. 2. We write the predicted mean and covariance of the system state as

$$\begin{array}{}\text{(22)}& \begin{array}{rl}& \langle \mathit{\phi }{\rangle }^{\mathrm{f}}={\left[\begin{array}{lll}\langle \mathit{h}{\rangle }^{\mathrm{f}}& \langle {\mathit{h}}^{\mathrm{w}}{\rangle }^{\mathrm{f}}& \langle \mathit{Y}\rangle \end{array}\right]}^{\mathrm{T}};\\ & {\mathbf{P}}^{\mathrm{f}}=\left[\begin{array}{lll}{\mathbf{C}}_{\mathrm{h}}^{\mathrm{f}}& {\left[{\mathbf{C}}_{{\mathrm{h}}^{\mathrm{w}}{\mathrm{h}}^{\mathrm{f}}}\right]}^{\mathrm{T}}& {\left[{\mathbf{C}}_{Y\mathrm{h}}^{\mathrm{f}}\right]}^{\mathrm{T}}\\ {\mathbf{C}}_{{\mathrm{h}}^{\mathrm{w}}\mathrm{h}}^{\mathrm{f}}& {\mathbf{C}}_{{\mathrm{h}}^{\mathrm{w}}}^{\mathrm{f}}& {\left[{\mathbf{C}}_{Y{\mathrm{h}}^{\mathrm{w}}}^{\mathrm{f}}\right]}^{\mathrm{T}}\\ {\mathbf{C}}_{Y\mathrm{h}}^{\mathrm{f}}& {\mathbf{C}}_{Y{\mathrm{h}}^{\mathrm{w}}}^{\mathrm{f}}& {\mathbf{C}}_{\mathrm{Y}}^{\mathrm{f}}\end{array}\right].\end{array}\end{array}$$

Here, the superscript “f” represents predicted quantities obtained in the forecast step, 〈h〉^f is the predicted mean head (Eq. 5), 〈h^w〉^f is the predicted mean water level at monitoring borehole (Eqs. 11 and 12), 〈Y〉 is the updated natural logarithm of conductivity obtained at the previous data assimilation cycle, ${\mathbf{C}}_{\mathrm{h}}^{\mathrm{f}}$ is the predicted N×N-dimensional head covariance matrix, (Eq. 9), ${\mathbf{C}}_{{\mathrm{h}}^{\mathrm{w}}\mathrm{h}}^{\mathrm{f}}$ is the N_w×N-dimensional predicted cross-covariance between well and aquifer head (Eq. 14), ${\mathbf{C}}_{{\mathrm{h}}^{\mathrm{w}}}^{\mathrm{f}}$ is the predicted N_w×N_w-dimensional covariance of well head (Eq. 13), ${\mathbf{C}}_{Y\mathrm{h}}^{\mathrm{f}}$ is the predicted N×N-dimensional cross-covariance between Y and aquifer head h (Eq. 8), ${\mathbf{C}}_{Y{\mathrm{h}}^{\mathrm{w}}}^{\mathrm{f}}$ is the predicted N×N_w-dimensional cross-covariance between Y and well head h^w (Eq. 15), and ${\mathbf{C}}_Y^{\mathrm{f}}$ is N×N-dimensional Y covariance and is equal to its updated counterpart associated with the previous updating step.

The equations used to evaluate the state updated vector 〈φ〉^up and the updated covariance matrix P^up are

$$\begin{array}{}\text{(23)}& \langle \mathit{\phi }{\rangle }^{\mathrm{up}}=\langle \mathit{\phi }{\rangle }^{\mathrm{f}}+{\mathbf{C}}_{\mathit{\phi }\mathrm{d}}{\left({\mathbf{C}}_{\mathrm{dd}}+\mathit{\alpha }{\mathbf{C}}_{\mathrm{D}}\right)}^{-\mathrm{1}}\left(\langle {\mathit{\phi }}_{\mathrm{so}}\rangle -{\mathit{d}}_{\mathrm{obs}}\right)\end{array}$$

and

$$\begin{array}{}\text{(24)}& {\mathbf{P}}^{\mathrm{up}}=\left(\mathbf{I}-{\mathbf{C}}_{\mathit{\phi }\mathrm{d}}{\left({\mathbf{C}}_{\mathrm{dd}}+\mathit{\alpha }{\mathbf{C}}_{\mathrm{D}}\right)}^{-\mathrm{1}}\mathbf{H}\right){\mathbf{P}}^{\mathrm{f}},\end{array}$$

where C_φd (=P^fH^T) is the $(\mathrm{2}N+N_{\mathrm{w}})\times d$-dimensional cross-covariance between the system state and the simulated observations, matrix H of dimension $d\times (\mathrm{2}N+N_{\mathrm{w}})$ is the observation operator that describes the relationship between the system state and the observations, C_dd (=HP^fH^T) is the d×d-dimensional covariance matrix of the simulated observations, C_D  is the d×d-dimensional covariance of observation errors, I is the identity matrix, α is a constant inflation factor, α=1 corresponding to the uninflated ensemble Kalman filter, 〈φ_so〉 is the mean vector of the simulated observations, and d_obs is the d-dimensional observation vector.

After the update step, 〈φ〉^up and P^up are expressed as

$$\begin{array}{}\text{(25)}& \begin{array}{rl}& \langle \mathit{\phi }{\rangle }^{\mathrm{up}}={\left[\begin{array}{lll}\langle \mathit{h}{\rangle }^{\mathrm{up}}& \langle {\mathit{h}}^{\mathrm{w}}{\rangle }^{\mathrm{up}}& \langle \mathbf{Y}{\rangle }^{\mathrm{up}}\end{array}\right]}^{\mathrm{T}};\\ & {\mathbf{P}}^{\mathrm{up}}=\left[\begin{array}{lll}{\mathbf{C}}_{\mathrm{h}}^{\mathrm{up}}& {\left[{\mathbf{C}}_{{\mathrm{h}}^{\mathrm{w}}\mathrm{h}}^{\mathrm{up}}\right]}^{\mathrm{T}}& {\left[{\mathbf{C}}_{Y\mathrm{h}}^{\mathrm{up}}\right]}^{\mathrm{T}}\\ {\mathbf{C}}_{{\mathrm{h}}^{\mathrm{w}}\mathrm{h}}^{\mathrm{up}}& {\mathbf{C}}_{{\mathrm{h}}^{\mathrm{w}}}^{\mathrm{up}}& {\left[{\mathbf{C}}_{Y{\mathrm{h}}^{\mathrm{w}}}^{\mathrm{up}}\right]}^{\mathrm{T}}\\ {\mathbf{C}}_{Y\mathrm{h}}^{\mathrm{up}}& {\mathbf{C}}_{Y{\mathrm{h}}^{\mathrm{w}}}^{\mathrm{up}}& {\mathbf{C}}_Y^{\mathrm{up}}\end{array}\right],\end{array}\end{array}$$

where all symbols have the same meaning (yet updated) as in Eq. (22).

When moving to a subsequent time interval during the assimilation process, we follow Panzeri et al. (2013) and (i) use the updated mean head vector 〈h〉^up as the initial condition of the governing equation for the zero-order mean head, i.e., Eq. (6); (ii) making use of 〈Y〉^up, evaluate the updated geometric mean N-vector ${\mathbf{K}}_{\mathrm{G}}^{\mathrm{up}}$; (iii) obtain the initial condition of Eq. (8) through the product ${\mathbf{K}}_{\mathrm{G}}^{\mathrm{up}}{\mathbf{C}}_{Y\mathrm{h}}^{\mathrm{up}}$; (iv) use ${\mathbf{C}}_{\mathrm{h}}^{\mathrm{up}}$ as the initial condition of Eq. (9); and (v) use ${\mathbf{K}}_{\mathrm{G}}^{\mathrm{up}}$ and ${\mathbf{C}}_Y^{\mathrm{up}}$ as inputs to Eqs. (6)–(9) and Eqs. (11)–(20).

It should be noted that if one neglects flux exchanges between the aquifer and Type B and/or C monitoring wells (or a Type A well is considered), moments including water level at well (i.e., 〈h^w〉, ${\mathbf{C}}_{{\mathrm{h}}^{\mathrm{w}}\mathrm{h}}$, ${\mathbf{C}}_{{\mathrm{h}}^{\mathrm{w}}}$, ${\mathbf{C}}_{Y{\mathrm{h}}^{\mathrm{w}}}$) should be omitted in Eqs. (21)–(25).

4 Illustrative examples

We consider a three-dimensional domain (Fig. 3a) of size $\mathrm{600}\times \mathrm{600}\times \mathrm{60}$ (hereafter, all quantities are considered in consistent units, following notation associated with prior studies, including, e.g., Panzeri et al. (2013) and references therein), the system being discretized onto a numerical mesh comprising $\mathrm{25}\times \mathrm{25}\times \mathrm{13}$ nodes, for a total of 34 560 tetrahedrons. A partially penetrating pumping well pumps at a constant rate of 1000 for $\mathrm{0}\le t\le \mathrm{30}$, after which water withdrawal stops and a recovering process takes place for $\mathrm{30}<t\le \mathrm{60}$. We subdivide the overall simulation time according to 20 uniform intervals, which can potentially be used for assimilation of head observations. The well pumping rate is uniformly distributed across the central nodes of layers no. 1 and 2 in the numerical mesh (numbering is from top to bottom of the domain). The left and right sides of the system are set as Dirichlet boundaries, where a deterministic head H=100 is fixed, the remaining boundaries being considered impervious. Initial head and storativity are deterministic and set equal to 100 and 10⁻³, respectively. The natural logarithm of conductivity, Y, is modeled as a spatially correlated second-order stationary random field with covariance given by

$$\begin{array}{}\text{(26)}& C_Y={\mathit{\sigma }}_Y^{\mathrm{2}}\exp \left(-\left[{\displaystyle \frac{{\mathit{\delta }}_{\mathrm{1}}}{{\mathit{\lambda }}_{\mathrm{1}}}}+{\displaystyle \frac{{\mathit{\delta }}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{2}}}}+{\displaystyle \frac{{\mathit{\delta }}_{\mathrm{3}}}{{\mathit{\lambda }}_{\mathrm{3}}}}\right]\right).\end{array}$$

Here, δ_i and λ_i denote, respectively, the lag and correlation scale between two points along direction x_i (with i=1, 2, 3).

Figure 3 Reference log-conductivity (Y) field (a) and location of monitoring wells within the flow domain (b).

Twenty virtual monitoring wells are regularly distributed across the domain (Fig. 3b). Type A boreholes are mimicked by considering three packers positioned at three distinct depths, corresponding to layers no. 4, 7, and 10. Type B wells are equipped with three screens (i.e., n=3) whose barycenter is set at the same depths of the packers in Type A wells (see Fig. 1) and b_i=5 with i=1, 2, 3. Type C wells are completely penetrating across the 13 layers of the domain (i.e., n=13) with $b_{\mathrm{1}}=b_{\mathrm{13}}=\mathrm{2.5}$ and b_i=5 for i=2, …, 12.

Reference hydraulic head values that are collected at Type B and C wells and employed in the data assimilation procedure are evaluated upon solving the flow problem on the reference hydraulic conductivity fields described in the following. Flux exchanges between the aquifer and monitoring wells are evaluated according to the procedure described in Sect. 2.2 and 2.3 setting the convergence criterion $\mathit{\epsilon }={\mathrm{10}}^{-\mathrm{6}}$.

The effective radius of the monitoring wells is evaluated as (Chen and Zhang, 2009) $r_{\mathrm{0}}={\mathrm{2}}^{\mathrm{0.25}}{e}^{-\mathrm{0.75}\mathit{\pi }}\mathrm{\Delta }x\approx \mathrm{0.113}\mathrm{\Delta }x=\mathrm{2.81}$ (Δx=25 being the horizontal size of a given element in the computational mesh). For the purpose of our illustrative example, we set r_w=0.1 (i.e., a=1.88).

We organize our exemplary settings according to the following four groups (for a total of 26 test cases, TCs) collected in Table 1.

Table 1 Summary of the test cases analyzed.

1. Group 1. This group includes seven TCs (TC1–TC7) that allow exploring the way conductivity estimates can be influenced by relying on the assimilation of head data collected at diverse types of virtual observation boreholes, while considering a simplified modeling approach where flux exchanges between the aquifer and Type B (or C) monitoring wells are neglected during the data assimilation procedure, head observations (considered in the data assimilation procedure) corresponding to depth-averaged values along the corresponding screens. We note that relying on this approach is tantamount to considering an imperfect flow model and would possibly oversimplify the mathematical representation of the system behavior when compared with the one employed for constructing the reference head field. Nevertheless, it has the advantage of requiring a straightforward numerical implementation.

   A zero-mean reference Y field is generated at the nodes of the computational mesh upon relying on the widely tested and used SGSIM code (e.g., Deutsch and Journel, 1998) by setting a unit variance, ${\mathit{\lambda }}_{\mathrm{1}}={\mathit{\lambda }}_{\mathrm{2}}=\mathrm{100}$ and λ₃=20. While the correlation scale values are considered as perfectly known, we aim to estimate the mean and variance of Y. The initial guesses employed for the variance and mean of Y during data assimilation are 1.0 and 0.2, respectively. Test cases 1–3 are designed in a way that all 20 observation boreholes are of Type A, B, or C, respectively. Test cases 4–7 enable one to explore the way conductivity estimates can depend on the use of Type A boreholes (i.e., equipped with packers) within different zones, while considering Type B or C wells in the remaining regions. Test case 4 comprises Type A wells within zone 1 in Fig. 3b (i.e., within distances shorter than λ₁ from the pumping well) and Type B boreholes being installed in zones 2 and 3 (at distances larger than λ₁ from the pumping well). Test case 5 is characterized by the presence of Type A wells in zones 1 and 2, and Type B wells in zone 3. Test cases 6 and 7 are characterized by the presence of Type A wells in zone 1 and in zones 1–2 respectively, and Type C wells being located in the remaining zones. We note that the spatial arrangement of the observation boreholes is designed to allow these to be spaced by a distance approximately corresponding to a correlation scale of Y, thus encompassing strong to low degrees of correlation with respect to the pumping well location.

2. Group 2. This group includes six TCs (TC#2–TC#7) that are variants of those of Group 1 and consider the solution of the data assimilation procedure without neglecting flux exchanges between virtual monitoring boreholes of Type B and/or C and the aquifer, i.e., data assimilation is performed by considering perfect knowledge of the groundwater flow model, which includes all of the processes underpinning the reference head fields.

3. Group 3. This group includes seven TCs designed to explore (i) the impacts of the mean and variance of the Y reference field on the data assimilation results associated with Type B and C boreholes (TC2#^*1–TC2#^*2, TC3#^*1–TC3#^*2) and (ii) settings where head data are assimilated solely from one depth (instead of all three locations) where a packer is installed along Type A wells (TC1^*1–TC1^*3).

   To elaborate, here we consider (i) a nearly uniform (while random) zero-mean Y reference field with variance equal to 0.01 (TC2#^*1 and TC3#^*1), the initial guesses employed for the variance and mean of Y during data assimilation being 0.09 and 0.2, respectively; (ii) a Y reference field with mean and variance equal to 0.2 and 1.70, respectively (TC2#^*2 and TC3#^*2), the initial guesses employed for the variance and mean of Y during data assimilation being 1 and 0, respectively; and (iii) three variants of TC1: TC1^*1 considers assimilating head information only from the upper packer (i.e., the one positioned at layer 4), and TC1^*2 and TC1^*3 are designed to assimilate head data only from the intermediate (positioned at layer 7) and bottom (positioned at layer 13) packer, respectively.

4. Group 4. This group includes six TCs where we explore the effect of inflating the measurement-error covariance matrix on the data assimilation when the latter is performed in a way similar to the corresponding TC2 and TC3 of Group 1. As such, data assimilation is based on an imperfect flow model (where flux exchanges between the aquifer and monitoring boreholes are disregarded). To cope with this, inflation on measurement-error covariance matrix is considered during data assimilation, the inflation factor being set to α=5 (TC2α₁ and TC3α₁), 10 (TC2α₂ and TC3α₂), and 100 (TC2α₃ and TC3α₃; note that TC2 and TC3 correspond to α=1).

Initial input quantities required to solve moment equations and spatial fields of K_G and C_Y are obtained through the generation of 10 000 realizations of Y. The latter form the collection of realizations upon which the traditional Monte Carlo (MC)-based ensemble Kalman filter (MC-EnKF) is also applied. Results based on the MEs-EnKF are then compared with those obtained through the MC-EnKF. Head observations in all TCs are considered to be noisy and are obtained by adding a Gaussian white noise with a standard deviation of 0.01 to the reference heads collected at the virtual boreholes and used in the data assimilation procedure. The strength of the noise is selected on the basis of the calculated reference head fields and considering the level of accuracy that is related to measuring devices commonly employed in practical settings (e.g., when considering water loggers, accuracy of pressure head observations is commonly comprised between $\sim \pm \mathrm{0.005}$ and $\sim \pm \mathrm{0.05}$ m).

We rely on the criteria illustrated in the following to appraise the quality of the data assimilation performance. These are (i) the average absolute difference between the estimated (or updated) Y field and its reference counterpart, E_Y; (ii) the square root of the average estimation variance, S_Y; and (iii) the average absolute difference between the estimated (or updated) aquifer head and its reference counterpart, E_h, evaluated as

$$\begin{array}{}\text{(27)}& {\displaystyle }& {\displaystyle }{E}_Y={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N\left|\langle Y_i{\rangle }^{\mathrm{up}}-Y_i^{\mathrm{r}}\right|\text{(28)}& {\displaystyle }& {\displaystyle }{S}_Y=\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N{\left({\mathit{\sigma }}_{Y,i}^{\mathrm{2}}\right)}^{\mathrm{up}}}\text{(29)}& {\displaystyle }& {\displaystyle }{E}_{\mathrm{h}}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N\left|\langle h_i{\rangle }^{\mathrm{up}}-h_i^{\mathrm{r}}\right|,\end{array}$$

where 〈Y_i〉^up, $({\mathit{\sigma }}_{Y,i}^{\mathrm{2}}{)}^{\mathrm{up}}$, and $Y_i^{\mathrm{r}}$ indicate the estimated mean, variance, and reference Y values at the ith node of the computational mesh, respectively; 〈h_i〉^up and $h_i^{\mathrm{r}}$ represent the estimated mean and reference aquifer head value at the ith node of the grid. We note that S_Y is a metric quantifying the uncertainty associated with the estimated Y field conditional on the data assimilated (see, e.g., Panzeri et al., 2014; Nowak, 2010).

Table 2 Comparison of results obtained by the MEs-EnKF and MC-EnKF (based on 100, 500, 1000, and 10 000 realizations) for TC1 and TC2#.

5 Results and discussion

5.1 Comparison between the MC-EnKF and MEs-EnKF

In this section we compare the results obtained with our MEs-EnKF approach and a standard MC-EnKF for two selected test cases, TC1 and TC2#. Table 2 summarizes the outcomes computed via the MEs-EnKF and MC-EnKF (increasing the number of MC simulations from 100 to 10 000) at the end of the assimilation process in terms of E_Y, S_Y, and E_h. These results suggest that the overall quality of conductivity estimates grounded on the MEs-EnKF is similar to what one can obtain upon relying on a MC-EnKF based on 10 000 realizations, which is also consistent with the results illustrated by Panzeri et al. (2014) in two-dimensional settings. Table 2 also includes the computational cost (CPU in seconds) needed for each test case and approach using the processor Inter Xeon(R) CPU E5-2650 v3 @ 2.30 GHz with 128 GB RAM. The CPU time required by the MEs-EnKF is 20 times lower than the one required by a standard Monte Carlo-based EnKF relying on 10 000 realizations. When compared with the findings of Panzeri et al. (2014) in their two-dimensional settings, our results further support the computational appeal and feasibility of relying on a MEs-EnKF approach also in a three-dimensional setting. Note also that the CPU time required by TC2# is significantly larger (about six times) than the one needed for TC1, due to the implementation of flux exchanges between the aquifer system and the boreholes.

As an additional term of comparison, Fig. 4 depicts the spatial distributions of the estimated values of the mean and variance of the log-conductivity field computed with the MEs-EnKF and MC-EnKF relying on 100, 500, 1000, and 10 000 realizations at the end of the data assimilation window at layers 4, 7, and 10 (where the packers are located) in TC1. The reference Y field is also depicted (see the left column of Fig. 4). Analogous outcomes are reported in Fig. 5 for TC2#. Visual inspection of these results provides further support to the ability of the MEs-EnKF to yield conductivity distributions consistent with the corresponding reference values in these settings. As expected, the degree of spatial variability of the estimated mean and variance values of Y tends to stabilize as the number of realizations increases, being very similar to their MEs-based counterparts when 10 000 realizations are employed.

Figure 4 Reference Y field (left column) and estimates of the mean and variance of Y at the end of the assimilation process across layers 4, 7, and 10 for TC1. Results computed via the MC-EnKF (with 100, 500, 1000, and 10 000 realizations) and MEs-EnKF are included.

Figure 5 Reference Y field (left column) and estimates of the mean and variance of Y at the end of the assimilation process at layers 4, 7, and 10 for TC2#. Results computed via the MC-EnKF (with 100, 500, 1000, and 10 000 realizations) and MEs-EnKF are included.

5.2 Effect of neglecting flux exchanges between boreholes and aquifer during data assimilation (Group 1)

Figure 6 shows the temporal behavior of E_Y (Fig. 6a), S_Y (Fig. 6b), and E_h (Fig. 6c) for TCs1–7 (i.e., Group 1 in Table 1) obtained through the MEs-EnKF. The lowest values of E_Y are associated with TC1, where packers are set for all observation wells. Very similar results are also obtained for TC5 and TC7, where Type B or C wells are installed only at the farthest locations (i.e., zone 3 in Fig. 3b) from the pumping well, respectively; Type A boreholes are installed within the regions (i.e., zones 1–2 in Fig. 3b) closest to the well.

Figure 6 Temporal evolution of E_Y (a), S_Y (b), and E_h (c) for TC1–TC7.

The highest values of E_Y correspond to TC3, where fully screened monitoring wells (Type C wells) are located in the entire domain. These are closely followed by the results associated with TC2, where observation wells screened across multiple (Type B wells) levels are considered. Considering the trend displayed by the results in Fig. 6a, one can then conclude that relying on point values of head contributes to increasing the overall quality of the data assimilation procedure, as expressed in terms of E_Y, when compared with considering vertically averaged head information while relying on a simplified groundwater flow model, which neglects flux exchanges between screened intervals of boreholes and the surrounding aquifer.

A comparison between the values of E_Y related to TC1 and TCs 5, 7 further suggests that the use of packers at locations far away (in terms of the horizontal correlation scale of Y) from the pumping well does not add additional information with respect to fully or partially screened wells, because vertical variability of head at such locations is modest. The importance of capturing vertical head variations during data assimilation is also manifest when comparing results related to TC4 and TC6. The values of E_Y related to the former are consistently lower than those associated with the latter, a result that is consistent with the use of fully screened (TC6) compared to partially screened (TC4) observation boreholes in most of the domain. The behavior of E_h is very similar to the one displayed by E_Y, thus strengthening the above conclusions.

One can note that the scenarios characterized by a dominance of Type C boreholes (i.e., TC3 and TC6) are characterized by the lowest values of S_Y (Fig. 6b). This result is related to the observation that depth-averaged well head information is here employed during data assimilation. Doing so tends to introduce a corresponding homogenization of the conductivity field resulting from the data assimilation procedure, which is reflected by the lower values of S_Y. As such, while the results associated with S_Y would suggest that the estimation variance associated with Y is low, the overall accuracy, as given in terms of E_Y, is also low when relying mostly on vertically averaged data. Otherwise, temporal values of S_Y are virtually indistinguishable for the other configurations considered.

The lowest E_h values are visually indistinguishable and are related to TC1, TC5, and TC7, which is generally consistent with the behavior of E_Y. The highest E_h values are associated with TC3, where fully screened monitoring wells (Type C wells) are located in the entire domain. These are followed by TC2, where partially screened monitoring wells (Type B wells) are distributed across the domain.

Figure 7 Temporal evolution of E_Y (a), S_Y (b), and E_h (c) for TC1 and TC2#–TC7#. Temporal evolution of the relative difference between E_Y (d), S_Y (e), and E_h (f) evaluated in TC2#–TC7# and their counterparts related to TC2–TC7.

5.3 Effect of including flux exchanges between boreholes and aquifer during data assimilation (Groups 2 and 3)

Figure 7 depicts the temporal behavior of E_Y (Fig. 7a), S_Y (Fig. 7b), and E_h (Fig. 7c) for TC1 and TC2#–TC7# (i.e., Group 2 in Table 1) obtained through the MEs-EnKF. The lowest values of E_Y are again associated with TC1. These are very closely matched by those obtained for TC5# and TC7#, where Type B or C wells are installed only at the farthest locations (i.e., zone 3 in Fig. 3b) from the pumping well. It is worth noting that the values of E_Y in TCs 5# and 7# are virtually identical. The highest values of E_Y correspond to TC3#, where only fully screened monitoring wells are distributed across the domain. These are very closely followed by the results associated with TC2#, where observation wells screened across multiple levels (Type B wells) are considered. Note that the temporal evolution of E_Y in these two TCs (TC2# and TC3#) is almost identical and tends to the same value at the end of the assimilation window. The comparison between the values of E_Y for TC4# and TC6# is also consistent with this finding.

Different from the results of Group 1, the temporal behavior of S_Y (Fig. 7b) is very similar to one displayed by E_Y. These results are in line with (i) the observation that each Type A borehole provides three head observations, while only one well head observation is essentially linked to Type B or C boreholes and (ii) the intuition that constraining the system with an increased number of observations would yield conductivity estimates characterized by an increased accuracy (in terms of lower E_Y and possibly S_Y values).

The lowest E_h values are related to TC1, TC5#, and TC7# and are visually indistinguishable (see Fig. 7c), a finding which is consistent with the behavior of E_Y (see Fig. 7a). The highest E_h values are associated with TC3#, closely followed by TC2#, TC6#, and then TC4#. Otherwise, the values of E_h become virtually indistinguishable.

Figure 7d depicts relative (percentage) differences between E_Y evaluated for TCs2#–7# and TCs2–7, considering the values of TCs2–7 as references (negative values correspond to lower values of E_Y in TCs2#–7# as compared with TCs2–7). Analogous results are reported for S_Y (Fig. 7e) and E_h (Fig. 7f). The largest accuracy improvement of Y estimates, as suggested by Fig. 7d, corresponds to TC3# (where the inclusion of flux exchanges between boreholes and aquifer is particularly relevant, since all monitoring wells fully penetrate the aquifer), followed by TC2#, TC6#, and TC7#, with respect to their counterparts in Group 1.

Figure 8 Temporal evolution of E_Y (a, d), S_Y (b, e), and E_h (c, f) for TC2#^*1 and TC3#^*1 (a–c) and for TC2#^*2 and TC3#^*2 (d–f).

Figure 9 Temporal evolution of E_Y (a), S_Y (b), and E_h (c) for TC1^*1, TC1^*2, TC1^*3, TC1, and TC2#.

Uncertainty associated with conductivity estimates increases in TCs2#–7# as compared with their counterparts in Group 1 (see Fig. 7e, where relative differences of S_Y are all positive), this result being related to the vertical variability of flux exchanges along Type B and C boreholes, which is embedded in Group 2 TCs. The lowest (negative) relative differences for E_h correspond to TC3#, a result that is consistent with the depiction offered in Fig. 7d.

The temporal evolution of E_Y (Fig. 8a and d), S_Y (Fig. 8b and e), and E_h (Fig. 8c and f) for TC2#^*1, TC3#^*1, TC2#^*2, and TC3#^*2 is displayed in Fig. 8. These results show that head observations collected at Type B or C screened wells yield conductivity and head estimates of similar quality (in terms of E_Y and S_Y, or E_h, respectively) for the degrees of heterogeneity analyzed. Mean absolute differences between E_Y values associated with TC2#^*1 and TC3#^*1 is 0.008, while being virtually null when considering TC2#^*2 and TC3#^*2. Results of corresponding quality are also obtained when comparing S_Y and E_h values related to TC2#^*1 and TC3#^*1, or TC2#^*2 and TC3#^*2. We further note that values of E_Y in Fig. 8d (or E_h in Fig. 8f) are always higher than their counterparts depicted in Fig. 8a (or Fig. 8c), consistent with the observation that the accuracy of conductivity (and head) estimates tends to deteriorate with increasing degree of spatial heterogeneity of conductivities.

Figure 9 juxtaposes the temporal variability of E_Y (Fig. 9a), S_Y (Fig. 9b), and E_h (Fig. 9c) values for TCs 1^*1, 1^*2, 1^*3 (see Group 3 in Table 1), 1 (Group 1), and 2# (Group 2). Values of E_Y for TC2# are close to those of TC1^*2 (where only data at the central layer of the system are assimilated). Otherwise, values of E_Y at the end of the assimilation window for TC1^*3 (where only data at the bottom of the system are assimilated) are lowest when considering the TCs TC1^*1, TC1^*2, TC1^*3, and TC2#. Values of S_Y for all test cases but TC1 are visually indistinguishable. Finally, the temporal behavior of E_h for each test case is consistent with the one displayed by E_Y. These results seem to suggest that the benefit (in terms of E_Y and E_h) of collecting head observations from packers installed along the borehole depends on the observation depth and on the duration of the assimilation period.

Figure 10 Temporal evolution of E_Y (a), S_Y (b), and E_h (c) with α=1 (TC2), 5 (TC2α₁), 10 (TC2α₂), and 100 (TC2α₃).

Figure 11 Temporal evolution of E_Y (a), S_Y (b), and E_h (c) with α=1 (TC3), 5 (TC3α₁), 10 (TC3α₂), and 100 (TC3α₃).

5.4 Effect of inflation on measurement-error covariance matrix (Group 4)

Figure 10 depicts the temporal evolution of E_Y (Fig. 10a), S_Y (Fig. 10b), and E_h (Fig. 10c) for TC2α₁, TC2α₂, and TC2α₃ (see Group 4 in Table 1), where partially screened monitoring wells (Type B wells) are located across the entire domain. The lowest values of E_Y are mainly associated with α=10 (i.e., TC2α₂) while the highest values correspond to α=1 or 100 (i.e., TC2 or TC2α₃, respectively). The magnitude of S_Y is seen to increase with α, in line with Eq. (24) according to which resorting to an inflation factor tends to decrease the strength of the dependence of conductivity estimates on head data. The highest E_h values are generally linked to TC2 at observation times shorter than 10 (i.e., corresponding to the stop of pumping), all other E_h values being otherwise visually indistinguishable.

Based on these results, we conclude that the accuracy of conductivity and head estimates is generally improved when inflating the measurement-error covariance matrix. As stated above, these results are consistent with the observation that inflating the measurement-error covariance matrix results in a reduced weight of the mismatch between modeled and observed values during data assimilation. We recall that using inflation (i.e., setting α>1) may be useful to compensate for relying on an imperfect mathematical model, a scenario which is consistent with TC2–7 in Group 1. For instance, the iterative ensemble smoothers (Chen et al., 2013; Luo et al., 2015) and the ensemble smoother with multiple data assimilation (ES-MDA; Emerick and Reynolds, 2013) rely on the action of an inflation factor on the measurement-error covariance matrix to cope with highly nonlinear systems.

Figure 11 shows the temporal behavior of E_Y (Fig. 11a), S_Y (Fig. 11b), and E_h (Fig. 11c) for TC3α₁, TC3α₂, and TC3α₃ (see Group 4 in Table 1), where fully screened monitoring wells (Type C wells) are located in the entire domain. The lowest E_Y values are mainly associated with α=100 (i.e., TC3α₃). The value of E_Y at the end of the assimilation period is largest for TC3. The magnitude of S_Y values also tends to increase with α in these cases. The highest and lowest E_h values are generally linked to TC3 and TC3α₃, respectively, for the early assimilation times and become virtually independent of α as time progresses. These results suggest that the highest inflation factors required to increase the quality of the data assimilation process (as measured through E_Y and E_h) are associated with the scenarios where head data are collected in fully screened boreholes (see, e.g., TC3α₃ as compared to TC2α₂).

6 Conclusions

We draw the following main conclusions based on this study:

• The use of packers to collect pointwise head data (Type A wells) yields higher accuracy of conductivity estimates than can be obtained by relying on partially or fully penetrating wells. The lowest values of E_Y (average absolute difference between the estimated mean of the logarithm of the conductivity field, Y, and its reference counterpart) are associated with the scenario where Type A wells are set across the domain. It is additionally worth noting that the benefit of installing Type A wells as opposed to partially (Type B) or fully screened (Type C) monitoring wells is mainly associated with regions (e.g., zone 1 in this study) where strong variations of head along the vertical can take place.

• Using depth-averaged head data from Type B and C wells leads to comparable results in our settings, in terms of E_Y, S_Y (square root of the average estimation variance of Y), and E_h (average absolute difference between the estimated aquifer heads and their reference counterparts).

• Neglecting flux exchanges between the aquifer and partiallyor fully screened monitoring wells in the groundwater flow model can significantly deteriorate the accuracy of conductivity estimates. Considering the application of an inflation technique to measurement-error covariance matrix can improve conductivity estimates when an imperfect flow model is applied.

• The computational feasibility and accuracy of the moment equations-based ensemble Kalman filter (MEs-EnKF) are explored. The MEs-EnKF is as accurate as a typical Monte Carlo (MC)-based ensemble Kalman filter, which relies on a large number (on the order of 10 000) of MC realizations. Otherwise, the MEs-EnKF is more efficient than its MC-EnKF counterpart, the latter requiring about 20 times the central processing unit (CPU) time of the former, on the basis of our examples.

Appendix A: Cross-covariance between water levels in partially and or fully screened monitoring boreholes

The water level at well I, $h_I^{\mathrm{w}}$ (with I=1, …, N_w), can be written as $h_I^{\mathrm{w}}=\langle h_I^{\mathrm{w}}\rangle +h_I^{{}^{\prime }\mathrm{w}}$. Making use of Eq. (10) one can obtain the following expression for the water level fluctuation $h_I^{{}^{\prime }\mathrm{w}}$

$$\begin{array}{}\text{(A1)}& h_I^{{}^{\prime }\mathrm{w}}\sum _{i=\mathrm{1}}^n{b}_i{K}_i=\sum _{i=\mathrm{1}}^n{b}_i{K}_i\left(h_i-\langle h_I^{\mathrm{w}}\rangle \right).\end{array}$$

n is the total number of cells according to which the screen of borehole I is discretized. In a similar way, the water level fluctuation at well J, $h_J^{{}^{\prime }\mathrm{w}}$, is given by

$$\begin{array}{}\text{(A2)}& h_J^{{}^{\prime }\mathrm{w}}\sum _{i=j}^m{b}_j{K}_j=\sum _{j=\mathrm{1}}^m{b}_j{K}_j\left(h_j-\langle h_J^{\mathrm{w}}\rangle \right),\end{array}$$

where m corresponds to the total number of cells according to which the screen of borehole J is discretized. Multiplying Eq. (A1) by Eq. (A2), taking expectation and disregarding moments of order larger than two yields

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}& \langle h_I^{{}^{\prime }\mathrm{w}}{h}_J^{{}^{\prime }\mathrm{w}}\rangle \sum _{i=\mathrm{1}}^n{b}_i\langle K_i\rangle \sum _{i=j}^m{b}_j\langle K_j\rangle =\sum _{j=\mathrm{1}}^m\sum _{i=\mathrm{1}}^n{b}_i{b}_j\\ & \left\{\begin{array}{l}\left[\langle K_i\rangle \langle K_j\rangle +\langle K_i^{\prime }{K}_j^{\prime }\rangle \right]\left[\langle h_j\rangle \langle h_i\rangle -\langle h_J^{\mathrm{w}}\rangle \langle h_i\rangle \right.\\ \left.-\langle h_I^{\mathrm{w}}\rangle \langle h_j\rangle +\langle h_I^{\mathrm{w}}\rangle \langle h_J^{\mathrm{w}}\rangle \right]\\ +\left[\langle K_i\rangle \langle K_j^{\prime }{h}_j^{\prime }\rangle +\langle K_j\rangle \langle K_i^{\prime }{h}_j^{\prime }\rangle \right]\langle h_i\rangle \\ +\left[\langle K_i\rangle \langle K_j^{\prime }{h}_i^{\prime }\rangle +\langle K_j\rangle \langle K_i^{\prime }{h}_i^{\prime }\rangle \right]\langle h_j\rangle \\ +\left[\langle K_i\rangle \langle K_j\rangle \langle h_i^{\prime }{h}_j^{\prime }\rangle \right]-\left[\langle K_i\rangle \langle K_j^{\prime }{h}_i^{\prime }\rangle \right.\\ \left.+\langle K_j\rangle \langle K_i^{\prime }{h}_{ii}^{\prime }\rangle \right]\langle h_J^{\mathrm{w}}\rangle \\ -\left[\langle K_i\rangle \langle K_j^{\prime }{h}_j^{\prime }\rangle +\langle K_j\rangle \langle K_i^{\prime }{h}_j^{\prime }\rangle \right]\langle h_I^{\mathrm{w}}\rangle \end{array}\right\}.\end{array}\end{array}$$

Evaluation of Eq. (A3) at second order yields Eq. (13).