Application of the Multimodel Ensemble Kalman Filter Method in Groundwater System

With the development of in-situ monitoring techniques, the ensemble Kalman filter (EnKF) has become a popular data assimilation method due to its capability to jointly update model parameters and state variables in a sequential way, and to assess the uncertainty associated with estimation and prediction. To take the conceptual model uncertainty into account during the data assimilation process, a novel multimodel ensemble Kalman filter method has been proposed by incorporating the standard EnKF with Bayesian model averaging framework. In this paper, this method is applied to analyze the dataset obtained from the Hailiutu River Basin located in the northwest part of China. Multiple conceptual models are created by considering two important factors that control groundwater dynamics in semi-arid areas: the zonation pattern of the hydraulic conductivity field and the relationship between evapotranspiration and groundwater level. The results show that the posterior model weights of the postulated models can be dynamically adjusted according to the mismatch between the measurements and the ensemble predictions, and the multimodel ensemble estimation and the corresponding uncertainty can be quantified.


Introduction
Groundwater resources are crucial to the development of human society and the sustainability of the environment.The rational management and effective protection of groundwater require accurate characterization of hydrogeological conditions.Due to the complexity of the groundwater system, various uncertainties exist and limit our ability to wisely manage the groundwater resources [1,2].The uncertainty associated with the groundwater system mainly originates from the spatial heterogeneity of the hydrogeological properties.It can be quantitatively described by using geostatistical methods, with the capability to constrain the spatial distribution of hydrogeological properties on direct measurements [3].To further improve the characterization of hydrogeological conditions, inverse modeling methods have been introduced to constrain the hydrogeological characterization on the available measurements of the state variables, such as hydraulic head, flow rate and solute concentration, and reduce the corresponding uncertainty associated with the characterization [4][5][6][7][8][9].
There has been a growing tendency to take conceptual hydrological model uncertainty into account in the past decade.The existence of conceptual uncertainty has been demonstrated in geoscience [10].If the conceptual uncertainty is neglected, two types of errors may occur: first, the prediction can be statistically biased due to the adoption of an invalid conceptual model; and second, the uncertainty and risk associated with the prediction can be underestimated due to insufficient sampling of the conceptual model space [11].The Generalized Likelihood Uncertainty Estimation (GLUE) pioneered the methods that can account for conceptual model uncertainty explicitly [12][13][14].A more statistically coherent multimodel analysis method is the Bayesian Model Averaging method [15][16][17][18].The GLUE and BMA (Bayesian Model Averaging) methods are combined together by [19,20].The Maximum Likelihood version of BMA (MLBMA) is proposed to incorporate the groundwater model calibration methods based on maximum likelihood estimation into the multimodel analysis framework more conveniently [11,[21][22][23][24].
With the development of in-situ monitoring technology in groundwater management, the ensemble Kalman filter method (EnKF) has gained significant attention due to its capability to update the model parameter and state simultaneously in a sequential manner [25][26][27].Recently, we developed a multimodel data assimilation method by embedding the traditional EnKF into the BMA framework [28].The novel multimodel EnKF method offers all of the advantages of the traditional EnKF method, and can also take conceptual model uncertainty into consideration during the data assimilation process.The methodology of the multimodel EnKF has been thoroughly described and the accuracy of the proposed method is quantitatively justified using the statistical criteria in [28].This work focuses on demonstration of the multimodel EnKF analysis procedures when applying them to a more practical case study.In this case, the meteorological, topographic and geological data are practically obtained from the Hailiutu River area in the Ordos Basin, China.Limited by the dataset of the detailed hydrogeological characterization and dynamical hydraulic head observation, we adopt a strategy similar to [29], which utilizes a synthetic hydraulic head observation data computed from a generated full heterogeneous reference field conditional on the available direct hydraulic conductivity measurements; therefore, this study is essentially synthetic.However, the benefits of this strategy are that: (1) the application of the proposed method can be easily extended to other practical cases by imitating the example in this research; (2) it can reflect available field conditions, such as meteorological, topographic and geological conditions, to the maximum extent; and (3) the performance of the proposed method to estimate the hydraulic conductivity field can be evaluated by comparing it against the generated reference hydraulic conductivity field.
The rest of this paper is arranged as follows.In Section 2, the multimodel EnKF method is briefly presented, and the study area, the formulation of the multiple models, and the method implementation are demonstrated in detail.In Section 3, the analysis results are discussed.The main conclusions are summarized in Section 4.

Theoretical Background
The general multimodel ensemble Kalman filter method is detailed in [28].Here, the implementation of the multimodel EnKF method is briefly introduced in a specialized manner for the groundwater modeling problem, in which the dynamic hydraulic head observations are assimilated to estimate the hydraulic conductivity field.
Suppose that a set M of Nk conceptual models, , are postulated to simulate the groundwater system.For a given model Mk, the data assimilation process consists of two steps: the forecast step and the update step.
For the convenience of the method presentation, Ne realizations of state vector, s, are used here, which includes hydraulic conductivity, hydraulic head, and predicted values at the measurement locations: , , , where ln K = Y is the natural logarithm of hydraulic conductivity K ; h is the hydraulic head; and m is the predicted values at the measurement locations.Usually, the hydraulic conductivity and hydraulic head in the state vector are updated simultaneously in the EnKF method.However, the updated hydraulic head and hydraulic conductivity are not consistent with each other if the nonlinearity of the groundwater model is strong.To obviate this issue, only the hydraulic conductivity is updated during each assimilation step, and the corresponding hydraulic head is predicted by: where t ξ is the model prediction error, with the assumption that ( ) where the superscript obs represents the observed value; and the superscript true represents the true (usually unknown) status; t ε is the Gaussian-distributed measurement error with ( ) ( ) , ζ is the added perturbation error [30], which is usually assumed to have the same distribution as the measurement error, i.e., ( ) , , , ; T is the transpose operator; and H is a matrix operator which only contains 0 and 1 elements, and can be written as: where I is the identity matrix.
The ensemble mean and covariance in the forecast step are: The Kalman gain at time t can be written as: The updated ensemble member of the state vector is: The ensemble mean and covariance of the updated state vector, conditional on all of the available observation data up to time t and the given model, are: ( ) where 1: is the observed hydraulic head data up to time t.
Bayesian model averaging framework provides a coherent way to deal with the situation in which multiple conceptual models are postulated.It shows that the posterior distribution of the predictions, conditional on all of the available data, is: This equation indicates that the posterior distribution of the predictions after averaging the model space is the weighted summation of the posterior distribution based on individual models.Now the key issue is to evaluate the posterior model weight.
The posterior model weight can be computed through Bayes' theorem as: The marginal likelihood can be evaluated through the Monte Carlo integration method as: Several types of likelihood functions can be found in [19].Here, we use the form: | , where E is a user-defined parameter (we choose E = 1 in this work).
The posterior mean and covariance of the estimated hydraulic conductivity, after model averaging, are, respectively: The posterior mean and covariance of the predicted hydraulic head, after model averaging, can be obtained simply by replacing , 15) and ( 16).

