Groundwater flow inverse modeling in non-MultiGaussian media: performance assessment of the normal-score Ensemble Kalman Filter

The normal-score ensemble Kalman filter (NSEnKF) is tested on a synthetic aquifer characterized by the presence of channels with a bimodal distribution of its hydraulic conductivities. This is a clear example of an aquifer that cannot be characterized by a multiGaussian distribution. Fourteen scenarios are analyzed which differ among them in one or various of the following aspects: the prior random function model, the boundary conditions of the flow problem, the number of piezometers used in the assimilation process, or the use of covariance localization in the implementation of the Kalman filter. The performance of the NS-EnKF is evaluated through the ensemble mean and variance maps, the connectivity patterns of the individual conductivity realizations and the degree of reproduction of the piezometric heads. The results show that (i) the localized NS-EnKF can characterize the non-multiGaussian underlying hydraulic distribution even when an erroneous prior random function model is used, (ii) localization plays an important role to prevent filter inbreeding and results in a better logconductivity characterization, and (iii) the NS-EnKF works equally well under very different flow configurations.


Introduction
Accurate characterization of the spatial variability of hydrogeologic properties and its corresponding uncertainty is a key issue for environmental risk assessment, site remediation and restoration engineering, and the design of underground repositories for radioactive material.
In a Monte Carlo framework, heterogeneity of hydrogeologic properties is commonly characterized by the following two steps: (i) on the basis of a limited amount of direct measurements (i.e., hard data), multiple representations of aquifer properties are generated by means of the geostatistical techniques such as sequential Gaussian simulation (Deutsch and Journel, 1998), sequential indicator simulation (Gómez-Hernández and Srivastava, 1990), multiplepoint geostatistical approach (Strebelle, 2002;Mariethoz et al., 2010b) or other related methods; and then (ii) on the basis of indirect measurements such as piezometric head and concentration data, inverse modeling is utilized to reduce the uncertainty by integrating these data to better characterize the spatial variability of hydrogeologic properties (e.g. for an overview see Yeh, 1986;McLaughlin and Townley, 1996;Zimmerman et al., 1998;Carrera et al., 2005;Hendricks Franssen et al., 2009;Zhou et al., 2012b).
Commonly used Monte Carlo type inverse algorithms (i.e. which generate many equally likely solutions to the inverse problems) include the self-calibration method (Sahuquillo et al., 1992;Gómez-Hernández et al., 1997;Capilla et al., 1999;Wen et al., 2002;Hendricks Franssen et al., 2003), the pilot point method (Ramarao et al., 1995;LaVenue et al., 1995;Alcolea et al., 2006), the Markov chain Monte Carlo method (Oliver et al., 1997;Fu and Gómez-Hernández, 2009;Alcolea and Renard, 2010), and the gradual deformation method (Hu, 2000;Capilla and Llopis-Albert, 2009), among others.These methods have in common that a multi-part objective function is minimized.The objective function is generally composed of the sum of squared differences (SSD) between simulated and observed L. Li et al.: Performance assessment of the normal-score EnKF state values plus the SSD of prior and calibrated parameters values.In order to minimize this objective function the hydrogeologic parameters are modified using derivative-based methods or by sampling the posterior distribution.The main difference between the various methods is how to solve the optimization problem.
The Ensemble Kalman Filter (EnKF) (Burgers et al., 1998;Anderson, 2001;Reichle et al., 2002;Evensen, 2003) is a further alternative to generate multiple equally likely solutions to the inverse problem.The EnKF is increasingly studied in hydrogeology as well as in petroleum engineering (e.g.Wen and Chen, 2005;Chen and Zhang, 2006;Hendricks Franssen and Kinzelbach, 2008;Sun et al., 2009;Nowak, 2009;Nan and Wu, 2010;Li et al., 2012b).The attractive characteristics of the EnKF are: (i) the efficiency in computation (e.g.Hendricks Franssen and Kinzelbach, 2009 conducted a synthetic example to demonstrate that the needed CPU-time was reduced by a factor of 80 compared with the self-calibration method, to achieve comparable results); (ii) the flexibility of handling multiple sources of uncertainty, for instance, Hendricks Franssen and Kinzelbach (2008) successfully used EnKF to account for the uncertainty of both recharge and hydraulic conductivity together; Li et al. (2012a) used the EnKF to jointly calibrate porosity and hydraulic conductivity by assimilating dynamic head and multiple concentration data; (iii) real-time data assimilation without the need to store all previous states, for instance, Hendricks Franssen et al. (2011) operationally implemented the EnKF to calibrate the conductivity and leakage coefficient together in real-time in a real-world case study.On the contrary, the traditional inverse approaches mentioned above (self-calibration, pilot point, etc.) are CPU-intensive, need recalibration when new data are available and handling multiple sources of uncertainty is less straightforward with these methods.
The EnKF provides an optimal solution when the state vector follows a multiGaussian distribution and the state function is linear (Arulampalam et al., 2002).In the literature of hydrogeology, most of studies applying the EnKF technique assume that the hydraulic logconductivities follow a multiGaussian distribution (e.g.Chen and Zhang, 2006;Hendricks Franssen and Kinzelbach, 2008;Huang et al., 2008;Li et al., 2012a).The significance of accounting for non-multiGaussian distributions of hydraulic conductivity for flow and transport predictions has been stressed in many studies (e.g.Gómez-Hernández and Wen, 1998;Zinn and Harvey, 2003;Knudby and Carrera, 2005).Hence, extending the EnKF to deal with non-Gaussian state vectors would facilitate more extensive applications.Sun et al. (2009) and Zhou et al. (2011Zhou et al. ( , 2012c)), developed variants of the EnKF which are better accommodated to handle non-Gaussianity of parameter distributions.Sun et al. (2009) resort to couple the EnKF with a Gaussian mixture model to update the parameters of a multi-modal distribution by assimilating head data.Zhou et al. (2011Zhou et al. ( , 2012c) ) transformed the augmented state vector (containing both the hydraulic conductivities and piezometric heads) with marginal multi-modal distributions into a new vector with marginal Gaussian distributions, perform the EnKF on the transformed state vector and backtransform the updated vector to the original space.Both studies show that these variants of EnKF result in significant improvements in the characterization of the aquifer properties and in the predictions of flow and transport for non-multiGaussian hydraulic conductivity fields.It is worth noting that, in both studies, it was assumed that the prior random function is known, i.e. the reference field and the initial set of non-multiGaussian conductivity fields are both generated with the same geostatistical approach and the same random function model, which, in this particular case, amounted to using the same training image in the multiple-point geostatistical generator.However, in practice, it is not trivial to decide about a reasonable training image given the limited amount of information available.Partly for this reason, traditional two-point variogram-based geostatistical methods are still employed in practice.It is worth recall that the normal-score transformations used in the NS-EnKF only ensure marginal Gaussianity, with the higher-order moments not necessarily closer to the multiGaussian distribution; it is also worth recall, that the EnKF is applied on the normalscore transformed variables what amounts to having a more non-linear forecast function that before the transformations.One of the motivations of this work is to explore the capacity of the normal score EnKF proposed by Zhou et al. (2011Zhou et al. ( , 2012c) ) to characterize non-multiGaussian media such as a highly channelized aquifer when there are no (hard) conductivity data available and the prior random function model is not the one used to generate the reference field.
Several studies (e.g.Gómez-Hernández and Wen, 1998;Western et al., 2001;Zinn and Harvey, 2003;Knudby and Carrera, 2005;Feyen and Caers, 2005;Lee et al., 2007) showed that flow and transport predictions strongly differ between multiGaussian and non-multiGaussian logconductivity fields, even when both types of fields share the same histogram and the same covariance function.These results demonstrated that the reproduction of the connectivity is of significance in practice.For the inverse conditioning, Kerrou et al. (2008) applied the self-calibration (Gómez-Hernández et al., 1997) method on a fluvial-sediment aquifer with a wrong prior model (i.e.multiGasussian instead of non-multiGaussian) and concluded that the channel structures cannot be retrieved, even when a large number of direct and indirect data are used for conditioning.
In this paper, we apply the normal-score EnKF (NS-EnKF) to a channelized aquifer characterized by a bimodal conductivity distribution, and analyze the performance of the method for fourteen scenarios which differ among them in one or several of the following aspects: the prior random function model, the boundary conditions of the flow problem, the number of piezometers used in the assimilation process, or the use of covariance localization in the implementation of the Kalman filter.None of the scenarios uses any conductivity conditioning data.The analysis focus on the ensemble mean and variance maps, the connectivity patterns of the individual conductivity realizations and the degree of reproduction of the piezometric heads.
This paper revisits the mathematical framework of the NS-EnKF as applied in the synthetic experiment and briefly presents its numerical implementation (Sects. 2 and 3); the impact of the choice of the prior random function model is discussed in Sect.4.1, the effect of using localization functions for the covariance while computing the Kalman gain is discussed in Sect.4.2, the performance under different boundary conditions inducing very different flow patterns in the aquifer is discussed in Sect.4.3, and, finally the impact of reducing the number of conditioning piezometers is discussed in Sect.4.4.The paper ends with some conclusions (Sect.5).

