Journal of Hydrology 594 (2021) 125729

Contents lists available at ScienceDirect

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

Research papers

Comparison of ensemble data assimilation methods for the estimation oftime-varying soil hydraulic parameters
Kun Liu a, b, Guanhua Huang a, b, *, Jiří Šimůnek c, Xu Xu a, b, Yunwu Xiong a, b,
Quanzhong Huang a, babcCenter 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, China
Department of Environmental Sciences, University of California Riverside, Riverside, CA 92 521, USA

A R T I C L E I N F O

A B S T R A C T

This manuscript was handled by Huaming Guo,
Editor-in-Chief

The hydraulic properties of the soil top layer may change during the growth period due to various factors such as
wetting and drying cycles, tillage practices, and crop root growth. In this study, the potential of the assimilation
method to estimate time-varying soil hydraulic parameters is explored. Four assimilation schemes, including thesimultaneous update state augmentation method, the partitioned update state augmentation method, thesimultaneous update parameter correction method, and the partitioned update parameter correction method, arecompared. The performance of four assimilation schemes on parameter estimations and soil moisture simulations
is tested first using the synthetic case. The influence of initial parameter values and a parameter update order on
assimilation is also analyzed. Finally, the partitioned update parameter correction method is applied to a realcase involving a field drip irrigation experiment. The results show that when the analyzed parameter has either aperiodical or a linear variation, there is a time lag between the assimilation value and the true value. The
assimilation method can respond immediately to an abrupt change of the parameter value. Using the simulta­
neous update method leads to an obvious parameter correlation problem. In contrast, the partitioned update
method can relieve the parameter correlation problems, thereby improving the accuracy of parameter estima­
tions and pressure head simulations. However, when the initial values of the parameters deviate from their true
values to a certain extent, the partitioned update method cannot obtain accurate parameter estimations.
Compared with the partitioned update state augmentation method, the partitioned update parameter correction
method is not sensitive to the parameter update order. The partitioned update parameter correction method has
higher computational efficiency and assimilation stability, and it can obtain more accurate parameter estima­
tions and soil moisture predictions in comparison with the traditional state augmentation method. The parti­
tioned update parameter correction method provides an assimilation tool for improving the predictions of soil
moisture by considering the time-varying parameters.

Keywords:
Data assimilation
Parameter estimation
Partitioned update
Parameter correction

1. Introduction
Influenced by wetting and drying cycles (Mubarak et al., 2009),
tillage practices (Mohanty et al., 1994; Cameira et al., 2003), and crop
root growth (Rasse et al., 2000; Iqbal et al., 2005), the hydraulic prop­
erties of the soil top layer can change during the growth period. Using
time-varying soil hydraulic parameters can improve the accuracy of soil
moisture simulations (Schwen et al., 2011). However, only a small
number of studies have focused on soil moisture modeling for timevarying parameters (Or et al., 2000; Xu and Mermoud, 2003; Schwen

et al., 2011).
In order to obtain time-varying soil hydraulic parameters, the com­mon practice is to conduct field infiltration experiments or take mea­
surements on soil samples in the laboratory at intervals during thesimulating period (Xu and Mermoud, 2003; Mubarak et al., 2009;
Schwen et al., 2011), which is labor-intensive and time-consuming.
Also, the parameters obtained in the laboratory may not be suitablefor field applications due to the scale effect. Inverse modeling provides
an alternative method and has been widely used to estimate soil hy­
draulic parameters (Pan and Wu, 1998; Abbaspour et al., 2001; Lambot* Corresponding author.
E-mail address: ghuang@cau.edu.cn (G. Huang).
https://doi.org/10.1016/j.jhydrol.2020.125729
Received 17 November 2019; Received in revised form 28 October 2020; Accepted 30 October 2020
Available online 6 November 20200022-1694/© 2020 Elsevier B.V. All rights reserved.

K. Liu et al.

Journal of Hydrology 594 (2021) 125729

Table 1
Soil hydraulic parameters for the synthetic case.
Parameter3− 3

θs (cm cm )θr (cm3 cm− 3)
KS (cm min− 1)
α (cm− 1)n (–)

1st layer

2nd layer

0.399
0.00010.02070.01741.376

0.399
0.00010.03150.01391.602

Note: adapted from Šimůnek et al. (2012).

et al., 2002; Li and Ren, 2011; Brandhorst et al., 2017).
There are mainly three kinds of inverse methods for estimating timevarying parameters. The first approach considers a functional form of
the time-varying parameter (Westra et al., 2014; Jeremiah et al., 2013).
The second approach divides observation data into several consecutive
time intervals and calibrates the model parameters for each interval
separately (Thirel et al., 2015; Merz et al., 2011; Gharari et al., 2013).
However, artificially selecting the functional form or dividing the timeintervals still requires scientific judgment, which is often difficult(Pathiraja et al., 2016). Finally, the third approach, i.e., assimilation
method, updates parameters automatically when observation data isavailable. Several studies have investigated the potential of data
assimilation methods to estimate time-varying parameters. Pathiraja
et al. (2016) proposed a multi-layer ensemble Kalman filter (EnKF)
method and a local linearization method to estimate the parameters of aconceptual hydrological model PDM (Moore, 2007). However, themulti-layer EnKF method requires estimation of the hyperparameter in
addition to the model parameters, thus increasing the number of pa­
rameters to be estimated. The local linearization method assumes that
there is a linear relationship between the parameters at the current time
and the previous time. Hence, it cannot deal with the situation when theparameter value has abrupt changes. Deng et al. (2016) and Xiong et al.
(2019) used the constrained EnKF method to estimate the parameters of
the water balance model TWBM (Xiong and Guo, 1999). The results
show that the constrained EnKF method can obtain accurate parameter
estimations, but there are time lags between the estimated value and itstrue value when the parameter has a periodical variation. Smith et al.
(2008) proposed a technique based on a particle filter to evaluate themodel’s structural inadequacy. In this method, the parameters weretreated as evolving in time. Vrugt et al. (2013) and Salamon and Feyen(2009) also tested the temporal evolution of model parameters based on
the particle filter method. However, the particle filter method iscomputationally very intensive.
Most of the existing studies estimating time-varying parameters
using the assimilation method involve the conceptual hydrologicalmodel. Since soil hydrodynamic models usually have strong nonlinear
behavior, the assimilation method’s potential to estimate time-varying
soil hydraulic parameters needs further testing. In addition, the stateaugmentation technique is usually used to estimate soil hydraulicproperties within the EnKF method (Li and Ren, 2011; Shi et al., 2015;
Bauser et al., 2016; Brandhorst et al., 2017; Wu and Margulis, 2011), in
which parameters and state variables are included in a state vector andupdated simultaneously. However, this technique’s disadvantage is thatwhen the number of unknown model states and parameters is large, thedegree of freedom for high-dimensional vectors increases, making theestimation unstable, especially in nonlinear dynamic models (Mor­adkhani et al., 2005; Xie and Zhang, 2013). Xie and Zhang (2013)
proposed a partitioned update scheme for state-parameter estimation
based on the EnKF method, in which the parameter set is partitioned
into several types according to their sensitivities, and each parameter
type is updated in an individual loop of assimilation. Lei et al. (2019)
proposed a parameter correction method based on the EnKF framework,in which the state vector consists of parameters to be estimated.
To the best of our knowledge, there exists no research on estimating(a) Computational domain(b) The initial pressure head profile and positions of observation points
Fig. 1. The computational domain for the synthetic case (a) and the initialpressure head profile and positions of observation points (b). Filled circles
represent observation points. .
Adapted from Šimůnek et al. (2012)

