Abstract

The localized normal-score ensemble Kalman filter is shown to work for the characterization of non-multi-Gaussian distributed hydraulic conductivities by assimilating state observation data. The influence of type of flow regime, number of observation piezometers, and the prior model structure are evaluated in a synthetic aquifer. Steady-state observation data are not sufficient to identify the conductivity channels. Transient-state data are necessary for a good characterization of the hydraulic conductivity curvilinear patterns. Such characterization is very good with a dense network of observation data, and it deteriorates as the number of observation piezometers decreases. It is also remarkable that, even when the prior model structure is wrong, the localized normal-score ensemble Kalman filter can produce acceptable results for a sufficiently dense observation network.

1. Introduction

Parameter identification is a critical step in constructing a reliable model. The process of recognizing model parameters is also referred to as inverse problem or data assimilation; that is, assimilate the system state data to identify the model parameters. In hydrogeology, the parameters describing the movement of groundwater can vary over several degrees of magnitude within the same aquifer, and the characterization of their spatial variability is very important. For this reason, inverse modeling has been a subject of active research, aiming to take the maximum advantage of the state observation data for the characterization of model parameters, such as hydraulic conductivities. One of the inverse methods that has attracted more attention in hydrogeology lately is the ensemble Kalman filter on its different implementations.

In this paper we focus on building a groundwater flow model, in which the model parameters are the hydraulic conductivities (also referred to as permeability) and the model responses are the piezometric heads (representing the elevation at which water rises in a perforation drilled in the aquifer). The groundwater flow state equation is given by [1] subject to initial and boundary conditions, where is the divergence operator (), is the gradient operator , is hydraulic conductivity [], is hydraulic head [], is specific storage [], is time [], and is source or sink []. The equation is discretized in space and time and solved numerically.

In a practical situation, the modeler will have geological information about the main features of the aquifer, which allows to have a first idea of the range and type of variability of hydraulic conductivity and to establish the aquifer boundary conditions; the modeler may also have a few local measurements of hydraulic conductivity (these measurements are expensive and this is why they are generally scarce) and will have some measurements of the piezometric heads (these measurements are less expensive since they only require inserting a small pipe in the aquifer). To capture the large heterogeneity of the aquifer, its discretization must be very fine, ending up with numerical models with tens of thousands of cells and, on occasion, millions of cells. Each model cell must be informed with a hydraulic conductivity value, a specific storage coefficient, and an initial piezometric head. Specific storage is not so heterogeneous as the hydraulic conductivity, so we will concentrate only on the heterogeneity of hydraulic conductivity. Piezometric heads vary smoothly in space, and they can be easily interpolated from their values at the available piezometers. The characterization problem resides in how to assign the hydraulic conductivity values to the model cells so that the heterogeneity observed in the field is respected by the model, and model predictions are consistent with observations. This characterization problem becomes acute when hydraulic conductivities (or its logarithms) cannot be statistically described by a multiGaussian distribution.

Inverse modeling is commonly used to solve this characterization problem. Hydraulic conductivity values (model parameters) are derived from the observed (in space and time) piezometric heads (system state). We will demonstrate a new implementation of the ensemble Kalman filter in a synthetic aquifer displaying complex curvilinear features, which cannot be statistically characterized by a multiGaussian distribution. We will show how the complex heterogeneity of the reference aquifer can be retrieved from the transient information about the state of the system, state that is recorded in a number of piezometers distributed over the aquifer domain.

The Kalman filter is an assimilation procedure based on the repetitive application in time of two steps: forecast and analysis. In the forecast step, a prediction of the spatial distribution of the piezometric heads is made based on the current characterization of hydraulic conductivity, then, in the analysis step, a correction of the conductivity values and of the piezometric heads is made based on the differences between piezometric head predictions and observations at the measurement locations. This analysis step relies in the calculation of a nonstationary covariance function involving parameter values and system states on the entire system domain. As time progresses, and more piezometric head data are assimilated, the characterization of the hydraulic conductivity distribution improves.

