Tong-chao Nan,Ji-chun Wu,*

aState Key Laboratory of Pollution Control and Resources Reuse,Nanjing University,Nanjing 210093,China

bDepartment of Hydrosciences,School of Earth Sciences and Engineering,Nanjing University,Nanjing 210093,China

Application of ensemble H-infinity filter in aquifer characterization and comparison to ensemble Kalman filter

Tong-chao Nana,b,Ji-chun Wua,b,*

aState Key Laboratory of Pollution Control and Resources Reuse,Nanjing University,Nanjing 210093,China

bDepartment of Hydrosciences,School of Earth Sciences and Engineering,Nanjing University,Nanjing 210093,China

Though the ensemble Kalman filter(EnKF)has been successfully applied in many areas,it requires explicit and accurate model and measurement error information,leading to difficulties in practice when only limited information on error mechanisms of observational instruments for subsurface systems is accessible.To handle the uncertain errors,we applied a robust data assimilation algorithm,the ensemble H-infinity filter(EnHF),to estimation of aquifer hydraulic heads and conductivities in a flow model with uncertain/correlated observational errors.The impacts of spatial and temporal correlations in measurements were analyzed,and the performance of EnHF was compared with that of the EnKF.The results show that both EnHF and EnKF are able to estimate hydraulic conductivities properly when observations are free of error;EnHF can provide robust estimates of hydraulic conductivities even when no observational error information is provided.In contrast,the estimates of EnKF seem noticeably undermined because of correlated errors and inaccurate error statistics,and filter divergence was observed.It is concluded that EnHF is an efficient assimilation algorithm when observational errors are unknown or error statistics are inaccurate.

Data assimilation;Hydraulic parameter estimation;Ensemble H-Infinity filter;Ensemble Kalman filter;Hydraulic conductivity;Robustness

1.Introduction

Water shortage and pollution are common problems worldwide.For instance,in China a recent national water quality survey has shown that only 38.7%of groundwater from wells meets quality criteria for source water supplies(Yu et al.,2014),and toxic organic chemicals threaten groundwater supplies and human health in the U.S.(National Research Council,2013).Aquifer characterization is always the first step to recognizing,managing,and protecting aquifers.Though a lot of pore-and fracture-scale studies manage to measure or characterize the exact geometries of porous or fractured media at small scales and simulate the detailed dynamics of flow and transport behaviors at such scales(Dou and Zhou,2014;Chen et al.,2014),to date,small-scale techniques are not applicable to field-scale problems in practice,like aquifer characterization and modeling,because of a lack of sufficient measurements.Unfortunately,due to considerable spatial variability in geology and a lack of information, affordable and accurate characterization of aquifer properties is still a challenge.Efficient and accurate estimation of parameters in groundwater models is always one of the most significant focuses of hydrology,since reliability and predictability of models are greatly dependent on model parameters.The parameter estimation or inverse problems insubsurface modeling are commonly considered ill-posed,and the“non-uniqueness of inverse problems should be addressed as uncertainty in the solutions”(Yeh et al.,2015)and may sometimes be handled via optimization approaches(Carrera and Neuman,1986).

In recent decades,filter methods used in data assimilation,which can combine and integrate different kinds of data,have drawn growing attention from hydrogeologists. These methods are capable of handling uncertainty with insufficient information in the framework of Bayes'theorem and improving the sensitivity of hydraulic variables(head, concentration,etc.)to hydraulic parameters(conductivity, dispersion coefficient,etc.)by introducing new kinds of data rather than simply adding more measurements of the same variable(McLaughlin and Townley,1996);these are the fundamental reasons for the popularity of filter methods.As one of the most popular data assimilation techniques,the ensemble Kalman filter(EnKF)has received broad recognition,and has been verified,improved,and applied in many areas,including hydrogeology(Evensen,2003;Chen and Zhang,2006;Franssen and Kinzelbach,2008;Sun et al., 2009;Nan and Wu,2011;Panzeri et al.,2013; Assumaning and Chang,2016).Applications show its capability of incorporating observations of various types from different sources to update the state of the model but expose its significant deficiencies as well.The performance of EnKF relies on accurate dynamic models of the system and on the premise that all the noises are Gaussian white noise whose statistical properties are clearly known.If the premise of accurate models or statistical properties of noises fails,the capability of the algorithm will deteriorate,and the method may even result in ridiculous system states(Sun et al.,2009;Nan and Wu,2011).In practice,many measurement noises are not white noise.They are often correlative in time or space or even change with time.Thus, a uniform error model is prohibitive in many situations, let alone in the cases when their statistical properties completely unknown(Wang et al.,2004).EnKF encounters big problems in practice when only limited information on error mechanisms of many observational instruments for subsurface systems is accessible.This is especially true for many newly developed hydrogeophysical techniques(Chung and Lin,2009;Tian et al.,2016).When water quality and security issues draw greater attention of the public,characterization of aquifer properties with a reliable risk control becomes a significant challenge,since variability of hydraulic properties and accuracy of parameter estimation are crucial for groundwater flow,solute transport,and incomplete mixing problems(Essaid et al.,2015;Tong et al., 2015).

