Coupled hydrogeophysical parameter estimation using a sequential Bayesian approach

Coupled hydrogeophysical methods infer hydrological and petrophysical parameters directly from geophysical measurements. Widespread methods do not explicitly recognize uncertainty in parameter estimates. Therefore, we apply a sequential Bayesian framework that provides updates of state, parameters and their uncertainty whenever measurements become available. We have coupled a hydrological and an electrical resistivity tomography (ERT) forward code in a particle filtering framework. First, we analyze a synthetic data set of lysimeter infiltration monitored with ERT. In a second step, we apply the approach to field data measured during an infiltration event on a full-scale dike model. For the synthetic data, the water content distribution and the hydraulic conductivity are accurately estimated after a few time steps. For the field data, hydraulic parameters are successfully estimated from water content measurements made with spatial time domain reflectometry and ERT, and the development of their posterior distributions is shown.


Introduction
In recent years, the worth of geophysical methods in hydrological applications has increasingly become apparent (e.g.Hubbard et al., 1999;Rubin and Hubbard, 2005;Linde et al., 2006;Vereecken et al., 2006).Typically, hydrological investigations rely on methods that disturb the soil, like soil coring, tensiometry or Time Domain Reflectometry (TDR).In contrast to that, non-invasive geophysical measurements provide the possibility to eavesdrop on subsurface flow and transport processes without disturbing them.
Correspondence to: J. Rings (j.rings@fz-juelich.de)This way, spatio-temporal patterns of hydrological states can be retrieved.In hydrogeophysical applications, hydrological states and parameters are estimated through observations of geophysical states (e.g.electrical resistivity).The possibly unknown parameters encompass both hydrological parameters determining flow and transport and petrophysical parameters that relate geophysical and hydrological state variables.
Geophysical surveys conducted to estimate subsurface states or parameters have been approached in a three-step manner (see e.g.Binley et al., 2002;Kemna et al., 2002;Chen et al., 2004;Cassiani and Binley, 2005;Koestel et al., 2008;Müller et al., 2010).First, a geophysical survey retrieves the geophysical state variables of the subsurface (e.g. the electrical resistivity distribution).Next, a petrophysical relationship is applied to transfer these into hydrological state variables (e.g.water content or tracer concentration).Finally, these state variables are used to drive a hydrological model or they are used in a parameter estimation process.
The numerical inversion methods to obtain the geophysical state variables from the survey data have to solve a highly nonlinear and mixed-determined problem.Typically, the resolution of the inversion result differs spatially, so that some regions may be well resolved while others are prone to exhibit interpretation errors (Day-Lewis et al., 2005).Furthermore, available knowledge about factors driving hydrologic processes (e.g.amount of infiltration or precipitation) does not enter into this conversion although it largely influences the observations.Rings and Hauck (2009) have studied the varying resolution for surface-based Electrical Resistivity Tomography (ERT).They found that unresolved contrasts and inversion artefacts may lead to quantification errors that would produce unphysical hydrological properties.Consequently, an improved inversion paradigm has become necessary.
Published by Copernicus Publications on behalf of the European Geosciences Union.
Integrated or coupled inversion approaches aim at inferring hydrological properties directly from the geophysical measurements.The parameters of the hydrological model and of the local-scale petrophysical relationship between geophysical and hydrological states are perturbed to minimize the difference between modelled and measured data.Several authors have investigated coupled hydrogeophysical inversion for hydrological modelling (Kowalsky et al., 2004;Lambot et al., 2006;Jadoon et al., 2008;Rucker, 2009;Hinnell et al., 2010).These contributions used synthetic data to evaluate the usefulness and applicability of coupled inversion for simplified hydrological systems.Additionally, newer studies are increasingly applying these methods to actual field data, where the structural inadequacies of the models, measurement errors and the assumptions inherent to the petrophysical relationship complicate the inversion problem (Kowalsky et al., 2005;Deiana et al., 2008;Looms et al., 2008).In a recent contribution, Huisman et al. (2010) have formulated the coupled inversion problem in a Bayesian framework to explicitly recognize uncertainty.They argued that the posterior uncertainty of the model parameters and predictions is useful to explore the value of different measurement types.
The coupled hydrogeophysical approach is especially suited for problems where hydrological states are observed by geophysical measurements on several occasions (i.e.timelapse geophysical surveys).When time-lapse geophysical measurements are used to monitor hydrological processes and to parameterize hydrological models, an inevitable question is whether additional measurements still provide additional information to improve the hydrological model parameterization.Although Huisman et al. (2010) showed that a Bayesian framework is useful to address this question, their framework only allows an answer a posteriori.For long-term monitoring, it is more attractive to provide regular updates of hydrological states, parameters and their uncertainty, using the information that is available so far.
Typically, filtering techniques, such as the popular Kalman filter (Kalman, 1960;Chen, 2003) are applied to update states whenever new measurements become available (e.g.Seppänen et al., 2001).The classical Kalman filter relies on linear behaviour in time, but additional filters have been developed that alleviate this problem, e.g. the extended Kalman filter (see e.g.Kaipio and Somersalo, 2004) or the ensemble Kalman filter (Evensen, 1994).Lehikoinen et al. (2009) have applied the extended Kalman filter to the problem of dynamical ERT inversion with explicit recognition of state uncertainty.However, all Kalman-type filters rely on Gaussian error distributions in the prior distribution.
A related family of filters are the particle filters.These filters are based on a sequential Bayesian framework (Gordon, 1993;Doucet and de Freitas, 2001;Arulampalam et al., 2002;Storvik, 2002;Chen, 2003;Cappe et al., 2007).They have been designed to cope with arbitrary prior distributions.The posterior probability distributions are approximated by sampling at discrete supporting points carried by weighted particles.This method is highly attractive for continuous state monitoring and also has potential for parameter estimation.A major advantage in a hydrogeophysical monitoring application is that the filter provides updated posterior distributions of states and parameters immediately after each measurement.Particle filters have recently been introduced into hydrology (see e.g.Weerts and El Serafy, 2006;Zhou et al., 2006;Hsu et al., 2009).Although filtering techniques are traditionally used to update state variables, they have also been applied in conjunction with hydrological models for parameter estimation (Kivman, 2003;Moradkhani et al., 2005;Vrugt et al., 2005;Hendricks Franssen and Kinzelbach, 2008;Salamon and Feyen, 2009).
In this study, we apply particle filtering to the problem of estimating hydraulic and petrophysical parameters from ERT measurements made during infiltration in the unsaturated zone.The rest of this study is organized as follows.First, we introduce the concept and implementation of the particle filter.Next, we apply the particle filter to estimate the water content distribution and the saturated hydraulic conductivity from a simple synthetic data set consisting of ERT measurements made during infiltration into a lysimeter.Finally, we apply the particle filter to estimate hydraulic and petrophysical parameters from ERT measurements made during infiltration into a full-scale dike model.As a reference, we also determine these parameters from detailled spatial time domain reflectometry (TDR) measurements.