time-varying parameters of a soil hydrodynamic model using the data
assimilation method. The purpose of this study is to explore the potentialof the assimilation method to estimate time-varying soil hydraulic pa­
rameters. We compare four assimilation schemes: the simultaneous
update state augmentation method, the partitioned update stateaugmentation method, the simultaneous update parameter correction
method, and the partitioned update parameter correction method. The
remainder of this paper is organized as follows. In Section 2, the soilhydrodynamic model and the assimilation algorithm are introduced.
Section 3 illustrates a synthetic case and a real case involving a drip
irrigation experiment. Results and discussion are presented in Section 4.
Finally, conclusions are summarized in Section 5.
2. Methodology
2.1. Soil hydrodynamic model
The soil moisture dynamics can be described using the Richards
equation:
[ (
)]
∂θ
∂
∂h
K KijA + KizA
− S
(1)
=
∂t ∂xi∂xj
where θ is the volumetric moisture content [L3 L− 3], h is the pressure
2

K. Liu et al.

Journal of Hydrology 594 (2021) 125729

Table 2
Soil physical properties for the real case.
Soil
depths(cm)

Soil particle fraction (%)
Sand >
0.05 mmSilt
0.002–0.05mm

Clay<0.002
mm

0–20

12.18

60.98

26.84

20–40

21.84

56.87

21.29

40–60

18.46

54.81

26.73

Soil
textureBulkdensity (gcm− 3)

Silt
loam
Silt
loam
Silt
loam

1.47

Table 5
Parameters of the van Genuchten–Mualem model estimated using theLevenberg-Marquardt optimization algorithm.
Soil depths(cm)

θr (cm3
cm− 3)

θS (cm3
cm− 3)αn (–)

KS (cmday− 1)

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

(cm− 1)

1.53
1.49

where KS is the saturated hydraulic conductivity [L T− 1], Se is theeffective water saturation [–], θr and θs are the residual and saturated
water contents, respectively [L3 L− 3], α is the inverse of the air-entry
value [L− 1], and n is the pore-size distribution index [–]. In this study,
unsaturated water flow was simulated using the CHAIN_2D code, a DOS
predecessor of the HYDRUS model (Šimůnek et al., 2008).

Table 3Growth stages of processing tomato.
Growth stage

InitialDevelopment

Mid-season

Late season

Day

1–26

27–43

44–83

84–1042.2. Parameter estimation based on the assimilation method
The assimilation method for parameter estimation consists of amodel operator, an observation operator, and an ensemble Kalman filter
algorithm. The model operator is used to pass the state vector to the next
moment, which can be expressed as follows:(5)

Xi+1 = Mi Xi + ηi+1

where Xi and Xi+1 are state vectors at current and next time, respec­tively; i and i + 1 denote time indicators; Mi is the model operator; ηi+1 isthe model error vector that is independent white noise for the model
operator, which is drawn from a multinormal distribution with zeromean and specified covariance Qi+1.
The observation operator constructs the mapping between the statevector and the observation vector, which can be written as:(6)

Yi+1 = Hi+1 Xi+1 + εi+1

in which Yi+1 is the observation vector; Hi+1 is the observation operator;
εi+1 is the observation error vector, which is also assumed to be inde­
pendent while noise is drawn from a multinormal distribution with zeromean and specified covariance Oi+1.
The EnKF algorithm is a Monte Carlo method, in which the initialensemble member can be generated by adding a random disturbance tothe initial state vector, which can be expressed as follows:

Fig. 2. The computational domain, boundary conditions, and positions ofobservation points for the real case.

head [L], K is the unsaturated hydraulic conductivity [L T− 1], KAij areAcomponents of a dimensionless anisotropy tensor K , xi and xj are thespatial coordinates [L], and S is the sink term accounting for root water
uptake.
The unsaturated hydraulic conductivity function is determined using
the van Genuchten-Mualem model (Mualem, 1976; van Genuchten,1980):[()m ]2
1
K(h) = Ks Se0.5 1 − 1 − Sem
(2)

Se =

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

where j indicates the jth member in the state vector; ξj is independent
white noise, which follows a multi-normal distribution with zero meanand specified covariance R.
In the implementation of EnKF, two steps, namely forecast andanalysis steps, are included (Evensen, 2003). First, the forecast statevector at the i + 1 moment is generated using the model operator and theanalysis state vector at the i moment, which can be written as:(8)

Xi+1,f= Mi Xi,ajjin which X i+1, f is the forecast state vector, and Xi+1, a is the analysis
state vector.
When the observation data is available, the analysis state vector is
updated using the analysis equation, which can be written as:

(3)

1
n

(7)

X0j = X0 + ξj

(4)

Table 4
Drip flux during the growth period.
Day− 1

Drip flux (cm day

)

1

5

37

43

52

61

68

77

88

22.00

5.12

5.57

3.26

6.08

5.28

4.22

5.98

5.34