Flow equation
The flow equation of an incompressible fluid in saturated porous media in a Cartesian coordinate system can be obtained by combining the continuity equation and Darcy's law (Bear, 1972): where h[L] is the piezometric head; ] is the volumetric source flow per unit volume; S s [L −1 ] is the specific storage coefficient; t[T ] is the time; ∇• = (∂/∂x + ∂/∂y + ∂/∂z) is the divergence operator of a vector field, and ∇ = (∂/∂x, ∂/∂y, ∂/∂z) T is the gradient operator of a scalar field.

The normal-score ensemble kalman filter with localization
The NS-EnKF method aims at generating equally-likely realizations of parameters and state variables, conditioned to real-time measurements such as piezometric head data following non-Gaussian marginal distributions.The basic algorithm is described in Zhou et al. (2011) and is summarized next for the case of assimilating observed piezometric head data at time t in a bimodal aquifer.
1. Initialization step.Generate the initial ensemble of logconductivity fields (in case that there are logconductivity measurements, these fields must be conditional to them); this is achieved by any of many geostatistical simulation algorithms, such as sequential indicator simulation (e.g.Deutsch and Journel, 1998) or the multiplepoint method (e.g.Strebelle, 2002).The choice of the algorithm will depend on the random function model adopted to describe the spatial variability of logconductivity.Here, it is assumed that the scale of hydraulic conductivity measurements have a support coinciding with that of the numerical model discretization gridblocks.If there were a discrepancy between the measurement scale and the gridblock scale, an upscaling technique would have been needed to reconcile these two scales (e.g.Zhou et al., 2010;Li et al., 2011a,b).
2. Normal-score function step.Build the logconductivity local conditional distribution functions at each grid cell from the ensemble of logconductivity realizations, and the corresponding local normal-score transfer functions (φ) (for the details see: Goovaerts, 1997).
3. Forecast step.For each realization, the transient groundwater flow Eq. ( 1) is solved from time t − 1 to t using standard block-centered finite differences (e.g.Harbaugh et al., 2000;Li et al., 2010).The solution can be schematically represented by: where X k−1 and Y k−1 denote the hydraulic conductivity and piezometric head estimates at time t k−1 respectively, and Y k are the forecasted piezometric heads at time step t k ; f represents the groundwater flow model including the boundary conditions, stresses and other known parameters.
4. Normal-score function step.Build the piezometric head local conditional distribution functions from the ensemble of forecasted realizations, and the corresponding local normal-score transfer functions (ϕ).
5. Filter step.Update the state variables by the NS-EnKF based on the observed head data at time t.
a. Normal-score transform each conductivity, each forecasted head, and also the observed heads b. Build the augmented normal-score transformed vector ˆ , which includes both the transformed conductivities ( X) and the transformed forecasted heads ( Ŷ).
where ˆ k,j is the j -th ensemble member of the augmented state vector at time t k .c. Calculate the localized Kalman gain (G k ) (e.g.Gaspari and Cohn, 1999;Hamill et al., 2001;Chen and Oliver, 2010).