State models
We want to describe the state x of a system, in our case e.g. a vector containing water content or pressure head for each element of a discretized representation of the subsurface.We can apply Bayesian analysis at each time step T =1,...,t,t+1,...At time t, this analysis can estimate past states x 1 ,...,x t−1 (smoothing), the current state x t (filtering) or the future state x t+1 (forecasting).x can be seen as a random variable in a stochastic model that translates the state in time: where the operator f is the time propagation function that translates the state in time from t to t+1 driven by some external forcing u t and modified by some static parameters ξ .α t+1 is Gaussian noise term that adds a stochastic diffusion process to the translation.From a state x t+1 , we can deduce an observation y t+1 from the observation process Hydrol.Earth Syst.Sci., 14, 545-556, 2010 www.hydrol-earth-syst-sci.net/14/545/2010/ where g is an operator to determine the system response from state and parameters, and β t+1 again is a Gaussian term describing the measurement noise.The noise term α t+1 and β t+1 are assumed to be independent random vectors.

Parameter estimation
The state-space has been formulated in the previous section, where we assumed the parameters ξ to be static.However, in this study we want to focus on parameter estimation rather than only filtering or forecasting of states.Therefore, we use a series of parameters ξ t , that shall converge to estimates of the true static parameters ξ for t → ∞.We assume these ξ t to be pseudo-static, meaning that they are not modified by external forcing, but only perturbed by a small Gaussian noise process as where ζ t+1 is an independent Gaussian random vector implemented as a static noise term as suggested by Liu and West (2000).For simplification, we abbreviate the notation as z=(x,ξ ), where z is also known as the augmented state variable.
We approximate this pdf by a discrete set of n=1...N particles z i by assigning a weight ŵi to each particle i.Then, we can write: where δ() are Dirac delta functions.Initially, all particles are assigned the same weight ŵ0 i = 1/N.Direct sampling from this posterior is generally very demanding or impossible, so we rather sample from a known distribution q(z 0:t+1 |y 1:t+1 ) called a proposal distribution.This way, we can define new weights w i as The choice of a proposal distribution is a critical decision (Arulampalam et al., 2002;Moradkhani et al., 2005), and is conveniently taken to be equal to the prior distribution The last step towards a sequential formulation is the assumption that the proposal distribution is chosen so that it factorizes in a recursive way: q(z 0:t+1 ) = q(z 0:t |y 0:t )q(z t+1 |z 0:t ,y t+1 ) (8) By combining Eqs. ( 7) and ( 8), we arrive at the simple form for weight updating Due to the Markovian formulation of the time propagation model in Eqs. ( 1) and ( 3) and the observation model in Eq. ( 2), this is a valid formulation and can be used to implement a particle filter.

