the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# Inflation method for ensemble Kalman filter in soil hydrology

### Hannes H. Bauser

### Daniel Berg

### Ole Klein

### Kurt Roth

The ensemble Kalman filter (EnKF) is a popular data assimilation method in soil hydrology. In this context, it is used to estimate states and parameters simultaneously. Due to unrepresented model errors and a limited ensemble size, state and parameter uncertainties can become too small during assimilation. Inflation methods are capable of increasing state uncertainties, but typically struggle with soil hydrologic applications. We propose a multiplicative inflation method specifically designed for the needs in soil hydrology. It employs a Kalman filter within the EnKF to estimate inflation factors based on the difference between measurements and mean forecast state within the EnKF. We demonstrate its capabilities on a small soil hydrologic test case. The method is capable of adjusting inflation factors to spatiotemporally varying model errors. It successfully transfers the inflation to parameters in the augmented state, which leads to an improved estimation.

Data assimilation combines information from models and measurements into an optimal estimate of a geophysical field of interest (Reichle, 2008). It has applications in all branches of the geosciences, with weather forecasting as the driving force behind many recent advances (van Leeuwen et al., 2015). The advantage of data assimilation methods (in contrast to, e.g., inverse modeling) is the possibility of considering model errors, which are characteristic of geophysical systems.

The ensemble Kalman filter (EnKF) (Evensen, 1994; Burgers et al., 1998) is a popular data assimilation method due to its simple conceptional formulation and ease of implementation (Evensen, 2003). It is an extension of the Kalman filter (Kalman, 1960) for nonlinear models.

In hydrology, the EnKF was used for soil moisture estimation from satellite data (e.g., Reichle et al., 2002; Crow and Van Loon, 2006) or from local measurements (e.g., De Lannoy et al., 2007, 2009; Camporese et al., 2009). However, the largest uncertainties in hydrology are associated with soil hydraulic material properties. They can neither be measured directly nor can they be transferred from the lab to the field, and are typically parameterized. Thus, including material properties in the estimation can be crucial in hydrology. Liu and Gupta (2007) called for an integrated assimilation framework including not only states, but also parameters and even model structure.

The joint estimation of states and parameters in data assimilation might be one possibility to reduce the influence of model errors on parameter estimation (Liu et al., 2012). Such a joint estimation in the EnKF with an augmented state was already demonstrated by Anderson (2001) for an atmospheric model. In hydrology Vrugt et al. (2005) combined an EnKF and the shuffled complex evolution Metropolis algorithm, while Moradkhani et al. (2005) used a dual EnKF approach to estimate states and parameters for a rainfall–runoff model. The joint assimilation of states and parameters in an augmented state was successfully performed for example in groundwater research (e.g., Chen and Zhang, 2006; Hendricks Franssen and Kinzelbach, 2008; Kurtz et al., 2012, 2014; Erdal and Cirpka, 2016), but also in soil hydrology for land surface models (e.g., Bateni and Entekhabi, 2012; Han et al., 2014; Zhang et al., 2017) and on smaller scales based on the Richards equation (e.g., Li and Ren, 2011; Wu and Margulis, 2011, 2013; Song et al., 2014; Erdal et al., 2014, 2015; Shi et al., 2015; Bauser et al., 2016; Brandhorst et al., 2017; Botto et al., 2018).

Due to unrepresented model errors and due to a limited ensemble size, the EnKF underestimates model errors, which can lead to filter inbreeding. Systematic model errors are common for example in land surface models (Vereecken et al., 2015). Additionally, in soil hydrology spatially and temporally varying model errors occur due to un- or ill-represented processes like preferential flow or hysteresis. Underestimated errors cause an insufficient ensemble spread in the augmented state. This is especially severe for parameters, which are typically not changed through a forward propagation and consequently cannot increase their uncertainty again. Due to the convergent dynamics in soil hydrology, the uncertainty in the state depends strongly on the parameter spread and becomes too small as well.

Covariance inflation can counteract filter inbreeding. Different methods have been proposed: (i) additive inflation, which adds a model error after the forward propagation. This method is especially useful if prior knowledge about the model error exists. In atmospheric sciences additive inflation has been successfully applied by, e.g., using reanalysis of historical weather prediction errors (Whitaker et al., 2008). (ii) Relaxation methods, which relax the analysis back to a prior perturbation or spread, have been proposed with tuning parameters (Zhang et al., 2004; Whitaker and Hamill, 2012) or based on deviations to measurements (Ying and Zhang, 2015). (iii) Multiplicative covariance inflation, which inflates the complete state with a scalar factor, where the inflation factor is either chosen manually (Anderson and Anderson, 1999) or is estimated based on deviations from measurements (e.g., Wang and Bishop, 2003; Anderson, 2007; Li et al., 2009). This method has been further extended to inflate each state component individually (Anderson, 2009).

All these inflation methods are developed in an atmospheric sciences context. Their transfer to soil hydrology is limited, due to the spatiotemporally varying model errors and the typically employed augmented state. For groundwater research, Kurtz et al. (2012) reported improved results by employing the inflation method by Anderson (2007), and Kurtz et al. (2014) used the constant inflation by Anderson and Anderson (1999). In soil hydrology, however, adjusted methods have been used: for example, Han et al. (2014) and Zhang et al. (2017) apply a special case of the inflation method by Whitaker and Hamill (2012) and keep the parameter spread constant to ensure a sufficient ensemble spread. Bauser et al. (2016) used the method by Anderson (2009), but adjusted the inflation of parameters.

Alternatively, no inflation method is reported (e.g., Li and Ren, 2011; Shi et al., 2015), but instead a damping factor (Hendricks Franssen and Kinzelbach, 2008), which can alleviate the issue, is employed. This is done by, e.g., Wu and Margulis (2011), Song et al. (2014), Erdal et al. (2014), Brandhorst et al. (2017), and Botto et al. (2018), where Erdal et al. (2014) and Brandhorst et al. (2017) combined this method with additive inflation.

