Particle filters are becoming increasingly popular for state and parameter estimation in hydrology. One of their crucial parts is the resampling after the assimilation step. We introduce a resampling method that uses the full weighted covariance information calculated from the ensemble to generate new particles and effectively avoid filter degeneracy. The ensemble covariance contains information between observed and unobserved dimensions and is used to fill the gaps between them. The covariance resampling approximately conserves the first two statistical moments and partly maintains the structure of the estimated distribution in the retained ensemble. The effectiveness of this method is demonstrated with a synthetic case – an unsaturated soil consisting of two homogeneous layers – by assimilating time-domain reflectometry-like (TDR-like) measurements. Using this approach we can estimate state and parameters for a rough initial guess with 100 particles. The estimated states and parameters are tested with a forecast after the assimilation, which is found to be in good agreement with the synthetic truth.

Mathematical models represent hydrological and other geophysical systems. They aim to describe the dynamics and the future development of system states. These models need the current state and certain system parameters (e.g., material properties and forcing) for state prediction. Both state and system parameters are typically ill known and have to be estimated.

Data assimilation methods, originally used for state estimation only, are adapted to also
estimate parameters and other model components like the boundary condition.
The ensemble Kalman filter

As the EnKF is based on Bayes' theorem, all uncertainties have to be
represented correctly; otherwise the method has a poorer performance

The particle filter has already been used for state and parameter estimation
for various hydrological systems. Since parameters do not have their own
model dynamics like the state, the resampling step is crucial.

Further development of the resampling for parameter estimation was done by

In this paper we introduce the covariance resampling, a resampling method that generates new particles using the ensemble covariance. This method conserves the first two statistical moments in the limit of large numbers while partly maintaining the structure of the estimated distribution in the retained ensemble. With the covariance, the unobserved parameters of the new particles are correlated to the observed state dimensions. The particle filter with covariance resampling is able to estimate state and parameters in case of a difficult initial condition without additional model evaluations, which are necessary for MCMC methods.

The particle filter is an ensemble-based sequential data assimilation method
that consists of a forecast and an analysis step. The ensemble members are
called particles. It is used to combine information from observation and
the model for a posterior estimate. For a detailed review, consider, for example,

If new information in the form of observations becomes available, the system
is propagated forward to the time the observation is taken (forecast). This
results in a prior probability density function. The prior is updated with
the information of the observation to get the posterior. This is accomplished
using Bayes' theorem,

This describes the assimilation process for one observation. For a set of
observations

The particle filter is a Monte Carlo approach, which directly approximates the probability density functions by a weighted ensemble of realizations (particles). This direct sampling allows the particle filter to have non-Gaussian probability density functions. This is in contrast to, for example, the EnKF, which is also based on Bayes' theorem and Monte Carlo sampling but assumes Gaussian distributions.

The posterior distribution of the previous analysis