Particle filter implementation
We implement a Sampling Importance Resampling (SIR) filter, which basically follows a three-step approach.The first step is the time propagation of state and parameters which samples the evolution density p(z t+1 |z t ) and follows from the models in Eqs. ( 1) and ( 3).
The second step is filtering, which assigns weights to the particles based on the probability for an observation.The observation is deduced from Eq. ( 2), then the weight is updated according to Eq. ( 9).In the case that measurements from j = 1,...,J different sources are available, we have vectors Y j of measurement data (e.g.measured transfer resistances or water content) and the modeled observations y j , we evaluate the weight w i,j as the inverse of the distance |Y t i,j −y t i,j |.If we assume to know the normalized, relative measurement error r j of each method, we can obtain w i = j r j w i,j .However, the filters implemented in this study have been restricted to one measurement method.
After assigning and normalizing the weights w i , resampling is used to prevent a concentration of weights in only one or a few particles.To find out whether this is necessary, an effective particle number N e is determined as If N e is larger than a threshold value τ (e.g.τ =0.5), the particle set remains unchanged.Otherwise, this is an indication that the particle set has become impoverished, as few particles contribute with considerable weight.To remedy this, a new set of particles is determined using a resampling algorithm.Douc et al. (2005)  The weight w i for all new particles is again set to w i =1/N.For T →∞ and N→∞, the posterior is found so that the weighted average ξ represents an appropriate parameter estimate: Figure 1 illustrates the particle filter method.We start with a particle cloud that has been propagated to time t.In the observation step, the particles are weighted according to how well they fit with the data.Then, they are redistributed according to their weights.A particle that has gained weight during observation may split up into two or more new particles, and particles that had small weight after observation may be removed from the cloud.Finally, the propagation carries the particle states to time t+1, which in our case would mean a forward model run from time t to t+1.These steps are repeated as long as new measurements become available.