In this paper, we propose a novel multiplicative inflation method, specifically designed for the needs in the soil hydrology community. The inflation method can vary rapidly in space and time to cope with the typically varying model errors and it is capable of a transfer of the inflation in the state to the parameters in the augmented state. The remainder of this paper is organized as follows: Sect. 2 describes (i) the EnKF, (ii) our proposed inflation method and (iii) a soil hydrologic test case. Section 3 shows the results of our method applied to the test case, followed by discussion and conclusion in Sects. 4 and 5.

## 2.1 Ensemble Kalman filter

The EnKF (Evensen, 1994; Burgers et al., 1998) is the Monte Carlo extension of the Kalman filter (Kalman, 1960) for nonlinear models and assumes unbiased Gaussian error distributions to combine model and measurement information. The filter is a sequential method and alternates between a forecast step and an analysis step. The forecast propagates a state including its uncertainty forward in time. The analysis combines uncertain model information with uncertain measurements at this time into an optimal estimate of the state. These two steps are now explained in more detail.

The forecast propagates an ensemble of states *φ*^{n}
forward from time *k*−1 to time *k* with a model *M*,

where the superscripts “f” and “a” denote forecast and analysis respectively,
while *n* denotes the ensemble members with *n*=1, …, *N*. The uncertainty in
the state is directly represented through the ensemble
${\mathit{\phi}}_{k-\mathrm{1}}^{\mathrm{a},n}$ and then propagated nonlinearly
with the model. Unrepresented model errors can be added through the unbiased
Gaussian process noise ** β**. This is also called additive
inflation. However, the details of the model error are typically unknown and
thus not represented adequately. The propagated uncertainties are directly
represented through the new forecast ensemble ${\mathit{\phi}}_{k}^{\mathrm{f},n}$.

The state can be extended by, e.g., model parameters ** ϕ** to
an augmented state

**=[**

*u***,**

*φ***]. This requires a forecast for each augmented state component. Parameters are typically assumed to be constant in time:**

*ϕ*The forecast of the state ${\mathit{\phi}}_{k}^{\mathrm{f},n}$ now also depends on the corresponding parameter set ${\mathit{\varphi}}_{k-\mathrm{1}}^{\mathrm{a},n}$. This way, uncertainties in the parameters are propagated as well and can be reduced jointly in the analysis.

Assuming unbiased Gaussian distributions, the ensemble of augmented states is
characterized through the forecast error covariance matrix **P**^{f},

where the overline denotes the expectation value and $\stackrel{\mathrm{\u203e}}{{\mathit{u}}_{k}^{\mathrm{f}}}$ is the ensemble mean.

The analysis combines model and measurement information based on the Gaussian
error assumption. The measurement error covariance matrix **R** of the
measurements ** d** is defined analogously as

where ** ϵ** is the measurement error. The measurements are
linked to the state through the linear measurement operator

**H**, which maps from the state space to the measurement space:

The Kalman gain **K** weighs the forecast error covariance matrix with
the measurement error covariance matrix and maps from the measurement space
back to the state space, based on the covariances in the forecast error
covariance matrix:

Based on the measurements, the Kalman gain updates the forecast ensemble to the analysis ensemble:

This update to the ensemble ${\mathit{u}}_{k}^{\mathrm{a},n}$ minimizes the analysis error covariance ${\mathbf{P}}_{k}^{\mathrm{a}}$, which fulfills

for infinite ensemble sizes.

Through spurious correlations and non-Gaussian distributions, ${\mathbf{P}}_{k}^{\mathrm{a}}$ will become too small, which can lead to filter inbreeding and ultimately filter divergence (e.g., Hamill et al., 2001). This is intensified if the model error required in Eq. (1) is unknown.

A common way to alleviate this issue in hydrology is the use of a damping
factor $\mathit{\gamma}\in [\mathrm{0},\phantom{\rule{0.125em}{0ex}}\mathrm{1}]$ (Hendricks Franssen and Kinzelbach, 2008), which is multiplied
to the correction vector in Eq. (7) and consequently
lessens the uncertainty reduction. The damping factor can be extended to a
vector ** γ** (and an entrywise multiplication) to treat augmented
state components differently (Wu and Margulis, 2011). Typically, parameters are
multiplied by a smaller factor than the state. However, the damping factor
can only alleviate and not completely prevent the inbreeding problem.

## 2.2 Multiplicative inflation for soil hydrology

Multiplicative inflation is another heuristic way to avoid filter inbreeding.
Anderson and Anderson (1999) proposed to increase the distance of each ensemble
member to the ensemble mean by multiplying this distance by $\sqrt{\mathit{\lambda}}$
for the inflation factor *λ*≥1:

This inflation factor is applied to the complete augmented state and has to
be adjusted to the specific problem. By construction, it does not alter the
mean value:
$\stackrel{\mathrm{\u203e}}{{\mathit{u}}_{\mathrm{inf}}^{\mathrm{f}}}=\stackrel{\mathrm{\u203e}}{{\mathit{u}}^{\mathrm{f}}}$.
A temporally varying inflation factor can be estimated by comparing
uncertainties with the distance of measurement and forecast
(e.g., Wang and Bishop, 2003; Anderson, 2007; Li et al., 2009). A spatiotemporally adaptive
inflation has been achieved by estimating a vector ** λ** for the
complete augmented state (Anderson, 2009). The author uses the
correlation between measurements and augmented state dimensions and asks this
question:

*How much inflation is required in each dimension to explain the observed differences to the measurements?*Anderson (2009) showed that this works excellently for the actual state. However, we experienced possible over-inflation in parameters (which do not have any dynamics to compensate for this), which can lead to filter collapses.

We propose a more conservative inflation method and ask this question:
*How much of the required change in the inflation are we allowed to transfer to the state dimensions based on the correlation information?* This
can be achieved by applying a Kalman filter for the inflation within the
EnKF.

In this Kalman filter, the inflation vector is treated as the state variable. As for parameters, we choose a constant model for the forecast in time:

