Ensemble Kalman filter versus ensemble smoother for assessing hydraulic conductivity via tracer test data assimilation

Estimating the spatial variability of hydraulic conductivity K in natural aquifers is important for predicting the transport of dissolved compounds. Especially in the nonreactive case, the plume evolution is mainly controlled by the heterogeneity of K. At the local scale, the spatial distribution of K can be inferred by combining the Lagrangian formulation of the transport with a Kalman-filter-based technique and assimilating a sequence of time-lapse concentration C measurements, which, for example, can be evaluated on site through the application of a geophysical method. The objective of this work is to compare the ensemble Kalman filter (EnKF) and the ensemble smoother (ES) capabilities to retrieve the hydraulic conductivity spatial distribution in a groundwater flow and transport modeling framework. The application refers to a two-dimensional synthetic aquifer in which a tracer test is simulated. Moreover, since Kalmanfilter-based methods are optimal only if each of the involved variables fit to a Gaussian probability density function (pdf) and since this condition may not be met by some of the flow and transport state variables, issues related to the nonGaussianity of the variables are analyzed and different transformation of the pdfs are considered in order to evaluate their influence on the performance of the methods. The results show that the EnKF reproduces with good accuracy the hydraulic conductivity field, outperforming the ES regardless of the pdf of the concentrations.


Introduction
One of the most challenging tasks in groundwater flow and transport modeling is the assessment of hydraulic properties of the porous media such as the hydraulic conductivity (K), whose spatial heterogeneity is typically very pronounced in natural aquifers and thus particularly difficult to characterize.For this reason, many efforts have been made for the development of inverse models capable of estimating aquifer hydraulic properties at various scales.Extensive reviews of these methods have been provided by Carrera et al. (1993); McLaughlin and Townley (1996); Zimmerman et al. (1998) and Vrugt et al. (2008).Hydraulic conductivity and head data have been typically used to constrain inverse models (e.g., Chen and Zhang, 2006;Hendricks Franssen and Kinzelbach, 2008;Rubin et al., 2010).More recently, the use of information derived by geophysical methods has been increasingly proposed and investigated for the estimation of parameters in groundwater modeling (e.g., Camporese et al., 2011;Pollock and Cirpka, 2012).
Among the several mathematical tools available for the inversion of hydrologic data, approximate Bayesian methods such as the ensemble Kalman filter (EnKF) and its variations (Evensen, 2009b) allow for seeking an ensemble of independent samples conditional to the measurements, all representing equally likely realizations of the actual variability of the hydraulic parameters.These methods have been applied, for example, by the following authors: Chen and Zhang (2006), who used the EnKF to estimate the hydraulic conductivity both in two-and three-dimensional domains assimilating hydraulic head data and K measurements; Liu et al. (2008), who obtained the K distribution at the MADE site (e.g., Boggs et al., 1992) by assimilating, in two subsequent steps, hydraulic head and concentration data collected during a tracer test; and Bailey and Baú (2010), who used the ensemble smoother (ES) by assimilating hydraulic head and E. Crestani et al.: EnKF vs. ES for tracer test data assimilation groundwater return flow volume measurements to estimate the hydraulic conductivity distribution in a synthetic twodimensional case.
More recently, Hendricks Franssen et al. (2011) applied the EnKF to jointly calibrate the hydraulic conductivity and leakage coefficient in real time in an unconfined aquifer.In Li et al. (2012), the EnKF was used to map the hydraulic conductivity and porosity fields by assimilating dynamic piezometric data and multiple concentration data.In Bailey and Baú (2012), the ES was iteratively applied to estimate the parameters of a geostatistical model through assimilation of water table elevation data.Tong et al. (2012) used the EnKF in a synthetic two-dimensional aquifer to estimate the hydraulic conductivity by assimilating solute concentration data measured in a large number of observation wells.
A fundamental hypothesis for the application of Kalmanfilter-based methods is that all the variables must be distributed as a joint multivariate Gaussian probability density function (pdf).In many applications, including groundwater flow and transport in heterogeneous aquifers, the "Gaussian approximation" is not honored and few efforts have been made recently trying to enforce this hypothesis by means of different data transformations.Béal et al. (2010), for instance, applied the Gaussian anamorphosis transformation to the variables of a three-dimensional coupled physical-biogeochemical model of the North Atlantic, while Schoeniger et al. (2012) used the same technique to estimate hydraulic conductivity through the assimilation of 3-D hydraulic tomography data.Another example includes Zhou et al. (2011), who applied the normal score transform to piezometric data used for the estimation of hydraulic conductivity.
A limitation of the above methods to enforce Gaussianity is that they operate univariate (marginal) transformations to multivariate problems.As a result, the multivariate dependence structure between the variables cannot be ensured (Zhou et al., 2011).A possible alternative to data transformation is the use of particle filter (PF) methods, which, like EnKF, are also based on a Monte Carlo approach, but do not require the variables to be normally distributed.On the other hand, a very large number of particles is needed for adequate sampling of high-dimensional state and parameter spaces, resulting in an excessive amount of computational time (Montzka et al., 2012).
Given the importance that the EnKF and the ES are acquiring as parameter estimation modeling tools in groundwater hydrology, there is the need to investigate in more detail their capabilities and the theoretical implications related to their use in a context that is much different from the one for which they were originally developed, i.e., the optimal estimation of system states only (Nowak, 2009).The objective of this work is thus to compare the EnKF and the ES capabilities to retrieve the hydraulic conductivity spatial distribution in a groundwater flow and transport modeling framework.In addition, the issues related to the non-Gaussianity of concentration data are analyzed by considering different transformations of the relevant pdfs.
The EnKF and the ES are here implemented in the same Lagrangian transport modeling framework proposed by Crestani et al. (2010) and Camporese et al. (2011) in order to estimate the hydraulic conductivity field by assimilating concentration measurements derived from a tracer test in a two-dimensional synthetic aquifer.We hypothesize to follow the full spatio-temporal evolution of the solute plume, which would be observed in an electrical resistivity tomography experiment (e.g., Perri et al., 2012).