Time propagation
For each particle i, the function f in Eq. ( 1) is evaluated by a hydrological forward model run from time t to t+1.The state x t i enters as the initial water content or pressure head at time t.
The hydrological models used in this study are HYDRUS (e.g.Simunek et al., 2008) for 2-D modelling domains and PARSWMS (Hardelauf et al., 2007) for 3-D modelling domains.They solve the Richards equation for unsaturated water flow using a linear Galerkin approach in a finite element scheme.We assume that the parametric model by van Genuchten (1980) can describe soil water retention and hydraulic conductivity as: where h is matrix potential (m), θ is the volumetric water content (m 3 /m 3 ), θ r is the residual water content for h→−∞, θ S is the saturated volumetric water content, α is the inverse of the air-entry value (m −1 ), n is an empirical shape factor (-), m is a factor that can be connected to n via m=1−1/n and k S is the hydraulic conductivity at θ=θ S (m/s).

Observations
To evaluate particle weights, the modelled observation y has to be determined.In this study, we concentrate on ERT as the geophysical observation method.To model the transfer resistances that are measured during an ERT survey, we have to determine the distribution of bulk resistivity ρ b from the soil water content distribution.Therefore, the petrophysical relationship by Archie (1942) is used: where ρ w is the pore water resistivity ( m), which is the inverse of the pore water conductivity σ w , is the porosity of the soil (m 3 /m 3 ), m A is the cementation exponent (-) and n A is the saturation exponent (-).Ulrich and Slater (2004) determined saturation exponents ranging from 1.0 to 2.7 for unconsolidated sands.
We use an ERT forward code by Ruecker et al. (2006) for 3-D modelling domains and the CRMOD code (Kemna et al., 2000) for 2-D modelling domains.Both ERT forward models were coupled to hydrological models following the approach and considerations presented in Huisman et al. (2010).

Numerical experiment
As a proof of concept, we use a synthetic data set of simulated lysimeter infiltration monitored by cross-borehole ERT.The model domain was set up as a cube of 1 m edge length filled with a homogeneous material (with θ r =0.001, θ S =0.27, k S =0.0015 m/s, α=4 and n=4.56).On each vertical face of the cube, two electrode arrays were installed at one third and two thirds of the width, with electrodes positioned below the surface from 0.1 m to 0.9 m every 0.1 m.Electrode arrays were created for all possible combinations of injection dipoles in one borehole and potential dipoles in the borehole immediately across the injection borehole at equal or greater depth.Model grids were created for the hydrological simulation with 1000 equally sized hexagonal cells that were further subdivided into 6000 tetrahedral cells.The parameter mesh for the ERT model was of equal size, but internally an irregular grid with tetrahedral cells was used (see Ruecker et al., 2006).
Water was infiltrated from all surface nodes into an initially dry medium (pressure head h=−1 m).After t=500 s, the wetting front reached the bottom of the volume.Artificial measurements were created every 100 s with the ERT forward model.For this simulation, a coarser internal ERT forward grid was used as recommended by Kaipio and Somersalo (2004).The measurements were perturbed with 3% Gaussian noise.

Implementation of the particle filter
For this 3-D simulation, the hydrological model PARSWMS (Hardelauf et al., 2007) was coupled to the ERT forward code by Ruecker et al. (2006).The state vector was initialized uniformly for each particle with 1331 pressure head values (one for each grid node) and a variable k S .The initial state and parameter distribution were sampled from uniform distributions for N=1000 particles.Distribution bounds and true values can be found in Table 1.The particle filter has been realized in a wrapper software code "Particle Filtering Inversion-Friendly Framework" (PFIFF) implemented in Python that connects the coupled models with the propagation, observation and resampling schemes.
For time propagation from t to t + 1, the hydrological model was run conditional to the particle's parameter k S .Then, α=ζ =2% Gaussian noise were added to the augmented state to simulate a stochastic diffusion process.
We assumed that the petrophysical relationship was known, and used it to transfer soil water content to bulk resistivity for each cell.The resistivity distribution thus obtained was used in the filtering step to simulate an ERT measurement at time t+1.Then, particle weights were assigned and normalized.Finally, the particles were resampled (each time step, so that the resampling threshold was τ =1) and their weights reset to 1/N.