3

K. Liu et al.

Journal of Hydrology 594 (2021) 125729(a) KS has a periodic variation(b) KS has a linear variation(c) KS has an abrupt change
Fig. 3. Parameter evolutions under the condition of a single time-varying parameter [(a) KS has a periodic variation, (b) KS has a linear variation, (c) KS has anabrupt change].

Xi+1,a= Xi+1,f+ Ki+1 (Yi+1 + εi+1 − Hi+1 Xi+1,f
)
j
j
j

(9)

where β is the damping factor, with values between 0 and 1. In thisstudy, β is considered with a value of 0.1.
In this study, two assimilation strategies are considered. The maindifference between the two strategies is the selection of the model
operator and the observation operator. In strategy 1 (S1), the soil hy­
drodynamic model is used as the model operator. Because the obser­
vation data includes soil moisture or pressure head, the observation
operator is a unit matrix I. The state variables and model parameters are
updated using the state augmentation technique. In strategy 2 (S2), thestate vector consists of soil hydraulic parameters, and the model oper­
ator is a unit matrix I. Since parameters are treated as state variables, S2
is called the parameter correction method. The soil hydrodynamic
model provides a direct link between parameters and observation ele­ments (soil moisture or pressure head), and hence the Richards equationis the nonlinear observation operator. In each assimilation strategy, twoparameter update methods are further analyzed. In method 1 (M1),parameters are updated simultaneously. In method 2 (M2), parameters

where Ki+1 is the Kalman gain, which is defined as:
TT

Ki+1 = Pi+1 Hi+1 (Hi+1 Pi+1 Hi+1 + Oi+1 )− 1

(10)

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

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

in which X

i+1,f

i+1,f

is the mean of Xj(11). The damping factor (Hendricks

Franssen and Kinzelbach, 2008) was used to reduce the inbreedingproblem, hence the Eq. (9) is revised as:
Xi+1,a= Xi+1,f+ βKi+1 (Yi+1 + εi+1 − Hi+1 Xi+1,f
)
j
j
j

(12)

4

K. Liu et al.

Journal of Hydrology 594 (2021) 125729

are updated separately, i.e., the partitioned update method (Xie andZhang, 2013). Therefore, a total of four assimilation schemes, namely
the simultaneous update state augmentation method (S1M1), the par­titioned update state augmentation method (S1M2), the simultaneous
update parameter correction method (S2M1), and the partitioned up­
date parameter correction method (S2M2) are compared in this study.

Table 6RMSE of the estimated parameters for four assimilation schemes.

S1M1
S1M2
S2M1
S2M2

Periodic variation

Linear variation

Abrupt change

KS (cmmin− 1)α(cm− 1)

KS (cmmin− 1)α(cm− 1)

KS (cmmin− 1)α0.00250.00260.00260.00260.00070.00030.00040.00020.00210.00170.00220.00190.00050.00020.00040.00020.00330.00310.00360.00300.00090.00040.00070.0003(cm− 1)

3. Case study
3.1. Synthetic caseIn the synthetic case, water infiltration from a single-ring infil­
trometer into a layered soil was simulated. The soil profile consists of
two layers: 40 cm thick A-horizon and an underlying B/C-horizon. The
hydraulic parameters of two soil layers are shown in Table 1 (Císlerová,1987; Hopmans and Stricker, 1989). Fig. 1 presents the axisymmetriccomputational domain. The computational domain was discretized into
684 finite elements. All sides of the computational domain are imper­
vious except for a small portion around the origin at the soil surface
where a constant pressure head was imposed and the lower right corner
where the groundwater level was kept constant. The initial pressure

Table 7
Correlation coefficient of the estimated parameters for a single time-varying
parameter.
S1M1
S1M2
S2M1
S2M2

Periodic variation

Linear variation

Abrupt change

0.99
0.61
0.99
0.45

0.99
0.54
0.99
0.34

0.99
0.65
0.99
0.50

Fig. 4. Absolute errors for simulated pressure head distributions (cm) for four assimilation schemes [(a) KS has a periodic variation, (b) KS has a linear variation, (c)KS has an abrupt change].
5

K. Liu et al.

Journal of Hydrology 594 (2021) 125729(a) Scenario 1(b) Scenario 2
Fig. 5. Parameter evolutions for different parameter variation scenarios [(a) Scenario 1, (b) Scenario 2].

obtained by adding an observation error to the pressure heads simulatedusing the true parameters. The standard deviation of the observation
error is 0.01 m. The parameters to be estimated are KS and α of the firstsoil layer because they are the most variable parameters (Carsel andParrish, 1988). Estimated parameters were log-transformed because
their probability distributions are known to approximate lognormal
distributions (Carsel and Parrish, 1988). The standard deviation of theinitial parameter ensemble is 0.1 for both log10(KS) and log10(α).
To test the ability of assimilation methods on identifying the dy­
namics of parameters, multiple scenarios, in which parameters havedifferent time variation trends, were designed. First, only a single timevarying parameter was considered, i.e., KS changed with time, and α wasconstant. Different variation trends of KS were developed, i.e., a periodic
variation, a linear variation, and an abrupt change. Similar assumptionsof parameter variations were also used by Deng et al. (2016) and Xiong
et al. (2019). Second, to check the parameter correlation issue that may
arise during the assimilation, both KS and α were considered to changewith time. Two scenarios were designed in this case: (1) KS with a pe­
riodic variation and α with a linear trend. The effect of the initial
parameter value and the order of parameter updates on the assimilation
results was also analyzed. In this case, (2) KS has a periodic variation,and α has an abrupt change. The results of four assimilation schemes
were compared with each other.

Table 8Correlation coefficient of the estimated parameters for multiple time-varying
parameters.
S1M1
S1M2
S2M1
S2M2

Scenario 1

Scenario 2

0.99
0.89
0.99
0.84

0.99
0.87
0.99
0.36

Table 9RMSE of the estimated parameters for different scenarios.
Scenario 1

Scenario 2

KS (cm min− 1)

α (cm− 1)KS (cm min− 1)

α (cm− 1)

S1M1

0.00440.00220.00340.0017

S1M2
S2M1
S2M2