The H-infinity(orH∞)filter,which is used in signal processing and control theory to achieve stabilization with guaranteed performance,has good robustness against uncertain system noises(Shaked,1990;Hassibi et al.,2000).It treats uncertain inputs and noises as random perturbations with limited energy and tries to minimize the H-infinity(H∞)norm of the transfer function from perturbations to estimation errors or to make it less than a given positive number.A lot of work has been performed to develop and investigate theories and algorithms of H-infinity in linear spaces or spaces which can be transformed into linear spaces(e.g.,Deng,2013; Yoneyama,2013;Zhang et al.,2014).Lu¨et al.(2010)used an H-infinity filter to estimate root zone water content by assimilating soil moisture data in a one-dimensional linearized Richards'equation.

To apply the H-infinity filter in nonlinear atmospheric systems,Han et al.(2009)combined the Monte Carlo method with the H-infinity filter and investigated the capability of the ensemble H-infinity filter(EnHF)in two synthetic data assimilation experiments.Luo and Hoteit(2011)proposed a time-local version of EnHF which utilized only the current state and observations of the system rather than the entire available history and found the equivalency between their algorithm and EnKF with covariance inflation.To the best of our knowledge,En HF has not been studied in hydrogeology or parameter estimation fields.Accounting for its robust performance in the case of insufficient information on model or measurement errors,this algorithm may turn out to be a powerful and practical tool for integrating hydrogeophysical data into groundwater models and be worthy of more attention from the hydrogeology community.In this paper,we first formulate groundwater flow models and the algorithms of En KF and EnHF,respectively.Afterwards,we examine the ability of EnHF to estimate aquifer properties with observations subject to uncertain(spatially and temporally correlated)errors,and compare it with that of the well-known En KF.

2.Methods

2.1.Groundwater flow simulation

Transient groundwater flow in confined and isotropic aquifers is considered to satisfy the following basic equation:

subject to the initial and boundary conditions:

where x is the spatial vector,tis time,K(x)is the hydraulic conductivity,h(x,t)is the pressure head,w(x,t)is the exchange between the element and outer space(source/sink term),Ssis the specific storage,H0(x)is the initial head distribution in the domain,HD(x,t)is the prescribed head on the Dirichlet boundary segments ΓD,Q(x,t)is the prescribed flux across the Neumann boundary segments ΓN,and n(x)is the outwardvector normal to the boundary ΓN.K(x)is treated as a spatial stochastic function and will be estimated by the two filter algorithms below.All flow models are simulated using MODFLOW-2000(Harbaugh et al.,2000),a widely used finite-difference groundwater flow simulator.

2.2.Ensemble Kalman filter

EnKF is a Monte Carlo method-based sequential assimilation method.It can modify both model variables and parameters included in the state vector according to observations. It can be briefly presented as follows:

2.3.Ensemble H-infinity filter

Considering the same model system handled by En KF,the state vector estimate can be written as follows,using observation vectors y from time step 0 to stepi:

The estimation error is

The transfer function matrixmaps the perturbation sequencesand the initial errorsto the estimation error sequenceThe key point of the H-infinity filter is to find a suitable estimator functionfthat minimizes the H-infinity norm of Ω(f),expressed as

