In this study, we discuss the EnKF and its variants in application of vadose zone soil parameters inversion.The inconsisitency of original EnKF and Confirming EnKF in sloving strong nonlinear saturated-unsaturated flow problems are investigated. A modified Restart EnKF method is put forward with lower computional expense and comparable results to Restart EnKF. The sensitivities of different algorithms with respect to factors such as observation error variance, initial guess, ensemble size and observation type are analyzed. Results indicate that the defect of lower accuracy cannot be make up by longer series of observations in fixed points, the priori estimate of first moment of initial realizations is much more important than the priori estimate of second moment. And using less sensitive variables in rapidly changing region (like topsoil) can avoid over correction.
Introduction
Vadose soil parameters are spatially heterogeneous in natural media. The measurements of the soil hydraulic properties based on laboratory experiments are very expensive and have been shown an inability to accurately characterize water flow processes at larger spatial scale, so inverse methods are commonly used to calibrate the spatial distribution of the parameters(Vrugt et al., 2008). Among all the inverse approaches, data assimilation method has become more and more popular due to its flexibility to various problems (Tong et al., 2011). Recently, data assimilation method has been applied to the estimation of vadose soil parameters.(Li and Ren, 2011; Lü et al., 2011; Montzka et al., 2011).
Data assimilation algorithms in recursive Bayesian scheme first emerged with the Kalman filter(Kalman, 1960)，and then，its nonlinear version Extended Kalman Filter (EKF). However, the number of simulations needed at each step in EKF is at the same order as the number of state variables, which will generally be very time-consuming for large systems (Geir et al., 2002). Another shortcoming of EKF is that it may not handle strongly nonlinear problems well because of the uses of a linearized equation for the error covariance evolution(Evensen, 2009).The Ensemble Kalman Filter (EnKF), proposed by (Evensen, 1994) and clarified by (Burgers et al., 1998), is another sequentially based data assimilation method. As a Monte Carlo method, EnKF avoids evolving the covariance matrix of the distribution of the state vector(Johns and Mandel, 2008). EnKF does not require the calculation of the objective function(as in traditional history matching) or the evaluation of the tangent linear operator(as in EKF) or adjoint equations(as in variational data assimilation),which make it easily adopted in different problems. Furthermore, EnKF is inherently compatible with parallel computation technique, since each ensemble member can be independently run on a different processor (Chen and Zhang, 2006).
Different from atmospheric and oceanic data assimilation systems which mainly concern estimation of system state variables，data assimilation in reservior models via EnKF and its variants are not limited to updating system state variables (such as hydraulic head in transient groundwater flow modeling) ,but also allow updating both state variables and model parameters (such as hydraulic conductivity in transient groundwater flow modeling) to yield more accurate model predictions(Tong et al., 2010). In this sense, EnKF and its variants are in effect used as easily performed inverse methods for solving history matching problem (Aanonsen et al., 2009). By assembling model parameters into an augmented state vector, (Geir et al., 2002) first used EnKF to update static parameters in near-well reservoir models. (Chen and Zhang, 2006) have successfully applied data assimilation for transient groundwater flow in geologic formations via EnKF.
The application of EnKF and other data assimilation methods in vadose zone modeling are also receiving more and more attention. Most current studies focuse on the prediction of soil moisure and other variables via assimilating remote sensing data (Reichle and Rolf H., 2008) like brightness temperature observations(Crow and Wood, 2003; Dunne and Entekhabi, 2006; Loew, 2008) into kinds of catchment land surface model(Reichle et al., 2007; Xie and Zhang, 2010; Chen et al., 2011) ,and studies exploring how to update soil hydraulic parameters are few. Up to now. (Montzka et al., 2011) assessed the performance of the Sampling-Importance-Resampling Particle Filter (SIR-PF) for a 1D mechanistic soil water model (hydrus 1D)，their results indicate that updating of both hydraulic parameters and state variables allowed better predictions of top soil moisture contents as compared with updating of states only. (Lü et al., 2011) used the OP-EKF(a method coupling optimal parameters and extended Kalman filter) to calibrate the soil hydraulic parameters and the root zone soil profile. (Li and Ren, 2011) extended the use of EnKF to parameter estimation in vadose zone hydrology. However ,all the studies above focused on the homogeneity soil and based on one-dimensional (1-D) Richards' equation，the performance of data assimilation in heterogeneous soil and in higher-dimensional situation (which means much more complicate flow pattern and a great increase of state vector space) needs to be further investigated.
Driven by the fierce change of top surface moisture, and more importantly, the nonlinear relationship between soil moisture and pressure head, the unsaturated water flow has a strong physical nonlinearity, which may limit the applicable range of original EnKF. The nonlinearity may lead to two major difficulties, one is the unphysical updates of the state variable (eg. water saturations greater than 1.0); another is the possible inconsistency between updated parameter(eg. hydraulic conductivity) and the updated state (eg. pressure head). To be specific, the consistency discussed here refers to that the updated state vectors at any time and the state vectors at the same time obtained by rerunning the simulator using final updated model parameters from time zero should obey a same probability distribution(Wang et al., 2010). To alleviate the negative effect from nonlinearity, some iterative EnKF(IEnKF) methods have been proposed. According to the starting time of each iteration, most IEnkF algorithms can be classified into two basic forms: 1) Rerun from previous assimilating step(Moradkhani et al., 2005; Wen and Chen, 2006; Li and Reynolds, 2007; Lorentzen and Nævdal, 2011); and 2) Rerun from time zero (Reynolds et al., 2006; Gu and Oliver, 2007; Hendricks Franssen and Kinzelbach, 2008; Zhao et al., 2008; Wang et al., 2010; Chen, 2012).The IEnKF algorithm reruning from time zero requires extremely high computational expense. And IEnKF algorithm reruning from previous assimilating step may still produce considerable inconsistency. As proved by (Zafari and Reynolds, 2005) theoretically, the popular IEnKF proposed by (Wen and Chen, 2006) actually cannot retain consistency even for simple cases under linear and Gaussian assumptions.
Although data assimilation of unsaturated flow has gained attention, to the best of our knowledge, the performance of iterative EnKF for strongly nonlinear unsaturated flow in heterogeneous soil has not been carefully investigated. Furthermore, although the so-called Confirming EnKF (Wen and Chen, 2006) has almost been accepted as a common practice in many studies and achieved some sound results (Gu and Oliver, 2006; Krymskaya et al., 2008; Li et al., 2010; Zagayevskiy et al., 2012; Zhang et al., 2012) , as we will show that it may produce considerable inconsistency under some situations. It is also not clear how the observation error variance, initial guess and ensemble size affect the performance of iterative EnKF in non-linear variably unsaturated flow. Moreover, as revealed by many studies (Liu et al., 2008; Fu and Jaime Gómez-Hernández, 2009; Bailey and Baù, 2010), different types of observation may provide unique information concerning the parameter.For example,both hydraulic head and solute concentration data have been used to estimate spatially-variable hydraulic conductivity(Liu et al., 2008). More recently, (Bailey et al., 2012) used nitrate concentration and mass of nitrate entering the river to estimate denitrification rate constant. The data value of different measurements (eg.pressure head and water content) has not been explored in variably saturated flow.
In this paper, we investigate the ability of EnKF and IEnKF to characterize two-dimensional unsaturated flow in heterogeneous soil through incorporating transient hydraulic head and water content measurements. Firstly, traditional Confirming EnKF and Restart EnKF are re-examined in the background of strong nonlinear unsaturated flow. Contrary to Hendricks Franssen and Kinzelbach's finding (2008) in saturated flow, inconsistency produced by EnKF and Confirming EnKF make great damage in unsaturated flow. Then, a modified Restart EnKF is developed to alleviate the computational effort of original Restart EnKF and meanwhile to retain consistency. Thirdly, the influences of the observation error variance, initial guess, and ensemble size are investigated. Lastly, the effect of observation type (pressure head and water content) on IEnKF stability and efficiency is investigated to provide reference for measurement network design.
Methodology and Formulation
Governing equation
We consider a two-dimensional transient unsaturated flow system in heterogeneous soil. The flow equation is given by
\* MERGEFORMAT ()
Where is the soil volumetric water content[L3/L3], is time[T], is the soil water pressure head[L], is the unsaturated hydraulic conductivity[L/T], x is the horizontal coordinate[L] ,z is the vertical coordinate[L], and is the root water uptake or/and other source/sink term[T-1]. The soil hydraulic properties were described by van Genuchten-Mualem constitutive relationships(Van Genuchten, 1980):
\* MERGEFORMAT ()
\* MERGEFORMAT ()
Where
\* MERGEFORMAT ()
is the effective water saturation[-];and are the residual and saturated volumetric water contents[L3/L3], respectively;is the saturated hydraulic conductivity[L/T]; and are parameters related to the soil pore-sized distribution[-],and [-].No hysteresis is considered in this study. Equation is solved using standard Galerkin finite element method (Å imunek and Vogel, 1992).
Ensemble Kalman filter with augmented state vector
We defined an augmented state vector at
\* MERGEFORMAT ()
where is the joint-model-state-observation vector consisting of model parameter vector, ( eg. hydraulic conductivity), state variable vector (eg. water pressure), and theoretical observation vector (corresponding to observations); is the dimension of the augmented state vector ; refers to that the current state is (often nonlinear) function of previous model and state; refers to that the theoretical data is (often nonlinear) function of model and state.
Assume a set of observations (eg. water pressure) are available at time , then the relation between the observations and true augmented state vector can be expressed as:
\* MERGEFORMAT ()
where is the unknown true value of observations;is the observation error which is assumed to be Gaussian with mean and covariance ,and the observation error is assumed to be uncorrelated at different observation times. is the observation operator which represents the relationship between the augmented state vector and the observation vector. can be written as
\* MERGEFORMAT ()
where is a matrix with all 0s as entries,is a identity matrix.
According to Bayes theorem, we can write the posterior probability based on all the observations up to current time as
\* MERGEFORMAT ()
Where is the collection of all observations from time 0 to time.
With the Gaussian assumption of prior probability and the likelihood probability, Equation is given by
\* MERGEFORMAT ()
Then a sample of the posterior density function should minimize the following objective function:
\* MERGEFORMAT ()
where the subscript is the sample index; is the prior value of ,usually obtained by running the simulator from time to ;is the random perturbed observations with mean equal to the observed and covariance equal to ,the perturbation is essential to prevent updated state vector with too low variances(Burgers et al., 1998);is the covariance of observations, and it is a diagonal matrix if the measurement errors are uncorrelated with each other;is the covariance of augmented state vector and can be calculated from an ensemble of realizations in EnKF as:
\* MERGEFORMAT ()
By setting the gradient of equal to zero, the updated augmented vector is given by following equation:
\* MERGEFORMAT ()
Where is the ensemble size;denotes the assimilated value of ;is known as the Kalman gain, which is defined as
\* MERGEFORMAT ()
One important advantage of EnKF is we in practice never need to calculate the whole of dimension in, instead we just calculate part of it such as and in by , the dimension of resulting matrices is only and , respectively.
\* MERGEFORMAT ()
Iterative Procedure
1.3.1 Confirming EnKF
The basic idea of Confirming EnKF(Wen and Chen, 2006) is to implement an additional procedure called "Confirming" to ensure the updated state variables and model parameters are plausible consistent(Aanonsen et al., 2009). At each assimilation step, the Confirming EnKF still forecasts from previous time and computes the kalman gain as in ,but only update the model parameter vector in ,then use the updated parameters to compute the state variable vector at time tn by running the simulator again from previous time tn-1 as
\* MERGEFORMAT ()
The updated model parameters and rerun state variables with theory observations compose the initial state vector of next assimilation step. The flowcharts of Confirming EnKF and original EnKF are presented in Fig.1. Confirming EnKF generates plausible consistency(Aanonsen et al., 2009) and has been successfully used in reservoir research(Gu and Oliver, 2006; Krymskaya et al., 2008; Li et al., 2010; Zagayevskiy et al., 2012; Zhang et al., 2012). (Hendricks Franssen and Kinzelbach, 2008)showed that Confirming EnKF and original EnKF have similar performance in saturated groundwater flow. However,(Zafari and Reynolds, 2005) proved that Confirming EnKF is theoretically inconsistent.
1.3.2 Restart EnKF
Considering the fact that the state variable can be computed with the updated model parameter and initial condition from time zero:
\* MERGEFORMAT ()
then the state variables part of can be taken from augmented state vector :
\* MERGEFORMAT ()
with the corresponding Bayes formula:
\* MERGEFORMAT ()
and an update equation
\* MERGEFORMAT ()
Although the update equation of Restart EnKF is similar to of original EnKF, it is emphasized that augmented state vector in these two schemes contains different components. Restart EnKF requires running the simulator not from previous time but from time zero at each assimilation step. The flowcharts of Restart EnKF are also presented in Fig.1. Comparing to Confirming EnKF, Restart EnKF is intrinsically consistent. It has been used in saturated groundwater flow (Hendricks Franssen and Kinzelbach, 2008) and reservoir research(Thulin et al., 2007; Wang et al., 2010) .The basic idea of using an augmented state vector without state variables and rerunning simulator from time zero is also shared by EnRML(Ensemble Randomized Maximum likehood Filter)(Gu and Oliver, 2007) and its variants(Wang et al., 2010; Chen, 2012; Sakov et al., 2012).However, the computational expense of Restart EnKF is unbearable for long-term simulation since it has to rerun all the ensemble members from time zero at every assimilation time.
1.3.3 Modified Restart EnKF
Restart EnKF is very time-consuming since it has to rerun all the ensembles from the time zero to the current time. As an alternative, the modified Restart EnKF is still under the theoretical framework of Restart EnKF but only need to rerun the mean of the ensemble. The modified Restart EnKF is based the observation that the ensemble spread from original EnKFã€Confirming EnKF and Restart EnKF are all very similar at any assimilation step even the prediction of original EnKF and confirming EnKF may be not correct.
In each assimilation step of modified Restart EnKF, only the model parameters are updated using kalman gain and the mean of updated model parameters is approximated by
\* MERGEFORMAT ()
then we use the parameters ensemble mean to rerun the simulator from time zero to current time to get a new "mean" of current state variables:
\* MERGEFORMAT ()
and then we can construct a new current state variables ensemble by replace the mean part of each realization with the new one:
\* MERGEFORMAT ()
\* MERGEFORMAT ()
The flowcharts of modified Restart EnKF are also presented in Fig.1. It should be noted that modified EnKF is based on two assumptions: one is the mean state computed by parameter ensemble mean can be regarded as approximation of state variable ensemble mean, which generally holds if the variance of parameter field is not very strong; the other is the variances of state variables change little in each assimilation step.
Numerical Experiments
2.1 General Setup
In this section,a synthetic experiment is constructed to test different algorithms and explore the influnence of various factor, incluing the observation error variance, initial mean and variance，ensemble size and observation type. Fig.2 shows the sketch of a general test-pit experiment, the 200Ã-200[L2] profile is discretized into uniform 50Ã-50 gird cells with the size of 4Ã-4[L2]. The bottom and two lateral sides are impervious boundaries, the top side is imposed on steady rainfall of 0.2[L/T]. The initial pressure head is -200[L] at the top,and evenly changes to 0 at the bottom.The pressure head h and/or moisture content θ observations at 20 locations are drawn from the simulation with the reference field. The total simualtion time is 200 [T],the observations are assimilated every 1[T]. A reference log hydraulic conductivity field is treated as Gaussian random field with mean <Y>=<lnKS>= 2 and variance σY2=0.7. The horizontal and vertical correlation length of reference Y field and initial Y realizations are 50[L] and 20[L],respectively.The remaining parameters in and are assumed to be deterministic constants,with θr=0.0001,θs=0 .399,α=0. 0174,n=1.3757. The detailed specifications of cases tested in the following section are shown in Table 1.
Case
Observation error
variance
Mean of initial realizations
Variance of initial realizations
Observation
type
No. of realization
Reference
-
2.0
0.7
-
-
Case1
0.64
2.0
0.7
20h
1000
Case2
64
2.0
0.7
20h
1000
Case3
0.64
2.0
0.4
20h
1000
Case4
0.64
2.0
1
20h
1000
Case5
0.64
3.0
0.4
20h
1000
Case6
0.64
3.0
1
20h
1000
Case7
0.64
2.0
0.7
20h
100
Case8
0.0025
2.0
0.7
20θ
1000
Case9
0.64/0.0025
2.0
0.7
10h/10θ
1000
Table 1 Parameter setup for all cases
2.2 Comparison of different algorithms
In this article，the criterions RMSE(root mean square error) and Ensemble Spread are used to evaluate the accuracy of different algorithms:
\* MERGEFORMAT ()
\* MERGEFORMAT ()
where is the estimated mean, are the reference values, denotes the ensemble variance at each grid, and is the number of grids. The RMSE measures the deviation between the ensemble mean and reference field. And the Ensemble Spread represents the estimated uncertainty based on the ensemble.
The inconsistency of H is measured by comparing the estimated mean with the mean after rerunning the simulator:
\* MERGEFORMAT ()
where is the mean computed by rerunning the simulator from time zero.
In case 1, the mean and variance of initial realizations is set to be 2.0 and 0.7, same as the reference Y. A relatively large ensemble size of 1000 is used to restrain the spurious correlation and resulting filter divergence which may happen under a small ensemble size. Fig.3 compares the RMSE and Ensemble Spread of log saturated hydraulic conductivity Y and pressure head H for original EnKF, Confirming EnKF, Restart EnKF, modified Restart EnKF and unconditional run(without data assimilation).It is shown that the Restart EnKF and modified Restart EnKF produce similar and much better results than other methods. Just at t= 25[T]，the Y RMSE of original EnKF and Confirming EnKF began to increase dramaticly.Because a very large ensemble size is used here, it can not be simply explained to be the problem of filter divergence, which caused by the underestimation of the ensemble covariance under a small ensemble size(Evensen, 2009). One possible reason is over correction caused by the EnKF's linear update in nonlinear problems(Gao et al., 2005), which can be resloved by a additional damping factor(Hendricks Franssen et al., 2011; Wu and Margulis, 2011), or equivalent to an enlarged observation error variance(Rommelse, 2008). The effect of diferent observation error variance would be discussed in the following sections and it turnd out not to be the key factor of the crash of original EnKF and Confirming EnKF. In another perspective, if the overshooting of update process is indeed serious, the Restart EnKF should also cash because all these methods share the same form update formula. The most possible reason is the implicit inconsistency beween updated model parameters and updated state variables of original EnKF and Confirming EnKF gradually accumulates with the simulation. High inconsistency will exhibit a pseudo correlationship between theoretical data and model parameters(or state variables) in next assimilation step and bring deviation in update of model parameters(or state variables).The deviation gradually accumulates in each assimilation step and will finally destroy the covariance structure of state vector. We show the obtained IC from different algorithms in Fig.4. The inconsistency of pressure head from modified Restart EnKF is at a low level during the whole data assimilation process, while the inconsistency from original EnKF and Confirming EnKF grow rapidly. The spatial distribution of IC from different methods is plotted in Fig.5. It is seen that the inconsistency gradually spreads from the top to the bottom with the wet front moving downward. Fig.6 presents the reference Y field and estimated <Y> from various methods at t=200[T]. It is shown modified Restart EnKF has the similar performance as Restart EnKF and both two methods capture the main patterns of the reference field, while Original EnKF and Confirming EnKF give unacceptable results. Fig.7 presents the contours of σY from different method at t=200[T]. The variances in the topsoil and at the observation points is quite low, which indicate the transient head at the measurement point provides positive information to calibrate the hydrulic conductivity around it. One interesting phenomena is that all methods have very similar variance distribution despite of the dissimilar <Y> distribution.
2.3 Sensitivity Analysis
This part will focus on the sensitivity of original EnKF, Confirming EnKF, Restart EnKF and modified Restart EnKF with respect to several different factors including the observation error variance, initial mean and variance, ensemble size and observation type. Based on Case 1,the observation error variance increases from 0.64 to 64 (Case 2) to examine the performance of the four methods in presence of different observation error variances. Also on the basis of Case 1,the variance of initial realizations is changed from 0.7 to 0.4(Case 3) and 1.0 (Case 4) to test the effect of initial variance. Then on the basis of Case 3 and Case 4,the mean of initial realizations is changed from 2.0 to 3.0(Case 5 and Case 6), repectively,to test the effect of initial mean.To discuss the impact of ensemble size, we recompute case 1 using an ensemble size of 100(Case 7), and 5 different ensembles are recomputed to avoid errors caused by the sample selection. Finally,in case 8 and case 9 we change the observation type of case 6 from 20 pressure heads to 20 water contents and 10 upper water contents plus 10 lower pressure heads. The detailed input parameters for these cases are all listed in Table 1.
2.3.1 Observation error variance
From the derivation of kalman gain, it is seen that a larger prior error variance results in more weight given to the observation and larger correction to the state vector. In verse, a larger observation error variance can relieve the correction to state vector prior and may be beneficial in sloving the divengence problems of original EnKF and Confirming EnKF.The comparision of different observation error variances(Case 1 vs Case 2) is ploted in Fig.8. It is seen that a larger observation error variance(=64) indeed slows down the rapid increase of RMSE(Y) from original EnKF and Confirming EnKF,however,the two methods still diverge in later time.This result indicates the divergence can be accelerated/decelerated by larger/smaller corrections but the occurrence of divergence is not detemined by the amplitude of correction.As mentioned before, the deteriorations of original EnKF and Confirming EnKF are essentially caused by inconsisitency.
Comparing the perforance of Restart EnKF and modified Restart EnKF under different observation error variance, we can see that stable RMSE(Y) from high observation variance is higher than that from low variance. This finding may suggest that low-precision measurements in fixed points are not able to inverse the Y field accurately even by assimilating a very long-term obersations.For Case 2, the contours of the reference Y and the estimated <Y> at t=200[T] of the four algorithms are ploted in Fig.9.It is seen Restart EnKF and modified Restart EnKF produces less effective prediction than Case 1(Fig.6),while the results of orginal EnKF and Confirming EnKF are improved partly.
2.3.2 Initial guess
In practice, the initial relizations would be generated based on prior knowledge of study area , and initial guess may deviate from the reality. In this subsection ,we investigate performance of the four algorithms on impact of the poor prior knowledge.The effects of the first and second moments of prior are considered by shifting the mean and variance of Y initial realizations.
The RMSE comparisons of different initial variances are displayed in Fig.10( case 1 vs case 3 and case 4). It can be seen that a overestimated initial variance may relieve the correction of original EnKF and Confirming EnKF, and then produce better results than underestimated ones in early time .The results of Restart EnKF and modified Restart EnKF are also slightly improved by a larger initial variance, even it's overestimated. However, the guess of initial variance seems not to have significant impact on the finally assimilated results.
Effect of wrong initial mean(<Y>=3.0) are shown in Fig.11 with two kinds of intial variance(σY2=0.4 and σY2=1.0). It can be seen that all RMSE(Y) curves of four algorithms rise sharply in first ten assimilation steps under the impact of wrong intial mean,then the lines begin to drop.Restart EnKF and modified Restart EnKF still perform much better than original and Confirming EnKF ,although the gaps have been narrowed due to the wrong intial mean and the final results are all unsatisfactory.It also can be observed that the performance of modified Restart EnKF is not so close to Restart EnKF (Fig.11(c) vs Fig.11(d) ) as in right initial mean,which means that using modifed Restart EnKF as an alternative of Restart EnKF when lack of prior knowledge still need careful consideration. After the comparsion we also can draw the conclusion that the initial mean initial variance plays an important role in vadosen zone data assimilation.
2.3.3 Ensemble Size
As Monte Carlo methods, EnKF and its variants require large number of realizations to approach statistical convergence. However, the heavy computational expense of Monte Carlo methods limits the ensemble size in real applications. RMSE comparisons under different ensemble size (1000 vs 100) are shown in Fig. 12(Y) and Fig.13(H).All methods generally generate larger RMSE with an ensemble size of 100,especially for RMSE(Y). It is seen that small ensemble size also leads to considerable estimation uncertainty, i.e. different series of 100 realizations may lead to highly different results. Since restart EnKF and modified restart EnKF have the best performance, a particular inspection is given to Fig. 12(c,d) and Fig. 13(c,d). It is observed that the assimilated <Y> from 100 realizations has a much larger deviation to reference field than that from 1000 realizations. However, if one only focuses on the prediction of pressure head, it seems the ensemble size does not have great effect on the accuracy.
2.3.4 Observation type
In a regular test-pit experiment, two types of contentious observation information, i.e pressure head and water content are usually available. The pressure head information is very sensitive to soil water flow but usually with a lower precision, while the measurement technique for water content is more sophisticated and accurate but may not sensitive to the subtle movement of water flow. In Cases 1-7, the observations are 20 pressure heads in two vertical lines. In this subsection, we explore the effect of different observation types.
Based on case 7, Case 8 is designed to have 20 water content measurements instead of pressure head. The observation error variance is set as σθ2=0.0025. The comparison of case 7 and case 8 are given in Fig.14 .It is shown that with water content as observations the RMSE(Y) curves overlap for all the methods, while at the same time very limited information is extracted in the inverse modeling of Y. This is because water content have a certain range and therefore the measurements are far more insensitivity than pressure head to water flow movement. A slight change of water content usually corresponds to a dramatic change of pressure head, especially in dry soils. And the statistic-based update of kalman gain cannot distinguish the characteristics of different observation types correctly, so the assimilation of water contents is too weak to show significant correction.
In Case 9, a new observation strategy is employed. The 10 points at the upper pits are equipped with water content sensor, and the bottom 10 points remain observing pressure head. A surprisingly improved prediction of Y is obtained (Fig.15) during the first 100 [T] even with original EnKF and Confirming EnKF. The RMSE(Y) from original EnKF and Confirming EnKF increases at later time, while RMSE(Y) from Restart EnKF and modified Restart EnKF goes stable at a much lower value than in Case 7.The improvement is due to two aspects: using water content measures instead of pressure head measures in upper dry area can alleviates the impact of wrong initial guess in early time, at the mean time assimilating pressure heads in bottom area can capture more details of water movement. The contours of the reference Y and the estimated <Y> at t=200[T] of Case 8 are ploted in Fig.16. It is seen the assimilated Y of Restart EnKF and modified Restart EnKF( Fig.16(e,f)) give decent but pooper results compared to case 1(with correct initial guess)in Fig.6(e,f),and the contours of orginal EnKF and Confirming contains many abnormally high/low values.The wrong initial guess also brought especially high deviation in surface layer ( Fig.16(e,f)).Fig.17 plots the mean and variance of <Y> in whole simulated domain.In which we can clearly see how the Restart EnKF and modified Restart EnKF gradually changing totally wrong Y fields toward to the reality, superior to which of original EnKF and Confirming EnKF.
These results suggest us that in practical applications we may use insensitive measures in rapidly changing region (eg. topsoil) to avoid over correction, and use sensitive measures in stability region (eg. subsoil) to achieve more valuable information.
3 Conclusions and discussions
The Ensemble Kalman Filter (EnKF) provides a flexible tool for combining the information from prediction of model and observation to obtain better estimation of system state and parameters. However, the application of EnKF was limited by its linear update scheme. Various iterative variants of EnKF have been developed to overcome the inconsistency in nonlinear problems. In this study, we extended the study of iterative EnKF to unsaturated flow in vadose zone.We test the application of original EnKF and so-called Confirming EnKF and Restart EnKF,which represent two basic forms of complicated iterative EnKF methods, and both can continuously update the unsaturated flow model though assimilating dynamic observation of presure head.A modified Restart EnKF with lower computational expense and acceptable result is put forward.A synthetic 2-D steady rainfall example is designed to test the original EnKF, Confirming EnKF, Restart EnKF and modified Restart EnKF. The performance has been evaluated for different observation error variance, intial guess, ensemble size and observation type. This study has the following major conclusions:
In all cases, the Restart EnKF and modified Restart EnKF perform much better than other two methods.The implicit inconsistency exists in original EnKF and Confirming EnKF make them soon collapse even in this simple 2-D steady rainfall example. Thus the widely used confirming procedure in reservoir research may be no longer applicable in strong nonlinear problems in vadose zone. And the low computational expense and comparable results of modified Restart EnKF make it a of a preferable alternative Restart EnKF.
A larger observation error variance can slower the correction and convergence of ensembles.The collapse of original EnKF and Confirming EnKF would be relieved but can't be avoided.And stable results from high observation variance is obviously higher than that from low variance,which suggest that low-precision measurements in fixed points may not be able to inverse the Y field accurately even by assimilating a very long-term obersations.
The choice of initial realizations has a profound influence on data assimilation in vadose zone, poorer initial guess in generally leads to worse performance. It's worth noting that the first order moment shows up much more influence than second order moment.
Besides the superior performance in relative large ensembles(1000 in this study), Restart EnKF and modified Restart EnKF remain their advantages in small ensembles(100 in this study).Although small ensemblese get worse estimation parameters with considerable estimation uncertainty, the ensemble size seems not have great effect on the accuracy pressure head prediction.
Due to the rapid soil water movement driven by atmospheric boundary,using moisture content instead of pressure head observations in topsoil can contribute to the stability of the assimilation process. On the contrary,using more sensitive variable like pressure head in subsoil can better capture the details of water movement to achieve better assimilation effect. According to the characteristics of observations and combining different types of observation properly can achieve much better assimilation results.Although the analyses and discussions are preliminary in this article,it may be helpful in observation system design.