0.00400.00410.00320.00160.00230.00180.00370.00290.00240.00160.00230.0016head profile and locations of observation points are shown in Fig. 1b.
The ensemble size was set to 300, which has been tested to be large
enough to meet the accuracy requirements, i.e., the relative errors of
estimated parameters are less than 5%. Shi et al. (2015) reported thatthe type of observation data might affect assimilation results. To better
test the ability of the assimilation methods to identify time-varying
parameters, pressure heads were chosen as the observation data inthis case. The total simulation time was 10 days, and pressure heads at
observation points were recorded every 2 h. The measurements were

3.2. Real case
A field drip irrigation experiment with processing tomatoes was
conducted at the upper reaches of the Yellow River basin, China (lati­
tude 40◦ 41′ N, longitude 107◦ 18′ E, and 1041 m altitude). The soil profile
60 cm deep can be divided into three layers. The physical properties of
the three soil layers are given in Table 2. The growth period started from
6/2/2015 and lasted till 9/1/2015, i.e., 104 days in total. According to6

K. Liu et al.

Journal of Hydrology 594 (2021) 125729

Fig. 6. Absolute errors of simulated pressure head distributions (cm) for Scenario 2.

FAO-56 (Allen et al., 1998), the growth period of processing tomatoes
can be divided into four growth stages (see Table 3). Meteorological data
was measured by an automatic weather station (HOBO, Campbell Sci­
entific Inc., USA). Daily reference evapotranspiration (ETo) was calcu­
lated using the modified Penman-Monteith equation (Allen et al., 1998).
Soil water contents at two horizontal locations and three depths (Fig. 2)
were monitored using HydraProbe sensors. The computational domain,boundary conditions, and the position of observation points are shownin Fig. 2. The computational domain was discretized into 2527 finite
elements. A time-variable flux boundary condition, determined from
irrigation rates, was specified on the BC boundary (Fig. 2). The drip flux
during the growth period is given in Table 4. A no flux boundary con­dition was specified on the CD side because of a plastic mulch cover, andan atmospheric boundary condition was specified on the DE side. The
lower boundary was treated as a free drainage boundary because thegroundwater table was relatively deep. Due to symmetry of flow, a noflux boundary condition was specified on both vertical sides. Water
uptake was simulated using the model of Feddes et al. (1978) withdefault values of parameters of the stress response function provided by
HYDRUS-2D (Šimůnek et al., 2012). The initial water content was 0.4
cm3 cm− 3, 0.37 cm3 cm− 3, and 0.45 cm3 cm− 3 in depths of 0–20 cm,
20–40 cm, and 40–60 cm, respectively. An ensemble size of 300 and a 1day interval were used in the assimilation. A value of 0.01 was used as
the standard deviation for the soil water content observation error (Shiet al., 2015). According to Haverkamp et al. (1996), the parameters θS,
KS, and α depend on the soil structure, and hence these three parametersof the soil top layer were estimated. Initial values of soil hydraulic pa­
rameters of the three soil layers (shown in Table 5) were obtained using

the Levenberg-Marquardt optimization algorithm (Moré, 1978). The
standard deviation of the initial parameter ensemble was 0.1 forlog10(KS) and log10(α), and 0.02 for θS. The parameter estimations andwater content simulations obtained using the partitioned updateparameter correction method are compared with those of the traditional
state augmentation method.
3.3. Evaluation index
The root mean square error (RMSE) is used to evaluate parameter
estimations for the synthetic case and simulated soil water contents forthe real case. The absolute error (AE) is used to evaluate the simulatedpressure head distribution.
√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅1 ∑N
RMSE =(13)(Si − Oi )2
i=1
N
AE = Si − Oi

(14)where N is the total number of observations, Si and Oi are the ith modelsimulated and observed values, respectively.
The Pearson correlation coefficient (R) is used to analyze the linear
relationship between estimated parameters.
∑N
i=1 (xi − x)(yi − y)̅√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
̅R = √̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
(15)
∑N
∑N
2
2
i=1 (xi − x)i=1 (yi − y)where x and y denote different model parameters; and x and y are themean values of these parameters.
7

K. Liu et al.

Journal of Hydrology 594 (2021) 125729(a) Equal to the true value(b) Larger than the true value(c) Less than the true value
Fig. 7. Parameter evolutions for different initial parameter values [(a) Equal to the true value, (b) Larger than the true value, (c) Less than the true value].

compared with the true value variation. This is consistent with the re­sults of Clark et al. (2008), Samuel et al. (2014), Deng et al. (2016), andXiong et al. (2019). This is due to the defect in the EnKF method. The
states and parameters are updated based on current observations andprevious predictions, and hence there is a time lag between the assim­
ilation value estimated using the EnKF method and the true value,
especially when a peak value occurs (Clark et al. 2008). The estimatedKS peak value is lower compared to its true value, especially for thesimultaneous update method (M1). This may be attributed to the mutual
influence between estimated parameters. As can be seen from Fig. 3(a),the parameter α estimated using the simultaneous update method (M1)
has a large deviation from its true value, which in turn affects the esti­
mation of KS. The parameter α estimated using the partitioned update
method (M2) can eventually converge to its true value.
When KS has a linear variation, values of KS estimated by the four
assimilation schemes are lower than the true value. This can also be

Table 10RMSE of the estimated parameters for different initial parameter values.
Larger

S1M1
S1M2
S2M1
S2M2

Equal

Less

KS (cmmin− 1)α(cm− 1)

KS (cmmin− 1)α(cm− 1)

KS (cmmin− 1)α0.00440.004
0.00410.00320.00220.00160.00230.00180.00430.00340.00420.00320.00230.00140.00230.00160.00450.00440.00440.00430.00220.00220.00230.0023(cm− 1)

4. Result and discussion
4.1. Single time-varying parameter
Fig. 3 presents the evolutions of the estimated parameters. When KS
has a periodic variation, estimated KS values show significant time-lags
8

K. Liu et al.

Journal of Hydrology 594 (2021) 125729

Fig. 8. Absolute errors of simulated pressure head distributions (cm) for different initial parameter values [(a) Larger than the true value, (b) Equal to the true value,(c) Less than the true value].

Fig. 9. Parameter evolutions for different parameter update orders.