The ensemble Kalman filter and the ensemble smoother
According to Evensen (2009a), the combined parameter and state estimation problem for a dynamical model can be formulated as finding the joint pdf of the parameters and model state, given a set of measurements and a dynamical model with known uncertainties.Using Bayes' theorem, the problem can be written in the simplified form where f (y, α) is the joint pdf for the model state y (as function of space and time) and the parameters α, f (z|y, α) are the likelihood function of the measurements z, and γ is a normalization constant whose computation requires the evaluation of the integral of Eq. ( 1) over the multi-dimensional solution and parameter space.Note that in writing Eq. (1), we implicitly assume that the only uncertainty in the model formulation lies in the parameters, while boundary and initial conditions are perfectly known.If, as usual, we work with a model state that is discretized in time, we can represent y at fixed time intervals as y i = y(t i ), with i = 0, 1, . . ., k.
If we further assume that the model is a first-order Markov process, we can define the pdf for the model integration from time t i−1 to t i as f (y i |y i−1 , α).Let us now assume that also the measurements z can be divided into subsets of measurement vectors z i , collected at the same time steps of the model t i , and that the measurement errors are uncorrelated in time.
Under these hypotheses and from Bayes' theorem, Eq. ( 1) becomes (Evensen, 2009a) in which the pdf of the parameters f (α) is expressed explicitly.By rewriting Eq. ( 2) as a sequence of iterations and integrating out the state variables at all previous times, we obtain  The combined parameter and state estimation can thus be formulated sequentially using Bayesian statistics, under the condition that measurement errors are independent in time and the dynamical model is a Markov process.Equations ( 1) and ( 3) can be solved numerically by means of a Monte Carlo approach, which approximates the probabilistic information conveyed by the conditional pdfs of the state, parameters, and measurements with an ensemble of realizations of size NMC.Each of the NMC state vectors is propagated in time according to the forecast model, which can be expressed as a vector-valued discrete-time state equation: where y j (t) is the j -th (j = 1, . . ., NMC) state vector predicted by the model at the time t, α j is the j -th set of model parameters, A is the operator relating the system state at the current time t to the system state at the previous time τ , and y j 0 is the initial condition at time t 0 .We assume here that the model is error free and only the parameters are affected by uncertainty.
At time t i , m i measurements are available and the model describing the relation between these measurements and the system state is also expressed as the following vectorial equation: where H is the operator that maps the model state to the measurement locations, y true (t i ) is the true state, and z j (t i ) is the j -th vector of observations at the time t i , which is obtained by perturbing the m i measurements with a random noise w j (t i ), representing the measurement errors.Here we assume that w j (t i ) is normally distributed with expected value equal to zero and assigned variance σ 2 meas .Under the hypothesis that the pdfs for the model prediction as well as the likelihood are Gaussian, it is possible to update each realization of the system state according to the following equation, which is obtained by minimizing the model error covariance matrix (Evensen, 2009a): where y j upd is the j -th realization of the updated system state, P e is the prior estimate of the system state error covariance matrix, which is computed by sampling the ensemble statistics, and R e is the measurement error covariance matrix.
Within the formulation of the ES, Eq. (1) requires that the pdfs are approximated from an integration of the ensemble through the whole assimilation time period.In other words, all of the data are processed by Eq. ( 6) in one step, and the solution is updated as a function of space and time, using the space-time covariances estimated from the ensemble of model realizations.
On the other hand, within the EnKF formulation, Eq. ( 3) is a sequential expression that can be solved through incremental updates given by Eq. ( 6) with the data at the current time step only.
Note that if the model is linear and all states, errors, and parameters are multivariate Gaussian, then both the EnKF and the ES are exact in the Bayesian sense, and lead to identical results at the final time of data assimilation.The technical difference between the EnKF and the ES is analogous to the difference between kriging and sequential kriging as described by Vargas-Guzman and Yeh (1999), who proved the identity between step-wise and one-in-all conditioning for the linear/multi-Gaussian case.
Moreover, in both the EnKF and the ES, when the model pdfs and the likelihood function are Gaussian, the Bayesian formulation corresponds to the minimization of a quadratic cost function (Evensen, 2009a).The only difference is that in the ES the time dimension is included in the optimization, while in the EnKF the minimization is performed at each assimilation time.However, when the pdfs are not Gaussian, Eq. ( 6) is an approximation and no longer yields optimal updates in terms of a cost function minimization, although the general Bayesian formulation remains valid, i.e., sampling of a posterior pdf is still performed.