Study Area
The study site, the Hailiutu River Basin, is located in a semi-arid region of northwest China as shown in Figure 1, covering an area of about 2600 km 2 between 38°06' N and 38°50' N and 108°37' E and 109°15' E. Groundwater level fluctuation is mainly driven by seasonally varying recharge and evapotranspiration.Based on the geological investigation information, two major layers in the aquifer system are identified: the overlying unconfined layer consists of Quaternary Shala-Usun fine sand, with thickness between 0.5 and 217 m; and the underlying confined layer consists of Cretaceous Luohe sandstone, with thickness between 20 and 637 m (see the elevation of the studied aquifer system in Figure 2).No regional aquitard has been found between these two aquifer layers.The Hailiutu River flows mainly through the upper layer and cuts into the second layer near the end of the river system.

Model Setup
The groundwater simulation software, MODFLOW, is used to build the numerical model of the study site.The study area is discretized into a 200 × 140 grid blocks with two layers, and each block represents a 500 m × 500 m area.The boundary of the study area is delineated based on the watershed divide, and it is treated as a no-flow boundary.The cells out of the boundary are characterized as inactive cells, as shown by the gray area in Figure 3.The number of the active cells is 20,840.The upper layer is an unconfined aquifer with a specific yield of 0.1, and the bottom layer is a confined aquifer with a specific storage of 1 × 10 −5 m −1 .The discretized grid for each layer is shown in Figure 3.The yellow line in Figure 3 represents the location of the river.The riverbed conductance is taken to be 300 m 2 /d for the part of the river flowing through the upper layer and 150 m 2 /d for the part of the river flowing through the bottom layer.The initial condition of this groundwater model is taken as the simulation result at the end of year 2006, which is obtained from a previous base flow analysis model.The simulation period is 3 years (field measurements of precipitation and pan evaporation from the beginning of 2007 to the end of 2009 are used), which are divided into 36 stress periods with the time step as 1 day, i.e., each stress period corresponds to a specific month.The differences of stress periods lie in the recharge and evapotranspiration.The precipitation (measured by rain gauge) and evaporation (measured by evaporation pan with a diameter of 20 cm) rate during the simulation period are depicted in Figure 4.The pan evaporation is converted to the potential or maximum evapotranspiration by multiplying a conversion constant.Based on a previous study, the conversion constant is set to be 0.3.The recharge coefficient of this area is set to be 0.32, and the evapotranspiration extinction depth is 4.6 m in this area.