where C0is the background error covariance matrix at the initial time.However,it is difficult to find the minimum value of the H-infinity norm for a complicated system in practice. Instead,a positive small number γ(called the level factor)is given,and the problem is turned into finding a suitable suboptimal filter that satisfies the following condition:

The H-infinity norm shows the maximum energy transferred from the input errors of the system and observations to the estimation errors.Ifis guaranteed to be minimal, the maximum energy from those perturbations is guaranteed to be minimal.Furthermore,if the total energy of the perturbations is fixed,the energy of the estimation error will be minimal,no matter the type and characteristics of the perturbations.This means that the H-infinity filter has good robustness against the errors.As γ has smaller values,the method becomes more robust(Khargonekar and Nagpal,1991; Shaked and Theodor,1992).

Constructing a filter satisfying Eq.(13)is actually an optimization problem.Hassibi et al.(2000)proposed a matrix square root algorithm of the H-infinity filter based on matrix triangularization.According to this algorithm,the gain matrixKHcan be evaluated by triangularizing a core matrix B with hyperbolic transformation:

The matrix square root algorithm of EnHF has several strong advantages.First,it is computationally stable because of avoidance of a loss of positive definiteness that is often induced by rounded-off errors in numerical approximation of sample covariance.Second,the condition for successfully implementing the algorithm is exactly the condition required for the existence of the filter for a given level factor γ,which means that the filter solutions exist if and only if the algorithm can be executed for a given γ.The direct criterion is the existence of hyperbolic transformation of the core matrix B.Last but not least,it can be easily implemented due to straightforward array operations,especially with array-friendly tools like MATLAB andRlanguages,with which basic array operations can be programmed with simple codes and run at a higher speed.

Once the gain matrix KHis found,the state vector can be updated as follows:

It is noteworthy that Hassibi et al.(2000)found that when,one has

This is equal to the Kalman gain matrix associated with white Gaussian noise of unit variance.As mentioned above,a large γ implies low robustness.Hence,Eq.(18)suggests that, compared to the(ensemble)H-infinity filter,the(ensemble) Kalman filter is not robust against uncertain perturbations, which is numerically demonstrated further below.

3.Model setup and tools

To focus on the main differences between the two algorithms and avoid distractions by other factors(e.g.,model errors,singular covariance matrices,and so on),a simple and well-controlled dynamic simulator is preferable.For this purpose,a synthetic groundwater flow example was constructed. The flow model was a horizontal confined aquifer with a length of 100 m,which was uniformly discretized into 100 elements. The thickness of the aquifer was 10 m.The left(x=0)and right(x=100 m)ends were assumed to be Dirichlet boundaries with heads of 100 m and 90 m,respectively.A pumping well with a volumetric flow rate of 0.03 m3/d was located at 30 m,and an injection well with the same rate was located at 70 m from the left end.Specific storage was assumed to be a constant equal to 0.00001.The initial heads were assumed to be completely known.

The natural logarithm of the hydraulic conductivity fieldY(x)=lnK(x)was assumed to be a Gaussian random function with a mean of-4 and a standard deviation of 0.5,and to be second-order stationary and characterized by the exponential covariance function with a correlation length equal to 10 m. We used the sequential Gaussian simulation algorithm in GSLIB to generate a set of logarithmic hydraulic conductivity fields(Deutsch and Journel,1998)with the statistics mentioned above and arbitrarily selected one as the true reference field,shown in Fig.1.In this field,the maximum and minimumYvalues were equal to-3.05 and-5.55,respectively,and the mean was equal to-4.04.The initial randomYfields were generated using the same algorithm with the same statistics.We set nine points(x=10,20,…,90)in the domain as the measurement locations forYandh.

For this groundwater model,the flow nearly reached a steady state after about 1.0 d.The simulation period wast=3.0 d,and the time duration was evenly divided into 30 steps with each step of 0.1 d.As mentioned above,all the flow models were subject to the boundary conditions described above and simulated via the same simulator,MODFLOW-2000,so the data assimilating process was free of model error.EnHF and EnKF.Furthermore,this reveals that,through a wider 95%confidence interval,EnHF provides more conservative estimation results of the model state,but has much safer control of system uncertainty.