The inversion model
The Lagrangian approach described in Dagan ( 1989) is here adopted to describe the movement of a tracer plume through a saturated, spatially heterogeneous porous medium.As described in Crestani et al. (2010) and Camporese et al. (2011), the solute cloud is driven by the effective velocity field obtained at steady state solving the groundwater flow equation with appropriate boundary conditions.In natural sedimentary aquifers the spatial variability of the porosity φ can be considered negligible with respect to that of K (e.g., Gelhar, 1993).Consequently, the concentration C, related to the pdf of the Lagrangian trajectory of equation x = X t (t; a, t 0 ), where a is the initial position of the considered particle and t 0 the injection time, is fully controlled by the spatial distribution of the hydraulic conductivity.When dealing with realworld applications, our interest is related to the average concentration C over a finite volume V whose centroid is at x.By defining as M = φC 0 V 0 , the total mass of initial concentration C 0 uniformly injected in the volume V 0 , an estimation of the concentration value is given by where N is the number of particles released to simulate the solute and δ is the Dirac delta function.Here we assume that any diffusion or local scale dispersion is negligible.

E. Crestani et al.: EnKF vs. ES for tracer test data assimilation
The problem of estimating the hydraulic conductivity field by using concentration measurements is solved by considering a system state constituted only by the model parameters: where Y 1 , . . ., Y n are the log-transformed hydraulic conductivity values (Y = ln K) at the n nodes discretizing the domain.At time t 0 , NMC realizations of the log-transformed hydraulic conductivity field are generated from a prior pdf f (α) and the state vectors are built as in Eq. ( 8).
In the EnKF application, starting with the same initial concentration C 0 for each Y field, the solute plume is propagated forward in time to the first measurement time t 1 , using the Lagrangian transport model.At time t 1 the log-transformed hydraulic conductivities (Y 1 , . . ., Y n ) j are updated based on the m 1 measurements available by means of Eq. ( 6), which leverages the cross correlation between Y and C, expressed by the product P e H T .The process continues sequentially: first, a propagation step over each interval between t 0 and t i and, second, an update step of the log-transformed hydraulic conductivity values at each measurement time t i are performed.The process stops at the time t tm corresponding to the last measurement.It must be observed that, in order to ensure mass conservation and consistency between updated hydraulic conductivities and concentrations throughout the simulation, after each assimilation step the plume evolves in the updated Y fields starting from t 0 and with the same initial solute distribution.This recursive application of the EnKF, also known as re-start EnKF (Wen and Chen, 2006;Hendricks Franssen and Kinzelbach, 2008), differs from the classic sequential one commonly adopted since by restarting the plume evolution with the same initial concentration in the updated Y fields, we apply the predictive analysis expressed by the Kalman gain on a system state that is improved at each assimilation step.
In the ES application there is no need for sequential updates, and the parameters are estimated in a single offline step.At time t 0 each realization of the ensemble of Y fields is initialized with the same concentration distribution C 0 and the solute plume is propagated forward in time until the last assimilation time, and the concentration distributions for all measurement times are recorded.The measurement vectors are thus assembled as i.e., the perturbed measurements for all the measurement times (tm) are stored in a vector of dimension [ tm i=1 m i ], for each Monte Carlo simulation.Since the ES does not require sequential updates, the computational time is significantly lower than in the EnKF (about 10 times in the tests presented in the following sections) even though the dimension of the vectors and matrices in Eq. ( 6) is larger than in the EnKF.  3 Numerical experiments