576
where C X Ŷ is the cross-covariance between transformed logconductivity and head data; C Ŷ Ŷ is the covariance of transformed head data, C D is the diagonal matrix of expected measurement error variances; • indicates the Schur product; and ρ X Ŷ and ρ Ŷ Ŷ are localization functions used for the C X Ŷ and C Ŷ Ŷ , respectively; the specific localization functions used in this work will be given later.
where superscripts u and f denote updated and forecasted, respectively; is a random observation error vector, and Ŷf is a vector containing the transformed forecasted heads at sampling locations (a subset of Ŷ).
e. Back transform the state vector ˆ .
6. Go to step 3 for the next time step.
The NS-EnKF assures that the non-Gaussianity of the marginal distributions of states and parameters is preserved, through the use of the normal-score transform.The distancedependent localization functions are used to reduce the influence of spurious correlations for large separation distances that may appear in the covariance computed through the ensemble of realizations.(i) the SNESIM code, a pattern-based multiple-point geostatistical algorithm, is utilized to generate a facies realization (Strebelle, 2002).(In contrast with traditional twopoint variogram-based models, the multiple-point method is able to reproduce complex and realistic curvilinear structures characteristic of fluvial deposits.)We used the training image in Strebelle (2002), which is commonly used for benchmarking to compare algorithms (e.g.Wu et al., 2008;Mariethoz et al., 2010a;Zhou et al., 2012.a).This training image serves as a conceptual model for the channelized non-multiGaussian aquifer, where the channels have high conductivities representing preferential paths and are embedded in a floodplain fine-grid deposits with low conductivities (see Fig. 1); (ii) after the facies are generated, the facies are populated with realizations of conductivity fields generated using sequential gaussian simulation algorithm (Gómez-Hernández and Journel, 1993) with the parameters listed in Table 1.This procedure results in the spatial distribution of logconductivity values shown in Fig. 2a, which serves as the reference lnK realization.This realization displays wellconnected sand channels (approximately 30 % of the system) on a low-conductivity matrix.The histogram of lnK, shown in Fig. 2b, clearly shows a bimodal logconductivity distribution typical of fluvial deposits.To explore the role of boundary conditions on pattern identification, two sets of boundary conditions are considered: (i) Boundary conditions inducing parallel flow (Fig. 3a), consisting of a combination of prescribed head boundary (West), prescribed flux boundary (East) and impermeable boundaries (North and South); (ii) boundary conditions inducing radial flow (Fig. 3b), with impermeable boundaries and a set of injection/extraction wells in a typical configuration for contaminant remediation or for reservoir exploitation.For both boundary conditions, groundwater flow is simulated for a period of 500 days starting with an initial head equal to zero everywhere.The simulation time is discretized into 100 time steps, the sizes of which follow a geometric series with a ratio of 1.05.The specific storage is assumed constant and equal to 0.003 m −1 .The MODFLOW 2000 code (Harbaugh et al., 2000) is used to solve the transient flow Eq. ( 1).The head data are sampled at the first 60 time steps (approximately 67.7 days) and will be used as the conditioning data in the NS-EnKF.The spatial distribution of the sampling wells is shown in Fig. 3.Although the sampled head data are errorfree since they are taken from the "true" field, in the stochastic model the sampled head data should have noise added to it.In the following data assimilation procedures it was assumed that they have a normalized measurement error with a mean of zero and a standard deviation of 0.05 m.