Fig.1.True Y reference field.

Fig.2.Y estimates and 95%con fidence intervals from EnHF and EnKF in Experiment 1 and comparison of Y estimates with true Y field.

Fig.3 indicates thatRMSEandWvalues ofYfrom EnHF and EnKF in Experiment 1 decrease rapidly in the first several steps and then tend to stabilize.This also shows thatWevolutions in EnHF and EnKF are highly similar whereas the correspondingRMSEcurves look different.With the very first update,theRMSEvalue in EnKF plunges to the value very close to its final value and nearly remains unchanged after the third update.In contrast,the decrease of theRMSEvalue in EnHF seems persistent and less abrupt.This suggests that EnHF tends to update the model state in a conservative manner because it does nottrustthe observations completely.This also leads to the fact that theRMSEvalue of 0.24 forYat the last step from EnHF is slightly higher than the corresponding value of 0.21 from EnKF.EnKF obtains nothing new from data of a near steady state,while EnHF seems to confirm thepartly trustedold data using the steady-state data and update the model state slowly and persistently.

It is noteworthy that EnHF led to a remarkably higher computational cost than EnKF in this study.EnHF spent about 1300 s updating the ensemble of 200 realizations 30 times, whereas EnKF spent only 560 s in an identical environment.

Fig.4 shows minimum γ values in EnHF as solutions to Eq.(13)at all time steps in Experiment 1.It can be seen that the minimum γ value gets smaller and tends to stabilize as assimilation proceeds.

Fig.3.RMSE and W curves of Y from EnHF and EnKF in Experiments 1 through 6.

Fig.4.Minimum γ value in EnHF at each time step in Experiment 1.

4.1.Assimilating observations with spatially correlated errors

In this phase,observations were generated by adding errors to the trueYandhvalues,which constitute Gaussian stochastic processes with a mean of zero and standard deviations of 0.05 and 0.3 m forYandh,respectively.In order to explore the impact of spatial correlations in observational errors,we carried out five experiments(Experiments 2 through 6)with EnHF and EnKF,in which the correlation lengths of the observational errors were different.We assumed that errors inYobservations were uncorrelated to those inhobservations, and that the autocovariance of errors inYandhobservations could be characterized by the exponential covariance function:

wheredis the separation distance,indicating the distance between two spatial points being considered;and σeis the standard deviation of errors.Again,GSLIB(Deutsch and Journel,1998)was used to generate a set of the observational errors in Experiments 2 through 6,witha=0 m,10 m, 50 m,100 m,and 500 m,respectively,as shown in Table 1. When these noisy observations were used,EnKF treated observational errors as white noise with the standard deviations of 0.01 forYand 0.1 m forh,representing the situation in which only inaccurate error statistics are available.No information about observational errors was given to EnHF.

Fig.3 shows evolutions ofRMSEandWofYin EnHF and En KF over all time steps.Four interesting facts are noted:(1) comparing Fig.3(a)with(b),we can see that theRMSEcurves from EnHF stay close to that of the control experiment(with perfect observations),while theRMSEcurves from EnKF oscillate wildly after the first several steps;(2) Fig.3(c)and(d)show that theWcurves from EnHF remain almost the same as that of the control experiment,while theWcurves from En KF consistently decrease and go far below that in the control experiment;(3)Fig.3(b)indicates thatRMSEvalues in En KF generally get higher as the correlation lengthaincreases,suggesting a stronger impact of errors on the EnKF performance;and(4)Wvalues in EnHF are always larger thanRMSEvalues,but in EnKF,Wvalues may get much smaller thanRMSEvalues and tend to diminish,a phenomenon typically referred to as filter divergence.Note that the evolution ofWin EnKF relies on error statistics recognized or used by En KF,instead of the actual error.That is why theWcurves in Experiments 2 through 6 are more or less the same but different from that in Experiment 1 in Fig.3(d).