regarded as a time lag for the parameter estimation. Values of parameter
α estimated using the simultaneous update method (M1) show large
deviations from its true value. On the other hand, the parameter αestimated using the partitioned update method (M2) can converge to itstrue value.
When KS has an abrupt change, all four assimilation schemes could

respond to the abrupt change of the parameter value immediately.
However, multiple assimilation operations (about 50 assimilations inthis case) are required for the parameter to converge to its true value.
The parameter α was overestimated using the simultaneous update
method (M1). The partitioned update method (M2) could obtain accu­
rate α estimation.
9K. Liu et al.

Journal of Hydrology 594 (2021) 125729

more accurate than those using parameters obtained by the simulta­
neous update method because of better parameter estimations. This isconsistent with the results of Xie and Zhang (2013). For the case with aperiodically varying KS, due to the time lags of parameter estimations,the accuracy of simulated pressure head distributions is similar for four
assimilation schemes.

Table 11
RMSE of the estimated parameters for different parameter update orders.
S1 – KS α
S1 – α KS
S2 – KS α
S2 – α KS

KS (cm min− 1)

α (cm− 1)

0.004
0.00280.00320.00320.00160.00180.00180.00194.2. Multiple time-varying parameters

Table 6 gives the RMSE values for the estimated parameter for four
assimilation schemes. It can be found that the partitioned update
method (M2) obtains more accurate results except for the case of aperiodically varying KS. Low RMSE values for the case of a periodic
varying KS may be attributed to the time lags of parameter estimation.
Standard deviations of the parameters are similar for four assimilation
schemes (see Fig. 3). Standard deviations have dropped to a very low
level after 50 assimilations.
As can be seen from Fig. 3, there is an obvious parameter correlation
issue, i.e., the evolutions of α display similar trends as KS when the pa­
rameters are updated simultaneously, which results in low estimation
accuracy of α. Compared with the simultaneous update method,parameter correlation coefficients decrease significantly for the parti­
tioned update method, especially when the parameter correction
method is used (see Table 7). This indicates that the partitioned update
method can effectively resolve parameter correlation issues, thereby
improving the accuracy of parameter estimations.
Simulated pressure head error distributions for four assimilation
schemes are shown in Fig. 4. Simulation errors in the 30–70 cm soil
depth are large, which is caused by the parameter estimation error.
Except for a periodically varying KS, simulated pressure head distribu­
tions using parameters obtained by the partitioned update method are

Fig. 5 presents the parameter evolutions for different parametervariation scenarios, i.e., for Scenario 1 with a periodic variation of KS
and a linear trend in α, and Scenario 2 with a periodic variation of KS andan abrupt change in α. It can be found that there are obvious parametercorrelations in the parameter estimations using the simultaneous update
method. The estimated KS has the same increasing trend as the αparameter in Scenario 1. The evolution of KS showed a significant in­flection point with an abrupt change of α in Scenario 2. The correlation
coefficients of the estimated parameters are given in Table 8. The R
values for the partitioned update method are smaller, resulting in ahigher accuracy of parameter estimation (see Table 9). Deng et al.
(2016) and Xiong et al. (2019) successfully estimated time-varying pa­
rameters of a two-parameter water balance model (Xiong and Guo,1999) using the simultaneous update state augmentation method.
However, four assimilation schemes cannot obtain the ideal estimationof time-varying soil hydraulic parameters in our study. This may be
attributed to the fact that the Richards equation has stronger nonlinear
properties compared to a two-parameter water balance model. Since
simulated pressure head error distributions for Scenarios 1 and 2 aresimilar, Fig. 6 only presents the results for Scenario 2. The accuracy ofpressure head simulations using parameters obtained by the partitioned

Fig. 10. Absolute errors of simulated pressure head distributions (cm) for different parameter updating orders.
10K. Liu et al.

Journal of Hydrology 594 (2021) 125729

Fig. 11. Parameter evolutions for the drip irrigation experiment.

update parameter correction method is the highest among the four
schemes because of good parameter estimations.

because the accuracy of simulated pressure heads is improved when thefirst parameter is updated, which is beneficial for updating the second
parameter. On the other hand, the parameter update order has only asmall influence on the estimations for the partitioned update parametercorrection method (Table 11), because the parameters and state vari­ables are not updated simultaneously in the parameter correction
method. Simulated pressure head error distributions are shown in
Fig. 10. For the partitioned update state augmentation method, updating
α first can improve the simulation accuracy, which indicates that in thisexample, the simulation is more sensitive to the parameter KS. Hence theparameter update order needs to be carefully selected when the parti­
tioned update state augmentation method is used. Xie and Zhang (2013)
recommended that the parameter updating order should be performedin the order of increasing sensitivity. In contrast, the partitioned updateparameter correction method can obtain good results of parameter es­
timations and simulations regardless of the parameter update order.

4.3. Effect of the initial parameter value
Since there is uncertainty in the selection of initial values of esti­
mated parameters in practical applications, the effects of the initial
parameter value on assimilation are analyzed next. Fig. 7 presents theparameter evolutions for different initial parameter values. It can befound that the time lags of the estimated parameter cannot be avoided
regardless of the initial parameter value. When the parameters’ initialvalues are greater than their true values, the trends in parameter esti­
mation are similar to when the initial parameter values are equal to theirtrue values. However, when the parameters’ initial values are smaller
than their true values, the parameter estimation accuracy is significantly
reduced (Table 10), which leads to inaccurate simulations of pressure
head distributions (see Fig. 8). This can be attributed to the fact that
smaller saturated hydraulic conductivity results in slower reactions toforcing. If soil water dynamics is slow, parameters cannot be estimatedwell. And if the van Genuchten α is low, the pressure entry head is high,and the soil is wet, and forcing will not have a strong effect. This makes
parameter estimation more difficult. The results indicate that when theinitial value of the optimized parameter deviates from the true value to acertain extent, even the partitioned update method cannot obtain ac­
curate parameter estimations and soil moisture simulations.