For convenience we will drop the time subscript *k* in the following.
Furthermore, we will use the same symbols as for the EnKF, but denote them
with the subscript *λ*. We approximate the forecast error covariance
matrix for lambda, ${\mathbf{P}}_{\mathit{\lambda}}^{\mathrm{f}}$, based on the covariance
matrix of the augmented state in the EnKF, **P**^{f}, as the
normalized absolute correlation matrix of the augmented state ensemble. The
matrix component *i**j* is determined as

where *σ*_{λ} denotes the uncertainty of the inflation factors. It
is a tuning parameter that is kept constant over time and is assigned to all
state dimensions. It influences how fast the inflation factors are adjusted.
This follows the idea by Anderson (2007, 2009) to avoid a closure
problem, where the inflation estimation would require its own inflation.
Instead, the uncertainty is kept constant. Furthermore, only the absolute
value of the correlation is considered, since the inflation is based on
differences between measurement and model, but ignores their direction. Note,
that this presumes that the correlations of the model state can be
transferred to the inflation. In the presence of unknown model errors this
assumption may or may not be correct. However, the estimation at measurement
locations will remain meaningful in any case.

For the analysis, the distance *d*_{λ} between mean forecast
and measurement is used as measurement for ** λ**:

The measurement error covariance matrix **R**_{λ}
of *d*_{λ} can be calculated based on the error covariance
matrices of ** d** and $\mathbf{H}\stackrel{\mathrm{\u203e}}{{\mathit{u}}_{\mathrm{inf}}^{\mathrm{f}}}$,

where the inflated forecast error covariance matrix ${\mathbf{P}}_{\mathrm{inf}}^{\mathrm{f}}$
can be calculated from the inflation vector
and the forecast error covariance matrix by combining Eq. (9) (with vector ** λ**
and entrywise multiplication) and Eq. (3):
${\mathbf{P}}_{\mathrm{inf}}^{\mathrm{f}}={\mathbf{P}}^{\mathrm{f}}\circ \left[\sqrt{{\mathit{\lambda}}^{\mathrm{f}}}{\sqrt{{\mathit{\lambda}}^{\mathrm{f}}}}^{T}\right]$.
The entrywise product is denoted by ∘ and the entrywise square root
of

**by $\sqrt{\mathit{\lambda}}$.**

*λ*The expected distance between measurement and mean forecast based on the current inflation is

which combines the uncertainties of ** d** and
$\mathbf{H}\stackrel{\mathrm{\u203e}}{{\mathit{u}}_{\mathrm{inf}}^{\mathrm{f}}}$. To be able to
determine the Kalman gain, we first calculate the Jacobian matrix

**H**

_{λ}of partial derivatives of

*h*_{λ}with respect to

**:**

*λ*
With this approximated measurement operator **H**_{λ}, the
Kalman gain **K**_{λ} and the analysis state *λ*^{a}
are obtained as

Note, that the matrices ${\mathbf{P}}_{\mathit{\lambda}}^{\mathrm{f}}$ and **R**_{λ}
can possibly become indefinite, due to the absolute
value in Eqs. (11) and (13). Consequently, the
inverse in Eq. (16) could become unfeasible. However, we never
encountered such a case. In a situation with uncorrelated measurements, the
issue can be resolved by reducing *σ*_{λ} just for that single time step.

With this Kalman filter, the inflation vector is updated at each time step
based on the difference of the mean forecast to the measurements. Following
Anderson (2007), we additionally prohibit a deflation by constraining
the inflation values to (** λ**)

_{i}≥1.

## 2.3 Model

We test the proposed inflation method on a small hydrologic test case. We constructed it specifically to require a strong inflation. This makes it possible to explore features of the inflation in detail on a rather short timescale. Due to a small ensemble size, the results vary depending on the seed of the random numbers. This however, is related to different performance of the EnKF itself. In simulations (results are not shown), we found that the behavior of the inflation remains consistent. We have also tested the inflation method with real-world data by reanalyzing the application by Bauser et al. (2016) with the main result shown in Appendix B.

The Richards equation describes the change of volumetric soil water
content *θ* (–) in a continuous porous medium,

where *K* (L T^{−1}) is the isotropic conductivity and *h*_{m} (L)
is the matric head. Both are related to the water content. This relation is
typically described through parameterized material properties. We choose the
Mualem–van Genuchten parameterization (Mualem, 1976; van Genuchten, 1980),

with the saturation Θ (–),

The parameterization is described by a set of six parameters:
*θ*_{s} (–), *θ*_{r} (–), *α* (L^{−1}),
*n* (–), *K*_{0} (L T^{−1}), and
*τ* (–).

We additionally consider small-scale heterogeneity through Miller scaling. It
assumes geometrical similarity. With this the microscopic geometry of the
pore space at a macroscopic position is parameterized by a single length
scale *ξ* and the macroscopic heterogeneity field can be generated with a
single scalar field of this length scale. Miller and Miller (1956) showed that the
hydraulic functions scale with this parameter according to

where the functions *K*(*θ*) and *h*_{m}(*θ*) are defined at a
reference point ^{*} with Miller scaling parameter *ξ*=1 and from
there are projected to all locations.

For the test, we choose a one-dimensional case with a depth of 50 cm for a
time of 6 days. We set a groundwater table as the lower boundary condition
throughout the whole time and start from equilibrium conditions. The upper
boundary condition is no flux, except for a rain event with
$\mathrm{2.0}\times {\mathrm{10}}^{-\mathrm{7}}$ (m s^{−1}) during the fourth day. As observations we
choose two water content measurements at depths of 9.5 and 19.5 cm as they
would be available from time domain reflectometry (TDR). We set the
measurement uncertainty to a standard deviation of 0.007
(e.g., Jaumann and Roth, 2017).

As material we choose sandy loam from Carsel and Parrish (1988):
*θ*_{s}=0.41, *θ*_{r}=0.065, $\mathit{\alpha}=-\mathrm{7.5}$ m^{−1},
*n*=1.89, ${K}_{\mathrm{0}}=\mathrm{1.23}\times {\mathrm{10}}^{-\mathrm{5}}$ m s^{−1}, and *τ*=0.5. For the
Miller scaling we choose *ξ*_{1}=0.32 at the upper measurement position and
*ξ*_{2}=3.2 at the lower measurement position. We reduce the description of
the heterogeneity to these two parameters. The full function of the scaling
factors is calculated by linearly interpolating between the measurement
positions and constantly extrapolating to the boundaries.