In Fig.5,theRMSEandWcurves ofYfrom EnHF and EnKF in Experiments 1 and 6 are compared.It shows clearly that,although EnKF converges with perfect observations and produces slightly better estimates than EnHF,EnKF still diverges significantly,with estimates getting worse but ensemble uncertainty diminishing,when perturbed by observation noises of long-distance correlation.In contrast,theRMSEcurves ofYfrom EnHF always change mildly and stay below or close to the correspondingWcurves,indicating the insensitivity of the EnHF performance to noises.This demonstrates that EnHF has good robustness against spatially correlated errors in observations and produces better estimates with the existence of observational errors.

Though in this study we focused on parameter estimation using the two algorithms,the variations in the model variablehreflected the behavior of the algorithms and were investigated as well.Fig.6 shows theRMSEandWcurves ofhfrom EnHF and EnKF.Note that all the initial values ofh(initial conditions)in Monte Carlo simulation were perfect and identical to those in the reference model.That is whyRMSEandWvalues ofhare zero at the first step for both algorithms. In addition to that,similar phenomena,e.g.,the low sensitivityof EnHF to noises and relatively strong impact of errors on EnKF results,can be observed in Fig.6.

Fig.5.Comparisons of RMS E and W curves of Y from EnHF and EnKF in Experiments 1 and 6.

Fig.6.RMSE and W curves of h from EnHF and EnKF in Experiments 1 through 6.

Fig.7.Y estimates and 95%con fidence intervals from EnHF and EnKF in Experiment 6 and comparison of Y estimates with true Y field.

The results above indicate that the existence of observational errors and the growing correlation length do not have significant influence onRMSEin EnHF.However,EnKF is prone to perturbations and may diverge and lead to misleading estimates of the model state.EnHF exhibits much better robustness against spatially correlated perturbations in observations than EnKF in estimating the hydraulic conductivity.

Again,we calculated the 95%confidence intervals forYestimates in all experiments.Fig.7 demonstratesYestimates (ensemble means)and the 95%confidence intervals from EnHF and EnKF in Experiment 6.The estimated mean and the 95% confidence interval evaluated by EnHF in Fig.7(a)look extremely similar to those of the control experiment in Fig.2(a), and only 4 outof101 trueYvalues are located outside ofthe 95% confidence interval.In contrast,the estimated mean and the 95% confidence interval evaluated by EnKF in Fig.7(b)are notably differentfromthose ofthe controlexperimentin Fig.2(b),and 43 out of 101 trueYvalues are outside of the 95%confidence interval.This suggests thatEnHF reduces parameter uncertainty in a safer and more conservative manner,and thus has no need for statistical information of observational errors.Unlike EnHF, EnKF was designed from the very beginning to eliminate parameter uncertainty as much as possible,and it may also foul up due to noisy observations and their inaccurate statistics, leading to misleading estimates.Results of other experiments demonstrated similar phenomena and were omitted here for the sake of brevity.

4.2.Assimilating observations with temporally correlated errors

For dynamic groundwater observational data that are sequentially collected via the same measurement instruments or methods,the observational errors are often temporally correlated.The influences of temporally correlated errors in observations on estimation results of EnHF and EnKF were investigated.

To generate the observational error vector eO(i)at theith time step,we used the equations below:

where e(i)is the independent and identically distributed Gaussian white noise with a mean of zero and variance equal toWe selectedto control the extent of the correlation of observational errors.As α increases,the temporal correlation of observational errors gets stronger.When α=0,observational errors are uncorrelated;when α=1, observation errors are time-invariant.

In order to investigate the impact of α on the performances of EnHF and EnKF,five parallel assimilation experiments (Experiments 7 through 11 in Table 1)were conducted,in which α=0.1,0.3,0.5,0.7,and 0.9,respectively.Fig.8 shows theRMSEandWcurves ofYfrom EnHF and EnKF in the five experiments as well as Experiment 1.For EnHF,theRMSEcurves in Fig.8(a)andWcurves in Fig.8(c)remain unchanged for different α values.In contrast,Fig.8(b)reveals that there is an α-dependent increase following an early drop inRMSEin EnKF.The larger the α value is,the higher theRMSEgrows, except for α=0.9,which corresponds to the second highestRMSEcurve rather than the highestRMSEcurve.This suggests that the temporal correlation in observational errors also impacts the performance of EnKF and exacerbates the filter divergence.Similar to Fig.3(d),Fig.8(d)demonstrates the consistency ofWvalues ofYfrom EnKF in Experiments 7 through 11 but the difference from that in the control experiment as well.As discussed above,the evolution ofWrelies on error statistics that are recognized or used by EnKF,rather than the actual errors or error statistics(including α).Hence,it is not a surprise to find thatWin Fig.8(d)does not vary with α.