4.5. Application to the drip irrigation experimentFig. 11 presents the parameter evolutions for the drip irrigation
experiment when the partitioned update parameter correction method
and the traditional state augmentation method are used. For compari­son, constant parameter values obtained using the Marquardt-Levenberg
optimization method are also shown. Fig. 11 shows that θS did notchange much during the growth period, and estimated θS values using
the partitioned update parameter correction method were similar to
those obtained using the traditional state augmentation method. There
are also significant differences in the parameter evolutions of KS and α
for the two methods. The parameters estimated using the traditional
state augmentation method showed larger fluctuations. The parametervalues changed significantly after each rainfall or irrigation. This isbecause the parameters and soil water contents are updated simulta­
neously in the state augmentation technique. When the soil water con­tent changed dramatically because of rainfall or irrigation, the model
parameters also changed significantly. The parameters estimated using

4.4. Effect of the parameter update order
For the partitioned update method, the parameter update order may
influence the parameter estimations (Xie and Zhang, 2013). Therefore,
the effect of the parameter update order is further investigated. Fig. 9
presents the parameter evolutions for different parameter update orders.
For the partitioned update state augmentation method, updating KS first
can result in better estimation accuracy of α. Conversely, updating α first
can improve the estimation accuracy of KS (Table 11). This may be
11

K. Liu et al.

Journal of Hydrology 594 (2021) 125729

out that crop rooting systems tend to increase water flow, resulting in ahigher KS.
The soil water content predicted using parameters obtained by thetwo assimilation methods and the Levenberg-Marquardt optimization
algorithm are shown in Fig. 12. Soil water contents for the 0–20 cm
depth predicted using parameters obtained by the two assimilation
methods and the Levenberg-Marquardt optimization algorithm weresimilar. This may be attributed to the fact that variations in the watercontent of the topsoil layer were mainly controlled by changing
boundary conditions. However, for deeper soils (20–60 cm), the accu­racy of the predictions of the partitioned update parameter correction
method and the Levenberg-Marquardt optimization algorithm washigher than that of the traditional state augmentation method (seeTable 12). After day 50, the soil water content predicted using the
traditional state augmentation method showed a significant deviationfrom observations, which may result from inaccurate estimations of theparameters. Compared with the Levenberg-Marquardt optimization al­
gorithm, the partitioned update parameter correction method could
obtain more accurate soil water content simulations, which indicatesthat using time-varying soil hydraulic parameters during the growth
period can improve soil moisture prediction. This is consistent with theresults of Mubarak et al. (2009).
To prove that the parameter variations estimated using the assimi­
lation methods were caused by the soil structure change rather than
assimilation errors, the assimilation methods were further used to esti­
mate constant parameters. The observations were generated by adding
an observation error to the soil water content simulated using the pa­
rameters estimated by the Levenberg-Marquardt optimization algo­rithm. The parameter evolutions of two assimilation methods are shownin Fig. 13. The parameter values estimated using the partitioned updateparameter correction method were very close to the values estimated by
the Levenberg-Marquardt optimization algorithm. There were no
obvious value changes during the growth period. This indicates that theparameter variations were caused by the soil structure change rather
than the assimilation errors in the previous simulations. In contrast, theparameters estimated using the traditional state augmentation method
deviated from their true values, which could be attributed to unstableassimilation resulting from an increasing degree of freedom in the sys­tem (Moradkhani et al., 2005).
Due to a large number of state variables in the simulations of multidimensional soil water flow, the size of the matrix to be updated in the
traditional state augmentation method is larger than in the partitioned
update parameter correction method, leading to much more computa­
tional time. The state augmentation method is prone to an inconsistency
problem during assimilations, resulting in non-convergence of the iter­
ative solution of the Richards equation, which makes the assimilation
unstable. This is because the unsaturated water flow is a highly
nonlinear problem, and the assimilation of state variables and parame­ters using the linear equation may not follow the nonlinearity of the
Richards equation (Song et al., 2014). Since only parameters are
assimilated in the partitioned update parameter correction method, theinconsistency problem may be avoided. Therefore, in multi-dimensional
unsaturated water flow problems, using the partitioned update param­
eter correction method can improve the calculation efficiency and theassimilation stability while obtaining more accurate parameter estima­
tions and soil moisture predictions.
It should be noted that the sequential data assimilation requires
consecutive data to achieve statistical convergence of parameter esti­
mation. If the observation data is not enough, the parameter may notconverge to its true value, especially for parameters with specific (suchas linear) trend. Hence, a small assimilation time interval or a long timeseries of data are required to improve the convergence of parameter
estimates for short-period simulation. For example, tillage practice may
result in the temporal variation of soil parameters. In this study, tillage
practice’s temporal variation of soil parameters was not consideredsince the tillage occurred before crop transplanting, and the soil(a) 0 - 20 cm(b) 20 - 40 cm(c) 40 - 60 cm
Fig. 12. Comparison of predicted and measured soil water contents [(a) 0–20
cm, (b) 20–40 cm, (c) 40–60 cm].

Table 12RMSE of the soil water contents predicted using parameters obtained by thepartitioned update parameter correction method, the traditional state augmen­
tation method and the Levenberg-Marquardt (L-M) optimization algorithm.
Soil depth(cm)

Partitioned (cm3
cm− 3)

Augment (cm3
cm− 3)

L-M (cm3
cm− 3)

0–20
20–40
40–60

0.039
0.007
0.019

0.038
0.015
0.028

0.039
0.010
0.021

the partitioned update parameter correction method varied aroundvalues estimated by the Levenberg-Marquardt optimization algorithm.
The parameters showed significant trend changes on days 5, 23, 27, and
68. The values of KS and α decreased after days 5, 23, and 68, which may
be attributed to soil compaction due to rainfall or irrigation. The
parameter values showed an increasing trend after day 27, likely
because of rapid cop growth (see Table 3). The KS value may increase
due to the root growth, which would be consistent with the results ofRasse et al. (2000) and Mubarak et al. (2009). Rasse et al. (2000) found12

K. Liu et al.

Journal of Hydrology 594 (2021) 125729

Fig. 13. Parameter evolutions of the partitioned update parameter correction method and the traditional state augmentation method under the condition of con­stant parameters.

