Journal of Hydrology 589 (2020) 125210

Contents lists available at ScienceDirect

Journal of Hydrology
journal homepage: www.elsevier.com/locate/jhydrol

Research papers

A gaussian process-based iterative Ensemble Kalman Filter for parameter
estimation of unsaturated ﬂow

T⁎

Kun Liua,b, Guanhua Huanga,b, , Zheng Jianga,b, Xu Xua,b, Yunwu Xionga,b, Quanzhong Huanga,b,Jiří Šimůnekc
aCenter for Agricultural Water Research, China Agricultural University, Beijing 100083, China
Chinese-Israeli International Center for Research and Training in Agriculture, China Agricultural University, Beijing 100083, Chinac
Department of Environmental Sciences, University of California Riverside, Riverside, CA 92 521, USAb

A R T I C LE I N FO

A B S T R A C T

Keywords:
Data assimilation
Gaussian process
Ensemble Kalman Filter
Surrogate model

Water ﬂow in the unsaturated zone is an important component of the water cycle. Accurate estimation of soil
hydraulic parameters ensures precise simulations of water ﬂow in the unsaturated zone. In this study, a GaussianProcess-based Iterative Ensemble Kalman Filter method (GPIEnKF) is proposed and applied for estimating soil
hydraulic parameters for two-dimensional soil water ﬂow. The method involves a surrogate model for twodimensional soil water ﬂow that is based on the Gaussian Process method. The accuracy and eﬃciency of the
GPIEnKF method are validated using two synthetic cases and a real dataset involving a ﬁeld experiment withdrip irrigation. The impact of the layout of observation points as well as the number of training base points andobservation points on the estimation of parameters are also analyzed. The results show that the surrogate model
can accurately predict pressure heads at observation points. The layout of observation points that precisely
describes inﬁltration water movement allows the surrogate model to better emulate the original model, thereby
improving the accuracy of parameter estimation. The number of training base points is found to have only asmall inﬂuence on parameter estimation. The accuracy of parameter estimation and pressure head predictions
can be further improved by increasing the number of observation points. However, the accuracy of predictions isaﬀected by the uncertainty of boundary conditions and the soil spatial heterogeneity. The GPIEnKF method caneﬀectively estimate multiple parameters characterizing water ﬂow in layered soils. As compared with thestandard Iterative Ensemble Kalman Filter, the GPIEnKF method can greatly improve computational eﬃciency
while obtaining comparable results. The GPIEnKF method is an eﬃcient tool for parameter estimation of multidimensional soil water ﬂow.

1. Introduction
Water ﬂow in the unsaturated zone is an important component of
the water cycle. From a hydrological point of view, soil moisture determines the separation of precipitation into runoﬀ, evaporation, andinﬁltration, and the division of net radiation into latent, sensible, andgeothermal heat ﬂuxes (Pauwels et al., 2007). As concern regarding theprotection of the environment and the quality of groundwater rises, theimportance of an accurate soil-physical description of porous systems
within the unsaturated zone is increasingly recognized in environmental engineering and groundwater hydrology (Hopmans andŠimůnek, 1999).
The quantiﬁcation of water ﬂux exchanges between the land surface, the unsaturated zone, and groundwater requires an accurate⁎

parameterization of soil hydraulic properties (Zhang et al., 2018). Many
inverse methods, such as the Levenberg-Marquardt algorithm (Koolet al., 1985), the shuﬄed complex evolution algorithm (Duan et al.,1992), and the AMALGAM algorithm (Vrugt and Robinson, 2007), havebeen used for estimating soil hydraulic parameters (Šimůnek et al.,
1998; Abbasi et al., 2003; Vrugt et al., 2003; Durner et al., 2008;
Wöhling et al., 2008). The main problem with respect to these methods
is that the uncertainty is attributed to model parameters, ignoring theadditional uncertainty of inputs, outputs, and model structures (Vrugt
et al., 2008). In the process of model construction, there are many errors, both in the model structure, due to the simpliﬁcation of physicalprocesses, and in observation data, due to the limitation of observation
methods. It is incorrect to attribute the uncertainty only to model
parameters.

Corresponding author at: Center for Agricultural Water Research, China Agricultural University, Beijing 100083, China.
E-mail address: ghuang@cau.edu.cn (G. Huang).

https://doi.org/10.1016/j.jhydrol.2020.125210
Received 6 October 2019; Received in revised form 16 June 2020; Accepted 17 June 2020
Available online 22 June 20200022-1694/ © 2020 Elsevier B.V. All rights reserved.

Journal of Hydrology 589 (2020) 125210

K. Liu, et al.

and apply it for the estimation of hydraulic parameters for two-dimensional unsaturated ﬂow problems. The paper is organized as follows: the methodology is introduced in Section 2, the accuracy andeﬃciency of the proposed method are presented using two synthetic
cases and a real dataset involving a ﬁeld experiment with drip irrigationin Section 3, and ﬁnally, conclusions are summarized in Section 4.

Data assimilation methods provide a strong framework for dealing
with the uncertainty problem. Due to its superiority in solving highdimensional and nonlinear problems, the Ensemble Kalman Filter
method (EnKF) (Evensen, 1994) has been a preferred method for assimilating state variables and for estimating parameters of unsaturated
zone (De Lannoy et al., 2007; Li and Ren, 2011; Pasetto et al., 2012;
Medina et al, 2014; Song et al., 2014; Shi et al., 2015; Brandhorst et al.,2017). However, applying the EnKF method to unsaturated zone problems may encounter an inconsistency problem: after being updated
simultaneously using a linear update function, the state variables andparameters may not follow the nonlinear Richards equation (Song et al.,2014). Since variance-minimizing analyses, rather than physical principles, are used to update state variables and parameters, these may
violate certain constraints, such as the water balance or nonnegativityof parameters and states (Pan and Wood, 2006). To resolve the inconsistency problem, several improved EnKF algorithms have been proposed in the literature. For example, Wang et al. (2009) proposed thenaive method, the projection method, and the accept/reject method toincorporate the constraints into the EnKF method. Moradkhani et al.
(2005) proposed the Dual-EnKF method, in which two parallel EnKF
structures are used to update state variables and parameters. Song et al.
(2014) proposed a modiﬁed Restart EnKF method and evaluated three
iterative EnKF algorithms. The results show that the Restart-EnKF and
its modiﬁcation are better than the original EnKF method. Whileiterative methods are an eﬀective way to resolve the inconsistency
problem, they are computationally expensive, which limits their applications in complicated cases. So far, most applications of the EnKF in
the unsaturated zone have considered only one-dimensional water ﬂow
problems. When water ﬂow is two- or three-dimensional, the computational cost increases dramatically. Therefore, it is necessary to propose a more eﬃcient data assimilation method for high-dimensional
water ﬂow in the unsaturated zone.
Surrogate models have been extensively used in hydrologicalmodeling to improve computational eﬃciency. An excellent review of
the use of surrogate modeling in water resources research was presented by Razavi et al. (2012). There are two types of surrogate modeling techniques, namely response-surface surrogates and lower-ﬁdelity
surrogates. The response surface is the relationship between the model
inputs and the model outputs. The response-surface surrogate model isa kind of data-driven technique for approximating the response surfaceof the original model based on a limited number of original model
evaluations. A lower-ﬁdelity surrogate is an alternative model of theoriginal model with a lower level of accuracy (Razavi et al., 2012).
When a problem is low-dimensional so that the response surface can becharacterized using few parameters, the response-surface surrogate ispreferred since a limited number of original model evaluations is required (Brunetti et al. 2017). In recent years, many scholars havecombined response-surface surrogate models with data assimilation
methods to improve assimilation eﬃciency. For example, Chang et al.
(2017) developed a surrogate model, based on an iterative ensemble
smoother algorithm, in which the polynomial chaos expansion method
and the interpolation method were used. Li and Xiu (2009) proposed a
generalized polynomial chaos method based on the EnKF algorithm.
Zeng et al. (2011) developed a Probabilistic Collocation-based Kalman
Filter to achieve an eﬃcient quantiﬁcation of uncertainty. Ju et al.
(2018) proposed an adaptive Gaussian Process-based Iterative Ensemble Kalman smoother. Compared with other surrogate methods, theGaussian Process method (GP), which can provide output uncertainty ofa surrogate model, has attracted much attention in hydrological modeling. Another advantage of the GP surrogate is that it is an exact
emulator, i.e., the GP surrogate can precisely predict values at all design
sites used in training of the surrogate, while an inexact emulator may
introduce bias at design sites (Brunetti et al., 2017).
However, to our knowledge, there are no studies that combine theGP surrogate with the EnKF method for parameter estimation.
Therefore, in this study, we propose a GP-based iterative EnKF method