Model setup
A two-dimensional reference Y field is taken into account to compare the proposed inversion models in a number of numerical experiments.Following the recommendation by Ababou et al. (1989), the domain has dimensions of 8 L × 8 L, where L is an arbitrary and consistent unit length, and is discretized along each direction into L/4 sided cells, for a total of 33 × 33 = 3297 corresponding nodes.The multivariate normal distribution of the reference field is the result of a single unconditional generation fitting to an isotropic exponential covariance model with spatial mean Y = 0.35, variance of σ 2 Y = 0.42, and correlation scale λ = 1 L. The random function Y is generated by an improved sequential Gaussian simulation algorithm (Baú and Mayer, 2008).The flow field is simulated using a standard finite volume solver at steady state with appropriate boundary conditions that ensure a constant mean gradient.The resulting Eulerian velocity field is used for the computation of the trajectories of N particles suitable for the simulation of a contaminant release.Dirichlet boundary conditions are applied at x = 0 L (h = 100.0L) and at x = 8 L (h = 95.2L), while Neumann noflow boundary conditions are imposed along the remaining sides of the domain.A graphic representation of the reference field is given in Fig. 1.A tracer test is simulated by assuming an instantaneous solute injection with initial transversal size of 6 L (from y = 1 L to y = 7 L) and longitudinal size of 0.5 L (centered in x = 0.875 L).The solute is simulated by 16 983 particles uniformly distributed, with the particle trajectories computed by means of the Pollock's particle tracking post-processing algorithm (Pollock, 1988) and the concentration computed according to Eq. ( 7) every 0.5 T , where T is any consistent time unit.The total simulation time is t = 4 T .Figure 2 shows the plume evolution in the reference field at t = 2 T and t = 4 T .In this work the inversion is carried out by assimilating concentration data.In particular, we include in the measurement vectors all the concentration data greater than zero in order to use all the available information on the plume evolution.For the EnKF, the resulting dimension of the measurement vector (number of concentration data assimilated) is 79,90,92,102,111,116,123, and 122 at times 0.5, 1, 1.5, 2, 2.5, 3, 3.5, and 4 T , respectively.Consequently, the measurement vector in the ES has dimension 79 + 90 + 92 + 102 + 111 + 116 + 123 + 122 = 835.The performances of the EnKF and the ES in retrieving the Y fields are compared to one another in a number of different scenarios described in Table 1.In these scenarios we also analyze the implications of the concentration non-Gaussianity through various marginal transformations of the pdf.The prior geostatistical parameters and the measurement uncertainty are kept constant in all scenarios as the objective is to study the issues related to the effect of the pdf of the model variables on the EnKF and the ES performances in terms of parameter estimation.Note that the spatial mean and the variance of the prior distributions of Y are significantly different from the corresponding parameters for the reference system (see Table 1).These assumptions serve to represent typical conditions in which both < Y > and σ 2 Y of the true system are unknown.In scenario 1, concentration measurements are assimilated without any transformation, i.e., with their original pdf, whereas in scenarios 2, 3, and 4, different pdf transformations are used to evaluate how they affect the parameter estimation by approaching the Gaussian requirement.
Based on preliminary sensitivity analyses, all scenarios are simulated using an ensemble size of 2000.This number of realizations is computationally affordable and guarantees a proper description of the dispersion process for the range of log-transformed hydraulic conductivity variance used here (Bellin et al., 1992;Salandin and Fiorotto, 1998).
It should be noted that the ensemble of prior Y fields is synthetically generated by the same algorithm used to create lines show the direction along which the cross-correlation structure between Y and C is evaluated.the true field and that the application of the Lagrangian transport model in the true field yields exactly the true concentration distribution.The knowledge of the true state allows us to select the concentration measurements used during the assimilation and to evaluate the performance of the EnKF and the ES with respect to a known reference system.
The estimate of the hydraulic conductivity fields in the various simulated scenarios is assessed by means of the rootmean-square error (RMSE), computed as where n is the total number of nodes of the discretized domain, Y sim,i is the ensemble mean of the Y values estimated at the i-th node, and Y true,i is the true Y value at the i-th node.Finally, we also assess the plume evolution as simulated in the reconstructed Y fields by means of the concentration RMSE, whose formulation is analogous to that of Y .