The ensemble Kalman filter is a variant of the Kalman filter designed to work with nonlinear prediction models. The ensemble Kalman filter (EnKF) is widely applied in multiple disciplines such as meteorology, oceanography, petroleum engineering, and hydrology (e.g., [26]). In the ensemble Kalman filter, multiple realizations of the parameters are generated, all of which considered as a plausible representation of reality; these realizations are progressively (as time passes and state observation data is acquired) modified so that they are consistent with the observations. The final outcome should be an ensemble of realizations, all with the same patterns of variability, but not with exactly the same values, all of them coherent with the state observations. The variability of parameters and state variables through the ensemble of realizations will serve to make predictions with an uncertainty measurement attached.

The popularity of the EnKF can be attributed to its computational efficiency compared with other Monte Carlo inverse method [7], to its easy implementation, and to the simplicity to adapt it to different forward models. The EnKF is optimal when both model parameters and state variables are multiGaussian, and they are related by a linear state equation [8]. The EnKF is still very effective when the state equation is not linear, but it fails when parameters and/or state variables do not follow a multiGaussian distribution [9]. Such is the case of aquifers from fluvial deposits, in which neither the hydraulic conductivities nor the piezometric heads can be characterized by multiGaussian distributions. We propose the use of the normal-score ensemble Kalman filter (NS-EnKF) proposed by Zhou et al. [10]. They introduced the normal-score transformation of both parameters and state variables into the standard EnKF. This transformation does not ensure that the resulting transforms are multiGaussian, only their marginal distributions are Gaussian thus the EnKF is applied on variables that follow, at least, marginal Gaussian distributions. As will be shown, this transformation, which makes parameters and state variables one step closer to multi-Gaussianity, has been proven quite effective in the characterization of very heterogeneous, non-multiGaussian conductivities.

The NS-EnKF has already been demonstrated to perform well in characterizing the hydraulic conductivities of bimodal aquifers [1012]. Here, we present a localized version of the NS-EnKF. Localization eliminates spurious correlations in the calculation of the different covariances as will be seen later. The objective of this paper is to evaluate, in a synthetic aquifer, the localized NS-EnK against the number of state observation locations (piezometers), to analyze its performance under different types of flow conditions (steady and transient states), and also to see what happens when the prior structure of the model departs from the true one. (Since we are using a perfectly known synthetic aquifer, we can analyze the impact of prior model, i.e., using the correct one from the reference or not.)

2. The Localized Normal-Score Ensemble Kalman Filter

The standard EnKF algorithm is proposed by Evensen [13] and Burgers et al. [14] and described in detail by Evensen [15]. The NS-EnKF algorithm is explained by Zhou et al. [10]. The difference between the NS-EnKF and the standard EnKF resides in the introduction of pre- and postprocessing steps carried out on the joint state vector ensemble. The localized NS-EnKF method mainly consists of the following steps.

(1) Generate the Initial Model
The localized NS-EnKF starts with an ensemble of realizations that have been generated following a given prior random function model. In our example, hydraulic conductivities are the unknown parameters to be identified. Multiple hydraulic conductivity realizations are generated serving as the starting point of the data assimilation process. There are many algorithms for the generation of the spatially correlated parameters that can be looked up in the geostatistics or in the spatial statistics literature.

(2) Normal Score Transformation of the Hydraulic Conductivities
Compute the local cumulative distribution functions (CDFs) at each cell through the ensemble of realizations [16]. These local CDFs serve to construct normal-score transform functions that are used to transform each log-conductivity realization into a new realization , that is, , with being a vectorial function with a different normal-score transform function for each cell. When analyzing the new ensemble of transformed realizations, at each grid cell the components of the new ensemble follow a marginal Gaussian distribution with zero mean and unit variance.

(3) Forward Simulation (Forecast Step)
For each realization of the ensemble and given boundary conditions, the state vector at time is computed based on the state vector at time using a state transfer equation while the parameter vector is assumed unchanged during the forward simulation. In our case, and represent hydraulic log conductivities and piezometric heads at time , respectively, and is the solution to the groundwater flow model, that is, the numerical solution of (1.1) numerically in a discretized aquifer.

(4) Normal Score Transformation of the Piezometric Heads
Similarly as it was done for the hydraulic conductivities, from the ensemble of predicted piezometric heads, the local CDFs for time at each grid are built, noted as . Then the CDFs are used to transform both simulated piezometric heads and observation heads into the univariate Gaussian distribution space, that is, and , where is the observation head at time . (Recall that is a vectorial function with a different expression for each grid cell, for simplicity we assume that observation data has been measured at grid cell centers, and therefore the same normal-score transform can be applied to predictions and observations; otherwise it would be necessary to construct specific normal-score transform functions for the observations based on the interpolation of the nearby predicted piezometric heads.)