2. Methodology
2.1. Governing equation
Two-dimensional water ﬂow in the unsaturated zone can be described using the Richards equation:∂θ∂ ⎡ ⎛ A ∂h⎤=
K ⎜K ij+ K izA⎞⎟ ⎥ − S∂t∂x i ⎢ ⎝
∂x j
⎠⎦
⎣

(1)

where θ is the volumetric moisture content [-], h is the pressure head[L], K is the unsaturated hydraulic conductivity [LT−1], K ijA are components of a dimensionless anisotropy tensor KA, xi and xj are thespatial coordinates [L], and S is the sink term [T−1] accounting for root
water uptake.
The highly nonlinear relationships between the unsaturated soil
hydraulic properties and the pressure head can be described by theusing the van Genuchten-Mualem model (Mualem, 1976; van
Genuchten, 1980):(

1

K (h) = Ks Se0.5 ⎡1 − 1 − Sem
⎣
Se =

m 2

) ⎤⎦

(2)

θ − θr= [1 + (αh)n]−m
θs − θr

(3)

1
n

m=1−

(4)
−1

where Ks is the saturated hydraulic conductivity [LT ], Se is the effective water saturation [-], θr and θs are the residual and saturated
water contents, respectively [-], α is the inverse of the air-entry value[L−1], and n is the pore-size distribution index [-]. In this study, unsaturated water ﬂow was simulated using the CHAIN_2D code (Šimůnek
et al., 2008), i.e., an open source version of HYDRUS-2D.
2.2. Parameter estimation based on the iterative EnKF framework
The iterative EnKF framework for parameter estimation consists of amodel operator, an observation operator, and an iterative ensemble
Kalman ﬁlter algorithm. The model operator is used to pass the statevector to the next time level, which can be expressed as follows:

Xi + 1 = MXi + η
i

(5)
i+1

are the state vectors at the current and next time
where X and X
levels, respectively, i and i + 1 denote time indicators, M is the model
operator, and η is the model error vector that is independent white
noise for the model operator, which is drawn from a multinormal distribution with a zero mean and a speciﬁed covariance Q. In this study,
the soil hydraulic parameters to be estimated constitute the state vector.
It is assumed that the true values of these parameters will not change
during the simulation period and their estimates will eventually converge to the true values during the assimilation process. Therefore, theM operator is the unit matrix I.
The observation operator constructs the mapping between the statevector and the observation vector, which can be written as:

Y i + 1 = H (Xi + 1) + ε
i+1

(6)

is the observation vector, H is the observation operator,
in which Y
and ε is the observation error vector, which is also assumed to be independent white noise drawn from a multinormal distribution with a2

Journal of Hydrology 589 (2020) 125210

K. Liu, et al.

zero mean and speciﬁed covariance O. The two-dimensional soil waterﬂow model establishes the relationship between the soil hydraulicparameters and the observation vector (e.g., pressure heads) and thus
the two-dimensional Richards equation is the observation operator inthis study. The Richards equation errors were included in the observation errors, which are assumed to be independent white noise.
Grooms et al. (2015) and Berardi et al. (2016) provided alternative
methods to deal with the model errors in ensemble Kalman ﬁlters.
Unlike previous studies of assimilating state variables and parametersof unsaturated zone (Li and Ren, 2011; Shi et al., 2015; Brandhorst
et al., 2017), the Richards equation was used as the nonlinear observation operator rather than the model operator. The EnKF framework is used to deal with the nonlinear observation operator. Similarly,
Lei et al. (2019) also used the EnKF to deal with the nonlinear observation operator, i.e., the one-dimensional hydrodynamic model in
their study.
The EnKF method is a Monte Carlo method, in which an initialensemble member can be generated by adding a random disturbance toinitial estimates of the soil hydraulic parameters, and is expressed asfollows:

X 0j = X 0 + ξ

2.3. Gaussian process surrogate
In this section, we brieﬂy introduce the GP method and then propose a GP-based iterative EnKF method. The GP surrogate model approximates the output of the original model F(m) using a Gaussianstochastic process. The surrogate model G(m) is speciﬁed using a meanfunction μ(m) and a covariance function k(m, m'):G(m) N(μ (m), k (m, m'))Suppose we have N training base points of input parameters
B = {m1B, ⋯, mNB } and the outputs of the original model
F(B) = {F(m1B), ⋯, F(mNB )} , where mB are the soil hydraulic parameters
to be estimated. Then an approximate of the original model can beobtained by conditioning G(m) on these base points. For an arbitrary
parameter point m, the output of the surrogate model is represented by
the mean, and the uncertainty of the output is expressed by the variance:

μ|B (m) = μ (m) + k (m, B) k (B, B)−1 [F(B) − μ (B)]σ|2B (m) = k (m, m) − k (m, B) k (B, B)−1k (B, m)nμ (m) = a + ∑i=1

where K

(9)

(10)

J

− i + 1, f
− i + 1, f
1
∑ (Xij+1,f − X )(Xij+1,f − X )J − 1 j=1

(11)

− i + 1, f

is the mean of Xij+ 1, f . The damping factor (Hendricks
in which X
Franssen and Kinzelbach, 2008) was used to reduce the inbreedingproblem, hence the Eq. (9) is revised as:

Xij+ 1, a = Xij+ 1, f + β Ki + 1(Y i + 1 + ε i + 1 − H (Xij+ 1, f ))(17)(1) First, the parameter ensemble Mke is generated randomly using theLatin-Hypercube method for a given mean and standard deviation.
(2) At each assimilation time tk, the training base point parameter MkB
is obtained by being randomly drawn from Mke.
(3) Pressure heads at observation points HB are simulated by running
the water ﬂow model from the initial time to the current time tk
using the training base point parameter MkB.
(4) The GP surrogate is established using MkB as input and HB as
output.
(5) Pressure heads at observation points He are predicted by the GP
surrogate using the parameter ensemble Mke.
(6) The parameter ensemble Mke is updated using Eq. (12), in which
the error term εn+1 = εO + εGP accounts for the observation error
and the GP surrogate error to avoid an over-conﬁdent estimation
result (Zhang et al., 2016);(7) Steps (2) to (6) are repeated until the simulation time reaches the
maximum assimilation time T.

is the Kalman gain, which is deﬁned as:

where Pi+1 is the predictive error covariance matrix, which can be
expressed as follows:

Pi + 1 =(16)

in which a and bi are hyperparameters of the mean function, and σf andλ are hyperparameters of the covariance function. The hyperparametersof GP surrogate model are inferred by minimizing the negative log
marginal likelihood using the exact inference method. More details
about this method can be found in GPML help ﬁles of Matlab toolbox(Rasmussen and Nickisch, 2010).
To improve the eﬃciency of assimilation, we propose a GP-based
Iterative EnKF method (GPIEnKF). The ﬂow chart of the calculation
procedure of the GPIEnKF method is shown in Fig. 1. The steps aredescribed as follows:(8)

Ki + 1 = Pi + 1HT (HPi + 1HT + Oi + 1)−1

bi mi