Synthetic Data Generation
Owing to the lack of sufficient and reliable real dynamic head observation data, a strategy similar to [29] is adopted, i.e., a series of simulated head observation values are used for the analysis.The head observation values are computed based on the synthetically generated high-resolution log hydraulic conductivity reference field.
To capture the geological characterization of the field when generating the reference log hydraulic conductivity field by using the geostatistical method, the proper variogram model needs to be selected and the variogram parameters should be calibrated against the field measurements of the log hydraulic conductivity.The hydraulic conductivity measurements are listed in Table 1.It can be found that 10 hydraulic conductivity measurements are available for the upper layer.For a specific potential variogram model, the adjoint state maximum likelihood cross validation (ASMLCV) method proposed by [31] is employed to estimate the variogram model parameters of the upper layer of the aquifer.The ASMLCV minimizes the cross validation errors between the measured values and the Kriging estimated values, which is expressed in the likelihood function form: where M is the number of the calibration data; ei is the cross validation error, which is the difference between the ith measurement and Kriged value using the rest of the measurements; and σ is the Kriging variance.Three common variogram models are calibrated against the hydraulic conductivity measurements in the upper layer.They are exponential, Gaussian, and spherical variogram models, respectively.The proper model to describe the geostatistical property of the upper layer is selected based on the Kashyap information criterion (KIC) [22]: where ( ) p θ is the prior probability of the parameters θ , in this case, the parameters are the sill and integral scale or range of the variogram models; Nk is the number of the parameters (here, it is 2); and F is the determinant of the normalized Fisher information matrix.The lesser the KIC value is, the higher the model ranks.The estimated parameters, NLL, and KIC values are listed in Table 2.The exponential model is selected due to its lower KIC value.There are only 3 measurements for the bottom layer because the wells were not drilled deep enough; thus, the exponential variogram parameters for the bottom layer are inferred from a similar field in the Ordos Basin around the study area.Here, sill and integral scale of the bottom layer are set to be 0.15 and 20,000 m, respectively.The conditional (on all of the hydraulic conductivity data listed in Table 1) reference log conductivity field is generated using GSLIB [32].To obtain the corresponding groundwater dynamic measurements, 25 synthetically designed head observation locations are used in this case.The generated log hydraulic conductivity field, head observation locations, and initial head distribution for each layer can be found in Figure 5.