Scenarios with untransformed concentration pdfs
In scenarios 1a and 1b, the concentration values are assimilated in the update procedure without modifying their original pdf, for both the EnKF and the ES.The spatial distributions of Y resulting from the inversions are reported in Fig. 3.The comparison between the retrieved fields and the reference field (Fig. 1) shows that the EnKF is quite effective in estimating the lnK field, whereas the ES performs rather poorly.The same figure shows also the ensemble variance maps of Y estimation for both the EnKF and the ES.As expected, the final variance is much smaller for the EnKF because of the numerous updates that bring the ensemble to converge toward a unique solution.Obviously, the ES, relying only on one update, shows a larger final variance.This indicates that the ES solution could probably be improved if an iterative procedure was used.However, we elected not to follow such an approach since applying the ES iteratively would lead to losing one of the main advantages of the scheme, i.e., its offline, one-time application at the end of the simulation period, and thus its lower computational demand    (computational time for the ES in these simulations is about 1/10 of the time required by the EnKF).
In these scenarios, both the EnKF and the ES are affected by the approximations related to the Gaussian assumption of Eq. ( 6) as the tracer concentration pdf at early travel times departs significantly from the normal pdf (e.g., Salandin and Fiorotto, 1998).Nevertheless, the EnKF can handle these approximations better than the ES as, at each update, the realizations are steered toward the true solution and the Gaussian increments of the ensemble members lead to an approximately Gaussian ensemble distributed around the true solution.This property of the sequential updating is not exploited in the ES, where realizations evolve freely until the end of the simulation, exacerbating the effects related to non-Gaussian ensemble distributions.This is demonstrated clearly by Fig. 4, which shows, for both the EnKF and the ES, the bivariate pdf of concentration and log-transformed hydraulic conductivity computed at times 1 T and 4 T in the vertical line at the expected value of the plume centroid.In the ES scenario, where the prior plumes are left free to evolve, the pdf tends to develop into a bimodal distribution.Instead, in the EnKF scenario, this behavior is prevented by the sequential updates and the final distribution, although still non-Gaussian, can be better approximated by a multivariate normal pdf.
A further explanation for the good performance of the EnKF can be found in its recursive application, through which, at each assimilation step, the plume is restarted with the same initial concentration in the updated Y fields.This is a known property of the restart EnKF: the state update is performed by invoking the governing equations with the updated parameters.Therefore, the parameter-state relation is always backed by the physics of the transport process (Wen and Chen, 2006;Hendricks Franssen and Kinzelbach, 2008;Nowak, 2009).This procedure is analogous to the advanced first-order second-moment (AFOSM) method adopted in risk analysis (Yen et al., 1986).As in AFOSM, the lack of Gaussianity is overcome by approaching recursively the solution with improved values of the estimator, here represented by the Kalman gain, which leads to improved estimates of the solution itself.

Scenarios with log-transformed concentration pdfs
In scenarios 2a and 2b, we assimilate log-transformed concentration values in order to evaluate if a proper transformation of the concentration pdf can improve the effectiveness of the assimilation algorithms.Although previous analyses highlight that beta-type pdfs can be effectively used to reproduce the concentration distribution (Caroni and Fiorotto,   2005), the log-normal pdf is easier to handle and often represents a reasonable assumption (e.g., Bellin et al., 1994;Zhang et al., 2000).
Figure 5 shows the Y fields estimated by the EnKF and the ES when assimilating log-transformed concentration measurements.The EnKF produces again a good reproduction of the true hydraulic conductivity field, with only small differences with the results of scenario 1a, while the ES confirms its difficulty in retrieving the Y spatial distribution.Nonetheless, the comparison between scenarios 1b and 2b shows a slight improvement of the ES solution, indicating that the ES is more sensitive than the EnKF to the pdf of the assimilated variable.

Scenarios with normal-score-transformed concentration pdfs
Since the log-transformation of concentration values does not ensure a normal pdf in all cases, in scenarios 3a and 3b another type of transformation is applied to the concentration data.Here a Gaussian pdf is obtained by applying a normal score transform (NST) (Zhou et al., 2011)  The results obtained by assimilating normal-scoretransformed concentration data are reported in Fig. 6 for both the EnKF and the ES.Despite the Gaussian distribution of the model variables, the retrieved fields are unsatisfactory for both techniques, and when compared to the results of scenarios 1a and 2a, even the EnKF performs poorly.2a) .
Since in scenario 3 the approximations related to the Gaussian assumption are removed, another reason must explain the poor performance of both EnKF and ES.The results seem to suggest that the NST might corrupt the cross-correlation structure between Y and C in the measurement locations.Indeed, when the C values are log-transformed, the relation that maps the original pdf of C to the transformed one is the same for all the different positions in space and thus the correlation structure is conserved.This is not the case for scenario 3, where a different NST is independently applied to the concentration ensemble at each node.
In order to illustrate this point, the Y − C cross-correlation structures, evaluated in the measurement nodes by means of the product P H T in Eq. ( 6), are compared for scenarios 1 and 3.The cross-correlation structures for scenario 1 at time t = 2 T and t = 4 T are reported in panels a and b of Figs. 7  and 8, respectively, while panels c and d of the same figures refer to scenario 3. The analysis considers only the longitudinal behavior of the cross correlation by reporting in the figures only the results calculated on the nodes aligned on the dotted line shown in Fig. 2. Each of the profiles shown in Figs.7 and 8 provides the cross correlation between Y at the node of interest, to which the profile corresponds, and the concentration values simulated at any given lag distance from that node.
In Figs. 7 and 8, we note that there are relevant differences between the EnKF and the ES cross-correlation structures.The EnKF cross-correlation behavior shows always higher values at the origin (zero lag) and a rapid decay with increasing lag, i.e., moving away from the measurement location.This confirms that the effectiveness of the EnKF is limited in a portion of the domain around the measurement location (Camporese et al., 2011).With the ES the peak of correlation at the measurement location is usually smaller than with the EnKF, and the cross-correlation structure is more spread out.As the cross correlation between Y and C is usually significant only for a limited lag distance, proportional to the product of the Y correlation scale and the length of the area covered by the plume, the relatively high correlation values characterizing the ES results at large lag values are probably spurious.For the same reason, the ES tends to produce a Y field more smoothed, and thus less accurate, than that estimated with the EnKF.Second, when the NST is applied (compare panels c and d with a and b in Figs.7 and 8), there is an overall decrease of the peaks, and the cross-correlation structures are even more smoothed, showing the significant alterations operated by the transformation.
In addition to the cross-correlation corruption described above, it must also be noted that the NST "fixes" only the marginal distributions, but not the multivariate dependence.Thus, there remains also a condition of non-multi-Gaussian dependency among all states and parameters.