Fig.9 shows theRMSEandWcurves ofhfrom EnHF and EnKF in Experiment 1 and Experiments 7 through 11,which are similar to what we see in Fig.6.Different α values make significant differences in the results of EnKF whereas the impact of α on the performance of En HF is limited and consistent in bothRMSEandW.Furthermore,EnKF suffers from filter divergence while EnHF works properly in all experiments for any α value.

Fig.10 presentsYestimates(ensemble means)and 95% confidence intervals from En HF and EnKF in Experiment 11. Only 4 out of 101 trueYvalues are located outside of the 95% confidence interval evaluated by EnHF(Fig.10(a)),but 37 out of 101 trueYvalues are located outside of the 95%confidence interval evaluated by EnKF(Fig.10(b)).The narrow conf idence interval,the large number of outliers,and the deviationof the black solid curve from red dots in Fig.10(b)demonstrate filter divergence occurring in EnKF.

Fig.8.RMSE and W curves of Y from EnHF and EnKF in Experiment 1 and Experiments 7 through 11.

Fig.9.RMSE and W curves of h from EnHF and EnKF in Experiment 1 and Experiments 7 through 11.

Fig.10.Y estimates and 95%con fidence intervals from EnHF and EnKF in Experiment 11 and comparison of Y estimates with true Y field.

5.Conclusions

In this study,a data assimilation method called EnHF based on the H-infinity filter and ensemble forecasting was introduced into aquifer parameter estimation.The ability of EnHF to estimate the aquifer hydraulic conductivity with observations subject to uncertain(inaccurate/spatially/temporally correlated)errors was numerically investigated,and compared with that of EnKF,a well-known approach.The main conclusions are as follows:

(1)The Kalman filter(and EnKF)was developed to minimize the resulting variance of system states through a variance weighting strategy(Evensen,2003),and is thus prone to perturbations or inaccurate error information,which was verified in our numerical experiments.

(2)EnHF minimizes the maximum energy transferred from the input errors of the system and observations to the estimation errors,which causes it to update the model state in a conservative and persistent manner.

(3)When accurate and perfect observations are used,EnKF and EnHF yield similar estimation results;the computation cost seems to be the main advantage of EnKF in such situation.

(4)When observations are subject to spatially/temporally correlated noises,EnHF maintains proper performance, demonstrating good robustness against noises.On the otherhand,the performance of EnKF deteriorates with filter divergence occurring,and stronger correlation in observational errors leads to worse EnKF estimates and more serious filter divergence.

(5)EnHF is a promising tool for data assimilation,especially where there is insufficient,inaccurate,or incomplete observational error information.Hopefully,its capability and applicability to assimilation of hydrogeophysical data will be explored and recognized by the hydrogeology community in the future.This will require further investigation and deeper thought.A study of two-dimensional/three-dimensional synthetic models and an application to the field will be our next step.

Acknowledgements

We would like to thank Dr.Yue-qi Han from the Institute of Meteorology,PLA University of Science and Technology, for sharing with us the details of his work in Han et al. (2009).

Assumaning,G.A.,Chang,S.Y.,2016.Application of sequential dataassimilation techniques in groundwater contaminant transport modeling. J.Environ.Eng.142(2),01015073.http://dx.doi.org/10.1061/(ASCE) EE.1943-7870.0001034.

Carrera,J.,Neuman,S.P.,1986.Estimation of aquifer parameters under transient and steady state conditions:2.Uniqueness,stability,and solution algorithms.Water Resour.Res.22(2),211-227.http://dx.doi.org/10.1029/ WR022i002p00211.

Chen,L.,Kang,Q.J.,Mu,Y.T.,He,Y.L.,Tao,W.Q.,2014.A critical review of the pseudopotential multiphase lattice Boltzmann model:Methods and applications.Int.J.Heat Mass Transf.76,210-236.http://dx.doi.org/ 10.1016/j.ijheatmasstransfer.2014.04.032.