Multiple Conceptual Models
To construct multiple conceptual models to simulate this hydrogeological system, two uncertainty sources are considered.Geological model uncertainty due to different interpretations of geological and geophysical data is a well-known model uncertainty, and various zonations can be achieved by varying the number and distribution of hydraulic conductivity zones [33,34].In particular, the hydraulic conductivity field in the two-layer aquifer system is partitioned into 6 and 8 zones as alternative geological models to the reference one with high resolution.The effect evapotranspiration is crucial to groundwater management in semi-arid areas.The relationship between evapotranspiration rate and groundwater level is, however, uncertain.Since MODFLOW-2000, the evapotranspiration segments package (ETS) allows simulation of evapotranspiration with a user-defined relation [35].A two-segment relation model (with the parameter PXDP = 0.4 and PETM = 0.3 at the intermediate segment endpoint) as the reference one is utilized to generate the head observations.As alternative models, a single-segment (linear, no intermediate segment endpoint) relation model and a three-segment model (with the parameters PXDP = 0.2 and PETM = 0.4 at the first intermediate segment endpoint and PXDP = 0.5 and PETM = 0.1 at the second one) are postulated in this case.The relationships used to describe the evapotranspiration rate and hydraulic head are depicted in Figure 6.The combination of two different zonations (with 6 and 8 zones, respectively) and two different segmented evapotranspiration relations (with 1 and 3 segments, respectively) gives 4 alternatively postulated simulation models.Note that the true two-segment evapotranspiration relation model is intentionally excluded from the postulated model set to reflect that the reference model is essentially unknown in practice.The multimodel EnKF analysis procedures are the same as those presented in Section 2. The key parameters to implement the analysis are listed in Table 3.

Results and Discussion
In this study, the assimilation step coincides with the stress period, i.e., the hydraulic head measurements are assimilated to update the log hydraulic conductivity field at the end of each stress period.Figure 7 shows the dynamic change of posterior model weights in each assimilation step.In the first few steps, owing to the fact that less head observation information is available to be assimilated and the postulated geological models are both highly simplified through zonation (relative to the generated reference heterogeneous hydraulic conductivity field), the posterior weights of all the models are nearly uniformly assigned, e.g., each weight is nearly 25%.As more head observations become available, the posterior model weight of model Z8S3 tends to continue increasing and that of model Z6S1 exhibits the opposite tendency.For model Z6S3, the posterior model weight increases in the first 15 steps and starts decreasing after that.The change of the posterior model weight of model Z8S1 is oscillated without a general trend.After about 30 assimilation steps, these weights stabilize and the most sophisticatedly postulated model in this case, Z8S3, is ranked as the best model to describe the reference system, with a posterior model weight of 79.87%.Figures 8a,b and 9a,b depict the multi-model ensemble mean for each layer before and after 36 steps of data assimilation in the studied aquifer system.It can be found that the ensemble mean after 36 steps of data assimilation captures the general pattern of the log hydraulic conductivity distribution by comparing them with the generated reference log hydraulic conductivity field shown in Figure 5.The corresponding multimodel ensemble estimation uncertainties before and after 36 steps of data assimilation for both layers are depicted in Figures 8c,d and 9c,d.As expected, after the data assimilation process, the estimation uncertainty has been dramatically decreased.This quantified multimodel uncertainty is fundamentally more accurate than its counterpart based on the individual model because the conceptual model uncertainty can be taken into consideration, which can be an important uncertainty source.This type of uncertainty can propagate into the prediction of groundwater dynamics through the simulation of the groundwater system, and significantly affect the consequent decision-making or risk assessment based on the prediction.The conceptual model uncertainty has been ignored in most of the previous studies because no suitable algorithm can account for it.Here, with the proposed multimodel EnKF method, the overall uncertainty considering the conceptual model uncertainty can be thoroughly evaluated in each data assimilation step.
where Nj is the number of grid nodes; and j is the node index.It can be observed from Equation ( 16) that the total variance, , .The dynamic change of the spatial averaged variances, including total variance, within-model variance, and between-model variance, are plotted in Figure 10.If a full heterogeneous model is calibrated against the monitored hydraulic head data, it requires estimating 20,840 log hydraulic conductivity values, i.e., log hydraulic conductivity values in all of the active cells, which may bring significant uncertainties associated with the estimation.However, with the zonation patterns in this case, it only requires estimating 6 log hydraulic conductivity values for model Z6S1 and Z6S3, and 8 log hydraulic conductivity values for model Z8S1 and Z8S3.Owing to the fact that the zonation dramatically reduces the number of parameters that need to be estimated, the within-model uncertainty decreases very rapidly with the assimilation of head observation data, and the between-model uncertainty increases in the beginning steps and then decreases when a certain amount of head data has been assimilated.The total uncertainty is dominated by the within-model uncertainty in the beginning and by the between-model uncertainty in the end.This highlights the importance of taking model uncertainty into account; otherwise, the risk induced by the uncertainty will be underestimated.