Scenarios with modified normal score-transformed concentration pdfs
In order to maintain the original Y − C cross-correlation structure and, at the same time, to operate within the      assumption of Gaussianity required by both the EnKF and the ES, a modified application of the NST is proposed.At every time step, only one cumulative distribution function is built by using the concentration values simulated in all the measurement nodes.We underline that this is different from the previous application of the NST, in which a relation between C and its cdf is defined in each node independently.Now the cdf F (C), estimated at time t i , is , where m i is the number of measurement locations and i = 1, . . ., NMC × m i is the rank of the concentration values after sorting all data in ascending order.With this modified application of the NST we obtain a satisfactory reproduction of the Gaussian distribution in each node, and by using an invariant transformation, we manage not to alter the Y − C cross-correlation structure as much as in the previous scenario using the classic NST.This is shown in panels e and f of Figs. 7 and 8, where the cross-correlation structures are more similar to the original ones in scenario 1, especially at late times.
The results of the inversions, reported in Fig. 9 for both the EnKF and the ES, demonstrate the effectiveness of the modified NST as the estimated Y fields are now showing an improvement (with respect to the prior fields) comparable to that in scenario 2 and look much better than those of scenario 3 (Fig. 6).As in the previous scenarios, the EnKF outperforms the ES.The latter, however, shows significant improvements and seems to benefit more from the application of the modified NST.

Discussion
In Fig. 10, the results of all the scenarios are summarized and compared in terms of root mean square errors of both Y and C vs. time.All of the observations made in the previous sections, which were based merely on visual comparison, are confirmed by the RMSE profiles.The EnKF consistently outperforms the ES, regardless of the adopted concentration pdfs, except for scenario 3a, in which the NST deteriorates the EnKF solution due to alterations of the Y − C cross-correlation structure.This result, which reveals the inadequacy of the NST for this application, is in accordance with the conclusions drawn by Schoeniger et al. (2012).In their work, Schoeniger et al. (2012) applied the NST in conjunction with the EnKF to assimilate aquifer drawdown measurements and showed that the dependence between these data and the parameters is higher than the one between concentrations and parameters.The ES is more sensitive to the transformations operated to the assimilated data, as indicated by the RMSE of Y in the various scenarios, even though the same sensitivity is not reflected by the RMSE of C (Fig. 10).For all scenarios, with the exception of 2b (ES with logtransformed concentration data), the RMSE of C is consistently lower than in the open loop, i.e., a Monte Carlo simulation carried out with the prior ensemble of Y fields and no data assimilation.
Overall, the problem of retrieving the hydraulic conductivity field through the assimilation of concentration measurements is better handled by the EnKF due to the violation of Gaussianity investigated earlier with the different scenarios and to the high nonlinearity of the problem under consideration.With the EnKF, the Y fields are progressively updated and the simulated plumes gradually converge toward the true one.With the ES, the plumes evolve freely in the prior fields until the end of the simulation and, consequently, their evolution is very different from the true one, especially at late times.To highlight this point, the true plume distributions at t = 2 T and t = 4 T (Fig. 2) are compared with the corresponding ensemble means of the plume distribution simulated in the prior fields (Fig. 11a and b). Figure 11c-f also shows the evolution of the ensemble average of the plumes as simulated in the Y fields estimated at t = 2 T and t = 4 T , respectively, in scenario 1a (EnKF with original concentration pdfs).In other words, in Fig. 11c and d the EnKF is applied only until t = 2 T , and thereafter the plumes are left to evolve in the estimated Y fields.The progressive correction of the mean simulated plume is evident, and at t = 2 T it is already very similar to the true one.With the ES, instead, there are no recursive updates that modify the ensemble to resemble the true Y field, and the nonlinearity of the problem cannot be captured.
To further demonstrate our main findings, we simulated two additional scenarios, in which we changed the variance of the true Y field to 0.05 and 1.50 in order to decrease or increase, respectively, the non-Gaussianity and nonlinearity of the numerical experiments.We then carried out the inversions with both the EnKF and the ES with the same prior geostatistical parameters as in scenarios 1a, 1b, 4a, and 4b (Table 1).Y of up to 0.50, while some difficulties arise only for the strongly heterogeneous case, probably due to the lower variability in the initial ensemble (σ 2 Y = 0.75) that could not capture all the variability of the highly heterogeneous reference field.However, a RMSE decrease of 11 % is still achieved.These reductions are always larger than the corresponding values for the ES scenarios, confirming that the EnKF can handle conditions of nonlinearity and non-Gaussianity better than the ES does.The results also confirm that the ES is more sensitive than the EnKF to the transformation of the pdfs, with better performances consistently achieved by using the modified NST.Finally, the ES performance progressively worsens as the heterogeneity of the reference field increases, both with assimilation of untransformed and transformed data.

