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

We employ an approach based on ensemble Kalman filter coupled with stochastic moment equations (MEs-EnKF) 10 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), (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 15 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 20 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 an incomplete/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.

Abstract. 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) multinode 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.
The EnKF can assimilate data sequentially through a realtime 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 variabledensity 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 deeplearning-based model as a surrogate of a solute transport model to reduce the CPU time associated with ensemblebased 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(Elci et al., , 2003Konikow et al., 2009;Zhang et al., 2019). Elci et al. (2001Elci et al. ( , 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. (2001Elci et al. ( , 2003, Post et al. (2007), Konikow et al. (2009, 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 multinode 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., , 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 steadystate 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., 2013Panzeri et al., , 2014. Panzeri et al. (2013Panzeri et al. ( , 2014Panzeri et al. ( , 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 logtransmissivity 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.
subject to initial and boundary conditions 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 0 (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 0 (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 σ 2 Y . Tartakovsky and Neuman (1998a, b) derive integrodifferential 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 integrodifferential 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.

Mean head and flux
We start by expressing a given random quantity, I, as the sum of its (ensemble) mean, I , and a zero-mean random fluctuation, I . Here, and in the following, mean head and flux are approximated up to second order in σ Y as and the following equations hold 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 geo-metric mean of K(x), and r(x, t) is the second-order residual flux evaluated as r(x, t) = lim y→x [−∇ x u(y, x, t)] (e.g., Xia et al. (2019) and references therein), where u(y, x, t) = K (y)h (x, t) (2) is the second-order approximation of the cross-covariance between hydraulic head and conductivity, computed as detailed in Sect. 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(y, x, t) is the covariance of the logconductivity field, C fK (y, x, t) = K (y)f (x, t) (2) is the second-order cross-covariance between conductivity K(y) and forcing term f (x, t), and U 0 (y, x) = K (y)H 0 (x) (2) is the second-order approximation of the cross-covariance between K and initial hydraulic head. Note that U 0 (y, x) vanishes when H 0 is deterministic.

Head covariance
The equation governing the evolution of the (secondorder) head covariance between space-time locations (y, τ ) and (x, t), C h (y, x, τ, t) = h (y, τ )h (x, t) (2) , is given by (Zhang, 2002) where u(x, y, τ ) = K (x)h (y, τ ) (2) is given by Eq. (8) and 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., Tartakovsky (2000, 2002), . 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. (2002Winter et al. ( , 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.

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. 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 w I , (Type B and/or C) at a given time t (omitted in the following equations for brevity) can be evaluated through 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 multinode 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 i K 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.

Moments for hydraulic head at observation wells
where, starting from Eq. (10), one can obtain the zero-, h Here, i , K G,i , and σ 2 Y,i 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; is the cross-covariance between conductivity and head at the ith cell; K i h w I (2) 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 w I ) and (2) , can be evaluated as (see also Appendix A, Eqs. A1-A3) where T G,i = b i K G,i , T G,j = b j K 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 logconductivity covariance between the ith and j th cells of boreholes I and J , respectively; u j i (or u ij ) is the crosscovariance between the conductivity of the j th cell of well J (or the ith cell of well I ) and head of the ith cell of well I (or the j th cell of well J ); and C h,ij is the head covariance between the ith cell of well I and the j th cell of well J . Terms u j i (or u ij ) and C h,ij can be readily obtained by solving Eqs. (8) and (9). The cross-covariance between h w I at a given time t and aquifer head h at time τ (at a given location; omitted for brevity), i.e., C h w Here, u i,τ = K i h (τ ) (2) represents the cross-covariance between conductivity at the ith cell of well I and aquifer head (at a given location) at time τ . Likewise, C h,i,τ = h i (t)h (τ ) (2) 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 w I at time t and conductivity at a given location in the aquifer, C h w I K = h w I (t)K (2) , can be expressed as is the covariance between logconductivity at the ith cell along the monitoring borehole I and log-conductivity at a given location in the domain, and u i = K h i (t) (2) 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 crosscovariances 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.

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 w I , h i , and Q i can then be obtained by relying on the Thiem (1906) formulation as where r 0 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 The cross-covariance between Q i at time t and aquifer head at time τ , , can be expressed as where . Finally, one can obtain the following expression for the crosscovariance between Q i and K, C Q i K = Q i K (2) ,

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 (0) , u, and h (2) are sequentially obtained by solving Eqs. (6)-(8), respectively. Details associated with the evaluation of C h , which requires h (0) 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., . We then recall that Zhang (2002) is also unknown, t being a constant temporal step size) by solving for C h (y, x, τ = t, t ) from t = 0 to t = t. While this procedure can be computationally heavy for long times, Zhang (2002) points out that when flow changes only mildly, C h (x, y, τ = t, t − t) ≈ C h (x, y, τ = t − t, t − 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 h (y, x, τ = t, t) from C h (y, x, τ = t − t, t − t) via (i) computing C h (y, x, τ = t, t − 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 h (y, x, τ = t − t, t − t) as initial condition and then (ii) assessing C h (y, x, τ = t, t) by solving Eq. (9) using C h (y, x, τ = t, t − 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 (0) , f (2) , 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 − t, t], the main computational cost stems from the evaluation of u(y, x, t), C h (y, x, τ = t, t − t), and C h (y, x, τ = t, t), each of these requiring N times the computational cost (hereafter denoted as C MEs c ) 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 MEs c . Note that the computational effort is reduced to 2 N C MEs c 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 i , C fK = C Q i K , and C fh = C Q i h τ . The evaluation of the (ensemble) moments of interest across time interval [t − 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, h w I iter+1 , and iter, h w I iter , is lower than a preset value ε. The main computational effort required for these evaluations corresponds to 3 (iter + 1)N C MEs c for each time step.

Ensemble Kalman filter coupled with moment equations
We start by introducing the mean system state vector ϕ as 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 − 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 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, C f h is the predicted N × N-dimensional head covariance matrix, (Eq. 9), C f h w h is the N w × N -dimensional predicted cross-covariance between well and aquifer head (Eq. 14), C f h w is the predicted N w × N w -dimensional covariance of well head (Eq. 13), C f Y h is the predicted N × Ndimensional cross-covariance between Y and aquifer head h (Eq. 8), C f Y h w is the predicted N × N w -dimensional crosscovariance between Y and well head h w (Eq. 15), and C f Y 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 and where C ϕd (= P f H T ) is the (2N +N w )×d-dimensional crosscovariance between the system state and the simulated observations, matrix H of dimension d × (2N +N w ) is the observation operator that describes the relationship between the system state and the observations, C dd (= HP f H T ) is the d × ddimensional 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 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)  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 , C h w h , C h w , C Y h w ) should be omitted in Eqs. (21)-(25).

Illustrative examples
We consider a three-dimensional domain (Fig. 3a) of size 600×600×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 25 × 25 × 13 nodes, for a total of 34 560 tetrahedrons. A partially penetrating pumping well pumps at a constant rate of 1000 for 0 ≤ t ≤ 30, after which water withdrawal stops and a recovering process takes place for 30 < t ≤ 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 −3 , respectively. The natural logarithm of conductivity, Y , is modeled as a spatially correlated secondorder stationary random field with covariance given by Here, δ i and λ i denote, respectively, the lag and correlation scale between two points along direction x i (with i = 1, 2, 3).
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 1 = b 13 = 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 ε = 10 −6 .
The effective radius of the monitoring wells is evaluated as (Chen and Zhang, 2009) r 0 = 2 0.25 e −0.75π x ≈ 0.113 x = 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. 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, λ 1 = λ 2 = 100 and λ 3 = 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 λ 1 from the pumping well) and Type B boreholes being installed in zones 2 and 3 (at distances larger than λ 1 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α 1 and TC3α 1 ), 10 (TC2α 2 and TC3α 2 ), and 100 (TC2α 3 and TC3α 3 ; 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 ∼ ±0.005 and ∼ ±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 where Y i up , (σ 2 Y,i ) up , and Y r i indicate the estimated mean, variance, and reference Y values at the ith node of the computational mesh, respectively; h i up and h r i 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).

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.

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.
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.

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.
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 variabil- , and E h (c, f) for TC2# * 1 and TC3# * 1 (a-c) and for TC2# * 2 and TC3# * 2 (d-f). ity 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.

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α 1 , TC2α 2 , and  TC2α 3 (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α 2 ) while the highest values correspond to α = 1 or 100 (i.e., TC2 or TC2α 3 , 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α 1 , TC3α 2 , and TC3α 3 (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α 3 ). 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α 3 , 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α 3 as compared to TC2α 2 ).

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.