The forward simulations are performed using MuPhi (Ippisch et al., 2006) with a spatial resolution of 1 cm. This corresponds to a state with 50 dimensions.

To test the inflation method, we perform a perfect model experiment. With the
EnKF we estimate the water content state and four parameters (*ξ*_{1},
*ξ*_{2}, *K*_{0}, and *τ*) through the augmented state
** u**=[

**, log**

*θ*_{10}(

*ξ*

_{1}), log

_{10}(

*ξ*

_{2}), log

_{10}(

*K*

_{0}),

*τ*]. We choose to include the logarithm of

*ξ*

_{1},

*ξ*

_{2}, and

*K*

_{0}, because we expect a more linear relation to the water content state than for the actual parameters. For the water content state, we use the correct initial condition as the mean with an uncertainty of 0.005. The uncertainty is spatially correlated using the fifth-order piecewise rational function by Gaspari and Cohn (1999) with the length scale

*c*=5 cm. As an initial guess for the parameters, we start with unknown heterogeneity ${\mathrm{log}}_{\mathrm{10}}\left({\mathit{\xi}}_{\mathrm{1}}\right)={\mathrm{log}}_{\mathrm{10}}\left({\mathit{\xi}}_{\mathrm{2}}\right)=\mathrm{0.0}\pm \mathrm{0.25}$, corresponding to 2 SD (standard deviations) away from the true values of ${\mathrm{log}}_{\mathrm{10}}\left({\mathit{\xi}}_{\mathrm{1}}\right)=-\mathrm{0.5}$ and log

_{10}(

*ξ*

_{2})=0.5. For the saturated hydraulic conductivity, we choose a too small value of ${\mathrm{log}}_{\mathrm{10}}\left({K}_{\mathrm{0}}\right)=-\mathrm{5.5}\pm \mathrm{0.5}$,

*K*

_{0}in (m s

^{−1}), which is about 1 SD away from the true value of ${\mathrm{log}}_{\mathrm{10}}\left({K}_{\mathrm{0}}\right)=-\mathrm{4.9}$. For the tortuosity $\mathit{\tau}=\mathrm{0.5}\pm \mathrm{0.5}$ we start from the true value.

Through the unrepresented heterogeneity, we can mimic a model error, leading
to a bias towards smaller values for the estimation of *K*_{0} during times
without dynamics, which may necessitate inflation. The parameter *τ* is
expected to have a smaller influence, since the uncertainty is chosen small
and it is already at the true value. This way it can act as an indicator
parameter for the inflation as it does not require inflation.

The EnKF is set up with a total of 25 ensemble members and a damping vector
of ** γ**=[

**1.0**, 0.3, 0.3, 0.3, 0.3], which we also apply to the inflation. The damping factor of 0.3 is applied to the parameters to alleviate issues of nonlinear relations between observations and parameters. For the uncertainty of the inflation factors we choose

*σ*

_{λ}=1.0.

We estimate the water content state together with the four parameters *ξ*_{1},
*ξ*_{2}, *K*_{0} and *τ* with the EnKF as described in Sect. 2.3.
The development of the water content at the two measurement
locations at a depth of 9.5 and 19.5 cm is shown together with the
inflation factor at these locations in Fig. 1. The inflation
factor is applied to the forecast ensemble before the analysis. The standard
deviation of the inflated ensemble should describe the distance of the
estimated mean to the synthetic truth. Note, that the inflation factor is not
based on this distance and relies on the noisy measurements. Therefore, it is
only an indicator.

During the first 3 days without any dynamics, the uncertainty for the upper measurement is slightly underestimated, while the uncertainty in the lower measurement is slightly overestimated. This leads to an inflation factor of basically 1 for the lower measurement (factors smaller than 1 are not allowed), while the inflation factor for the upper measurement is larger. However, due to correlations between the measurement locations a stronger inflation to fully explain the difference to the truth is prevented.

The deviation from the synthetic truth is induced through the initial guess of no heterogeneity and can also be seen in the systematic deviation of the inflated mean (which is equal to the forecast mean) from the analysis mean. When the infiltration front reaches the measurements, the deviations from the truth, underestimation of the uncertainty, and inflation factors increase rapidly. All of them are more pronounced for the upper measurement location. After the main peak, the differences and also the inflation factors decrease rapidly again.

The inflation factor for the state is shown in Fig. 2. It shows the strong increase in the inflation factor during the infiltration and its fast decrease afterwards. The inflation is strongest at the measurement location at a depth of 9.5 cm. The inflation factor is transferred to the other state locations through the correlations, which decrease with distance. Directly below the measurement locations the inflation factors are increased less than above. This is due to the chosen interpolation of the Miller scaling factors. Through the interpolation between the measurement locations and extrapolation to the boundaries, the dynamics changes at the measurement locations. During the infiltration the dynamics is mainly influenced by the water content above and the correlations with these locations are stronger.

The development of the Miller scaling factors *ξ*_{1} and *ξ*_{2} at the
two measurement positions (9.5 and 19.5 cm depth) is shown in
Fig. 3a and b together with the estimated inflation factor
for these parameters. Both initial conditions assume no heterogeneity and
start at ${\mathrm{log}}_{\mathrm{10}}\left({\mathit{\xi}}_{\mathrm{1}}\right)={\mathrm{log}}_{\mathrm{10}}\left({\mathit{\xi}}_{\mathrm{2}}\right)=\mathrm{0.0}\pm \mathrm{0.25}$, corresponding to
2 SD away from the true value. At the upper location the
true value of ${\mathrm{log}}_{\mathrm{10}}\left({\mathit{\xi}}_{\mathrm{1}}\right)=-\mathrm{0.5}$ corresponds to a finer material.
Consequently, the water content drops, as seen in Fig. 1,
leading to a strong correlation with the scaling factor, and
log_{10}(*ξ*_{1}) is adjusted rapidly to lower values. Accordingly, the
inflation factor is increased quickly in the beginning and then reduced back
to 1 when the estimation of log_{10}(*ξ*_{1}) reaches and eventually
underestimates the true value. The underestimation of the scaling factor
corresponds to a too fine material, which leads to slower changes in the
water content state and therefore smaller correlations. The scaling factor is
corrected during the rain event on the fourth day, which also leads to an
inflation.