Summary and conclusions
The present work investigated the capabilities of the EnKF and the ES to retrieve the hydraulic conductivity spatial distribution through the assimilation of concentration measurements.The objective was to compare the performance of the two techniques and to analyze the effects of the lack of Gaussianity in the system variables.A tracer injection test was simulated in a two-dimensional domain representing a heterogeneous aquifer.Different scenarios were analyzed to determine how different transformations of the concentration probability distribution impact the inversion results.In the first scenario, concentration data were assimilated in the model without any manipulation, while in other three scenarios we considered a log-transformation of the data and two variants of the normal score transform (NST).In addition, we analyze also the cross-correlation structure between logtransformed hydraulic conductivity Y and concentration C.
The main conclusion of our study is that the EnKF can reproduce with good accuracy the hydraulic conductivity field and consistently outperforms the ES, regardless of the pdf of the concentrations.This is due to two reasons: (i) with the EnKF the lack of Gaussianity is overcome via the recursive Gaussian increments given by the EnKF updates, which eventually lead to an ensemble of members normally distributed around the true solution; (ii) the same recursive procedure progressively steers the ensemble of realizations toward the true solution, thus easing the inversion problem of the strong inherent nonlinearity of the dispersion process.The only case in which the EnKF does not work properly is when the NST is applied to ensure that the concentration pdf is Gaussian in each node of the domain.This  structure, which instead must be correctly evaluated in order to assure the effectiveness of the EnKF inversion procedure.This suggests that the NST must be applied with caution in any Kalman-filter-based inversion scheme by checking for possible corruptions of the cross correlation between parameters and assimilation variables.The ES performs always worse than the EnKF as it does not involve recursive updates of the Y fields.This has two consequences: (i) the solute plumes are free to evolve in the prior fields without corrections, eventually leading to significant differences from the true plume evolution; (ii) non-Gaussian contributions in the concentration pdf are not kept under control.

Fig. 1 .
Fig. 1.Spatial distribution of log-transformed hydraulic conductivity in the reference field.The color bar indicates log-transformed hydraulic conductivity Y in ln(L/T ).

Fig. 1 .
Fig. 1.Spatial distribution of log-transformed hydraulic conductivity in the reference field.The color bar indicates log-transformed hydraulic conductivity Y in ln(L/T ).

Fig. 2 .
Fig. 2. Plume evolution in the reference field at a) t = 2 T and b) t = 4 T .The color bar denotes dimensionless concentration C. Dotted lines show the direction along which the cross-correlation structure between Y and C is evaluated.

Fig. 2 .
Fig. 2. Plume evolution in the reference field at (a) t = 2 T and (b) t = 4 T .The color bar denotes dimensionless concentration C. Dotted lines show the direction along which the cross-correlation structure between Y and C is evaluated.

Fig. 3 .
Fig. 3. Scenario 1: log-transformed hydraulic conductivity field retrieved by a) the EnKF and b) ES using untransformed concentration data.The color bar indicates log-transformed hydraulic conductivity in ln(L/T ).Panels c) and d) show the estimation variance maps after the inversion for the EnKF and ES, respectively.

Fig. 3 .
Fig. 3. Scenario 1: log-transformed hydraulic conductivity field retrieved by (a) the EnKF and (b) ES using untransformed concentration data.The color bar indicates log-transformed hydraulic conductivity in ln(L/T ).(c) and (d) show the estimation variance maps after the inversion for the EnKF and ES, respectively.

12
Crestani et al.: EnKF vs ES for tracer test data assimilation

Fig. 4 .
Fig. 4. Scenario 1: empirical bivariate pdf of concentration and log-transformed hydraulic conductivity computed in the vertical line at the expected x value of the plume centroid for a) the EnKF at time 1 T , b) the EnKF at time 4 T , c) ES at time 1 T , and d) ES at time 4 T .

Fig. 5 .
Fig. 5. Scenario 2: log-transformed hydraulic conductivity field retrieved by a) the EnKF and b) ES using log-transformed concentration data.The color bar indicates log-transformed hydraulic conductivity in ln(L/T ).

Fig. 4 .
Fig. 4. Scenario 1: empirical bivariate pdf of concentration and log-transformed hydraulic conductivity computed in the vertical line at the expected x value of the plume centroid for (a) the EnKF at time 1 T , (b) the EnKF at time 4 T , (c) ES at time 1 T , and (d) ES at time 4 T .