Conclusions
The proposed multimodel ensemble Kalman filter method is applied to a three-dimensional case study in this work to estimate the hydraulic conductivity field by assimilating dynamic hydraulic head observation data.The main contribution of the proposed multimodel EnKF method is that it can take conceptual model uncertainty into account and also possesses all of the capabilities of the widely used single-model based EnKF method.The accuracy and superiority of this method have been investigated in detail in a previous work published by the same author using a two-dimensional synthetic example.The aim of this paper is to show how to implement this method step-by-step in a more practical case.Despite the fact that part of the dataset is synthetically generated, it does not affect the demonstration of the application of this method in a practical case study.On the contrary, it helps us to assess the performance of this method explicitly while considering all of the collected field information.
Multiple conceptual models are constructed by considering two common but crucial uncertainty sources that influence the groundwater dynamics in a semi-arid study area, such as the Hailiutu River Basin, in this study.One is the zonation pattern of the hydraulic conductivity field, and the other is the relationship between the evapotranspiration and the groundwater level.The posterior model weight can be adjusted dynamically in each assimilation step, according to the mismatch between the predicted values and the actual observations.The mismatch is quantified by the likelihood function, and the posterior model weights are obtained through Bayes' theorem.The obtained multimodel ensemble mean is a weighted average of the ensemble means based on individual models, and the multimodel ensemble uncertainty consists of within-model uncertainty and between-model uncertainty.All of these terms can be quantified in each assimilation step.The results show that the posterior model weights change dynamically and tend to be stabilized when sufficient observations are assimilated.With the increasing of the assimilated information, the within-model uncertainty continues decreasing but the between-model uncertainty tends to increase in the beginning and stabilize in the late stage, and the overall uncertainty tends to decrease in the beginning and stabilize in the late stage.

Figure 1 .
Figure 1.Location of the Hailiutu River Basin.

Figure 2 .
Figure 2. The elevation (in meters) contour of the studied aquifer.

Figure 3 .
Figure 3.The discretized study area: (a) the upper layer; (b) the bottom layer.The yellow line represents the location of the river.

Figure 4 .
Figure 4. Precipitation (measured by rain gauge) and evaporation (measured by evaporation pan with a diameter of 20 cm) rates for 36 stress periods.

Figure 6 .
Figure 6.The reference relationship (black line) and two postulated relationships (red and blue lines) used in the ETS package.

Figure 7 .
Figure 7. Posterior model weights changing with time for individual postulated models.

Figure 8 .
Figure 8.The multi-model ensemble mean and variance of log hydraulic conductivity for each layer before data assimilation.

Figure 9 .
Figure 9.The multi-model ensemble mean and variance of log hydraulic conductivity for each layer after 36 steps of data assimilation.

Figure 10 .
Figure 10.Dynamic changing of spatial averages of total multi-model ensemble variance, within-model variance, and between-model variance with time.
is treated as the log hydraulic conductivity at time t in the forecast step.Together with the hydraulic head at time t and the observed hydraulic head at time t, they form the state vector ,

Table 1 .
The hydraulic conductivity measurements in the study area.

Table 2 .
The parameter, NLL and KIC values, and model rank for the upper layer.

Table 3 .
Summary of the implementation parameters. )