Results
Figure 2 shows the distribution of particle weights immediately before the first resampling (t=100 s) and at the end of the simulation (t=500 s).The figure implicitly contains the variability in the state estimate, as can be seen from different weights of particles with similar conductivity.The initial distribution of the states was only h=4 cm wide, so it is evident that errors in the state have a large influence on the weight.While most particles have a low weight, some particles near the true value of k S =0.0015 m/s have a markedly higher weight.At the final time step, t=500 s, the filtering and resampling have caused most particles to converge to the correct value.The weighted average of k S is 0.00165 m/s, which is a slight overestimation.The state variable is shown for a 1-D vertical profile taken in the middle of the model domain for two time steps (Fig. 3).The initial distribution was sampled from an interval that underestimated the pressure head (see Table 1) for all particles.While at time t=300 s, the

Parameter
True value Initial range Initial pressure head −0.5 m −0.54...−0.50 m Saturated hydraulic k S 0.0015 m/s 0.0008...0.0042 m/s conductivity particle states show an even greater variability than in the initial distribution, the true state (blue line) is no longer on the boundary of the distribution, but the particles gather around it.This is true not only for the model parts near the surface, where changes due to the infiltrating water front have influenced the particle states, but also near the bottom, where the stochastic diffusion process that is added to the time propagation step combined with the resampling have influenced the particle states to move into the right direction.In the final time step, at t=500 s, the mean state value of all particles successfully retrieved the true state.Only in the lowest parts of the model (where ERT sensitivity is very low), the variability in pressure head is much higher.Finally, the results of the parameter estimation are shown as a comparison of the initial prior and posterior pdf of k S in the histograms in Fig. 4. The prior distribution was randomly sampled from a uniform distribution, but the parameter estimation successfully transformed the pdf into a posterior distribution approximating a Gaussian distribution with a mean near the true value.

Site description and measurements
Measurements have been taken on a full-scale dike model (see Fig. 5) located at the Federal Waterways Research Institute in Karlsruhe, Germany.It has a height of 3.6 m and a length of 22.4 m.The dike model was built homogeneously of a well-graded, uniform sand, which was covered by a thin soil layer overgrown by grass.Below the dike, there is a waterproof sealing that creates a hydrological no-flow boundary.The waterside can be flooded to just below the crest.At the foot of the landside slope, there is a drain that removes excess water and ensures the stability of the model dike.More details on the dike model can be found in Scheuermann et al. (2009) and Rings et al. (2008).
The dike model is equipped with 12 vertically installed TDR cables.By applying a reconstruction algorithm (Schlaeger, 2005), permittivity profiles can be obtained along these cables.These can be converted to spatially resolved soil water content profiles by applying a petrophysical relationship, in this case a site-specific calibration (see Scheuermann et al., 2009).In 2007, water content measurements were taken using TDR cables.During a similar experiment in 2005, ERT measurements were made using a SYSCAL Junior instrument (IRIS instruments) on a 8 m long line with 48 electrodes with a spatial spacing of 0.16 m down the landside slope of the dike (see Fig. 5).Each ERT measurement consisted of 348 discrete arrays in the Wenner-Schlumberger configuration with a fixed spacing of 0.16 m between potential electrodes and a current electrodes spacing ranging between 0.48 m and 4.32 m.From these 348 arrays, only the 120 arrays with the largest absolute change over time were selected for use in the particle filter to reduce the computational burden.
A flooding was simulated, in which the water level was raised to 2.4 m (1.2 m below the crest) within 48 h.Unfortunately, no TDR measurements were taken during the flooding experiment in 2005.On the other hand, fewer ERT measurements were taken down the landside slope in 2007.As both experiments were performed under very similar conditions (with the difference that the water level was kept at the highest level for an additional 48 h in 2007 with a subsequent more rapid lowering of the water level), we use ERT measurements taken at 12 different times over the course of 92 h from 2005 and TDR measurements taken at 14 different times over the course of 96 h in 2007 for two runs of the particle filter.Figure 6 shows the evolution of the water level and the measurement times for ERT and TDR.
At this point, we want to emphasize that the TDR data are far superior to the ERT data because of the large amount of invasive TDR data as compared to the non-invasive ERT data.This is confirmed by Huisman et al. (2010), which showed that the posterior uncertainty in the hydrologic parameters estimated from TDR is much smaller.Therefore, the TDR measurements should be considered as a reference to compare with the ERT parameter estimation results.This large difference in information content is also the reason that a simultaneous inversion of ERT and TDR data is not attempted in this study.