1
k (m, m') = σ f2 exp ⎛− 2 (m − m')T (m − m')⎞⎝ 2λ⎠in which Xf is the forecast state vector, and Xa is the analysis statevector.
When there are observation data, the analysis state vector is updated using the analysis equation as follows:

i+1

(15)

The performance of the GP surrogate is inﬂuenced by choice ofmean and covariance functions deﬁned in Eq. (13). In this study, theGaussian Process from a Machine Learning (GPML) toolbox developedby Rasmussen and Nickisch (2010) is adopted. A sum of a constant anda linear mean function is used as the mean function, and a diagonalsquared exponential covariance function is used as the covariancefunction of the GP surrogate:

where j indicates the jth member in the parameter ensemble, and ξ isindependent white noise, which follows a multi-normal distribution
with a zero mean and a speciﬁed covariance R.
The calculation of the iterative EnKF method is similar to the EnKF
method. The major diﬀerence between the EnKF and iterative EnKF
methods is the starting time level of rerunning the forward model.
Instead of rerunning the forward model from the previous time level,
the iterative EnKF method reruns it from time zero (initial time).
Here we brieﬂy introduce the EnKF method; more details can befound in Evensen (2003). In the implementation of the EnKF, two
steps–namely forecast and analysis steps–are included. First, the forecast state vector at the i + 1 time level is generated using the model
operator and the analysis state vector at the ith time level:

Xij+ 1, a = Xij+ 1, f + Ki + 1(Y i + 1 + ε i + 1 − H (Xij+ 1, f ))

(14)

and

(7)

Xij+ 1, f = IXij, a

(13)

(12)

where β is the damping factor, with values between 0 and 1. In thisstudy, a series of values e.g. 0.01, 0.05, 0.1, 0.2, 0.5 for β were considered and tested, it was found that β with a value of 0.1 is reasonable.
The same value of β was also used by Brandhorst et al., (2017) forparameter estimation of one-dimensional case.
3

Journal of Hydrology 589 (2020) 125210

K. Liu, et al.

Fig. 1. Flow chart of the calculation procedure of the GPIEnKF algorithm.

3. Case studies
In this section, the performance of the GPIEnKF method for estimating soil hydraulic parameters is evaluated using two synthetic casestudies.

3.1. Case 1: Inﬁltration into a homogeneous soil
In Case 1, inﬁltration into a homogeneous soil under a linear inﬁltration source was considered. The computational domain is
50 cm × 50 cm and was discretized into 200 ﬁnite elements. A constantpressure head with a value of 2 cm was considered as the upperboundary condition in the central 20 cm of transport domain, while freedrainage was considered as the bottom boundary condition (see Fig. 2).
No ﬂow was considered on all other boundaries. The initial condition isa uniformly distributed pressure head with a value of −100 cm. In thisstudy, the initial and boundary conditions were assumed to be known
without uncertainty. This is a common assumption in parameter estimation studies using the EnKF (Li and Ren, 2011; Chen and Zhang,2006; Song et al. 2014). A loamy soil was considered with soil hydraulicparameters according to Carsel and Parrish (1988), i.e.,
Ks = 0.0173 cm/min, α = 0.036 cm−1, n = 1.56, θs = 0.43, andθr = 0.078. The total simulation time is 400 min, and pressure heads at
observation points were recorded every 20 min.
The parameters to be estimated are Ks and α, which are the mostvariable parameters (Carsel and Parrish, 1988). The other soil hydraulicparameters were assumed to be known. Estimated parameters were logtransformed because their probability distributions are known to approximate lognormal distributions (Carsel and Parrish, 1988). To test
the accuracy of the GPIEnKF method for parameter estimation, theFig. 2. Computational domain for Case 1.

initial estimates of Ks and α are larger than their true values, i.e.,
Ks0 = 1.5Ks, α0 = 1.5α. The standard deviation of the initial parameter
ensemble is 0.1 for both log10(Ks) and log10(α). The number of trainingbase points was set to 30. To examine the eﬀect of ensemble size on
estimation accuracy, diﬀerent ensemble sizes, i.e., 50, 100, 200, 300,
4

Journal of Hydrology 589 (2020) 125210

K. Liu, et al.

Fig. 3. Relative errors of the saturated hydraulic conductivity versus an ensemble size.

GP surrogate are in good agreement with the measurements as thenumber of assimilations increases. This is because the estimated parameters converge to their true values during the data assimilation,
thereby resulting in good pressure head predictions by the GP surrogate. Since the measurements were generated by the original modelusing the true parameters, the results indicate that the GP surrogate canpredict the state variables of the original model well. We compared theprediction accuracy at observation point 3 for diﬀerent layout schemes,
as the location of this point was the same in both the horizontal andvertical layout schemes. Pressure heads predicted using the verticallayout scheme are slightly more precise than those predicted using thehorizontal layout scheme, with the RMSE values of 6.56 cm and6.59 cm, respectively.

400, 500, and 600 were considered. The mean values of the parameter
ensembles at the maximum assimilation time were the ﬁnal parameter
estimations. As shown in Fig. 3, the relative error of the saturated hydraulic conductivity decreases signiﬁcantly as the ensemble size increases from 50 to 300, and then it declines only slightly with a further
increasing ensemble size. Therefore, the ensemble size of 300, the sameas also used by Li and Ren (2011), was used in the following simulations.
The measurements were obtained using the following steps.
Pressure heads at observation points were ﬁrst simulated using the true
values of soil hydraulic parameters. The measurements were then obtained by adding observation errors to the simulated pressure heads.
The standard deviation of the observation error was 0.01.
The inﬂuence of a layout scheme of observation points, the numberof training base points, and the number of observation points onparameter estimation were analyzed. The relative error (RE) was used
as an indicator to evaluate the accuracy of parameter estimation. The
root mean square error (RMSE) was used to evaluate the accuracy ofpressure head predictions.

3.1.2. Inﬂuence of the number of training base points
To examine the inﬂuence of the number of training base points (NB)
on predicted results, diﬀerent numbers of training base points, i.e.,
NB = 10, 20, 30, … , 100, were tested. Fig. 7 shows that parameter
evolutions for diﬀerent NBs are similar. Fig. 8 presents the RE for different NBs. For the vertical layout scheme, the smallest RE (0.13%) was
obtained for NB = 40 and the largest RE (0.96%) for NB = 30. For thehorizontal layout scheme, the smallest RE (0.36%) was obtained for
NB = 40 and the largest RE (2.07%) for NB = 100. The results showthat increasing the number of training base points does not improve theaccuracy of parameter estimation. This may be due to the homogeneityof the soil proﬁle in this case study, in which a small number of basepoints can produce an accurate prediction of pressure heads at observation points. In contrast, Ju et al. (2018) reported that a larger
number of base points was required to obtain accurate simulations
when a strongly heterogeneous ﬁeld was considered.
The comparison of GP surrogate predictions with original model
simulations for diﬀerent NB is shown in Fig. 9. It can be seen thatpressure heads predicted by the GP surrogate match those simulated by
the original model reasonably well. However, the prediction accuracyof the GP surrogate is not improved by increasing NB. The RMSE values
for NB = 30, 60, and 100 are 2.65, 2.69, and 2.86 cm, respectively.
This may be due to the diﬀerence in the initial parameter ensemble andtraining base points for diﬀerent NBs. Since the initial parameter ensemble was randomly generated, even though the same mean andstandard deviation were used, there is no guarantee that every elementof the parameter ensemble is identical for diﬀerent NBs. Apart from
NB = 30, the parameter estimation accuracy of the vertical layout

3.1.1. Inﬂuence of the layout scheme of observation points
Two types of layout schemes of observation points, namely verticaland horizontal layout schemes, were considered. For the horizontal
layout the observation points were located in upper soil, whereas forvertical layout the observation points were located at various depths.
The locations of observation points are shown in Fig. 4. Fig. 5 presents
the evolutions of the mean of the parameter ensemble for diﬀerent
observation point layout schemes. Fig. 5 shows that the optimized
parameters converged to their true values for both vertical and horizontal layout schemes as the number of assimilations increased. The
standard deviation of the optimized parameters decreases gradually as
the number of assimilations increases, indicating that the uncertainty ofparameters decreases. The accuracy of estimated parameters for thehorizontal layout scheme is slightly higher for the vertical layout
scheme. The relative errors (RE) is less than 1% (0.96% and 0.87% for
the horizontal and vertical layout schemes are 0.96% and 0.87%, respectively), which meets the accuracy requirements. The results indicate that the soil hydraulic parameters can be well estimated by the
proposed method.
Fig. 6 shows pressure heads at observation points measured andpredicted by the GP surrogate for the horizontal and vertical layout
schemes. The ﬁgure demonstrates that pressure heads predicted by the5Journal of Hydrology 589 (2020) 125210

K. Liu, et al.

Fig. 4. Layout schemes [(a) a horizontal layout, (b) a vertical layout] of observation points (ﬁlled circles represent observation points, red numbers indicate
observation points discussed in the text).

Fig. 5. Parameter evolutions (Ks - top, α - bottom) for diﬀerent layout schemes (vertical and horizontal) of observation points (V-stdev and H-stdev are standard
deviations for the vertical and horizontal layouts, respectively).
6

Journal of Hydrology 589 (2020) 125210

K. Liu, et al.

Fig. 6. Comparison of pressure heads at observation points measured and predicted using the GP surrogate for diﬀerent layout schemes [(a) a horizontal layout, (b) avertical layout].

Fig. 7. Parameter evolutions (Ks - left, α - right) for a diﬀerent number of training base points (NB) and diﬀerent layout schemes [(a) a horizontal layout, (b) a verticallayout].
7

Journal of Hydrology 589 (2020) 125210

K. Liu, et al.

3.2. Case 2: Inﬁltration into a layered soil
In Case 2, water inﬁltration from a single-ring inﬁltrometer into alayered soil was considered. The soil proﬁle consists of two layers:
40 cm thick A-horizon and an underlying B/C-horizon. The hydraulicparameters of the two-layered soil are shown in Table 1 (Cislerova,
1987; Hopmans and Stricker, 1989). The axisymmetric computational
domain is shown in Fig. 12a. The computational domain is
125 cm × 130 cm and was discretized into 684 ﬁnite elements. Allsides of the computational domain are impervious except for a small
portion around the origin at the soil surface where a constant pressure
head was imposed, as well as the lower right corner where thegroundwater level was kept constant. The initial pressure head proﬁleand locations of observation points are shown in Fig. 12b. The ensemble
size was set to 300, and the number of training base points was set to
30. The total simulation time was 10 days, and pressure heads at observation points were recorded every 2 h. The measurements wereobtained by adding an observation error to the pressure heads simulated using the true parameters. The standard deviation of the observation error is 0.01. More hydraulic parameters were estimated inthis case, i.e., Ks, α, and n of each soil layer. Similarly, as in Case 1, theinitial estimates of parameters deviate from their true values, i.e.,
Ks0 = 1.5Ks, α0 = 1.5α, and n0 = 1.2n. The standard deviation of theinitial parameters ensemble is 0.1 for both log10(Ks) and log10(α) and
0.01 for n. To illustrate the accuracy and eﬃciency of the proposed
method, the results of the GPIEnKF method were compared with thoseof the original Iterative Ensemble Kalman Filter method (IEnKF).
Fig. 13 presents the parameter evolutions of the GPIEnKF and IEnKF
methods. It can be seen that Ks and α estimated by the GPIEnKF method
are closer to their true values at the beginning of the simulation incomparison with the IEnKF method. This may be caused by the difference in pressure heads predicted by the surrogate and originalmodels. The ﬁnal values of estimated parameters by the IEnKF methods
are more accurate than those by the GPIEnKF method, with the RE
value of 3.8% and 4.4% for the IEnKF and GPIEnKF methods, respectively. For n, the parameters estimated by the IEnKF converge to theirtrue values more quickly. This is because the simulation of water
movement is more sensitive to Ks and α than to n in this case. The Ks andα values estimated by the GPIEnKF are more accurate at the beginningof the simulation. Even when the n value deviates from its true value,
the simulated pressure heads can match the measurements at the beginning of the simulation. The ﬁnal estimates of n by the IEnKF are
more accurate. The RE value for the ﬁrst soil layer is 1.9% and 2.3% for
the IEnKF and GPIEnKF methods, respectively, whereas the RE value forthe second soil layer is 4.5% and 4.8% for the IEnKF and GPIEnKF
methods, respectively. The parameters estimated using IEnKF and
GPIEnKF did not completely converge to their true values. This can be
attributed to two factors. First, diﬀerent combinations of parameters

Fig. 8. Relative errors (RE) for a diﬀerent number of training base points anddiﬀerent layout schemes.

scheme is higher than that of the horizontal layout scheme (see Fig. 8).
This is likely because water movement during the inﬁltration process isbetter described using the vertical layout scheme of observation points,
resulting in better emulation of the original model by the GP surrogate.
The results are consistent with Montzka et al. (2011), who reported that
assimilation eﬀects are limited when only soil surface observations
were used. These results indicate that the layout scheme of observation
points has a greater inﬂuence on parameter estimation than the numberof training base points. Increasing the number of training base points
leads to longer training time of the surrogate model. The CPU timeincreased from 49 s to 168 s when NB was increased from 10 to 100.
Since the RE for NB = 30 is less than 1%, which meets the accuracy
requirement, in the subsequent experiments, NB was set to 30.
3.1.3. Inﬂuence of the number of observation points
To analyze the inﬂuence of the number of observation points (NO),a diﬀerent number of observation points, i.e., NO = 8, 12, and 16, were
tested. The locations of observation points are shown in Fig. 10. Fig. 11
presents parameter evolutions for diﬀerent NOs. The estimation accuracy can be improved for a higher NO, with the RE decreasing from
0.87% to 0.56%. For the same observation points, the RMSE of predicted pressure heads decreases from 5.95 to 5.82 cm. The results are
consistent with those of Brandhorst et al. (2017), who reported thatusing a higher number of observation points can improve the accuracyof parameter estimation. This is because the assimilation system canutilize more observation information, which allows the GP surrogate tobetter emulate the original model, thereby resulting in more accurate
parameter estimations and pressure head predictions. A higher NO does
not require signiﬁcantly longer calculation time. Since the number oforiginal model evaluations is the same for diﬀerent NOs, the training
time of the surrogate model does not grow much as NO increases.
Hence the total calculation time for diﬀerent NOs is almost the same.

Fig. 9. Comparison of pressure heads predicted using the GP surrogate model and simulated using the original model for a diﬀerent number of training base points(NB = 30, 60, and 100).
8

Journal of Hydrology 589 (2020) 125210

K. Liu, et al.

Fig. 10. Layout schemes for 8, 12, and 16 observation points (ﬁlled circles represent observation points).

Fig. 11. Parameter evolutions (Ks - left, α - right) for a diﬀerent number ofobservation points (NO).
Table 1
Soil hydraulic parameters for Case 2.
Parameter

1st layer

2nd layerθsθrKs (cm/min)α (cm−1)
n (-)

0.399
0.00010.02070.01741.376

0.399
0.00010.03150.01391.602

Fig. 12. Computational domain for Case 2 (a) and the initial pressure headproﬁle and locations of observation points (b). (Filled circles represent observation points, red numbers indicate observation points discussed in the text).
.
Adapted from Šimůnek et al. (2011)

methods, respectively. This is because parameters estimated by the
IEnKF method are more accurate than those estimated by the GPIEnKF
method.
It should be pointed out that the two most time-consuming steps ofassimilation with the IEnKF method are original model evaluations andmatrix calculations. In this study, the updated state vector includes soil
hydraulic parameters, and is thus low dimensional, requiring less time
for matrix calculations. Since the time needed for training of the surrogate model and its predictions is almost negligible, most computational time is consumed in original model evaluations. The computational eﬃciency of the GPIEnKF and IEnKF methods can be evaluatedby comparing the number of original model evaluations. In this case,Note: adapted from Šimůnek et al. (2011).

may lead to similar simulation results. Second, the errors of the numerical solution for Richards equation may aﬀect the results of parameter estimation. However, the RE values were less than 5%. The results indicate that the GPIEnKF method can obtain a comparable
accuracy of parameter estimations as the IEnKF method. Fig. 14 shows
the measured and predicted pressure heads at two observation points.
The pressure heads predicted by the GPIEnKF and IEnKF methods are ingood agreement with the measurements. The prediction accuracy of the
IEnKF method is slightly higher than that of the GPIEnKF method, with
the RMSE value of 0.66 and 1.00 cm for both the IEnKF and GPIEnKF
9Journal of Hydrology 589 (2020) 125210

K. Liu, et al.

Fig. 13. Comparison of parameter evolutions (Ks - top, α- middle, and n - bottom) for the GPIEnKF and IEnKF methods [(a) the ﬁrst layer, (b) the second layer].

Fig. 14. Comparison of pressure heads measured and predicted by the GPIEnKF and IEnKF methods.

10Journal of Hydrology 589 (2020) 125210

K. Liu, et al.

40–60 cm, respectively.
An ensemble size of 300, a number of training base points of 30, anda 2-day interval were used in the following assimilation. Unlike in theprevious two synthetic cases, soil water contents rather than pressure
heads were used as observation data. A value of 0.01 was used as thestandard deviation for the soil water content observation error (Shiet al., 2015). The initial values of soil hydraulic parameters to be estimated (θr, θs, Ks, α, and n) for each soil layer (see Table 3) were obtained using the Rosetta pedotransfer functions (Schaap et al., 2001)
with the bulk density and the percentage of sand, silt, and clay as input.
The standard deviation of the initial parameter ensemble was 0.3 forlog10(KS) and log10(α), and 0.01 for θr, θs, and n (Brandhorst et al.,2017). The initial condition error was also considered by adding a
disturbance σh = 0.005 m (Brandhorst et al., 2017) to the initialpressure head, which was calculated from the initial soil water content.
The results of the parameter estimation are shown in Fig. 17. Deviations can be found between the initial parameter values and theirestimated values, especially for the most variable parameters KS and α(Carsel and Parrish, 1988). Since the true values of parameters were notknown, the results of the GPIEnKF method were compared with those of
the Levenberg-Marquardt optimization algorithm (L-M) (Moré, 1978),which are shown in Table 4. It can be seen that the values of theparameters KS and n estimated using the GPIEnKF method (see Table 5)
are close to those obtained using the L-M method, contrary to those ofother parameters, which showed large deviations. This may be attributed to the fact that the L-M method estimates parameters by minimizing the objective function, i.e., the sum of squared deviations between simulated and observed soil water contents, while the GPIEnKF
method obtains parameters using the variance-minimizing analysis(Evensen, 2003).
Fig. 18 shows a comparison between the simulated and observed
soil water contents. Simulated soil water contents agree reasonably wellwith measurements at soil depths of 20–60 cm. RMSEs for the GPIEnKF
method are 0.012 cm3/cm3 for the 20–40 cm soil depth and 0.031 cm3/
cm3 for the 40–60 cm depth, while for the L-M method they are
0.007 cm3/cm3 for the 20–40 cm depth and 0.016 cm3/cm3 for the
40–60 cm depth. However, simulated soil water contents for the
0–20 cm soil depth show large deviations from measurements, with
RMSEs of 0.049 cm3/cm3 for the GPIEnKF method and 0.035 cm3/cm3
for the L-M method. Larger deviations between the simulated and observed soil water contents may be attributable to the following reasons.
First, there is uncertainty in the upper boundary condition, which can
not be quite well characterized due to frequent changes between evaporation and inﬁltration associated with multiple irrigation and precipitation events. Second, the simulations are aﬀected by an inherent
soil spatial heterogeneity. The classiﬁcation of the soil proﬁle into
multiple homogeneous layers in numerical simulations is a signiﬁcantsimpliﬁcation of the complex reality. Finally, the assumption thatparameter values are constant in time may be unreasonable. Soil
structure is continuously changing during the growing season due toirrigations, crop growth, and other factors, which may lead to timevariable soil hydraulic parameters (especially KS and α).
More accurate simulations can be obtained using parameters estimated by the GPIEnKF method compared to those estimated by Rosetta(initial parameter values) (Fig. 18). This further indicates that the assimilation algorithm can improve prediction accuracy. Similar results

the ensemble size was set to 300, and the number of training basepoints for the GP surrogate was set to 30. Therefore, the computational
cost of the GPIEnKF method is approximately 1/10th of that of the
IEnKF method. Compared with the IEnKF method, the GPIEnKF method
has much higher calculation eﬃciency with almost no loss of accuracy.
Simulated pressure heads vary smoothly without abnormal ﬂuctuations (see Fig. 14), which indicates that the iterative method is aneﬀective way to avoid the inconsistency problem. The results show that
the GP surrogate can accurately predict pressure heads at observation
points similarly to the original model. Hence, for high-dimensional soil
water ﬂow problems, when only state variables at critical locations
need to be predicted, the GP surrogate provides an eﬃcient forecastingtool.
Though the properties of the two soils used in the numerical simulation did not show big diﬀerences, the GPIEnKF method with
HYDRUS-2D could also be used to estimate hydraulic parameters forlayered soils with signiﬁcantly diﬀerent properties. This is because similar methods have been eﬀectively used to assimilate the hydraulicparameters. For example, Li and Ren (2011) and Shi et al. (2015) used
the EnKF method with one-dimensional HYDRUS to assimilate soil
hydraulic parameters of layered soils of which the hydraulic conductivities diﬀered by one order of magnitude. In addition, Ju et al.
(2018) used the GP surrogate model to approximate the soil hydrodynamic model, i.e., HYDRUS-2D, and demonstrated that the GP surrogate model based method can improve the assimilation eﬃciency ofheterogeneous hydraulic conductivity ﬁeld.
3.3. Case 3: A real case involving a ﬁeld experiment with drip irrigation
A ﬁeld experiment involving drip irrigation of processing tomato
was conducted at the upper reaches of the Yellow River basin, China(latitude 40◦41‘N and longitude 107◦18′E, 1041 m altitude). The soil
proﬁle (down to a depth of 0–60 cm) can be divided into three layers.
Physical properties of these soil layers are given in Table 2. Meteorological data measured by an automatic weather station (HOBO,
Campbell Scientiﬁc Inc., USA) were used to calculate daily reference
evapotranspiration (ETo) using the modiﬁed Penman-Monteith equation (Allen et al., 1998). Daily ETo and precipitation are shown in
Fig. 15. Soil water contents at two horizontal locations and three depths(Fig. 16) were monitored using HydraProbe sensors.
Soil water contents were simulated for a computational domain
40 cm wide and 60 cm deep during the growing season (6/2/2015–9/
1/2015, a total of 104 days). The computational domain, appliedboundary conditions (BC), and the position of observation points areshown in Fig. 16. The computational domain was discretized into 2,527
ﬁnite elements. A time-variable ﬂux boundary condition, determinedfrom irrigation rates, was speciﬁed on the BC boundary (Fig. 16). Whilea no ﬂux boundary condition was speciﬁed on the CD boundary, which
was covered by a plastic mulch, an atmospheric boundary condition
was speciﬁed on the DE boundary. Since the groundwater table is relatively deep, the lower boundary was treated as a free drainage
boundary. Due to symmetry of ﬂow, a no ﬂux boundary condition was
speciﬁed on both vertical sides. Water uptake was simulated using themodel of Feddes et al. (1978) with default values of parameters of thestress response function. The initial water content was 0.4 cm3/cm3,
0.42 cm3/cm3, and 0.43 cm3/cm3 in depths of 0–20 cm, 20–40 cm, andTable 2
Soil physical properties for the experimental site.
Soil depths (cm)

0–20
20–40
40–60

Soil particle fraction (%)
Sand > 0.05 mm

Silt 0.002–0.05 mm

Clay < 0.002 mm

12.18
21.84
18.46

60.98
56.87
54.81

26.84
21.29
26.73

11

Soil texture

Bulk density (g/cm3)

Silt loam
Silt loam
Silt loam

1.47
1.53
1.49

Journal of Hydrology 589 (2020) 125210

K. Liu, et al.

Fig. 15. Daily reference evapotranspiration (ETo) and precipitation during the growing season.

investigated in the follow up studies. For example, the uncertainties of
initial and boundary conditions might be integrated in surrogate
method with considering them as the input variables of the constructed
model. In addition, the present study focused on parameter estimationof layered soils, the capability of GPIEnKF for the parameter estimationof multi-dimensional cases with high heterogeneities needs further investigation.

4. Conclusions
We propose a Gaussian Process-based Iterative Ensemble Kalman
Filter method for parameter estimation. The proposed method combines the advantage of the high computational eﬃciency of the GP
surrogate with the advantage of dealing with the uncertainty problem
by the EnKF method. Since the number of original model evaluations isequal to the number of training base points, which is much smaller thanthe ensemble size, the calculation time of the GPIEnKF method can begreatly reduced in comparison with the standard Iterative Ensemble
Kalman Filter method. We veriﬁed the accuracy and eﬃciency of the
proposed method using two synthetic cases and a real dataset involving
a ﬁeld experiment with drip irrigation. For a homogeneous soil proﬁle,the results show that the number of training base points has little eﬀect
on parameter estimation. Soil water movement needs to be described
well by the layout scheme of observation points so that the surrogate
model can correctly emulate the original model, thereby leading tomore accurate parameter estimation and pressure head predictions.
Increasing the number of observation points provides more observation
information for the assimilation system and improves the prediction
accuracy of the surrogate model, thereby resulting in more accurate
parameter estimation. For the layered soils, the GPIEnKF method canobtain accurate results for multi-parameter conditions. However, theaccuracy of predictions is aﬀected by the uncertainty of the boundary
conditions and the soil spatial heterogeneity. In comparison with thestandard Iterative Ensemble Kalman Filter method, the GPIEnKF
method can obtain comparably accurate results with a greatly reduced
computational cost. The proposed method is an eﬀective tool forparameter estimation for high-dimensional soil water ﬂow problems. In
this study, we only considered unsaturated ﬂow model, but the proposed method can also be easily adapted for other hydrologic models.
The parameters were assumed to be constant over time in this study.
However, due to climatic variations and human activities, this assumption may result in large simulation errors. As a sequential assimilation method, the GPIEnKF method can update the state variables andparameters in real time according to observed data. In future research,

Fig. 16. The computational domain, boundary conditions, and observation
points for the ﬁeld case.
Table 3
Parameters of the van Genuchten–Mualem model for diﬀerent soil layers.
Soil depths (cm)θr (cm3/cm3)

θs (cm3/cm3)

α (cm−1)n (-)

Ks (cm/day)

0–20
20–40
40–60

0.08150.07030.0790

0.456
0.431
0.447

0.00660.00530.0066

1.569
1.619
1.562

15.37
18.25
15.82

were reported by Li and Ren (2011) and Shi et al. (2015). For the realexperimental case, the CPU time for the GPIEnKF method (about 15 h)
is about 10 times shorter than that for the standard iterative EnKF. Thisfurther indicates that the GPIEnKF method can improve the computational eﬃciency, especially for cases, with which the original model
simulations are time-consuming.
The initial and boundary conditions are very important to the simulation of soil water dynamics. However, their uncertainties were notconsidered in the surrogate model of the present study, due to that theiruncertainties mainly depend on the accuracy of observations. In general, the initial and boundary conditions are assumed to be deterministic and known in the parameter estimation using assimilation
method (Li and Ren, 2011; Song et al., 2014; Shi et al., 2015). It should
be pointed out that the inﬂuence of initial and boundary conditions onparameter estimation with the surrogate method is worth to be
12

Journal of Hydrology 589 (2020) 125210

K. Liu, et al.

Fig. 17. Parameter evolutions at depths 0–20 cm, 20–40 cm, and 40–60 cm for the real ﬁeld case (GP denotes the GPIEnKF method, Initial denotes an initial
parameter value).

we will explore the potential of the GPIEnKF method for estimatingtime-varying parameters.

Table 4
Parameters of the van Genuchten–Mualem model estimated using theLevenberg-Marquardt optimization algorithm.
Soil depths (cm)θr (cm3/cm3)

θs (cm3/cm3)

α (cm−1)n (–)

Ks (cm/day)

CRediT authorship contribution statement

0–20
20–40
40–60

0.06380.09100.0250

0.510
0.375
0.473

0.01290.00340.0053

1.523
1.655
1.463

13.83
15.41
7.57

Kun Liu: Conceptualization, Methodology, Writing - original draft,
Software. Guanhua Huang: Supervision, Writing - review & editing.
Zheng Jiang: Methodology, Software. Xu Xu: Conceptualization,
Methodology. Yunwu Xiong: Conceptualization, Formal analysis.
Quanzhong Huang: Conceptualization, Methodology. Jiří Šimůnek:Supervision, Writing - review & editing.

Table 5
Parameters of the van Genuchten–Mualem model estimated using the GPIEnKF.
Soil depths (cm)θr (cm3/cm3)

θs (cm3/cm3)

α (cm−1)n (-)

Ks (cm/day)

0–20
20–40
40–60

0.08150.07090.0796

0.45630.43160.44710.00690.00560.0074

1.575
1.624
1.568

9.96
11.83
10.25

Declaration of Competing Interest
The authors declare that they have no known competing ﬁnancial
interests or personal relationships that could have appeared to inﬂuence the work reported in this paper.
Acknowledgments
This research was jointly supported by the National Key R&D
Program of China (No. 2017YFC0403301) and National Natural Science
13

Journal of Hydrology 589 (2020) 125210

K. Liu, et al.

Fig. 18. Comparison of soil water contents measured and simulated using the Gaussian process based IEnKF method (GPIEnKF), the Levenberg-Marquardt method (LM) and the initial (Rosetta-estimated) parameter values (Original) [(a) 0–20 cm, (b) 20–40 cm, (c) 40–60 cm].

14

Journal of Hydrology 589 (2020) 125210

K. Liu, et al.

observations into the Richards equation - Part 2: A dual ﬁlter approach for simultaneous retrieval of states and parameters. Hydrol. Earth Syst. Sc. 18, 2521–2541.
Montzka, C., Moradkhani, H., Weihermüller, L., Franssen, H.J.H., Canty, M., Vereecken,
H., 2011. Hydraulic parameter estimation by remotely-sensed top soil moisture observations with the particle ﬁlter. J. Hydrol. 399 (3–4), 410–421.
Moradkhani, H., Sorooshian, S., Gupta, H.V., Houser, P.R., 2005. Dual state–parameter
estimation of hydrological models using ensemble Kalman ﬁlter. Adv. Water Resour.
28 (2), 135–147.
Moré, J.J., 1978. The Levenberg-Marquardt algorithm: implementation and theory. In:
Numerical analysis. Springer, Berlin, Heidelberg, pp. 105–116.
Mualem, Y., 1976. A new model for predicting the hydraulic conductivity of unsaturated
porous media. Water Resour. Res. 12 (3), 513–522.
Pan, M., Wood, E.F., 2006. Data assimilation for estimating the terrestrial water budget
using a constrained ensemble Kalman ﬁlter. J. Hydrometeorol. 7 (3), 534–547.
Pasetto, D., Camporese, M., Putti, M., 2012. Ensemble Kalman ﬁlter versus particle ﬁlter
for a physically-based coupled surface-subsurface model. Adv. Water Resour. 47,
1–13.
Pauwels, V.R., Verhoest, N.E., De Lannoy, G.J., Guissard, V., Lucau, C., Defourny, P.,
2007. Optimization of a coupled hydrology–crop growth model through the assimilation of observed soil moisture and leaf area index values using an ensemble Kalmanﬁlter. Water Resour. Res. 43, W04421. https://doi.org/10.1029/2006WR004942.
Rasmussen, C.E., Nickisch, H., 2010. Gaussian processes for machine learning (GPML)
toolbox. J. Mach. Learn. Res. 11 (Nov), 3011–3015.
Razavi, S., Tolson, B.A., Burn, D.H., 2012. Review of surrogate modeling in water resources. Water Resour. Res. 48, W07401. https://doi.org/10.1029/2011WR011527.
Schaap, M.G., Leij, F.J., van Genuchten, M.T., 2001. ROSETTA: a computer program forestimating soil hydraulic parameters with hierarchical pedotransfer functions. J.
Hydrol. 251 (3), 163–176.
Shi, L., Song, X., Tong, J., Zhu, Y., Zhang, Q., 2015. Impacts of diﬀerent types of measurements on estimating unsaturated ﬂow parameters. J. Hydrol. 524, 549–561.
Šimůnek, J., van Genuchten, M.T., Gribb, M.M., Hopmans, J.W., 1998. Parameter estimation of unsaturated soil hydraulic properties from transient ﬂow processes. Soil
Tillage Res. 47 (1–2), 27–36.
Šimůnek, J., van Genuchten, M.T., Šejna, M., 2008. Development and applications of the
HYDRUS and STANMOD software packages and related codes. Vadose Zone J. 7 (2),587–600.
Šimůnek, J., van Genuchten, M.T., and Šejna, M., 2011. The HYDRUS software package
for simulating the two-and three-dimensional movement of water, heat, and multiple
solutes in variably-saturated porous media. Technical Manual, Version 2, PC
Progress, Prague, Czech Republic, 258.
Song, X., Shi, L., Ye, M., Yang, J., Navon, I.M., 2014. Numerical comparison of iterative
ensemble Kalman ﬁlters for unsaturated ﬂow inverse modeling. Vadose Zone J. 13(2). https://doi.org/10.2136/VZJ2013.05.0083.
van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils 1. Soil Sci. Soc. Am. J. 44 (5), 892–898.
Vrugt, J.A., Gupta, H.V., Bouten, W., Sorooshian, S., 2003. A Shuﬄed Complex Evolution
Metropolis algorithm for optimization and uncertainty assessment of hydrologic
model parameters. Water Resour. Res. 39, 1201. https://doi.org/10.1029/
2002WR001642, 8.
Vrugt, J.A., Robinson, B.A., 2007. Improved evolutionary optimization from genetically
adaptive multimethod search. Proc. Nat. Acad. Sci. 104 (3), 708–711.
Vrugt, J.A., Stauﬀer, P.H., Wöhling, T., Robinson, B.A., Vesselinov, V.V., 2008. Inverse
modeling of subsurface ﬂow and transport properties: A review with new developments. Vadose Zone J. 7 (2), 843–864.
Wang, D., Chen, Y., Cai, X., 2009. State and parameter estimation of hydrologic models
using the constrained ensemble Kalman ﬁlter. Water Resour. Res. 45 (11).
Wöhling, T., Vrugt, J.A., Barkle, G.F., 2008. Comparison of three multiobjective optimization algorithms for inverse modeling of vadose zone hydraulic properties. Soil Sci.
Soc. Am. J. 72 (2), 305–319.
Zeng, L., Chang, H., Zhang, D., 2011. A probabilistic collocation-based Kalman ﬁlter for
history matching. SPE J. 16 (02), 294–306.
Zhang, J., Li, W., Zeng, L., Wu, L., 2016. An adaptive Gaussian process-based method foreﬃcient Bayesian experimental design in groundwater contaminant source identiﬁcation problems. Water Resour. Res. 52 (8), 5971–5984.
Zhang, Y., Schaap, M.G., Zha, Y., 2018. A High-Resolution Global Map of Soil Hydraulic
Properties Produced by a Hierarchical Parameterization of a Physically Based Water
Retention Model. Water Resour. Res. 54 (12), 9774–9790.

Foundation of China (Nos. 51639009, 51621061, 51125036).
References
Abbasi, F., Šimůnek, J., Feyen, J., van Genuchten, M.T., Shouse, P.J., 2003. Simultaneous
inverse estimation of soil hydraulic and solute transport parameters from transientﬁeld experiments: Homogeneous soil. Trans. ASAE. 46 (4), 1085.
Allen, R.G., Pereira, L.S., Raes, D., Smith, M., 1998. Crop evapotranspiration: guidelines
for computing crop water requirements. In: FAO Irrigation and Drainage Paper No.
56. FAO, Rome, Italy.
Berardi, M., Andrisani, A., Lopez, L., Vurro, M., 2016. A new data assimilation techniquebased on ensemble Kalman ﬁlter and Brownian bridges: An application to Richards’
equation. Comput. Phys. Commun. 208, 43–53.
Brandhorst, N., Erdal, D., Neuweiler, I., 2017. Soil moisture prediction with the ensemble
Kalman ﬁlter: Handling uncertainty of soil hydraulic parameters. Adv. Water Resour.
110, 360–370.
Brunetti, G., Šimůnek, J., Turco, M., Piro, P., 2017. On the use of surrogate-based modeling for the numerical analysis of Low Impact Development techniques. J. Hydrol.
548, 263–277.
Carsel, R.F., Parrish, R.S., 1988. Developing joint probability distributions of soil water
retention characteristics. Water Resour. Res. 24 (5), 755–769.
Chang, H., Liao, Q., Zhang, D., 2017. Surrogate model based iterative ensemble smootherfor subsurface ﬂow data assimilation. Adv. Water Resour. 100, 96–108.
Chen, Y., Zhang, D., 2006. Data assimilation for transient ﬂow in geologic formations via
ensemble Kalman ﬁlter. Adv. Water Resour. 29 (8), 1107–1122.
Cislerova, M., 1987. Comparison of simulated water balance for ordinary and scaled soil
hydraulic characteristics (No. Landbouwuniversiteit Wageningen, pp. 82)..
De Lannoy, G.J., Houser, P.R., Pauwels, V.R., Verhoest, N.E., 2007. State and bias estimation for soil moisture proﬁles by an ensemble Kalman ﬁlter: Eﬀect of assimilation
depth and frequency. Water Resour. Res. 43, W06401. https://doi.org/10.1029/
2006WR005100.
Duan, Q., Sorooshian, S., Gupta, V., 1992. Eﬀective and eﬃcient global optimization forconceptual rainfall-runoﬀ models. Water Resour. Res. 28 (4), 1015–1031.
Durner, W., Jansen, U., Iden, S.C., 2008. Eﬀective hydraulic properties of layered soils at
the lysimeter scale determined by inverse modelling. Eur. J. Soil Sci. 59 (1), 114–124.
Evensen, G., 1994. Sequential data assimilation with a nonlinear quasi-geostrophic modelusing Monte Carlo methods to forecast error statistics. J. Geophys. Res.-Oceans. 99(C5), 10143–10162.
Evensen, G., 2003. The ensemble Kalman ﬁlter: Theoretical formulation and practicalimplementation. Ocean Dyn. 53 (4), 343–367.
Feddes, R.A., Kowalik, P.J., Zaradny, H., 1978. Simulation of Field Water Use and Crop
Yield. Centre for Agricultural Publishing and Documentation, Wageningen, pp. 189.
Grooms, I., Lee, Y., Majda, A.J., 2015. Ensemble ﬁltering and low-resolution model error:
covariance inﬂation, stochastic parameterization, and model numerics. Mon Weather
Rev. 143 (10), 3912–3924.
Hendricks Franssen, H.J., Kinzelbach, W., 2008. Real-time groundwater ﬂow modelingwith the ensemble Kalman ﬁlter: Joint estimation of states and parameters and theﬁlter inbreeding problem. Water Resour. Res. 44, W09408. https://doi.org/10.1029/
2007WR006505.
Hopmans, J.W., Stricker, J.N.M., 1989. Stochastic analysis of soil water regime in awatershed. J. Hydrol. 105 (1–2), 57–84.
Hopmans, J. W., Šimůnek, J., 1999. Review of inverse estimation of soil hydraulicproperties. p. 643–659. In M. Th. van Genuchten et al. (ed.) Characterization andmeasurement of the hydraulic properties of unsaturated porous media. Proc. Int.
Workshop, Riverside, CA. 22–24 Oct. 1997. U.S. Salinity Lab., Riverside, CA.
Ju, L., Zhang, J., Meng, L., Wu, L., Zeng, L., 2018. An adaptive Gaussian process-based
iterative ensemble smoother for data assimilation. Adv. Water Resour. 115, 125–135.
Soil Science Society of America Journal 49 (6), 1348–1354. https://doi.org/10.2136/
sssaj1985.03615995004900060004x.
Li, C., Ren, L., 2011. Estimation of unsaturated soil hydraulic parameters using the ensemble Kalman ﬁlter. Vadose Zone J. 10 (4), 1205–1227.
Li, J., Xiu, D., 2009. A generalized polynomial chaos based ensemble Kalman ﬁlter withhigh accuracy. J. Comput. Phys. 228 (15), 5454–5469.
Lei, X., Tian, Y., Zhang, Z., Wang, L., Xiang, X., Wang, H., 2019. Correction of pumping
station parameters in a one-dimensional hydrodynamic model using the Ensemble
Kalman ﬁlter. J. Hydrol. 568, 108–118.
Medina, H., Romano, N., Chirico, G., 2014. Kalman ﬁlters for assimilating near-surface

15