(5) Update the State Vector (Analysis Step)
The analysis step is applied on the normal-score transforms of both the parameters and state variables.

(i) Build the Augmented Transformed State Vector Ensemble
This vector contains all realizations of the ensemble. The member of the augmented state vector ensemble can be expressed as where and are normal-score transformed log conductivities and piezometric heads, respectively. (It is called augmented, because it includes state variables but also model parameters.)

(ii) Compute the Localized Kalman Gain
For the sake of simplicity the subscript is ignored. The Kalman gain is computed as where is the model ensemble covariance, is the observation matrix, and is the observation error covariance. Assuming that the observation locations coincide with the grid nodes, the matrix will consist of only 0’s and 1’s. In this case we do not have to compute all the components of the covariance explicitly but instead compute directly the and which correspond to the covariances and , respectively (in the standard implementation of the EnKF, these covariances are computed on the untransformed augmented vector). Then, the expression of Kalman gain can be rewritten as where is the covariance between the normal-score transformed vector and the piezometric head (the latter at measurement locations only) and is the covariance between the normal-score transformed piezometric heads (at measurement locations). These covariances are computed from the ensemble of realizations as where is the number of realizations in the ensemble. These numerically computed covariances generally display larger than zero values between locations that are far apart; these spurious correlations, which are most noticeable when the ensemble size is small or the heterogeneity large, should be removed before updating the augmented vector. For this purpose, we use the localization of the covariance as recommend by Franssen and Kinzelbach [7]. Two localization matrices ( and ) modify the ensemble covariance matrix estimates using a Schur product, resulting in the following expression for the Kalman gain [17, 18]: A distance function is used for the localization matrices [19], with values that decrease from 1 to 0 as a function of the distance; beyond a certain distance the localization matrix coefficients equal zero. (The function distance used in the synthetic example that follows is radially symmetric; it is shown for any direction in Figure 1 for distances up to 80 m, and it has a value of zero for distances larger than 80 m.)

(iii) Update the State Vector (as in the Standard EnKF)
Consider where is the vector with the updated state variables at time and is the vector computed from the forward simulation and then normal-score transformed; is the normal-score transformed observation vector at time ; is an observation error characterized by a normal distribution with zero mean and a diagonal covariance .

(6) Back Transform
Each component of the updated state vector is back-transformed using the previously constructed local CDFs corresponding to the hydraulic log conductivity and piezometric head, that is, for and for .

Repeat the process from step 3 to step 6 until all the available head observations are assimilated.

3. Synthetic Experiment

3.1. Reference Field

A synthetic bimodal aquifer composed of sand and shale is used to evaluate the performance of the localized NS-EnKF in different scenarios. The domain size is 300 m × 240 m × 10 m discretized into 100 columns by 80 rows by 1 layer. The aquifer is assumed confined with a thickness of 10 m. The reference facies field is generated with a multiple-point geostatistical algorithm, using a training image shown in Figure 2. Each facies is then populated with continuous log-conductivity values generated by GCOSIM3D [21] with parameters listed in Table 1. The reference ln field and its histogram are presented in Figure 3. The histogram of log conductivity is bimodal, with modes coinciding with the means of the sand and shale distributions, 2.1 ln() and −1.4 ln(), respectively, and it has a global mean of −0.33 ln and a global standard deviation of 1.68 ln. This field is clearly non-multiGaussian, both in its univariate statistics as in its curvilinear patterns of spatial continuity.

The groundwater flow equation is solved for the reference field using the MODFLOW computer code [22] under the following boundary conditions: impermeable boundaries in the north and south, a river with water level at 2 m and heterogeneous river bed thickness and conductivity in the west, and a prescribed flow rate of 270.5 m3/d along the east boundary. The flow rates in the east boundary are distributed depending on their water supply capability; that is, large flow rates correspond to zones with big conductance. The initial head is 0 m over the domain. The total simulation period of 500 days is discretized into 100 steps with step length that increases by a ratio of 1.05. Specific storage is assumed constant and equal 0.003 m−1.