Implementation of the particle filter
We follow the model setup by Huisman et al. (2010), where a 2-D section coinciding with the ERT transect is modelled using the hydrological code HYDRUS (Simunek et al., 2008) coupled to the ERT forward code by Kemna et al. (2000).Different model grids were used for the hydrological and ERT model.For the hydrological model, a discretization with a total of 7603 elements was chosen with a denser distribution near the soil surface to account for steeper gradients of pressure head at the soil-atmosphere interface.In contrast, the ERT model grid needs only be refined near the electrode positions.A discretization with 5924 elements was used so that there was at least one node in between neighbouring electrodes.The domain was extended into the subsurface, as the hydrological no-flow boundary permits flow of electrical current.In Huisman et al. (2010), the electrical conductivity assigned to this domain extension was optimized as an additional parameter.In this study, we fixed the subsurface electrical conductivity to this optimized value.
The particle filter approach needs considerable computational resources, but has the benefit that propagation and observation can be run independently for each particle.Therefore, the PFIFF code was implemented in a parallelized ver- sion which uses a multi-threaded root process that distributes and starts model runs for different particles on other processors using MPI.
Two versions of the filter were run with different observation models.In the first version, the observation model simply compared modelled and measured soil water content derived from spatial TDR.The second version used the coupled electrical forward code to model ERT data at all times when field ERT measurements were taken; then particles were weighted by the inverse of the difference between measured and modelled apparent resistivity.The TDR run used N =3000 particles with 1% parameter noise, the ERT run used 4096 particles with 3.5% parameter noise to ensure sufficient exploration of the parameter space.Initial parameter estimates were sampled from uniform distributions (see Tables 2 and 3) by stratified sampling.In both filter runs, resampling was used if the effective particle number fell below a threshold of τ = 0.95.The TDR filter run took about 12 h on four quadcore CPUs, the ERT run took 50 h on 6 quadcore CPUs.