The initial guess for the scaling factor for the depth of 19.5 cm underestimates the scaling factor, which corresponds to a too fine material. Again, the correlations are small. The value increases slowly during the dry period in the beginning, but is inflated and adjusted strongly during the rain event.

The saturated hydraulic conductivity *K*_{0} (Fig. 3c) was
chosen to start a little more than 1 SD below the true value.
Due to the unrepresented heterogeneity in the beginning, the value decreases
even further. The inflation remains very small due to correlations with both
measurement locations. However, as soon as the infiltration event reaches the
first measurement location, the value is corrected towards the true value. At
the same time, the inflation factor is increased due to the too small
uncertainty. After the rain event the inflation factor drops rapidly back to
1. The hydraulic conductivity remains below the true value. Another rain
event would be required to improve the estimation further.

The tortuosity *τ* (Fig. 3d) also influences the
hydraulic conductivity function, but has in this case a much smaller impact
and consequently smaller correlations with the measurements than *K*_{0}. We
use it as an indicator parameter and start at the true value. During the
infiltration event the value is changed due to its correlation. The
corresponding inflation factor is increased as well, but remains small enough
and drops back to 1 quickly enough to not cause any over-inflation of the
parameter.

To emphasize the need for a fast-adapting inflation factor, we reduce the
uncertainty of the inflation factors to *σ*_{λ}=0.5 to slow down
their adjustment. The results are summarized in Fig. 4.
The inflation of the water content state (Fig. 4a) shows
that the inflation factor does not reach as high values as before (see
Fig. 2). To compensate for this, the inflation acts
over a longer period of time. The same effect is also observed in the
inflation of the parameters (Fig. 4b and c). This leads
to a smaller inflation during the rain event and consequently a too small
uncertainty. At later times, when the cause of the error is not active any
more, the correlations with measurement locations are reduced, leading to a
slower reduction of the inflation in the parameters. In the indicator
parameter *τ* the beginning of an over-inflation can be seen towards later
times. This necessitates a more rapid inflation when correlations are used to
update inflation information.

The results for the parameters *K*_{0} and *τ* of a run without inflation
(and only damping) are shown in Fig. 5. Again, *K*_{0} moves
further away from the true value due to the unrepresented heterogeneity and
comes closer to the true value during the infiltration event. However, since
the Miller scaling factor is not inflated in the beginning, it is adjusted
slower. Consequently, the *K*_{0} is corrected longer in the wrong direction.
The uncertainty eventually becomes too small and in the end the mean is more
than 5 SD away from the true value, since the uncertainty
cannot be increased any more.

The proposed inflation method uses a Kalman filter to estimate inflation factors within the EnKF. It is based on the difference between measurements and mean forecast state. It transfers correlations from the forecast of the augmented state to the inflation. Consequently, the performance will be limited if model errors are structurally not represented in the forecast error covariance matrix. The estimation of the inflation factors with a Kalman filter is, like the EnKF itself, based on a linearized analysis. The use of a damping factor can alleviate issues with estimating nonlinear dependent parameters. To keep the inflation consistent with the analysis in the EnKF, we apply the same damping factor for both.

We designed a small synthetic hydrologic test case for the inflation. This test case mimics a model error through initially unrepresented heterogeneity. We designed the test case so that a strong temporally varying inflation is necessary, as it can occur with real data. We choose a short time so that the details of the behavior of the method can be explored. The method showed that it is capable of inflating states and parameters. The inflation is adjusted quickly and differentiates between parameters with strong and not so strong correlations. No over-inflation of weakly correlated parameters occurred. In this specific test case the estimation with inflation is far superior to an estimation without inflation.

The fast adjustment speed of the inflation factor is important because of the
fast-changing model errors and correlations with parameters. The adjustment
speed is determined by the uncertainty of the inflation factor. This
uncertainty is set to a constant value and has to be adjusted. For all our
cases a value of *σ*_{λ}=1 was sufficient, but larger values were
possible too. The need for such a fast adjustment is shown by estimating the
same case with a reduced uncertainty of *σ*_{λ}=0.5, which leads
to a slower adaptation of the inflation factor. This leads to smaller
inflation factors, which is compensated for by maintaining them for a longer
period of time. In this test case this leads to inflation at times after the
infiltration front has passed the measurements already and the model error is
small again. This can cause over-inflation of weakly correlated parameters.
Too large uncertainties of the inflation (in our test case
*σ*_{λ}=4), where the uncertainty is larger than the typical range
for the values of lambda, can also lead to over-inflation of weakly
correlated parameters. Reasons for this can be the linearizations in the
analysis and the calculation of the Jacobian (Eq. 15). This
limits the adjustment speed of the inflation.

Fast-dropping correlations between measurements and parameters are a limit for the method. An example could be a parameter only acting on an infiltration boundary condition. After the infiltration is over, correlations with this parameter would drop to zero and the inflation factor for this parameter will not be changed any more. If the inflation factor is not equal to 1 at this time, the parameter spread will keep increasing. In such a case, when there is no correlation, the parameter should be excluded from the estimation and consequently also from the inflation.

The method is in principle capable of compensating unrepresented model errors. However, it relies on correlations calculated from the forecast ensemble of the augmented state. If parameters have correlations with measurement locations with underestimated forecast uncertainties, the inflation will keep increasing the parameter spread until the forecast uncertainties are increased sufficiently. Therefore the correlations have to contain useful information. This means that inflating the parameters based on their correlations with measurement locations has to increase the forecast spread at these measurement locations. If the parameters have an insufficient influence on the state uncertainty, an over-inflation of the parameters can occur. An example are measurements with underestimated measurement uncertainties and short times between measurements compared to the timescale of the dynamics. Then the parameters are not able to increase the state uncertainty in the short forecast time between measurements and the forecast dynamics is not able represent the measurement noise. If such errors occur intermittently, e.g., the closed-eye period as proposed by Bauser et al. (2016) could be used to bridge these times. A rather heuristic solution could be a decay of the inflation factor towards values of 1, as already proposed by Anderson (2009).