Chen,Y.,Zhang,D.X.,2006.Data assimilation for transient flow in geologic formations via ensemble Kalman filter.Adv.Water Resour.29(8), 1107-1122.http://dx.doi.org/10.1016/j.advwatres.2005.09.007.

Chung,C.C.,Lin,C.P.,2009.Apparent dielectric constant and effective frequency of TDR measurements:Influencing factors and comparison. Vadose Zone J.8(3),548-556.http://dx.doi.org/10.2136/vzj2008.0089.

Deng,Z.H.,2013.Robust finite-time H-infinity filtering for uncertain systems subject to missing measurements.J.Inequal.Appl.236.http://dx.doi.org/ 10.1186/1029-242X-2013-236.

Deutsch,C.,Journel,A.,1998.GSLIB:Geostatistical Software LIBrary and User's Guide,second ed.Oxford University Press,New York.

Dou,Z.,Zhou,Z.F.,2014.Lattice Boltzmann simulation of solute transport in a single rough fracture.Water Sci.Eng.7(3),277-287.http://dx.doi.org/ 10.3882/j.issn.1674-2370.2014.03.004.

Essaid,H.I.,Bekins,B.A.,Cozzarelli,I.M.,2015.Organic contaminant transport and fate in the subsurface:Evolution of knowledge and understanding.Water Resour.Res.51(7),4861-4902.http://dx.doi.org/10.1002/ 2015WR017121.

Evensen,G.,2003.The ensemble Kalman filter:Theoretical formulation and practical implementation.Ocean.Dyn.53(4),343-367.http://dx.doi.org/ 10.1007/s10236-003-0036-9.

Franssen,H.J.H.,Kinzelbach,W.,2008.Real-time groundwater flow modeling with the ensemble Kalman filter:Joint estimation of states and parameters and the filter inbreeding problem.Water Resour.Res.44(9),W09408. http://dx.doi.org/10.1029/2007WR006505.

Han,Y.Q.,Zhang,Y.C.,Wang,Y.F.,Ye,S.,Fang,H.X.,2009.A new sequential data assimilation method.Sci.China Ser.E Technol.Sci.52(4), 1027-1038.http://dx.doi.org/10.1007/s11431-008-0189-3.

Harbaugh,A.W.,Banta,E.R.,Hill,M.C.,McDonald,M.G.,2000.MODFLOW-2000,the U.S.Geological Survey Modular Ground-water Model: User Guide to Modularization Concepts and the Ground-water Flow Process,Open-File Rep 00-92.U.S.Geological Survey,Reston.

Hassibi,B.,Kailath,T.,Sayed,A.H.,2000.Array algorithms forH∞estimation.IEEE Trans.Autom.Control 45(4),702-706.http://dx.doi.org/ 10.1109/9.847105.

Khargonekar,P.P.,Nagpal,K.M.,1991.Filtering and smoothing in anH∞setting.IEEE Trans.Autom.Control 36(2),152-166.http://dx.doi.org/ 10.1109/9.67291.

Lu¨,H.S.,Li,X.L.,Yu,Z.B.,Horton,R.,Zhu,Y.H.,Hao,Z.C.,Xiang,L.,2010. Using a H-infinity filter assimilation procedure to estimate root zone soil water content.Hydrol.Process.24(25),3648-3660.http://dx.doi.org/ 10.1002/hyp.7778.

Luo,X.D.,Hoteit,I.,2011.Robust ensemble filtering and its relation to covariance inflation in the ensemble Kalmanfilter.Mon.Weather Rev.139, 3938-3953.http://dx.doi.org/10.1175/MWR-D-10-05068.1.

McLaughlin,D.,Townley,L.R.,1996.A reassessment of the groundwater inverse problem.Water Resour.Res.32(5),1131-1161.http://dx.doi.org/ 10.1029/2006WR005144.

Nan,T.C.,Wu,J.C.,2011.Groundwater parameter estimation using the ensemble Kalman filter with localization.Hydrogeol.J.19(3),547-561. http://dx.doi.org/10.1007/s10040-010-0679-9.