Results
Tables 2 and 3 show the results for the hydraulic and petrophysical parameters estimated either from water content measurements (TDR) or electrical measurements (ERT).For each estimation method, the weighted average (see Eq. 13), the median and the 90% confidence interval have been calculated from the posterior distribution.Figures 7 and 8 show, for the six parameter estimates, the posterior distributions.For each parameter, the development of the 98% and 50% confidence intervals, median and weighted average are shown.
The saturated hydraulic conductivity k S has been constrained by both methods.ERT has estimated a slightly higher k S , but both values are near the values determined in www.hydrol-earth-syst-sci.net/14/545/2010/ Hydrol.Earth Syst.Sci., 14, 545-556, 2010 The distributions of α and n in Fig. 7 have hardly been constrained.A strong correlation between these parameters was observed in Huisman et al. (2010), which leads to a very slow convergence as many possible pairs of α and n provide acceptable simulations.This is also apparent in Fig. 9. Here, the distribution of particle weights is shown as a function of the parameters n and α.The ridge of the weight landscape has approximately the same height for a large range of parameter combinations, leading to a slowly converging estimation.In the limited number of time steps of the filtering process, no convergence could be achieved, however, the median and weighted average approach the laboratory values.The petrophysical parameters n A and σ w have been well constrained as shown in Fig. 8.The value of n A agrees very well with the estimate of n A =1.16 by Rings et al. (2008).
A comparison of water content distributions using laboratory and ERT calibrated k S for the 2005 experiment at t=29 h is shown in Fig. 10.As the subsurface seen by ERT is limited to a depth of about 1.5 m on the landside (see Rings et al., 2008), ERT starts to detect changes in water content from t=20 h on.At t=29 h, the parameter k S as estimated by the coupled approach leads to a simulation with mostly established saturated area (Fig. 10, top row, right).The lower k S as determined in the laboratory, however, leads to a saturated area that would barely reach the region visible to ERT (Fig. 10,top row,left).After t=29 h, a saturated area is established and the measurements seem to carry few information to additionally constrain the posterior distributions.However, these time steps are important to allow the posterior distribution to converge through additional resampling steps.
The RMS error of the modelled ERT response in the last time step is 311 m with a variance of only 16 m.This is twice the RMSE as found by the MCMC optimization in Huisman et al. (2010), which had 12 free parameters, but a significant reduction from the initial RMS of 1966 m with a variance of 989 m.From the small variance of the RMSE in the final set of particles, it is concluded that the filter has approximately converged to the posterior distribution that reflects the uncertainty and information content of the measurements.
To evaluate the influence of the number of particles, we repeated the ERT filter run with only 400 particles.While the constraints on the parameter distributions look similar to those in Fig. 8 with only slightly wider bounds, the median and weighted averages differ considerably between runs with 400 and 4096 particles.Figure 11 shows the particle weight distributions in the k S parameter space at time t=78 h.The vertically striped structure is the effect of resampling and the Gaussian noise on the parameters.For the 4096 parameters, the density of particles is high enough so that it can be assumed that the particle distribution approximates the posterior distribution, although the space is not systematically explored.For 400 particles, due to resampling, gaps in the parameter distribution become apparent.In this case, the filter run still worked for 400 particles because of the lowdimensional parameter space.However, it is questionable whether our standard implementation of the particle filter can be successfully used to explore higher dimensional parameter spaces.A possible remedy could be a larger Gaussian noise on the parameters, e.g.expressed in terms of variance multiplied with a factor as suggested by e.g.?.The factor can also be systematically reduced over time.A combination of particle filters with optimization concepts borrowed from MCMC methods also seems promising to improve convergence and parameter space exploration issues.