In this work we propose a novel spatiotemporally adaptive inflation method, specifically designed for soil hydrology, which nevertheless is expected to work in similar systems as well. The inflation method is based on a Kalman filter acting within the EnKF. The method is capable of rapid adjustments of inflation factors, treating each augmented state dimension individually. This rapid adjustment is required due to temporally varying model errors, as they can appear through violation of the local equilibrium assumption of the Richards equation, hysteresis, or unrepresented heterogeneity.

We demonstrate the use of our inflation method in combination with a damping
factor on a small hydrologic example. We choose heterogeneity as a possible
model error, but allow the heterogeneity to be estimated along with the soil
hydrologic parameters *K*_{0} and *τ* of the Mualem–van Genuchten
parameterization. Our proposed inflation method proved to be stable in
combination with parameter estimation. The performance of the estimation
improved and parameter uncertainty remained consistent. The method requires
that the correlations from in the forecast ensemble contain useful
information for the inflation. However, we demonstrate that it even works for
only weakly correlated parameters. We expect the inflation method to
generally improve data assimilation with the EnKF and to thus lead to better
state and parameter estimations in soil hydrology.

The synthetic data are available upon request from the corresponding author.

We briefly show the derivation of the Jacobian matrix **H**_{λ}
for the inflation (Eq. 15). Again, the entrywise product is
denoted by ∘ and the entrywise square root of ** λ** by $\sqrt{\mathit{\lambda}}$:

We also applied the inflation method to reanalyze the case presented earlier by Bauser et al. (2016), where measurements from 11 TDR probes were assimilated with an EnKF. There, the inflation method confirmed the behavior observed in the small synthetic case presented in this paper. For the details of the real-world case as well as the concept of the closed-eye period please refer to Bauser et al. (2016) or Bauser (2018, chap. 5), the latter of which also includes the reanalysis of the case.

In this paper, we only show the inflation related to the closed-eye period (Fig. B1), which presents the major challenge to the inflation in that particular application. During this time, preferential flow occurs and the underlying local equilibrium assumption of the Richards equation is violated. With a standard approach, parameters become biased to compensate these errors. To avoid this, Bauser et al. (2016) introduced the closed-eye period, which pauses the parameter estimation and only guides the water content states through times, when assumptions are violated. Compared to the standard approach, this leads to a reduced bias in the parameters, but effectively increases the model errors during the closed-eye period. A strong inflation is required to compensate this error. The inflation method used in Bauser et al. (2016) was just able to accomplish this and the authors were concerned that a too slow adjustment speed of the inflation limits the applicability of the closed-eye period for cases with larger model errors.

Figure B1 confirms the fast adjustment speed of the new inflation method proposed in this paper for the real-world application. The strong required inflation stays well within the closed-eye period. This enables the EnKF to pick up the parameter estimation after the period from a water content state consistent with the TDR measurements and facilitates the use of the closed-eye period.

HHB designed, implemented, performed and analyzed the presented study. DB provided computational software. All authors participated in continuous discussions. HHB prepared the manuscript with contributions from all authors.

The authors declare that they have no conflict of interest.

We thank editor Insa Neuweiler and two anonymous reviewers for their comments, which helped to improve this paper.

This research is funded by the Deutsche Forschungsgemeinschaft (DFG) through
project RO 1080/12-1 and the Ministerium für Wissenschaft, Forschung und
Kunst Baden-Württemberg (Az 33-7533.-30-20/6/2). HGS MathComp provided
travel expenses for Hannes H. Bauser and Daniel Berg.

Edited by: Insa Neuweiler

Reviewed by: two
anonymous referees

Anderson, J. L.: An ensemble adjustment Kalman filter for data assimilation, Mon. Weather Rev., 129, 2884–2903, https://doi.org/10.1175/1520-0493(2001)129<2884:AEAKFF>2.0.CO;2, 2001. a

Anderson, J. L.: An adaptive covariance inflation error correction algorithm for ensemble filters, Tellus A, 59, 210–224, https://doi.org/10.1111/j.1600-0870.2006.00216.x, 2007. a, b, c, d, e

Anderson, J. L.: Spatially and temporally varying adaptive covariance inflation for ensemble filters, Tellus A, 61, 72–83, https://doi.org/10.1111/j.1600-0870.2008.00361.x, 2009. a, b, c, d, e, f

Anderson, J. L. and Anderson, S. L.: A Monte Carlo implementation of the nonlinear filtering problem to produce ensemble assimilations and forecasts, Mon. Weather Rev., 127, 2741–2758, https://doi.org/10.1175/1520-0493(1999)127<2741:AMCIOT>2.0.CO;2, 1999. a, b, c

Bateni, S. M. and Entekhabi, D.: Surface heat flux estimation with the ensemble Kalman smoother: Joint estimation of state and parameters, Water Resour. Res., 48, w08521, https://doi.org/10.1029/2011WR011542, 2012. a

Bauser, H. H.: Knowledge Fusion in Soil Hydrology, PhD thesis, Ruperto-Carola University Heidelberg, Heidelberg, Germany, https://doi.org/10.11588/heidok.00024713, 2018. a, b

Bauser, H. H., Jaumann, S., Berg, D., and Roth, K.: EnKF with closed-eye period – towards a consistent aggregation of information in soil hydrology, Hydrol. Earth Syst. Sci., 20, 4999–5014, https://doi.org/10.5194/hess-20-4999-2016, 2016. a, b, c, d, e, f, g, h

Botto, A., Belluco, E., and Camporese, M.: Multi-source data assimilation for physically based hydrological modeling of an experimental hillslope, Hydrol. Earth Syst. Sci., 22, 4251–4266, https://doi.org/10.5194/hess-22-4251-2018, 2018. a, b