National Research Council,2013.Alternatives for Managing the Nation's Complex Contaminated Groundwater Sites.The National Academies Press,Washington,D.C.

Panzeri,M.,Riva,M.,Guadagnini,A.,Neuman,S.P.,2013.Data assimilation and parameter estimation via ensemble Kalman filter coupled with stochastic moment equations of transient groundwater flow.Water Resour. Res.49(3),1334-1344.http://dx.doi.org/10.1002/wrcr.20113.

Rubin,Y.,Hubbard,S.S.,2007.Hydrogeophysics.Springer,Dordrecht.

Shaked,U.,1990.H∞minimum error state estimation of linear stationary processes.IEEE Trans.Autom.Control 35(5),554-558.http://dx.doi.org/ 10.1109/9.53521.

Shaked,U.,Theodor,Y.,1992.H∞optimal estimation:A tutorial.In:Proceedings of the 31st IEEE Conference on Decision and Control,IEEE, Tucson,pp.2278-2286.

Sun,A.Y.,Morris,A.P.,Mohanty,S.,2009.Comparison of deterministic ensemble Kalman filters for assimilating hydrogeological data.Adv.Water Resour.32(2),280-292.http://dx.doi.org/10.1016/j.advwatres.2008. 11.006.

Tian,Z.C.,Li,Z.Z.,Liu,G.,Li,B.G.,Ren,T.S.,2016.Soil water content determination with cosmic-ray neutron sensor:Correcting aboveground hydrogen effects with thermal/fast neutron ratio.J.Hydrol.540,923-933. http://dx.doi.org/10.1016/j.jhydrol.2016.07.004.

Tong,J.X.,Yang,J.Z.,Hu,B.X.,2015.Analysis of soluble chemical transfer from soil to surface runoff and incomplete mixing parameter identification. Water Sci.Eng.8(3),217-225.http://dx.doi.org/10.1016/j.wse.2015. 04.011.

Wang,Y.F.,Wang,B.,Han,Y.Q.,Zhu,M.,Hou,Z.M.,Zhou,Y.,Liu,Y.D., Kou,Z.,2004.Variational data assimilation experiments of Mei-yu front rainstorms in China.Adv.Atmos.Sci.21(4),587-596.http://dx.doi.org/ 10.1007/BF02915726.

Yeh,T.J.,Mao,D.Q.,Zha,Y.Y.,Wen,J.C.,Wan,L.,Hsu,K.C.,Lee,C.H., 2015.Uniqueness,scale,and resolution issues in groundwater model parameter identification.Water Sci.Eng.8(3),175-194.http://dx.doi.org/ 10.1016/j.wse.2015.08.002.

Yoneyama,J.,2013.Robust H-infinity filtering for sampled-data fuzzy systems.Fuzzy Sets Syst.217,110-129.http://dx.doi.org/10.1016/ j.fss.2012.08.014.

Yu,Z.B.,Yang,T.,Schwartz,F.W.,2014.Water issues and prospects for hydrological science in China.Water Sci.Eng.7(1),1-4.http://dx.doi.org/ 10.3882/j.issn.1674-2370.2014.01.001.

Zhang,W.A.,Dong,H.,Guo,G.,Yu,L.,2014.Distributed sampled-data H-infinity filtering for sensor networks with nonuniform sampling periods.IEEE Trans.Ind.Inf.10(2),871-881.http://dx.doi.org/10.1109/ TII.2014.2299897.

Received 11 October 2016;accepted 2 January 2017

Available online 14 March 2017

This work was supported by the National Natural Science Foundation of China(Grant No.41602250)and the Project of the China Geological Survey (Grant No.DD20160293).

*Corresponding author.

E-mail address:jcwu@nju.edu.cn(Ji-chun Wu).

Peer review under responsibility of Hohai University.

http://dx.doi.org/10.1016/j.wse.2017.03.009

1674-2370/?2017 Hohai University.Production and hosting by Elsevier B.V.This is an open access article under the CC BY-NC-ND license(http:// creativecommons.org/licenses/by-nc-nd/4.0/).

?2017 Hohai University.Production and hosting by Elsevier B.V.This is an open access article under the CC BY-NC-ND license(http:// creativecommons.org/licenses/by-nc-nd/4.0/).