3.2. Initial Model and Conditioning Data Scenarios

The assimilation of piezometric head data is performed by the localized NS-EnKF. Observed piezometric heads from the piezometers in 60 time steps (67.7 days) serve as conditioning data. In the synthetic experiment we will test the influence of the following conditions.(1)Flow Regime: Steady-State Flow versus Transient Flow. The amount of information contained in the transient piezometer observations is much larger than that contained in steady-state observations. This will imply a larger characterization power by the transient heads than by the steady-state ones, as will be shown below.(2)Influence of the Number of Piezometers. With the increase of the number of the piezometers, more observation data are collected and the structure of the conductivities is better characterized. The influence of the number of the piezometers is examined in this experiment; that is, the piezometers are reduced until a number below which the localized NS-EnKF does not work. Figure 4 shows the locations of the piezometers in four scenarios.(3)Prior Model Selection. Lack of data may lead to the use of a prior model that departs from the real one. This issue has already been addressed by Li et al. [11], where they show that even with a wrong prior model, the identification is acceptable if the density of the observation network is sufficient. Here, we focus on the influence of the number of piezometers on conductivity structure recognition with a correct and a wrong prior model. The wrong prior model is fully described by Li et al. [11]; we only mention here that the wrong prior model shares the same histogram and two-point covariance with the correct one, but all other higher-order moments approach a multiGaussian model (“Figure 4" in [11]). Figure 5 shows the ensemble mean and variance for both correct and wrong prior log conductivity model. The correct and wrong prior models are used to generate the first batch of conductivity realizations and also to build the local CDFs used for the normal-score transformations.

Nine scenarios (Table 2) are considered in this experiment.

4. Results and Discussions

Figure 6 displays the ensemble mean of the updated hydraulic log conductivities in different scenarios (listed in Table 2) along with the reference. The ensemble mean should reflect the main patterns of the reference aquifer. Figure 7 shows the ensemble variance of log conductivities corresponding to the ensemble mean in Figure 6. The variance field is used to evaluate the uncertainty of the estimates. Figure 8 is the evolution of the average absolute bias and the average ensemble spread as more observation data are assimilated. The and the quantitatively evaluate accuracy and uncertainty of the estimation, respectively. They are defined as where is the estimated log conductivity at time step , grid and realization , is the reference log-conductivity at node , is the number of grid cells, is the number of realizations, and is the variance over the ensemble at time step and location . Figure 9 shows the hydraulic head evolution with time at one of the observation piezometers (no. 21). Figure 10 is the hydraulic head evolution with time at another piezometer (no. 103), the data from which is used as observation only in scenarios 2 and 6; this latter piezometer is used for validation of the updated model for the remaining scenarios. The head evolution shown in Figures 9 and 10 past day 66.7 (the end of the assimilation period) serves to evaluate the ability of the updated conductivity fields to perform predictions.

4.1. Steady-State Flow versus Transient Flow

Ensemble mean and variance of the updated log conductivity for the steady-state flow are the images labelled as scenario 1 in both Figures 6 and 7, respectively. Comparing the ensemble mean with the reference we can hardly see any channel pattern, and comparing the ensemble variance with the prior variance, the uncertainty is reduced to a very small extent. We argue that, for hydraulic conductivity identification purposes, steady-state information is clearly insufficient, even, as is the case, a very dense network of observations is used (111 wells as shown in Figure 4(a)). On the contrary, transient data is extremely valuable for identification purposes as can be seen when comparing the ensemble means and variances of scenarios 1 with 2 in Figures 6 and 7. Scenarios 1 and 2 share the same prior model and the same configuration of piezometers, but one is at steady state and the other is at transient state. The ensemble mean of scenario 2 delineates almost perfectly the channels of the reference, and the ensemble variance is reduced quite drastically with respect to the prior variance, showing also some ridges that could be used to identify the channel boundaries.

4.2. Number of Piezometers

Here we will test the influence of the number of observation piezometers under transient flow conditions: the piezometers are reduced gradually from 111 to 63, 32 and 10 as shown in Figure 4. The prior ensemble mean of log conductivity is almost uniform and shows no channel before piezometer observations are assimilated (Figure 5). In Figure 6, we find, as expected, that the ability of the localized NS-EnKF for conductivity identification deteriorates as the number of piezometers is reduced, for both the correct prior model (scenarios 2 to 5) and the wrong prior model (scenarios 6 to 9). When 111 observation wells are available, the ensemble mean resembles closely the reference and the ensemble uncertainty is reduced substantially; when the number of piezometers is reduced to 63, the channels are not so connected; when the number of conditioning piezometers is further reduced to 32, the borders between the two facies are not distinct; when only 10 piezometers are located, the updated ensemble mean cannot capture the main channel structure no matter whether the correct model is used or not. In parallel, the uncertainty of the estimation (Figure 7) increases as the number of piezometers decreases. Through the comparison we find that keeping a dense observation piezometer is important to identifying the log conductivities since the influence area of each observation is limited. Also we argue that for parameter identification, it is very important that observation piezometers be located in both facies, the highly conductive channels and the less conductivity background shale; therefore, in designing a sampling network, whatever geological information available should be used to help in the configuration of the piezometers. In this experiment, when 32 piezometers are located, the longest channel feature in the reference aquifer is captured by the ensemble mean (scenarios 4 and 8 in Figure 6). If we overlap the piezometer distribution (Figure 4(c)) on the reference log-conductivity field we can find that in each stretch of the channel there exists a piezometer, which explains the fuzzy but recognizable ensemble mean. When the number of piezometers is reduced further to 10, the piezometer observation network is too scarce to capture the general structure of the underlying parameters (log conductivities).

We compute the and the to quantitatively evaluate the performance of the inverse method regarding different scenarios (Figure 8). We find that the values of both and decrease with time and with the number of data. Both values approach a plateau at the end of the 60 time steps.

Figures 9 and 10 show the hydraulic head evolution at two of the piezometers (no. 21 and no. 103), both of which serve as the conditioning observations wells when 111 piezometers are used and no. 103 is used only for model validation when 63, 32 or 10 conditioning wells are available. We find that for both wells the head prediction uncertainty is reduced in all updated models compared with the results by the prior model, in which no data assimilation is carried out. Furthermore, the uncertainty spread for the 111 piezometers cases (S2 and S6) is, as expected, narrower than that in the rest scenarios.

It is very important to notice that no hydraulic conductivity data have been used for the generation of the prior ensemble realizations. The only hydraulic conductivity information used is the prior random function model, that is, the prior bimodal histogram and the training image.

4.3. Correct Prior versus Wrong Prior

The influence of the prior model on the performance of the NS-EnKF has been discussed by Li et al. [11], and the main point here is to test the influence of the number of piezometers with regard to correct and wrong prior model. In Figure 6, the left column from the second to the last row shows the results with the correct prior model, and the right column shows the result with the wrong prior model. We find that the ensemble mean of the updated log-conductivity field exhibits always more noise in the case of the wrong prior model but it still is capable of identifying the main patterns of the reference. In Figure 8 we find that the values of the and the of the correct prior model are smaller than those of the wrong prior model. Moreover, the magnitudes of the and the are comparable indicating no apparent filter inbreeding in the application of the localization NS-EnKF. Note that the hydraulic head predictions (Figures 9 and 10) with the correct prior model are comparable to those with wrong prior model although the characterization of the conductivity structure by the former one is better than by the latter one. This can be attributed to the nonuniqueness of the inverse problem solution and the smoothing effect of the groundwater flow equation [23, 24].

5. Conclusion

The objective of this paper is to investigate the influence of different factors on the performance of the localized NS-EnKF in identifying ln patterns by assimilating hydraulic heads. The ln in the synthetic aquifer is characterized by a bimodal distribution. The main conclusions are as follows (1) information from steady-state flow conditions is insufficient to identify the channel structure; transient flow information is critical for this purpose; (2) the data assimilation algorithm is able to characterize the hydraulic conductivity when using transient piezometric data, with as little as 32 piezometers uniformly distributed in the aquifer. The geology of the aquifer should be taken into account for the design of the observation network, and it is the one that will mark the minimum number of piezometers needed; (3) the transient data contains enough information so that even when a wrong prior model is used, the identification of the conductivity channels is acceptable.

Acknowledgments

The authors gratefully acknowledge the financial support by the Spanish Ministry of Science and Innovation through project CGL2011-23295. The authors want to thank the reviewer for the comments which help improving the quality of the paper.