Brandhorst, N., Erdal, D., and Neuweiler, I.: Soil moisture prediction with the ensemble Kalman filter: Handling uncertainty of soil hydraulic parameters, Adv. Water Resour., 110, 360–370, https://doi.org/10.1016/j.advwatres.2017.10.022, 2017. a, b, c

Burgers, G., van Leeuwen, P. J., and Evensen, G.: Analysis scheme in the ensemble Kalman filter, Mon. Weather Rev., 126, 1719–1724, https://doi.org/10.1175/1520-0493(1998)126<1719:ASITEK>2.0.CO;2, 1998. a, b

Camporese, M., Paniconi, C., Putti, M., and Salandin, P.: Ensemble Kalman filter data assimilation for a process-based catchment scale model of surface and subsurface flow, Water Resour. Res., 45, W10421, https://doi.org/10.1029/2008WR007031, 2009. a

Carsel, R. F. and Parrish, R. S.: Developing joint probability distributions of soil water retention characteristics, Water Resour. Res., 24, 755–769, https://doi.org/10.1029/WR024i005p00755, 1988. a

Chen, Y. and Zhang, D.: Data assimilation for transient flow in geologic formations via ensemble Kalman filter, Adv. Water Resour., 29, 1107–1122, https://doi.org/10.1016/j.advwatres.2005.09.007, 2006. a

Crow, W. T. and Van Loon, E.: Impact of incorrect model error assumptions on the sequential assimilation of remotely sensed surface soil moisture, J. Hydrometeorol., 7, 421–432, https://doi.org/10.1175/JHM499.1, 2006. a

De Lannoy, G. J. M., Houser, P. R., Pauwels, V. R. N., and Verhoest, N. E. C.: State and bias estimation for soil moisture profiles by an ensemble Kalman filter: Effect of assimilation depth and frequency, Water Resour. Res., 43, w06401, https://doi.org/10.1029/2006WR005100, 2007. a

De Lannoy, G. J. M., Houser, P. R., Verhoest, N. E. C., and Pauwels, V. R. N.: Adaptive soil moisture profile filtering for horizontal information propagation in the independent column-based CLM2.0, J. Hydrometeorol., 10, 766–779, https://doi.org/10.1175/2008JHM1037.1, 2009. a

Erdal, D. and Cirpka, O. A.: Joint inference of groundwater-recharge and hydraulic-conductivity fields from head data using the ensemble Kalman filter, Hydrol. Earth Syst. Sci., 20, 555–569, https://doi.org/10.5194/hess-20-555-2016, 2016. a

Erdal, D., Neuweiler, I., and Wollschläger, U.: Using a bias aware EnKF to account for unresolved structure in an unsaturated zone model, Water Resour. Res., 50, 132–147, https://doi.org/10.1002/2012WR013443, 2014. a, b, c

Erdal, D., Rahman, M., and Neuweiler, I.: The importance of state transformations when using the ensemble Kalman filter for unsaturated flow modeling: Dealing with strong nonlinearities, Adv. Water Resour., 86, 354–365, https://doi.org/10.1016/j.advwatres.2015.09.008, 2015. a

Evensen, G.: Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics, J. Geophys. Res.-Oceans, 99, 10143–10162, https://doi.org/10.1029/94JC00572, 1994. a, b

Evensen, G.: The ensemble Kalman filter: theoretical formulation and practical implementation, Ocean Dynam., 53, 343–367, https://doi.org/10.1007/s10236-003-0036-9, 2003. a

Gaspari, G. and Cohn, S. E.: Construction of correlation functions in two and three dimensions, Q. J. Roy. Meteorol. Soc., 125, 723–757, https://doi.org/10.1002/qj.49712555417, 1999. a

Hamill, T. M., Whitaker, J. S., and Snyder, C.: Distance-dependent filtering of background error covariance estimates in an ensemble Kalman filter, Mon. Weather Rev., 129, 2776–2790, https://doi.org/10.1175/1520-0493(2001)129<2776:DDFOBE>2.0.CO;2, 2001. a

Han, X., Hendricks Franssen, H.-J., Montzka, C., and Vereecken, H.: Soil moisture and soil properties estimation in the Community Land Model with synthetic brightness temperature observations, Water Resour. Res., 50, 6081–6105, https://doi.org/10.1002/2013WR014586, 2014. a, b

Hendricks Franssen, H. J. and Kinzelbach, W.: Real-time groundwater flow modeling with the ensemble Kalman filter: Joint estimation of states and parameters and the filter inbreeding problem, Water Resour. Res., 44, w09408, https://doi.org/10.1029/2007WR006505, 2008. a, b, c

Ippisch, O., Vogel, H.-J., and Bastian, P.: Validity limits for the van Genuchten–Mualem model and implications for parameter estimation and numerical simulation, Adv. Water Resour., 29, 1780–1789, https://doi.org/10.1016/j.advwatres.2005.12.011, 2006. a

Jaumann, S. and Roth, K.: Effect of unrepresented model errors on estimated soil hydraulic material properties, Hydrol. Earth Syst. Sci., 21, 4301–4322, https://doi.org/10.5194/hess-21-4301-2017, 2017. a

Kalman, R. E.: A new approach to linear filtering and prediction problems, J. Basic Eng., 82, 35–45, 1960. a, b

Kurtz, W., Hendricks Franssen, H.-J., and Vereecken, H.: Identification of time-variant river bed properties with the ensemble Kalman filter, Water Resour. Res., 48, w10534, https://doi.org/10.1029/2011WR011743, 2012. a, b

Kurtz, W., Hendricks Franssen, H.-J., Kaiser, H.-P., and Vereecken, H.: Joint assimilation of piezometric heads and groundwater temperatures for improved modeling of river-aquifer interactions, Water Resour. Res., 50, 1665–1688, https://doi.org/10.1002/2013WR014823, 2014. a, b

Li, C. and Ren, L.: Estimation of unsaturated soil hydraulic parameters using the ensemble Kalman filter, Vadose Zone J., 10, 1205–1227, https://doi.org/10.2136/vzj2010.0159, 2011. a, b