Conclusions
We have presented a sequential Bayesian framework for estimation of hydrological states and parameters from hydrogeophysical measurements.This particle filter method approaches the pdfs of state and parameters by discrete sets of particles which each carry state and parameters sampled from initial distributions.Through time propagation and comparison to measurement data, the particles are assigned weights according to how well they describe the data.Over time and with the help of a resampling technique, the particles approach the true distributions.Furthermore, the filtering approach has the benefit of providing updated posterior distri- butions of states and parameters whenever new measurement data become available.
For a synthetic data set simulating lysimeter infiltration monitored by ERT, this method was shown to be able to retrieve the correct model states and parameters and provide a reasonable estimate of the remaining predictive uncertainty.However, it was also apparent that small errors in the state estimation can strongly affect parameter estimation, which necessitates the use of an increased number of particles.
For field data, the focus was on parameter estimation; therefore the states were not varied.The coupled hydrogeophysical approach has been applied to ERT and TDR data collected during infiltration experiments on a full-scale dike www.hydrol-earth-syst-sci.net/14/545/2010/ Hydrol.Earth Syst.Sci., 14, 545-556, 2010  model.This model dike allows experiments with temporal changes in water content up to full saturation.It was monitored by TDR cables in the dike body that measured soil water content, and by ERT perpendicular to the crest on the land side of the dike.Two particle filter runs were made.
In the first run, water content measured with TDR was used as the observation model, while in the second run only ERT measurements were analyzed using the coupled hydrogeophysical inversion approach.The first filter run estimated three hydrological parameters, while the second estimated k S and two petrophysical parameters.The parameter ranges were successfully constrained for k S and the petrophysical parameters.For all parameters, the weighted average of the parameter distributions were in good agreement with values obtained in laboratory measurements and previous studies.
The results of this study are encouraging and show that sequential Bayesian methods are an appropriate estimation technique for parameter estimation in hydrogeophysical surveys.Even more, as updates of the posterior distributions become available with new measurements, this can be used for on-line parameter estimation in a permanent monitoring installation, e.g. for the continuous monitoring of dike water content.The filter can also be used to forecast future states and to optimize experimental design.This may involve a spatial focussing in regions where change will probably occur or an optimization of the timing of measurements in a permanently installed system.
Filtering techniques are mostly applied for state estimation, where the propagation model places strict constraints on the variability.For parameter estimation, no such constraints exist.The presented particle filter implementation does not include a search strategy for the parameter space, but instead relies on stratified sampling of the initial distribution and random perturbations after each time step.This necessitates a large number of particles, which is computationally prohibitive for a larger number of parameters.The present technique provides a promising strategy for applications that focus on state estimation with few unknown parameters.To deal with higher-dimensional parameter spaces, the current particle filter implementations are most likely not suitable because of insufficient exploration of the parameter space.Future developments of particle filtering are required to overcome this limitation.In this respect, a combination of particle filters with state-of-the-art Monte Carlo methods that have been developed for (non-sequential) parameter space exploration seems especially promising.

Fig. 1 .
Fig. 1.Illustration of the filtering steps of observation and propagation and the resampling scheme.Adapted from Chen (2003).

Fig. 2 .
Fig. 2. Weights of particles before resampling steps at t=100 s and t=500 s.

Fig. 3 .
Fig.3.Vertical 1-D pressure head profiles.The blue line marks the true value, while each red line corresponds to one particle's state.Green lines mark the state pdf median and the 5% to 95% confidence interval.

Fig. 4 .
Fig. 4. Prior and posterior probability distribution of k S for the numerical case.

Fig. 5 .
Fig. 5. Full-scale dike model at the Federal Waterways Research Institute in Karlsruhe, Germany.

Fig. 6 .
Fig. 6.Heights of the water level in the 2005 and 2007 flooding experiment.Points mark the times where measurements were taken, either with ERT in 2005 or TDR in 2007.

Fig. 7 .
Fig. 7. Posterior probability obtained from TDR measurements.Distributions of k S , α and n are shown with the 98% percentile marked by light color, the 50% percentile marked in blue, the blue line marking the median and the green dots marking the weighted averages at each time step.Green diamonds on the axes mark the laboratory values.

Fig. 8 .
Fig. 8. Posterior probabilities obtained from ERT measurements.Distributions of k S , n A and σ w are shown with the 98% percentile marked by light color, the 50% percentile marked in blue, the blue line marking the median and the green dots marking the weighted averages at each time step.Green diamonds on the axes mark the laboratory values.

Fig. 9 .
Fig. 9. Particle weights, n and α for all particles in the filter run calibrated to SWC measurements at the final time step.

Fig. 10 .
Fig. 10.Water content distribution in the dike model at time t=29 h for the 2005 experiment using either the laboratory measured k S or the (final) k S estimated from ERT.

Fig. 11 .
Fig. 11.Particle weights and k S for ERT filter runs with 400 and 4096 particles at tme t=78 h.

Table 1 .
Initial state and parameter value and range.

Table 2 .
Initial and posterior parameter values for the dike and calibration to TDR.Units are m/s for k S and 1/m for α.

Table 3 .
Initial and posterior parameter values for the dike and calibration to ERT.Units are m/s for k S and ( m) −1 for σ w .