simulations were first tested using the synthetic case. The influence of
the initial parameter value and the parameter update order on assimi­
lation was also analyzed. The partitioned update parameter correction
method was then applied to a real case involving a field drip irrigation
experiment. The results showed that when the parameter had a peri­
odical or a linear variation, there were time lags between the assimila­tion value and the true value. The assimilation method could respond toan abrupt change of the parameter value immediately. There were
obvious parameter correlation issues when the simultaneous update
method was used. The partitioned update method could effectively
resolve the parameter correlation issue, thereby improving the accuracyof parameter estimations and pressure head simulations. However,when the initial values of the parameters deviated from their true values
to a certain extent, the partitioned update method could not obtain ac­
curate parameter estimations. The partitioned update parametercorrection method was not sensitive to the parameter update order.
Since only parameters are updated in the assimilations, the inconsis­
tency of the assimilation could be avoided. Compared with the tradi­tional state augmentation method, the partitioned update parametercorrection method has high computational efficiency and assimilation
stability, it provides an efficient tool for the assimilation of multidimensional unsaturated water flow problems by considering the timevarying parameters.

structure had stabilized when crop growth began. However, the pro­
posed assimilation method can also be applied to soil disturbed bytillage practice if sufficient observation data can be obtained.
In the real test case, on the one hand, the change of boundary con­ditions, e.g., time-varying irrigation flux and evapotranspiration, androot growth might cause the change of soil structures, leading to timevarying parameters. Compared with the traditional methods, the parti­
tioned update parameter correction method can improve soil moisturepredictions. On the other hand, soils are heterogeneous, which cannot
be well captured by available models. This will result in the predictions
deviating from the observations. The partitioned update parametercorrection method provides an efficient tool for assimilating multidimensional soil water flow by considering the time-varying parame­ters representing the influence of all relevant factors, e.g., soil hetero­geneity, change of boundary conditions, root growth and model errors,
etc. However, it should be mentioned that the temporal change of pa­
rameters is inherent in the assimilation method and that there can beother reasons why parameters change during the updates. Most of thesereasons are related to experimental and modeling errors. For example,
the measurement errors, simplifications of initial and boundary condi­
tions, inadequacy of governing transport equations to describe observed
processes fully (e.g., the presence of preferential flow during rainfalls),
etc., can also result in assimilation methods predicting temporal varia­
tions of optimized parameters.

Declaration of Competing Interest

5. Summary and conclusions

The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to influence
the work reported in this paper.

In this study, we explored the potential of the assimilation method toestimate time-varying soil hydraulic parameters. Four assimilation
schemes were compared, including the simultaneous update stateaugmentation method, the partitioned update state augmentation
method, the simultaneous update parameter correction method, and thepartitioned update parameter correction method. The four assimilation
schemes’ performance on parameter estimation and soil moisture

Acknowledgments‘This research was jointly supported by the National Key R&D Pro­
gram of China (No. 2017YFC0403301) and the National Natural Science13

K. Liu et al.

Journal of Hydrology 594 (2021) 125729

Foundation of China (Nos. 51639009, 51621061, 51125036).

Moradkhani, H., Sorooshian, S., Gupta, H.V., Houser, P.R., 2005. Dual state–parameter
estimation of hydrological models using ensemble Kalman filter. Adv. Water Resour.
28 (2), 135–147.
Moore, R.J., 2007. The PDM rainfall-runoff model. Hydrol. Earth Syst. Sci. Discuss. 11(1), 483–499.
Moré, 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.
Mubarak, I., Mailhol, J.C., Angulo-Jaramillo, R., Ruelle, P., Boivin, P., Khaledian, M.,
2009. Temporal variability in soil hydraulic properties under drip irrigation.
Geoderma 150 (1–2), 158–165.
Or, D., Leij, F.J., Snyder, V., Ghezzehei, T.A., 2000. Stochastic model for posttillage soil
pore space evolution. Water Resour. Res. 36 (7), 1641–1652.
Pan, L., Wu, L., 1998. A hybrid global optimization method for inverse estimation of
hydraulic parameters: annealing-simplex method. Water Resour. Res. 34 (9),2261–2269.
Pathiraja, S., Marshall, L., Sharma, A., Moradkhani, H., 2016. Hydrologic modeling in
dynamic catchments: a data assimilation approach. Water Resour. Res. 52 (5),3350–3372.
Rasse, D.P., Smucker, A.J., Santos, D., 2000. Alfalfa root and shoot mulching effects on
soil hydraulic properties and aggregation. Soil Sci. Soc. Am. J. 64 (2), 725–731.
Salamon, P., Feyen, L., 2009. Assessing parameter, precipitation, and predictive
uncertainty in a distributed hydrological model using sequential data assimilation
with the particle filter. J. Hydrol. 376 (3–4), 428–442.
Samuel, J., Coulibaly, P., Dumedah, G., Moradkhani, H., 2014. Assessing model state andforecasts variation in hydrologic data assimilation. J. Hydrol. 513, 127–141.
Schwen, A., Bodner, G., Loiskandl, W., 2011. Time-variable soil hydraulic properties innear-surface soil water simulations for different tillage methods. Agric. Water
Manag. 99 (1), 42–50.
Shi, L., Song, X., Tong, J., Zhu, Y., Zhang, Q., 2015. Impacts of different types ofmeasurements on estimating unsaturated flow parameters. J. Hydrol. 524, 549–561.
Šimůnek, J., van Genuchten, M.T., Šejna, M., 2008. Development and applications of the
HYDRUS and STANMOD software packages and related codes. Vadose Zone J. 7 (2),587–600.
Šimůnek, J., Van Genuchten, M.T., Šejna, M., 2012. The HYDRUS software package forsimulating the two-and three-dimensional movement of water, heat, and multiple
solutes in variably-saturated porous media. Technical manual, version, 2, 258 p.
Smith, P.J., Beven, K.J., Tawn, J.A., 2008. Detection of structural inadequacy in processbased hydrological models: A particle-filtering approach. Water Resour. Res. 44,
W01410. https://doi.org/10.1029/2006WR005205.
Song, X., Shi, L., Ye, M., Yang, J., Navon, I.M., 2014. Numerical comparison of iterative
ensemble Kalman filters for unsaturated flow inverse modeling. Vadose Zone J. 13(2) https://doi.org/10.2136/vzj2013.05.0083.
Thirel, G., Andréassian, V., Perrin, C., Audouy, J.N., Berthet, L., Edwards, P., Folton, N.,
Furusho, C., Kuentz, A., Lerat, J., Lindström, G., 2015. Hydrology under change: anevaluation protocol to investigate how hydrological models deal with changing
catchments. Hydrol. Sci. J. 60 (7–8), 184–1199.
van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulicconductivity of unsaturated soils 1. Soil Sci. Soc. Am. J. 44 (5), 892–898.
Vrugt, J.A., ter Braak, C.J., Diks, C.G., Schoups, G., 2013. Hydrologic data assimilation
using particle Markov chain Monte Carlo simulation: theory, concepts andapplications. Adv. Water Resour. 51, 457–478.
Westra, S., Thyer, M., Leonard, M., Kavetski, D., Lambert, M., 2014. A strategy for
diagnosing and interpreting hydrological model nonstationarity. Water Resour. Res.
50 (6), 5090–5113.
Wu, C.C., Margulis, S.A., 2011. Feasibility of real-time soil state and flux characterization
for wastewater reuse using an embedded sensor network data assimilation approach.
J. Hydrol. 399 (3–4), 313–325.
Xie, X., Zhang, D., 2013. A partitioned update scheme for state-parameter estimation of
distributed hydrologic models based on the ensemble Kalman filter. Water Resour.
Res. 49 (11), 7350–7365.
Xiong, L., Guo, S., 1999. A two-parameter monthly water balance model and its
application. J. Hydrol. 216 (1–2), 111–123.
Xiong, M., Liu, P., Cheng, L., Deng, C., Gui, Z., Zhang, X., Liu, Y., 2019. Identifying timevarying hydrological model parameters to improve simulation efficiency by theensemble Kalman filter: a joint assimilation of streamflow and actual
evapotranspiration. J. Hydrol. 568, 758–768.
Xu, D., Mermoud, A., 2003. Modeling the soil water balance based on time-dependent
hydraulic conductivity under different tillage practices. Agric. Water Manag. 63 (2),
139–151.