Scenarios
Fourteen scenarios are analyzed with the characteristics indicated in the Table 2.The impact of the prior model, the use of a localization function and the different boundary conditions is tested with help of fourteen scenarios (Table 2).
The scenarios are organized to show how the characterization of the hydraulic conductivities evolves starting from an unconditional ensemble of realizations as piezometric heads are incorporated, and then when localization is used.This evolution is analyzed for each of the two flow patterns induced by the two sets of boundary conditions, and also for the cases in which the correct and an erroneous prior random functions are used.Scenarios S1, S2, S7 and S8 are unconditional, they are the base cases containing the initial ensemble realizations to be used by the NS-EnKF in the other scenarios.They will produce the worst results in terms of conductivity characterization and piezometric head predictions and will serve as a starting point to assess the benefits of using the correct prior model, the assimilation of dynamic piezometric head information, and the use of localization.(Gómez-Hernández and Srivastava, 1990) is used for the second set.The parameters defining the two multiGaussian random functions are given in Table 1. Figure 4 shows the first realization for the correct and the incorrect model, their histograms and their variograms.While both realizations follow the same histogram and variograms, the variogram-based sequential indicator simulation cannot generate the curvilinear patterns for the sand facies observed in the training image, which are well reproduced in the multipoint sequential simulation.
For the scenario S3, S5, S9, S11, S13 and S14 head data are assimilated without using localization whereas for the scenarios S4, S6, S10 and S12 localization is employed.The EnKF tends to filter inbreeding if the ensemble size is small, and the heterogeneity is large (Hendricks Franssen and Kinzelbach, 2009).It is therefore expected that the NS-EnKF applied to a complex geological environment will suffer this problem, too.To analyze this aspect, a fifth-order distancedependent localization function (Hamill et al., 2001) is used to reduce the sampling errors in the computation of the ensemble covariances.Considering previous works (Chen and Oliver, 2010), the size of the domain, the separation between piezometers, and the distance between wells for the boundary conditions inducing radial flow, the localization function is set to zero at a distance of 80 m, implying that for distances larger than that, the sampled non-stationary covariances used to compute the Kalman gain are zero.More specifically, the distance function that is used for both ρ X Ŷ and ρ Ŷ Ŷ in Eq. ( 7) is: ( where d is the distance between observed head and lnK data (when applied to ρ X Ŷ ) or between observed heads (when ap- Here, the distance function is constant in time, isotropic, and used both for the cross-covariances between observed heads and parameters and the covariances of the head data.Scenarios 13 and 14 stand apart from the rest of the conditional scenarios in that the number of piezometers used for assimilation is reduced to one third.Recall that no conductivity conditioning data are used for any scenario, therefore, the ability to identify the conductivity patterns is solely thanks to the dynamic head data that are assimilated at each time step, and it is expected that the quality of the characterization of the conductivities will be affected by the number of conditioning piezometers.
For all the scenarios, 1000 unconditional realizations of logconductivity are generated.

Evaluation criteria
For every scenario, the following criteria are used to assess the results (Chen and Zhang, 2006;Zhou et al., 2012c): 1. Root mean square error (RMSE).RMSE measures the accuracy of the estimation.
where X i is either logconductivity lnK or hydraulic head h at location i, X i represents its ensemble mean value at node i, X ref,i is the reference value at location i, N b is the number of nodes.

ES(
where σ 2 X i is the ensemble variance at location i.The smaller the values for RMSE and ES, the better the prediction of variable X.As discussed in Chen and Zhang (2006), when RMSE and ES have a similar magnitude the resulting ensemble variance provides a more realistic measure of the uncertainty associated to the ensemble mean estimate.

Connectivity. The flow and transport behaviors are
strongly impacted by the presence of connected highconductivities and connected low-conductivities (e.g.Wen and Gómez-Hernández, 1998;Knudby and Carrera, 2005;Fernàndez-Garcia et al., 2010).Here, we adopt the connectivity function defined by Pardo-Igúzquiza and Dowd ( 2003) to evaluate the connectivity.In this work, two steps are needed to analyze the connectivity for the continuous variable: i.An indicator image is obtained as: where α is the threshold that classifies the logconductivity into sand and shale.In our case, α is set to zero since this value approximately splits the histogram of logconductivity in two (see Fig. 2).
ii.The code CONNEC3D (Pardo-Igúzquiza and Dowd, 2003), a computer program for connectivity analysis of 3-D random model, is used to calculate the probability that two points with logconductivity larger than α are connected.The result is a function of distance that gives the probability that two points located in the sand facies are connected by a continuous path in that same facies.

Results and discussion
Ensembles of conductivity realizations are generated for the fourteen scenarios according to the configurations described earlier.As mentioned earlier, scenarios S13 and S14 stand apart in that they are the only ones using a smaller number of conditioning piezometric heads, for this reason most of the figures group the results for scenarios S1 to S12.Scenarios S13 and S14 are analyzed later by themselves.

Effect of prior random function model
For the scenarios S1, S2, S7 and S8 no conditioning data were used.For those scenarios, the ensemble average of hydraulic logconductivity is approximately spatially uniform and equal to the prior mean (see Fig. 5) and, similarly, the ensemble variance map is approximately equal to the prior variance (see Fig. 6).As expected, the unconditional average maps do not capture any spatial pattern, although individual realizations do display the spatial variability according to their prior random function models (see the first realization in Fig. 4a and d).
For the scenarios where head data are used for assimilation, the ensemble means of hydraulic conductivity (estimated after 60 time steps of assimilation), depict the facies distribution well in visual comparison with the reference.
The ensemble variances of logconductivity are clearly reduced in comparison with the prior variances and display a feature with higher variance at the boundaries of the channels and at places far away from the head measurements such as the boundaries of the domain.Comparing the two columns on Figs. 5 and 6, it is quite striking to realize that, even for the wrong prior random function model, assimilating the piezometric head data with the NS-EnKF, yields ensembles of realizations which, on average, capture the location of the channels in the reference field and display the highest uncertainty at the boundary channels.
For a more quantitative analysis, Fig. 7 shows how both the RMSE and the ES evolve for the different conditional scenarios as a function of the updating step.It is not surprising that the wrong prior always gives higher values for the average departure between individual realizations and the reference and also for the average spread, being most noticeable for the case of parallel flow.
Regarding the reproduction of the connectivity function, Fig. 8 shows the individual connectivity functions for all realizations, their average and the connectivity of the reference field for the first 12 scenarios.The discrepancy between the results obtained with the correct and the wrong prior are more significant in this case, something that is quite understandable since the wrong prior only uses the indicator variogram of the facies distribution as the starting set of ensemble realizations.The average connectivity of the wrong prior is significantly smaller than the reference connectivity for the unconditional case; then, conditioning to piezometric heads helps increasing the average connectivity as well as the spread of the individual connectivity functions, however, the connectivity remains below that of the reference for all cases.The average connectivity of the correct prior starts slightly above the reference one with a very wide spread of the individual functions, and gets closer to the reference when the piezometric heads are assimilated.
Finally, regarding the reproduction of the piezometric heads, Fig. 9 shows the RMSE and ES as a function of the updating step for the first 12 scenarios.In these figures we have as a reference the metric values for the unconditional cases to show the important reduction in both metrics induced by conditioning to the piezometric heads.Like in Fig. 7, these metrics are average values computed over the entire aquifer, not just for the conditioning piezometers.The behavior of the wrong prior is very similar to that of the correct prior with the only noticeable difference that the localization tends to worsen the RMSE for the wrong model whereas for the correct prior, localization improves the RMSE.
In summary, the NS-EnKF displays a great potential to capture the patterns of variability of logconductivities on the sole basis of piezometric information for clearly non-multiGaussian fields even when the prior random function model is not fully consistent with the reference realization.

Effect of localization
It is well known that the standard EnKF has the problem of filter inbreeding, which leads to a final ensemble of realizations too similar among them resulting in the underestimation of uncertainty of flow and transport predictions (Hendricks Franssen and Kinzelbach, 2009;Devegowda et al., 2010;Zhou et al., 2011).The ensemble size we have used is large, 1000 realizations, yet, when analyzing the piezometric head covariances, and cross-covariances with logconductivity we observe the same type of spurious correlation values that, eventually, lead to filter inbreeding; therefore, we have decided to apply localization techniques and compare the resulting ensembles with and without localization.In general, localization should increase the ensemble spread and, at the same time, reduce the RMSE, that is, there is a gain in accuracy at the cost of precision, which also helps to avoid filter inbreeding.
Analyzing Figs. 5 and 6 we observe that the ensemble mean hydraulic conductivity maps produced with the localized NS-EnKF are smoother than the same mean maps without localization, and capture better the main channels of the reference field.All artifacts in the non-localized realizations are removed after localization; the local ensemble variances in the localized ensembles increase, particularly significant values appear in the northwestern corner of the parallel flow exercise an area without conditioning data and too close to a prescribed head boundary where the sensitivity of the piezometric head to a change in logconductivity is close to zero.These results may indicate that even in this case with a large number of realizations, the ensemble variance estimate obtained from a non-localized NS-EnKF may be too optimistic.
The effect of localization is also quite clear when analyzing the RMSE and ES of the calibrated conductivity fields (see Fig. 7).On one hand, localization induces a reduction of the RMSE in all cases, and an increase in the ES.In all cases, these changes result in values of RMSE and ES that are very similar in value at the end of the updating steps, which is indicative, according to Chen and Zhang (2006), of reduced filter inbreeding.Taking the ratio ES/RMSE as an index of filter inbreeding, this is worse for the radial flow scenarios without localization.In previous synthetic studies, Chen and Zhang (2006) and Hendricks Franssen and Kinzelbach (2008) have found that 200 realizations are enough to avoid filter inbreeding without localization for multiGaussian fields.Here we observe that filter inbreeding can be serious when EnKF is used without localization in a fluvial-type aquifer.
Regarding sand connectivity it is not clear whether localization improves its characterization.The spread of the connectivity functions is much larger for the localized ensemble than for the non localized one, and, at least for parallel flow, the average of the connectivity functions does not get closer to the real one after localization; the small spread exhibited by the non localized connectivity functions makes the reference connectivity one fall outside the cloud of individual functions for the radial flow case.After localization, the reference connectivity is within the cloud for the correct prior function model, while still falls partly outside for the wrong prior.The falling of the reference connectivity function outside the cloud of individual functions is also an indication of underestimation of the ensemble variability induced by the non localized NS-EnKF.
Finally, localization has a minimal impact in the reproduction of the piezometric heads as seen in Fig. 9, the RMSE is almost the same with and without, although localization serves to reduce the RMSE for the correct prior scenarios, but increases it for the wrong prior ones.Localization, as expected, increases the ES in all cases.
It is concluded that, for this particular case, 1000 realizations is not enough to avoid filter inbreeding in the application of the EnKF.For this reason we recommend the use of localization together with the NS-EnKF to improve the results.

Effect of boundary conditions
The NS-EnKF has been tested for two sets of boundary conditions which basically induce parallel flow and radial flow, the rest of the setup being the same.
The characterization of the hydraulic conductivities is good under both flow conditions.The mean of the final ensemble of realizations gets closer to the reference field for the parallel flow case than for the radial flow one, as can be observed by looking at the channels identified in Fig. 5, their width is more uniform and their extent go through the entire aquifer for parallel flow.As already mentioned, the local ensemble variance (Fig. 6) show the highest values for the cells close to the prescribed head boundary and extending through the shale area in the parallel flow case as it would be expected given the lack of sensitivity of head to conductivity changes at those locations.
There is not a significantly different behavior on the evolution of the RMSE and the ES (Fig. 7) for the two sets of boundary conditions.Just notice that the final values of both metrics tend to converge to the same value after localization for both flow conditions, while prior to localization they are different (larger RMSE and smaller ES for radial flow than for parallel flow).
Regarding the reproduction of the connectivity, it is apparently more difficult to reproduce it for radial flow case than for parallel flow; for radial flow, there seems to be an important filter inbreeding that leaves the reference connectivity function outside the cloud of individual functions; apart from that, there is no significant difference on the performance of the NS-EnKF in characterizing the connectivity under radial or parallel flow conditions.
Finally, regarding the reproduction of the piezometric heads (Fig. 9), the main difference between parallel and radial flow is that for radial flow the overall reduction of RMSE with regard to the unconditional realizations is smaller than for parallel flow, but this is easily understandable since the piezometric head fluctuations with respect to the initial state are stronger for radial flow than for parallel due to the several pumping and recharge wells used for the radial flow scenarios.
We conclude that the NS-EnKF performs similarly well for either radial flow or parallel flow conditions.However, from the analysis of the RMSE(lnK) and ES(lnK), the radial flow case apparently is more prone to filter inbreeding than the parallel flow case.Inbreeding that disappears after localizing the filter.

Effect of number of conditioning piezometers
Scenarios 1 to 12 have been analyzed considering that there are 111 piezometers (Fig. 3) providing dynamic hydraulic head measurements during the first 66.7 days (first 60 simulation time steps), which are assimilated through the NS-EnKF to update an ensemble of hydraulic conductivity realizations.For this number of piezometers, the NS-EnKF seems to perform reasonably well for a variety of scenarios, including the case in which an incorrect prior random function model is used.Since there are no logconductivity data, the information provided by the piezometer is crucial for the identification of the hydraulic conductivity spatial patterns.We decided to reduce the number of piezometers to one third of the original ones and keep only 32 of them (Fig. 10a) to analyze two more scenarios, under parallel flow and without localization; scenario 13 assumes the correct prior random function model and scenario 14 the wrong one.
Maps B and D in Fig. 10 show the ensemble mean of logconductivity and maps C and E the ensemble variance for the correct and wrong prior, respectively.As expected, the characterization of the hydraulic conductivity field is worse and the ensemble variance is larger when fewer piezometers are used; however, the characterization is not bad.The overall patterns are well captured and the highest local variances occur at channel borders, with a much better performance for the correct prior random function than for the wrong one.Although not shown, the connectivity in scenarios 13 and 14 is clearly deteriorated when the number of piezometers goes down to 32.
We conclude that the number of piezometers is important for the characterization of the logconductivity, particularly when the prior random function is wrong.

Discussion
This paper presented a performance assessment of the NS-EnKF, a real-time data assimilation algorithm, in a two dimensional bimodal aquifer, and focused on four aspects: (1) the impact of random function choice; (2) the effect of introducing a localization function (3) the flow regimes and (4) the number of piezometers used during assimilation.
The results show that assimilating a sufficiently large number of transient piezometric head data with the NS-EnKF allows a good characterization of non-multiGausian patterns such as those depicted by channelized aquifers, even if the prior geostatistical model is erroneous and based on twopoint indicator variograms.This result may have important implications for groundwater modeling studies, since it implies that a sufficiently large number of piezometers may offset the negative impact of an erroneous choice of the random function model used to generate the starting ensemble of conductivity realizations.For those who are not keen on using training images to characterize the spatial variability of the hydraulic conductivity, a variogram-based simulation algorithm may serve as the starting point, leaving to the NS-EnKF the task to identify the curvilinear features that the variogram-based algorithm cannot.Although it might be argued that we had access to the "correct" prior histogram in the wrong prior analysis, if the contrast in conductivities between the channel and non-channel facies is large enough, the really important a priori information that would be needed would be the proportions of each facies and the contrast between their average conductivities.
Our results are more optimistic than the ones from Kerrou et al. (2008), who concluded that an erroneous prior geostatitical model cannot be corrected by inverse modeling, even if many hydraulic head data are available.We attribute this impossibility to the inverse modeling method utilized and claim that the NS-EnKF can help to produce an acceptable characterization of the hydraulic conductivities even when the prior random function model is incorrect.
We also demonstrate the importance of coupling the NS-EnKF with a distance-dependent localization function, even for a case such as this in which the ensemble consisted of 1000 realizations.A future research topic is the improvement of the localization functions, allowing for anisotropy and non-stationarity (i.e.time and/or space dependence), and with the option to have different formulations for piezometric head covariances and piezometric-conductivity cross-covariances.
Apparently, the NS-EnKF works equally well for both parallel and radial flow conditions.And, it is clear that for a case such as this one in which the conductivity identification is based solely on the information conveyed by the dynamic piezometric head data, the number of piezometers is very important, more so, when the prior random function model is wrong.

Conclusions
The complexity of subsurface geological structures calls for a characterization technique, which can incorporate the conductivity measurement data as well as the responses of aquifers in an efficient and effective way.The NS-EnKF, a L. Li et al.: Performance assessment of the normal-score EnKF real-time data assimilation technique which couples EnKF with normal score transformation, presented by Zhou et al. (2011).The method avoids the classical multiGaussianity inherent to most simulation algorithms.
In this paper, we have analyzed 14 scenarios to present a detailed analysis of the impact of (i) prior model choice, (ii) use of localization functions, (iii) flow regime, and (iv) number of piezometers used, on the performance of the NS-EnKF algorithm in a synthetic 2-D bimodal aquifer.The following conclusions can be drawn from the simulation exercises: -The NS-EnKF performs relatively well when an incorrect prior geostatistical model is used.Results are worse than when the correct prior is used, but acceptable.
-Coupling the NS-EnKF with a distance-dependent localization function improves both the characterization of conductivity and the prediction of groundwater flow.
-The performance of NS-EnKF was similar for two different flow scenarios (parallel flow and radial flow), in both cases, the patterns of conductivity are identified very well.
-The quality of the characterization is directly related to the number of piezometers used, more so when a wrong a priori model is used.

Fig. 3 .
Fig. 3. Spatial distribution of sampled head data and flow configurations: (A (B) Radial flow (the unit of flow in wells: m 3 /d). 34

Fig. 3 .
Fig. 3. Spatial distribution of sampled head data and flow configurations: (A) Parallel flow (B) Radial flow (the unit of flow in wells: m 3 d −1 ).

Fig. 4 .
Fig. 4. (A-C) show the spatial distribution of lnK, histogram and variogram of the 1st realization from the correct prior model; (D-F) show the same characteristics of the 1st realization from wrong prior model.In (C) and (F), solid line and dotted line correspond to the experimental facies variogram in X and Y direction, respectively.

Fig. 7 .
Fig. 7. RMSE and ES of lnK for the different cases.

Fig. 9 .
Fig. 9. RMSE and ES of simulated heads for the different scenarios.

Fig. 10 .
Fig. 10.(A) flow configuration and spatial distribution of sampled head data.(B) the ensemble mean of logconductivity with correct prior model.(C) the ensemble variance of logconductivity with correct prior model.(D) the ensemble mean of logconductivity with wrong prior model.(E) the ensemble variance of logconductivity with wrong prior model.

L. Li et al.: Performance assessment of the normal-score EnKFTable 1 .
Parameters of the random functions describing the spatial continuity of the sand and shale logconductivities.Mean Standard deviation Variogram type λ x [m] λ y [m] sill

Table 2 .
Definition of scenarios.