Li, H., Kalnay, E., and Miyoshi, T.: Simultaneous estimation of covariance inflation and observation errors within an ensemble Kalman filter, Q. J. Roy. Meteorol. Soc., 135, 523–533, https://doi.org/10.1002/qj.371, 2009. a, b

Liu, Y. and Gupta, H. V.: Uncertainty in hydrologic modeling: Toward an integrated data assimilation framework, Water Resour. Res., 43, w07401, https://doi.org/10.1029/2006WR005756, 2007. a

Liu, Y., Weerts, A. H., Clark, M., Hendricks Franssen, H.-J., Kumar, S., Moradkhani, H., Seo, D.-J., Schwanenberg, D., Smith, P., van Dijk, A. I. J. M., van Velzen, N., He, M., Lee, H., Noh, S. J., Rakovec, O., and Restrepo, P.: Advancing data assimilation in operational hydrologic forecasting: progresses, challenges, and emerging opportunities, Hydrol. Earth Syst. Sci., 16, 3863–3887, https://doi.org/10.5194/hess-16-3863-2012, 2012. a

Miller, E. E. and Miller, R. D.: Physical theory for capillary flow phenomena, J. Appl. Phys., 27, 324–332, https://doi.org/10.1063/1.1722370, 1956. a

Moradkhani, H., Sorooshian, S., Gupta, H. V., and Houser, P. R.: Dual state–parameter estimation of hydrological models using ensemble Kalman filter, Adv. Water Resour., 28, 135–147, https://doi.org/10.1016/j.advwatres.2004.09.002, 2005. a

Mualem, Y.: A new model for predicting the hydraulic conductivity of unsaturated porous media, Water Resour. Res., 12, 513–522, https://doi.org/10.1029/WR012i003p00513, 1976. a

Reichle, R. H.: Data assimilation methods in the Earth sciences, Adv. Water Resour., 31, 1411–1418, https://doi.org/10.1016/j.advwatres.2008.01.001, 2008. a

Reichle, R. H., McLaughlin, D. B., and Entekhabi, D.: Hydrologic data assimilation with the ensemble Kalman filter, Mon. Weather Rev., 130, 103–114, https://doi.org/10.1175/1520-0493(2002)130<0103:HDAWTE>2.0.CO;2, 2002. a

Shi, L., Song, X., Tong, J., Zhu, Y., and Zhang, Q.: Impacts of different types of measurements on estimating unsaturated flow parameters, J. Hydrol., 524, 549–561, https://doi.org/10.1016/j.jhydrol.2015.01.078, 2015. a, b

Song, X., Shi, L., Ye, M., Yang, J., and Navon, I. M.: Numerical comparison of iterative ensemble Kalman filters for unsaturated flow inverse modeling, Vadose Zone J., 13, https://doi.org/10.2136/vzj2013.05.0083, 2014. a, b

van Genuchten, M. T.: A closed-form equation for predicting the hydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. J., 44, 892–898, https://doi.org/10.2136/sssaj1980.03615995004400050002x, 1980. a

van Leeuwen, P. J., Cheng, Y., and Reich, S.: Nonlinear data assimilation, in: vol. 2, Springer International Publishing, Switzerland, https://doi.org/10.1007/978-3-319-18347-3, 2015. a

Vereecken, H., Huisman, J. A., Hendricks Franssen, H. J., Brüggemann, N., Bogena, H. R., Kollet, S., Javaux, M., van der Kruk, J., and Vanderborght, J.: Soil hydrology: Recent methodological advances, challenges, and perspectives, Water Resour. Res., 51, 2616–2633, https://doi.org/10.1002/2014WR016852, 2015. a

Vrugt, J. A., Diks, C. G. H., Gupta, H. V., Bouten, W., and Verstraten, J. M.: Improved treatment of uncertainty in hydrologic modeling: Combining the strengths of global optimization and data assimilation, Water Resour. Res., 41, w01017, https://doi.org/10.1029/2004WR003059, 2005. a

Wang, X. and Bishop, C. H.: A Comparison of breeding and ensemble transform Kalman filter ensemble forecast schemes, J. Atmos. Sci., 60, 1140–1158, https://doi.org/10.1175/1520-0469(2003)060<1140:ACOBAE>2.0.CO;2, 2003. a, b

Whitaker, J. S. and Hamill, T. M.: Evaluating methods to account for system errors in ensemble data assimilation, Mon. Weather Rev., 140, 3078–3089, https://doi.org/10.1175/MWR-D-11-00276.1, 2012. a, b

Whitaker, J. S., Hamill, T. M., Wei, X., Song, Y., and Toth, Z.: Ensemble data assimilation with the NCEP global forecast system, Mon. Weather Rev., 136, 463–482, https://doi.org/10.1175/2007MWR2018.1, 2008. a

Wu, C.-C. and Margulis, S. A.: Feasibility of real-time soil state and flux characterization for wastewater reuse using an embedded sensor network data assimilation approach, J. Hydrol., 399, 313–325, https://doi.org/10.1016/j.jhydrol.2011.01.011, 2011. a, b, c

Wu, C.-C. and Margulis, S. A.: Real-time soil moisture and salinity profile estimation using assimilation of embedded sensor datastreams, Vadose Zone J., 12, https://doi.org/10.2136/vzj2011.0176, 2013. a

Ying, Y. and Zhang, F.: An adaptive covariance relaxation method for ensemble data assimilation, Q. J. Roy. Meteorol. Soc., 141, 2898–2906, https://doi.org/10.1002/qj.2576, 2015. a

Zhang, F., Snyder, C., and Sun, J.: Impacts of initial estimate and observation availability on convective-scale data assimilation with an ensemble Kalman filter, Mon. Weather Rev., 132, 1238–1253, https://doi.org/10.1175/1520-0493(2004)132<1238:IOIEAO>2.0.CO;2, 2004. a

Zhang, H., Hendricks Franssen, H.-J., Han, X., Vrugt, J. A., and Vereecken, H.: State and parameter estimation of two land surface models using the ensemble Kalman filter and the particle filter, Hydrol. Earth Syst. Sci., 21, 4927–4958, https://doi.org/10.5194/hess-21-4927-2017, 2017. a, b