Fig. 5 .
Fig. 5. Scenario 2: log-transformed hydraulic conductivity field retrieved by a) the EnKF and b) ES using log-transformed concentration data.The color bar indicates log-transformed hydraulic conductivity in ln(L/T ).

Fig. 6 .
Fig. 6.Scenario 3: log-transformed hydraulic conductivity field retrieved by (a) the EnKF and (b) ES using normal-score-transformed data.The color bar indicates log-transformed hydraulic conductivity in ln(L/T ).
to the C values.The NST is a tool through which any cumulative probability distribution function (cdf) F (x) is mirrored to the standard normal distribution cdf G(y).In other words, the generic variable x of the F (x) distribution can be transformed into the corresponding normally distributed variable y through the relation F (x) = G(y), i.e., y = G −1 [F (x)].In this case x = C and a cdf is built for each node of the domain with the ensemble of C values simulated by the model, using the Hazen formula F (C) = (i − 0.5)/NMC, where i = 1, . . ., NMC is the rank of the concentration values after sorting the data in ascending order.For each node, different values of C are univocally associated to the cdf values, which are always the same and depend only on NMC.

Fig. 8 .
Fig. 8. Y − C cross-correlation at t = 4 T for the nodes located at y = 5.75 L in scenario 1a (subpanel a), scenario 1b (subpanel b), scenario 3a (subpanel c), scenario 3b (subpanel d), scenario 4a (subpanel e), and scenario 4b (subpanel f).Each color corresponds to a correlation structure centered at a different node sampled by the plume (see Figure 2, subpanel b) Fig. 8. Y − C cross correlation at t = 4 T for the nodes located at y = 5.75 L in scenario 1a (a), scenario 1b (b), scenario 3a (c), scenario 3b (d), scenario 4a (e), and scenario 4b (f).Each color corresponds to a correlation structure centered at a different node sampled by the plume (see Fig. 2b).

Fig. 9 .
Fig. 9. Scenario 4: log-transformed hydraulic conductivity field retrieved by a) the EnKF and b) ES using concentration data transformed with a modified normal score transform.The color bar indicates log-transformed hydraulic conductivity in ln(L/T ).

Fig. 10 .
Fig. 10.Root mean square error of a) the retrieved log-transformed hydraulic conductivity field and b) the concentration distribution in the retrieved field for the all scenarios.

Fig. 9 .
Fig. 9. Scenario 4: log-transformed hydraulic conductivity field retrieved by (a) the EnKF and (b) ES using concentration data transformed with a modified normal score transform.The color bar indicates log-transformed hydraulic conductivity in ln(L/T ).

Fig. 9 .
Fig. 9. Scenario 4: log-transformed hydraulic conductivity field retrieved by a) the EnKF and b) ES using concentration data transformed with a modified normal score transform.The color bar indicates log-transformed hydraulic conductivity in ln(L/T ).

Fig. 10 .
Fig. 10.Root mean square error of a) the retrieved log-transformed hydraulic conductivity field and b) the concentration distribution in the retrieved field for the all scenarios.

Fig. 10 .
Fig. 10.Root mean square error of (a) the log-transformed hydraulic conductivity field and (b) the concentration distribution in the retrieved field for the all scenarios.

Fig. 11 .
Fig. 11.Ensemble mean of the plume evolution at a) t = 2 T and b) t = 4 T in the prior Y fields, at c) t = 2 T and d) t = 4 T in the Y field estimated at t = 2 T by the EnKF in scenario 1a and at e) t = 2 T and f) t = 4 T in the Y field estimated at t = 4 T by the EnKF in scenario 1a.The color bar indicates concentration values.

Fig. 11 .
Fig. 11.Ensemble mean of the plume evolution at (a t = 2 T and (b) t =4 T in the prior Y fields, at (c) t = 2 T and (d) t = 4 T in the Y field estimated at t = 2 T by the EnKF in scenario 1a and at (e) t = 2 T and (f) t = 4 T in the Y field estimated at t = 4 T by the EnKF in scenario 1a.The color bar indicates concentration values.

Table 1 .
Description, prior geostatistical parameters, and measurement errors of the numerical experiments.CV stays for the measurement coefficient of variation.In scenarios labeled with "a", the normal score transform is applied independently to the ensemble of concentration data in each assimilation node.In scenarios labeled with "b ", the normal score transform is applied to the ensemble of concentration data over all assimilation nodes.Crestani et al.: EnKF vs ES for tracer test data assimilation 11 Table 2 reports the inversion performance for all scenarios in terms of relative changes of Y RMSE computed as (RMSE t=4 − RMSE t=0 )/RMSE t=0 .With both assimilation strategies (assimilation of untransformed and modified NST-transformed C), the EnKF is able to decrease RMSE by almost 40 % for σ 2

Table 2 .
Inversion performance for scenarios 1a, 4a, 1b, and 4b (in terms of Y RMSE relative changes between prior and estimated solutions) as a function of the reference field Y variance.