Using Eq. (

In general,

To estimate state and parameters simultaneously we use an augmented state. In our case the augmented state

Particle filters tend to filter degeneracy, which is also referred to as filter impoverishment. After several analysis steps, one particle gets all statistical information as its weight becomes increasingly large, whereas the remaining particles only get a small weight such that the ensemble is effectively described by this one particle. In this case, the filter does not react to new observations and follows the particle with the large weight.

The stochastic universal resampling

Illustration of the universal resampling process. A random number

If the model does not have a stochastic model error, like we consider in this study, it is necessary to perturb the particles; otherwise they would be identical and the filter would still degenerate. Even in the presence of a model error it can be useful to perturb the particles after the resampling step. For example if the model error is ill known or structurally incorrect, it can help to guarantee a sufficient ensemble spread and diversity.

There are different suggestions for how to perturb. For example,

We neither perturb the duplicated states nor draw a complete new ensemble.
The covariance resampling we propose uses the universal resampling – other
resampling methods can be equally used – to choose the ensemble members
that are kept. Instead of duplicating the particles and setting the weights to

The total ensemble reduces to

Since the dropped particles are sampled from a Gaussian distribution, the mean and the covariance are conserved in the limit of large numbers. However, the structure of the non-Gaussian distribution is only partly conserved through the retained ensemble. In more difficult situations, where an increasing fraction of particles is resampled, the posterior is dominated by the approximated multivariate Gaussian distribution. However, the approximation allows the use of the covariance information in the ensemble, which facilitates the generation of meaningful new particles and improves the exploration of the state space. In less difficult situations, when only a few particles are resampled, the distribution remains close to that previously estimated, which includes the full structure of the estimated distribution.

Using the multivariate Gaussian distribution utilizes the information of the covariance but sacrifices the more accurate description of the univariate distribution that could be achieved by a kernel density estimation. However, it requires a much smaller sample size compared to a multivariate kernel density estimation.

The whole resampling process is
illustrated in Fig.

New particles are generated with a Cholesky decomposition of the covariance matrix.
The calculation of the covariance from the ensemble can result in small numerical errors
that have to be regularized; otherwise the decomposition would fail.
For details about the generation of new particles and regularization of
the covariance matrix see Appendix

Illustration of the particle filter with covariance resampling. The
green bars show the weight of each ensemble member (10 in this example), and
the black dots show the position of them.

The algorithm is tested using a synthetic case study of a one-dimensional
unsaturated porous medium with two homogeneous layers. The system has a
vertical extent of 1

The dynamics in an unsaturated porous medium can be described by
the Richards equation:

The macroscopic material properties represent the averaged subscale physics
with the functions

True Mualem–van Genuchten parameters and range of the uniformly distributed initial guess.

Fixed Mualem–van Genuchten parameters.

For the true trajectories and the observations, parameters by

Since state and parameters are estimated, we separate the model equation Eq. (

The forcing is reflected in the boundary condition of the simulation. For the
lower boundary, a Dirichlet condition with zero potential (groundwater table)
is used. The upper boundary condition is chosen as a flux boundary (Neumann),
representing four rain events with increasing intensity and dry periods in
between (see Fig.

Using infiltrations with an increasing intensity has the advantage that the system has more time to adjust the parameters. This way fewer observations are necessary for resolving the infiltration front and the information of the observations can be incorporated in the state and parameters. The stronger infiltration front in the end is used to test the robustness of the estimated state and parameters.

Upper boundary condition for the data assimilation case. Four rain events (blue) followed by a dry period (orange) are used for the test of the data assimilation run. After this run, two additional rain events and dry periods are used in a free run to test the assimilation results (grey background). The intensity and duration of these events is set to be equal to the first events of the data assimilation run. Note the different axes for infiltration and evaporation.

The last component of the system is the state. Initially, the system is in
equilibrium and will be forced by the boundary condition. The initial state
is depicted in Fig.

For the initial state of the data assimilation, the observations at time zero
are used. The measured water content is interpolated linearly between the
measurements and kept constant towards the boundary. Additionally, the
saturated water content for loamy sand, which is

The approximated state is perturbed by a correlated multivariate Gaussian distribution.
The main diagonal of the covariance matrix is

The use of the covariance increases the diversity of the ensemble and also helps to avoid degeneration. Using uncorrelated Gaussian random numbers with a zero mean would result in a fast degeneration of the particle filter caused by the dissipative nature of the system. The perturbation would simply dissipate, and the ensemble would collapse.

Initial state for the data assimilation run. Observations (purple)
at time zero are connected linearly and set constant towards the layer and
upper boundary. For the lower boundary, the saturated water content

The initial parameters for the ensemble are uniformly distributed. The ranges
of the uniform distributions are given in
Table

After the assimilation, the ensemble is used to run a forecast without data
assimilation, and the mean is calculated from the propagated ensemble using the weights of the last assimilation time.
In this run two additional infiltration events
and evaporation periods (see Fig.

The development of the parameters is depicted in
Fig.

If dynamic is induced in the system, the ensemble spread in the water content
space increases because of different parameter sets. This makes the particles
and their corresponding parameter sets distinguishable and parameter
estimation possible. The parameters

Estimation of all six parameters (

To see the effect of the parameters on the forward propagation, it is
necessary to have a closer look at the conductivity function Eq. (

Conductivity function

The final water content state estimated with the particle filter agrees with
the synthetic truth in a narrow band, while the mean of the ensemble
propagated without data assimilation is far-off (see Fig.

Final water content state after the assimilation run. The truth
(black dashed line) is almost congruent
with the estimated mean (orange line),
so only the 95 % quantile of the ensemble (light orange) is visible. The
final ensemble with the corresponding weights is used to start a free forward
run afterwards. The mean without data assimilation (green line) is calculated from
a forecast of the initial ensemble (see Figs.

To analyze the ensemble, we take a closer look at the effective sample size
and the number of particles that are resampled. The effective sample size is
defined as

Figure

The effective sample size is a crucial parameter for the covariance
resampling. If it drops to low values the filter can degenerate because,
effectively, too few particles contribute to the weighted covariance (Eq.

Amount of particles that are resampled (orange), and the effective sample size (green dots). The lines connecting the dots are for better visibility of the time-dependent behavior. The effective sample size increases while the number of resampled particles decreases. Every infiltration reduces the effective sample size and leads to more resampled particles.

For further analysis, the RMSE is calculated based on the difference of the
true water content and the weighted mean at every observation time. This
includes also the unobserved dimensions. The RMSE shows a similar behavior,
like the parameters and the effective sample size (see Fig.

The RMSE (red line) of the
water content calculated between the truth and the estimated mean including
all dimensions.
After 160 h the free run starts (grey background). The mean of the free run is calculated using
the propagated ensemble members with their corresponding weights at the last assimilation time.
During this time, the RMSE is about

After the data assimilation, the final ensemble including the weights are
used for a forecast. This forward run without data assimilation shows that
the RMSE oscillates in a narrow range. These oscillations are part of the
unobserved space next to the boundary and are mainly caused by the wrong
value of

For the presented case study, this section explores two issues relevant when applying the proposed covariance
resampling method: (i) the choice of the factor

To explore the convergence of the particle filter with
covariance resampling, the filter was run for

The RMSE and the standard deviation of the water content at the
last observation time (160

Figure

Figure

The tuning factor has similarities to multiplicative inflation for the EnKF

Model errors are omnipresent in real systems. They can have a structural or
stochastic nature and different intensities, and they can manifest, for example,
as biases in the results. For data assimilation of real measurements, the
consideration of model errors is an important part for the success of the
methods. Several extensions and modifications to sequential data assimilation
methods have been discussed (e.g.,

In the course of this paper, we briefly study the behavior of the particle
filter with covariance resampling for the case of a biased upper boundary
condition. Two cases are considered, one with a 10 % bias and one with a 20 % bias towards less water. This means that the amount of rain is reduced and the
evaporation is increased by these percentages. The observations are still
generated using the previous boundary condition (Fig.

Except for the ensemble size and the upper boundary condition the setup is
identical to that in Sect.

Final water content state after assimilation run using a bias for the upper boundary
condition of

Figure

Conductivity function

The conductivity function (see Fig.

We introduced a resampling method for
particle filters that uses the covariance information of the ensemble to
generate new particles and effectively avoids filter degeneracy. The method
was tested in a synthetic one-dimensional unsaturated porous medium with two
homogeneous layers. Even with just a rough initial guess, a broad parameter
range and only

The covariance connects information between observed and unobserved
dimensions. This has some similarity to the ensemble Kalman filter, but in our
approach, information from the non-Gaussian distribution is partly maintained in
the retained ensemble. Even though the RMSE of the water content includes the
unobserved state dimensions, it stays within a narrow range (RMSE is about

Transferring the information to the unobserved dimensions helps the filter in
not degenerating when only a rough initial guess is available. The states and
parameters are both altered actively. For the used initial condition,
perturbing the parameters only

The effective sample size (Eq.

The filter performance can be increased by a tuning parameter

Different parameter sets can approximately describe the same conductivity
function (Eq.

The covariance resampling connects observed with unobserved dimensions to effectively estimate parameters and prevent filter degeneracy. It conserves the first two statistical moments in the limit of large numbers, while partly maintaining the structure of the non-Gaussian distribution in the retained ensemble. The method is able to estimate state and parameters in case of a difficult initial condition without additional model evaluations and using a rather small ensemble size.

The data used for the figures are available online from heiDATA (

The following pseudocode describes the covariance resampling for a single time

Correlated random numbers are generated using the Cholesky decomposition. We
use the LDLT decomposition which is
part of the Eigen3 library

The calculation of the Cholesky decomposition (LDLT version) is only possible if
the matrix is not indefinite. Mathematically, a covariance matrix has to be positive semidefinite:

For this purpose, the identity matrix

In our experiments, the smallest variance in the main diagonal of the covariance matrix is
still 16 orders of magnitude larger than

The saturated conductivity in the second layer is analyzed in the same setup
as in Sect.

The mean saturated conductivity in the second layer after the data
assimilation run for

Figure

For all ensemble sizes, the filter either degenerates or finds the true
value. Increasing the ensemble size increases the number of successful runs.
The degeneration of the filter can directly be seen in the effective sample
size, which drops to

For the case

In Fig.

The apparent bias towards a larger saturated conductivity for

Increasing the factor to

The supplement related to this article is available online at:

DB designed, implemented, performed and analyzed the presented study. DB developed the algorithm with contributions from HHB. All authors participated in continuous discussions. DB prepared the paper with contributions from all authors.

The authors declare that they have no conflict of interest.

We thank the editor Harrie-Jan Hendricks Franssen and the two reviewers, Damiano Pasetto and Jasper Vrugt, for their comments, which helped to improve this paper.

This research is funded by the Deutsche Forschungsgemeinschaft (DFG) through project RO 1080/12-1. Hannes H. Bauser and Daniel Berg were supported in part by the Heidelberg Graduate School of Mathematical and Computational Methods for the Sciences (HGS MathComp), funded by DFG grant GSC 220 of the German Universities Excellence Initiative. Edited by: Harrie-Jan Hendricks Franssen Reviewed by: Jasper Vrugt and Damiano Pasetto