References
Abbaspour, K.C., Schulin, R., Van Genuchten, M.T., 2001. Estimating unsaturated soil
hydraulic parameters using ant colony optimization. Adv. Water Resour. 24 (8),827–841.
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.
Bauser, H.H., Jaumann, S., Berg, D., Roth, K., 2016. EnKF with closed-eye
period–towards a consistent aggregation of information in soil hydrology. Hydrol.
Earth Syst. Sci. 20 (12), 4999–5014.
Brandhorst, N., Erdal, D., Neuweiler, I., 2017. Soil moisture prediction with the ensemble
Kalman filter: Handling uncertainty of soil hydraulic parameters. Adv. Water Resour.
110, 360–370.
Cameira, M.R., Fernando, R.M., Pereira, L.S., 2003. Soil macropore dynamics affected bytillage and irrigation for a silty loam alluvial soil in southern Portugal. Soil Tillage
Res. 70 (2), 131–140.
Carsel, R.F., Parrish, R.S., 1988. Developing joint probability distributions of soil water
retention characteristics. Water Resour. Res. 24 (5), 755–769.
Císlerová, M., 1987. Comparison of simulated water balance for ordinary and scaled soil
hydraulic characteristics (No. 82). Landbouwuniversiteit Wageningen.
Clark, M.P., Rupp, D.E., Woods, R.A., Zheng, X., Ibbitt, R.P., Slater, A.G., Schmidt, J.,
Uddstrom, M.J., 2008. Hydrological data assimilation with the ensemble Kalmanfilter: use of streamflow observations to update states in a distributed hydrologicalmodel. Adv. Water Resour. 31 (10), 1309–1324.
Deng, C., Liu, P., Guo, S., Li, Z., Wang, D., 2016. Identification of hydrological model
parameter variation using ensemble Kalman filter. Hydrol. Earth Syst. Sci. 20 (12),4949–4961.
Evensen, G., 2003. The ensemble Kalman filter: theoretical formulation and practicalimplementation. Ocean Dyn. 53 (4), 343–367.
Feddes, R.A., Kowalik, P.J., Zaradny, H., 1978. In: Simulation of Field Water Use andCrop Yield. Centre for Agricultural Publishing and Documentation, Wageningen,p. 189.
Gharari, S., Hrachowitz, M., Fenicia, F., Savenije, H.H.G., 2013. An approach to identify
time consistent model parameters: sub-period calibration. Hydrol. Earth Syst. Sci. 17(1), 149–161.
Haverkamp, R., Arrúe, J.L., Vandervaere, J.P., Braud, I., Boulet, G., Laurent, J.P.,
Taha, A., Ross, P.J., Angulo-Jaramillo, R., 1996. Hydrological and Thermal
Behaviour of the Vadose Zone in the Area of Barrax and Tomelloso (Spain):Experimental Study, Analysis and Modeling. LTHE, France.
Hendricks Franssen, H.J., Kinzelbach, W., 2008. Real-time groundwater flow modelingwith the ensemble Kalman filter: joint estimation of states and parameters and thefilter inbreeding problem. Water Resour. Res. 44 (9), 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.
Iqbal, J., Thomasson, J.A., Jenkins, J.N., Owens, P.R., Whisler, F.D., 2005. Spatial
variability analysis of soil physical properties of alluvial soils. Soil Sci. Soc. Am. J. 69(4), 1338–1350.
Jeremiah, E., Marshall, L., Sisson, S.A., Sharma, A., 2013. Specifying a hierarchicalmixture of experts for hydrologic modeling: Gating function variable selection.
Water Resour. Res. 49 (5), 2926–2939.
Lambot, S., Javaux, M., Hupet, F. and Vanclooster, M., 2002. A global multilevel
coordinate search procedure for estimating the unsaturated soil hydraulic properties.
Water Resour. Res. 38(11), 1224. doi: 10.1029/2001WR001224, 2002.
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 filter. J. Hydrol. 568, 108–118.
Li, C., Ren, L., 2011. Estimation of unsaturated soil hydraulic parameters using theensemble Kalman filter. Vadose Zone J. 10 (4), 1205–1227.
Merz, R., Parajka, J., Blöschl, G., 2011. Time stability of catchment model parameters:
implications for climate impact analyses. Water Resour. Res. 47, W02531. https://doi.org/10.1029/2010WR009505.
Mohanty, B.P., Ankeny, M.D., Horton, R., Kanwar, R.S., 1994. Spatial analysis ofhydraulic conductivity measured using disc infiltrometers. Water Resour. Res. 30(9), 2489–2498.

14