UvA-DARE ( Digital Academic Repository ) Inverse modelling of in situ soil water dynamics : investigating the effect of different prior distributions of the soil hydraulic parameters

In situ observations of soil water state variables under natural boundary conditions are often used to estimate the soil hydraulic properties. However, many contributions to the soil hydrological literature have demonstrated that the information content of such data is insufficient to accurately and precisely estimate all the soil hydraulic parameters. In this case study, we explored to which degree prior information about the soil hydraulic parameters can help improve parameter identifiability in inverse modelling of in situ soil water dynamics under natural boundary conditions. We used percentages of sand, silt, and clay as input variables to the ROSETTA pedotransfer function that predicts the parameters in the van Genuchten-Mualem (VGM) model of the soil hydraulic functions. To derive additional information about the correlation structure of the predicted parameters, which is not readily provided by ROSETTA, we employed a Monte Carlo approach. We formulated three prior distributions that incorporate to different extents the prior information about the VGM parameters derived with ROSETTA. The inverse problem was posed in a formal Bayesian framework and solved using Markov chain Monte Carlo (MCMC) simulation with the DiffeRential Evolution Adaptive Metropolis (DREAM) algorithm. Synthetic and real-world soil water content data were used to illustrate the approach. The results of this study demonstrated that prior information about the soil hydraulic Correspondence to: B. Scharnagl (benedikt.scharnagl@tu-bs.de) parameters significantly improved parameter identifiability and that this approach was effective and robust, even in case of biased prior information. To be effective and robust, however, it was essential to use a prior distribution that incorporates information about parameter correlation.


Introduction
Simulation of soil water dynamics under transient conditions typically requires knowledge of the soil hydraulic properties, that is, the water retention function and the hydraulic conductivity function. A broad array of methods exists to determine these two constitutive relationships from laboratory or field experiments. An overview of these methods together with a discussion of their strengths and limitations can be found in Durner and Lipsius (2005), amongst others. With the ever increasing pace of computational power, availability of accurate and stable numerical solution schemes of the governing flow equations, and effective and efficient parameter optimization methods, the use of inverse modelling to determine soil hydraulic properties has become increasingly popular in the last few decades. A review of Vrugt et al. (2008a) discusses recent progress in inverse modelling of soil hydraulic properties. Laboratory methods such as the multistep outflow method (van Dam et al., 1994) have the advantage of being comparatively quick and precise. This allows the processing of a large number of samples, which opens up the possibility to study the spatial variability of the soil hydraulic properties. However, soil hydraulic properties derived from laboratory experiments on small soil cores are typically inadequate to simulate soil water dynamics at larger spatial scales (Ritter et al., 2003;Mertens et al., 2005;Guber et al., 2006;Wöhling et al., 2008;Baroni et al., 2010). There are multiple reasons for this discrepancy, most notably that the sample volume analysed in the laboratory is not a representative elementary volume (e.g. Mallants et al., 1997) or that the experimental conditions dictated in the laboratory are not representative for field conditions (e.g. Basile et al., 2003). Arguably, field methods such as the internal drainage method (Zhang et al., 2003) provide estimates of the soil hydraulic properties that are more representative for in situ soil water dynamics. However, such methods place a high demand on equipment and time and are labour intensive. Moreover, considerable difficulties arise when these local scale soil hydraulic properties are used to infer the effective retention and hydraulic conductivity function that characterize soil water dynamics at larger spatial scales (e.g. Smith and Diekkrüger, 1996;Zhu and Mohanty, 2002). Effective properties are defined here as the soil hydraulic properties of an equivalent homogeneous domain that produces the same response as the actual heterogeneous domain under some upscaled boundary conditions .
As an alternative to laboratory or field experiments, soil hydraulic properties can be derived from field measurements of soil water state variables under naturally occurring boundary conditions (Vereecken et al., 2008). One main advantage of this approach is that it allows for estimating effective soil hydraulic properties at larger spatial scales given observations at multiple locations within the considered soil domain. Several applications of this approach can be found in the soil hydrological literature. Jacques et al. (2002) estimated the soil hydraulic properties of a four-layer soil profile using pressure head and water content data collected at 12 different locations and 5 depths along a 5.5 m long trench. To minimize problems with nonuniqueness and to reduce the dimensionality of the inverse problem, they used a stepwise parameter estimation procedure that sequentially estimates the soil hydraulic parameters for each individual soil layer. Another study by Ritter et al. (2003) used soil water content measured at 6 locations and 3 depths within a 4800 m 2 field plot to infer the effective soil hydraulic properties of a three-layer soil profile. They also reported problems in finding welldefined values of the soil hydraulic parameters. To alleviate these problems, Ritter et al. (2003) fixed some of the soil hydraulic parameters at values derived from laboratory experiments. Wöhling et al. (2008) derived effective soil hydraulic properties of a three-layer soil profile using pressure head observations from 3 locations and 3 depths. They compared the efficiency of three multiobjective search algorithms in finding Pareto solutions of soil hydraulic parameters that characterize the trade-off in the fitting of pressure head data at different depths. Steenpass et al. (2011) estimated effective soil hydraulic properties of a two-layer profile from observed soil surface temperature data and measurements of spatially distributed soil water content at 36 locations and 2 depths within a 6 m × 6 m plot. They implemented a Bayesian inference scheme using Markov chain Monte Carlo (MCMC) simulation and reported considerable uncertainty in the estimated soil hydraulic parameters. Some of the parameters attained physically unrealistic values, which Steenpass et al. (2011) attributed to a lack of information in the wet range. Finally, Wöhling and Vrugt (2011) explored whether using observations of two different soil water state variables helps to better constrain the soil hydraulic parameters. They estimated effective soil hydraulic properties of a four-layer profile using pressure head and soil water content measurements from 3 locations and 5 depths. Depending on what kind of state variables were used for model calibration, the estimated parameters derived with MCMC simulation differed substantially, while the uncertainty in these estimates was generally small.
In summary, these studies suggest that in situ observations of soil water dynamics contain insufficient information to warrant accurate and precise estimation of the soil hydraulic properties. A general approach to improve parameter identifiability in case of data with limited information content is the inclusion of prior information about the parameters of interest. In this study, we define "parameter identifiability" as the antithesis of "parameter uncertainty" (Vrugt et al., 2003a). We will refer to a parameter as being identifiable if the uncertainty in its estimate is reasonably small. Note that our definition slightly differs from the classical definition used in inverse problem theory (e.g. Carrera and Neuman, 1986b). The use of prior information is well established in groundwater hydrology and hydrogeophysics (e.g. Carrera and Neuman, 1986a;Woodbury and Ulrich, 1993;Kowalsky et al., 2004). However, only few soil hydrological studies have investigated the effect of using prior information about the soil hydraulic parameters in inverse modelling of soil water dynamics. Wang et al. (2003) estimated soil hydraulic properties of four materials of a layered soil profile using neutron probe measurements of soil water content collected in 9 boreholes distributed evenly in a 50 m × 50 m plot during an extensive infiltration experiment. In their Bayesian analysis, they included prior information about the mean values and variances of the soil hydraulic parameters as provided by Carsel and Parrish (1988) to estimate the parameter values that have maximum posterior density. Mertens et al. (2004) estimated effective soil hydraulic properties for a two-layer profile using water content observations from 25 locations and 3 depths within a 80 m × 20 m hillslope plot. They used prior information derived from laboratory and field experiments to formulate a prior probability density function (pdf) of the soil hydraulic parameters. This prior distribution was subsequently used in the Generalized Likelihood Uncertainty Estimation (GLUE) method (Beven and Binley, 1992) to generate random samples of the soil hydraulic parameters. Mertens et al. (2004) found that incorporation of prior information improved the sampling efficiency of the GLUE method and that the resulting prediction uncertainty bounds better enclosed the observational data. Finally, Hou and Rubin (2005) employed the principle of minimum relative entropy (Woodbury and Ulrich, 1993) to infer prior distributions of the soil hydraulic parameters from prior information that incompletely characterizes these distributions. To illustrate this approach, they used prior information about the expected values, variances, and lower and upper bounds of the soil hydraulic parameters derived from the ROSETTA database (Schaap et al., 2001). This prior distribution was then applied to calibrate a onelayer vadose zone model against observations of soil water content profiles derived from neutron probe and time domain reflectometry (TDR) measurements and to investigate the effect of the length of the calibration period on model predictive uncertainty.
In this study, we investigated the effect of three different prior distributions of the soil hydraulic parameters in inverse modelling of in situ soil water dynamics. We used prior information on sand, silt, and clay percentages and translated this information into prior information on soil hydraulic parameters using the ROSETTA pedotransfer function (Schaap et al., 2001). Sand, silt, and clay percentages constitute basic soil information and are readily available in most vadose zone studies, making this approach widely applicable. In addition to the standard ROSETTA prediction that provides the mean values and standard deviations of the predicted parameters, we tested a Monte Carlo approach to derive the correlation structure of the predicted parameters. We formulated three prior pdfs that incorporate to different extents the prior information derived with ROSETTA. A Bayesian framework was used to combine the various prior pdfs with the information contained in observations of soil water content in a 50 m × 50 m bare soil plot exposed to natural boundary conditions. The resulting posterior distribution was explored by MCMC simulation using the DiffeRential Evolution Adaptive Metropolis (DREAM) algorithm (Vrugt et al., 2008b). The present study had two main objectives. The first objective was to explore the potential benefit of including additional information on parameter correlation in the prior distribution of the soil hydraulic parameters. This information helps to better constrain the parameter space, and we may expect that this additional constraint improves the identifiability of the estimated parameters. The second objective was to test the robustness of the Bayesian approach in case of biased prior information.

Measurements
We measured soil water content at 61 locations within a 50 m × 50 m bare soil plot (Fig. 1). The measurement site was located at the bottom of a gently sloping (<3 • )  Herbst et al., 2009). The measurement plot itself was plane. The soil in this part of the field had a silt loam texture and was classified as a Stagnic Luvisol (IUSS Working Group WRB, 2007). Soil texture within the 50 m × 50 m plot was fairly homogeneous. The fine earth fraction of the topsoil (0 to 30 cm) was composed of 14 ± 1 % sand, 70 ± 1 % silt, and 16 ± 1 % clay (mean ± standard deviation, N = 47). Due to clay accumulation, the subsoil had a slightly higher percentage of clay, with 14 ± 2 % sand, 66 ± 2 % silt, and 20 ± 2 % clay (N = 12). The soil was kept bare during the measurement campaign. Accumulation of weeds was prevented by occasional application of herbicides and manual removal.
Soil water content was measured using TDR. Two-rod probes (25 cm rod length, 2.3 cm rod spacing) were installed horizontally 6 cm below the soil surface. The waveforms were recorded manually using a TDR100 device (Campbell Scientific, Logan, UT, USA) and analysed using the algorithm described in Heimovaara and Bouten (1990). We used the empirical relationship of Topp et al. (1980) to convert the apparent dielectric permittivity to soil water content. Measurements were taken on 29 days between 19 March and 14 October 2009, comprising a measurement campaign of 210 days. Because measurements of soil temperature and carbon dioxide efflux were carried out simultaneously with the TDR measurements, a complete measurement cycle on all locations took about three hours. The soil temperature and carbon dioxide efflux data were not used in this study. The soil water content data did not show any significant temporal or spatial trend and could be well described with a normal distribution at any given measurement date. We arithmetically averaged the soil water content data from the 61 locations and used this time series of mean water content to calibrate an effective soil hydraulic model. The mean values were assigned the midpoints of each measurement cycles.

Governing flow equation
We simulated one-dimensional vertical water flow in a 100 cm deep, homogeneous profile using the Richards equation: where θ (cm 3 cm −3 ) denotes the soil water content, K (cm h −1 ) represents the soil hydraulic conductivity, h (cm) signifies the pressure head, t (h) is time, and z (cm, positive upward) defines the vertical coordinate. We used the HYDRUS-1D model (Šimůnek et al., 2008) to solve the Richards equation for given initial and boundary conditions.

Soil hydraulic properties
The soil hydraulic properties were parametrised using the van Genuchten-Mualem (VGM) model (van Genuchten, 1980). The water retention function θ(h), expressed in terms of effective saturation S (dimensionless), is given by: where θ r and θ s (cm 3 cm −3 ) represent the residual and saturated water content, respectively, and α (cm −1 ), n, and m = 1−1/n (both dimensionless) are shape parameters. The hydraulic conductivity function K(h) is given by: where K s (cm h −1 ) is the saturated hydraulic conductivity, and L (dimensionless) is an additional shape parameter.

Spatial discretisation
The 100 cm deep profile was discretised into 81 nonequidistant nodes. Nodal distance was shortest adjacent to the soil surface and gradually increased with depth, with a distance of 0.05 cm at the upper and 3.5 cm at the lower boundary.
We chose this particular discretization scheme to accommodate the large gradients in pressure head that occur close to the surface in response to atmospheric forcing. In general, if nodal spacing is too large, the numerical solution of the Richards equation becomes inaccurate due to linearisation errors of the pressure head and averaging errors of the hydraulic conductivity (van Dam and Feddes, 2000). In the present case, a further increase in the number of nodes did not significantly alter the simulation results, from which we concluded that the spatial discretization of the soil profile was adequate.

Boundary conditions
The water flux at the soil surface is controlled by potential evaporation E pot (cm h −1 ) and precipitation P (cm h −1 ). The numerical solution of the Richards equation was obtained by limiting the actual water flux across the upper boundary as follows: where h min UB = −100 000 cm and h max UB = 1 cm denote the minimum and maximum values of the pressure head allowed at the upper boundary, respectively. If the simulated pressure head reaches either of these two limits, the HYDRUS-1D model switches to a pressure head boundary condition to calculate the actual water flux.
Precipitation and other meteorological variables were recorded at a meteorological station located 100 m west of the measurement site. Potential evaporation was estimated using the FAO method (Allen et al., 1998). The FAO method consists of two steps. First, the potential evapotranspiration from a grass reference surface, ET ref , is calculated using a modified Penman-Monteith equation (Allen et al., 1998, p. 74). We used hourly averaged values of air temperature, relative humidity, wind speed, incoming shortwave radiation, and atmospheric pressure as input variables. In a second step, the reference evapotranspiration is scaled with an empirical coefficient, E pot = 1.15ET ref (Allen et al., 1998, p. 263). This coefficient reflects the increased evaporation potential of bare soils (as compared to the reference grass surface), which is mainly due to the lower albedo of wet soil surfaces.
In the absence of detailed information about the lower boundary of the considered soil domain, we tested different lower boundary conditions in HYDRUS-1D. From an inverse modelling point of view, a zero gradient pressure head boundary condition is most appealing because it does not require explicit information about soil water state variables or fluxes at the lower boundary. Readings from a nearby piezometer suggested that the ground water table was about 200 cm below the lower boundary of the simulated profile. Simulations with the zero gradient boundary condition indicated that a substantial amount of water was lost from the profile due to drainage, resulting in simulated water contents that considerably underestimated the actually observed soil water content data. Repeated model runs with different realizations of the VGM parameters demonstrated that this difference was persistent and could not be explained by an inappropriate selection of the soil hydraulic parameters. We therefore used a prescribed constant pressure head h LB as the lower boundary condition: Unfortunately, no measurements of pressure head or soil water content were made at the bottom of the simulation domain from which an appropriate value of h LB could be inferred.
We therefore treated h LB as an unknown parameter that was estimated jointly with the VGM parameters.

Initial condition
All model runs started from a uniform initial pressure head throughout the profile equal to the pressure head at the lower boundary h LB . A 75 day spin-up period was used to allow for the relaxation from the initial pressure head distribution and to reduce sensitivity of simulation results to soil water state initialization.

Inverse modelling
Let the pressure head at the lower boundary be stored in The difference between the soil water content observations y=[ỹ 1 ,...,ỹ N ] and the corresponding HYDRUS-1D predicted values y=[y 1 ,...,y N ] was computed using the following residual vector: where N is the number of observations. A common approach is to aggregate ε(x 1 ,x 2 ) into a single measure of model performance and, depending on its definition, minimize or maximize this criterion during model calibration. If the inverse problem is posed in a probabilistic framework, this criterion is called the likelihood. It gives the probability of observing the data given the model parameters. Under the assumption of independent, identically and normally distributed residuals, ε∼N (0,σ 2 ε ), the likelihood is given as: Note that the likelihood is not only a function of the model parameters, x 1 and x 2 , but also of the standard deviation of the residuals, σ ε . The value of σ ε is typically unknown a priori because it integrates over various error sources such as measurement errors, model structural errors, and errors in model input variables. In practice, σ ε should therefore be considered as an unknown parameter that needs to be inferred from the observational data.
Bayesian inference provides a formal way of combining information from observations with prior information about the system. This is achieved through Bayes' theorem: where p(x 1 ), p(x 2 ), and p(σ ε ) denote the prior probabilities, and p(x 1 ,x 2 ,σ ε |ỹ) represents the posterior probability, that is, the probability of the parameters after assimilating the observational data. Note that Eq. (8) implies independence between prior information of x 1 and x 2 . Based on the underlying physics of the system, we may expect that x 1 and x 2 are correlated, but in the absence of detailed prior information about this correlation, the assumption of independence is justifiable. It is common practice in Bayesian inference to eliminate the standard deviation of the residuals from the inference equations. This is expedient because we are not particularly interested in σ ε . In addition, elimination of σ ε has the advantage of reducing the number of estimated parameters. Assuming a Jeffreys prior for σ ε , p(σ ε )∝1/σ ε , the standard deviation of the residuals can be integrated out of the inference equations (e.g. Box and Tiao, 1992;Kavetski et al., 2006). The likelihood function (Eq. 7) and Bayes' theorem (Eq. 8) then reduce to: In many practical applications of Bayes' theorem, prior information about the estimated parameters x 1 and x 2 is vague. In this case, a uniform prior distribution is usually imposed. The only information provided by a uniform prior is that of the bounds of the feasible parameter space. Within these bounds, all parameter values have equal probability. In the present case, we had rather limited information about the pressure head at the lower boundary, x 1 =[h LB ]. We therefore specified a uniform prior distribution for this parameter defined as p(x 1 )∼U(a x 1 ,b x 1 ), where a x 1 and b x 1 denote the lower and upper bounds, respectively (Table 1). The bounds were selected based on the available information about the depth of the ground water table and assuming a hydrostatic equilibrium between the lower boundary of the simulated profile and the ground water table. For the VGM parameters, x 2 =[θ r θ s α n K s L], we tested three different prior distributions, one being multivariate uniform, p(x 2 )∼U(a x 2 ,b x 2 ), and the other two being multivariate normal, p(x 2 )∼N (µ x 2 , x 2 ). These prior distributions are defined in Sect. 2.4.
For models that are nonlinear in their parameters, such as HYDRUS-1D, the posterior distribution p(x 1 ,x 2 |ỹ) cannot be obtained by analytical means nor by analytical approximation. We therefore resort to iterative methods that approximate the posterior pdf by generating a large sample from this  , 1953). To simplify notation, we merge the two parameter vectors x 1 and x 2 into a single vector x.
The Metropolis algorithm generates a Markov chain, which has the property that the next position of the chain x i+1 only depends on its current position x i . The Markov chain is generated by alternating between two basic steps. First, a proposal x is generated. Second, the proposal is accepted with probability: If the proposal is accepted, the Markov chain moves to the proposal position, x i+1 = x . Otherwise, the current position is retained, x i+1 = x i . To generate samples from the posterior distribution, we used the DREAM framework of Vrugt et al. (2008bVrugt et al. ( , 2009). This MCMC scheme runs multiple chains simultaneously for global exploration of the parameter space and automatically tunes the scale and orientation of the proposal distribution during evolution of the chains to the target distribution. This scheme is an adaptation of the Shuffled Complex Evolution Metropolis algorithm (Vrugt et al., 2003b) and has the advantage of maintaining detailed balance and ergodicity while showing excellent efficiencies on complex, highly nonlinear, and multi-modal target distributions. The use of multiple chains protects against premature convergence and opens up a wide array of statistical tests to diagnose whether the chains have converged to a stationary distribution or not. We used the most recent variant of DREAM that uses sampling from past states and a mix of parallel direction and snooker updates to generate proposals in each individual chain. This algorithm, entitled DREAM (ZS) , has several desirable advantages. First, sampling from the past circumvents the requirement of using a large number of chains for posterior exploration. Just a few chains will suffice. This will speed up convergence to the target distribution, especially for highdimensional problems. Second, outlier chains do not need explicit consideration. By sampling historical states, aberrant trajectories can jump directly to the modal region at any time during the simulation. Third, the transition kernel defining the jumps in each of the chains does not require information about the current states of the chains. This is of great advantage in a multi-processor environment where the various candidate points can be generated simultaneously so that each chain can evolve most efficiently on a different processor. An implementation of sampling from past states and snooker updating is described in ter Braak and . An application of DREAM (ZS) appears in Schoups and Vrugt (2010).
In this study, we used three parallel chains to explore the parameter space and approximate the posterior distribution. The diagnostic of Brooks and Gelman (1998) was used to check when convergence to the target distribution had been achieved. After convergence, we continued to run DREAM (ZS) and generated an additional 50 000 samples, which were used to summarize the posterior distribution.

Predicting the soil hydraulic parameters
The ROSETTA program (Schaap et al., 2001) implements five hierarchical pedotransfer functions to estimate the VGM parameters from a varying degree of basic soil data, such as textural class, texture, bulk density, and soil water content at specific pressure head values. We estimated the soil hydraulic parameters from measured sand, silt, and clay percentages of the topsoil layer at the experimental site. The corresponding pedotransfer function is labelled with H2-C2 in Schaap et al. (2001). This pedotransfer function consists of an ensemble of artificial neural network models that were each calibrated to a different data set. These data sets were generated from the ROSETTA database using the bootstrap method (Efron, 1979). The use of multiple calibration data sets provides a simple way to address uncertainty in the predicted soil hydraulic parameters. Each artificial neural network model provides slightly different estimates of the VGM parameters, and the ensemble means p = [p 1 ... p 6 ] and standard deviations u = [u 1 ... u 6 ] (both size 1×6) constitute the ROSETTA output (Table 2). ROSETTA uses log 10 -transformed values of α, n, and K s , which induces an approximate normal distribution for each of the VGM parameters (Schaap et al., 2001). This transformation scheme was retained in the present study.

Deriving the covariance matrix of the predicted soil hydraulic parameters
The information provided by ROSETTA can readily be used to formulate an informative prior distribution that is multivariate normal with mean vector µ x 2 = p and diagonal covariance matrix x 2 = diag(u 2 1 ,...,u 2 6 ). This approach, however, is rather simplistic in that it ignores the correlation that is typically found between the soil hydraulic parameters. Parameter correlation is evident from statistical analysis of soil hydrological databases (e.g. Carsel and Parrish, 1988;de Rooij et al., 2004). From an inverse modelling point of view, it is apparent from first-order approximations of the parameter covariance matrix (e.g. Kool and Parker, 1988), objective function contour plots (e.g. Toormann et al., 1992), or scatter plots of the posterior sample (Vrugt et al., 2003a). Given these findings, it seems productive to include the correlation of the soil hydraulic parameters in the prior pdf. In general, if we neglect correlation, we assign too high prior probabilities to physically unrealistic combinations of the soil hydraulic parameters.
To formulate prior pdfs that consider parameter correlation, we need detailed information about the correlation structure of the predicted parameters, p. This information is not provided by ROSETTA but can be derived using the following Monte Carlo approach. The basic idea of this approach is that the correlation structure of p can be inferred from a random sample drawn in close vicinity of p. Using this approach, we can derive the full covariance matrix of the predicted parameters p in four subsequent steps: Step 1:. Draw a random sample of input variables. Let the measured percentages of sand, silt, and clay used as input variables to ROSETTA be denoted by f (size 1×3). Generate a random sample of input variables F (size 1000×3) from a multivariate normal distribution, N µ f , f , centred around f : Based on a preliminary analysis, we assigned the diagonal entries of f (variances) an arbitrary small value of 0.25 %. This value works well in practice, ensuring that the sampled percentages are in close vicinity of f , and hence, that the corresponding random sample of soil hydraulic parameters will be in close vicinity of p.
The negative values of the off-diagonal terms (covariances) were chosen such that they are consistent with the compositional nature of soil texture and that f is positive semidefinite.
Step 2:. Use ROSETTA to generate the corresponding random sample of soil hydraulic parameters. Propagate the random sample of input variables F through ROSETTA to obtain a random sample of soil hydraulic parameters P (size 1000 × 6): We used quantile-quantile plots and pairwise scatter plots to test whether P follows a multivariate normal distribution. These diagnostic plots (not shown) indicated that each of the parameters had a marginal distribution close to normal and that the pairwise correlations among these parameters were approximately linear, providing further justification for the use of the multivariate normal distribution as the prior model.
Step 3:. Calculate the correlation matrix of the random sample of soil hydraulic parameters. Calculate the covariance matrix C (size 6 × 6) of the random sample of soil hydraulic parameters P. Use C to calculate the corresponding correlation matrix R (size 6 × 6) defined as: This correlation matrix contains the additional information needed to derive the full covariance matrix of the predicted soil hydraulic parameters.
Step 4:. Derive the covariance matrix of the predicted soil hydraulic parameters. Derive the full covariance matrix of the predicted soil hydraulic parameters p by scaling the correlation matrix R with the corresponding standard deviations of the predicted parameters u: The scaling ensures that the diagonal entries of p correspond to the standard deviations of the predicted parameters while the offdiagonal terms reflect the correlation among the soil hydraulic parameters derived from the random sample.
Note that the resulting covariance matrix p only accounts for parameter uncertainty and correlation induced by the pedotransfer function. It does not account for uncertainty in the predicted parameters due to uncertainty in the input variables. The correlation matrix R derived with the Monte Carlo approach presented above is shown in Table 2.

Defining the prior distributions of the soil hydraulic parameters
We tested three different formulations of the prior distribution of the soil hydraulic parameters. These prior distributions incorporate to different extends the information derived from ROSETTA.
Prior 1: Multivariate uniform distribution. The first prior is multivariate uniform p(x 2 )∼U (a x 2 ,b x 2 ), with lower and upper bounds a x 2 and b x 2 , respectively, that jointly define the feasible parameter space. We calculated the lower and upper bounds as p ± 4 u ( Table 1).

Prior 3: Multivariate normal distribution with correlation.
The third prior is also multivariate normal but considers correlation among the soil hydraulic parameters using the full covariance matrix x 2 = p derived using the Monte Carlo approach.
Figure 2 illustrates how the three prior distributions affect the prior uncertainty of the water retention and hydraulic conductivity functions. Unsurprisingly, the largest prior uncertainties are found when using a uniform prior distribution (Prior 1, red line). This uncertainty is substantially reduced when prior information about the mean and standard deviation is considered (Prior 2, blue line). Even though parameter correlation is ignored in this prior distribution, the VGM parameters are much better constrained now, and this results in much smaller 95 % confidence intervals compared to the uniform prior case. When parameter correlation is considered in the prior distribution (Prior 3, grey area), the confidence interval around the inflection point of the water retention function shrinks symmetrically, but the uncertainty in the wet and dry ranges is hardly affected. A qualitatively similar picture results for the hydraulic conductivity function.

Generation of synthetic data
To test the effectiveness and robustness of the Bayesian approach, we used two synthetic data sets. In the first case, we created a time series of soil water content observations using the HYDRUS-1D model with ROSETTA predicted values of the soil hydraulic parameters (Table 2) and h LB = −150 cm. The upper boundary conditions as well as the observation dates were the same as for the real data set. A normally distributed error was added to the simulated soil water contents to represent the combined effect of model structural, boundary condition, and observational error. The magnitude of this random error was similar to the root mean square error of the best HYDRUS-1D fit to the real soil water content data (Sect. 3.2). This synthetic data set was then used with the three prior distributions to inversely estimate the VGM parameters and h LB , providing a test of the effectiveness of using prior information about the soil hydraulic parameters. We considered a second test case in which we have biased prior information. This test was specially designed to establish whether we can infer the appropriate values of the soil hydraulic parameters even though the prior distribution is biased, providing a test of the robustness of the Bayesian approach. There are essentially two ways in which we can generate this bias. One way is to maintain the synthetic time series of soil water content observations and to corrupt the three prior pdfs. A simpler approach followed herein is to leave the prior pdfs untouched, but to create a second synthetic time series of soil water content observations using soil hydraulic parameters that differ from the values predicted by ROSETTA and used to create the prior pdfs. These parameters were selected by drawing a large random sample from Prior 3 and then purposely picking a realization with very low prior probability. Model settings and input variables were the same as in the previous case, and the same random error was added to the synthetic observations. We then The blue lines represent the prior distributions. The red circles denote the 2.5, 10, 20, 30, 40, 50, 60, 70, 80, 90, and 97.5 % quantiles of the posterior distributions, respectively. The grey lines mark the parameter values used to generate the data. The x-axes cover the bounds listed in Table 1. The blue lines represent the prior distributions. The red circles denote the 2.5, 10, 20, 30, 40, 50, 60, 70, 80, 90, and 97.5 % quantiles of the posterior distributions, respectively. The grey lines mark the parameter values used to generate the data. The x-axes cover the bounds listed in Table 1.  Table 1. estimated the unknown model parameters from this second time series using the three (now biased) prior distributions.

Synthetic data
The results for the synthetic soil water content data using unbiased prior distributions are depicted in Fig. 3. This plot shows the cumulative prior (blue line) and posterior (red circles) distributions of the estimated parameters corresponding to each of the three prior distributions. When using the uniform prior (Prior 1, Fig. 3a), the water retention parameters (θ r , θ s , α, and n) were not identifiable, with posterior distributions that extended over the entire prior defined ranges. The posterior distributions of the two additional soil hydraulic parameters (K s and L) on the contrary differed markedly from their marginal prior distributions. Seemingly, the observational data contained sufficient information to constrain these two parameters. The pressure head at the lower boundary (h LB ) was not warranted by calibration against the synthetic soil water content data. Very similar findings were observed with the informative prior that neglects parameter correlation (Prior 2, Fig. 3b). Note that h LB was somewhat better constrained in this case but still demonstrated considerable uncertainty. This is a nice illustration of parameter interdependence. Even though the prior of h LB was the same for the three prior distributions, this parameter became better identifiable when the VGM parameters were more constrained in their prior pdf. Finally, the informative prior that considers parameter correlation (Prior 3, Fig. 3c) substantially improved the results. All the soil hydraulic parameters as well as the pressure head at the lower boundary suddenly became well identifiable by calibration against the observational data. Note that in this case the true parameter values (grey lines) were always located within the 95 % posterior confidence intervals and that these intervals were reasonably small. In this study, we refer to an estimate that satisfies both of these preconditions as an accurate and precise estimate.
The results for the biased prior distributions are shown in Fig. 4. They were qualitatively similar to those observed  10, 20, 30, 40, 50, 60, 70, 80, 90, and 97.5 % quantiles of the posterior distributions, respectively. The x-axes cover the bounds listed in Table 1. previously for the unbiased priors case. Again, when using prior distributions that neglect correlation among the soil hydraulic parameters (Prior 1 and Prior 2), accurate and precise estimation of the water retention parameters was not possible. Including prior information about parameter correlation (Prior 3) significantly improved parameter identifiability, enabling accurate and precise estimates of all parameters. Note that this was possible despite the bias in the prior information. This is illustrated in more detail in Fig. 5, which presents pairwise scatter plots of the posterior sample (grey dots) and associated 95 % confidence ellipses of the prior (blue line) and posterior (red line) distribution. These ellipses were calculated based on the assumption of a bivariate normal distribution of the respective parameters. Although the prior distribution of the VGM parameters was biased and therefore assigned very low prior probability to the actual parameter values (orange squares) used to generate the synthetic time series, the 95 % posterior confidence ellipses encompassed these values.

Real data
Posterior distributions of the estimated parameters for the real soil water content data are illustrated in Fig. 6. These results resemble the findings presented in the previous section for the synthetic data. The use of the uniform prior (Prior 1) and the multivariate normal prior that neglects parameter correlation (Prior 2) did not warrant accurate and precise identification of all estimated parameters. Consideration of parameter correlation in the prior distribution (Prior 3) substantially improved the outcome of the Bayesian inference scheme.
To verify whether the prior distribution of the parameters is consistent with the information from the soil water content data, please consider Fig. 7 which compares the 95 % confidence ellipses of the prior distribution (blue line) with those derived from the posterior sample (red line). This figure highlights several important observations. First, the confidence ellipses of the posterior distribution are much smaller than their respective counterparts of the prior distribution. Second, the prior distribution was apparently unbiased and consistent with the posterior sample. And finally, the prior and posterior distributions exhibit a very similar correlation structure, particularly for the highly correlated parameters. Figure 8 compares the observed and simulated soil water dynamics. The dark grey region represents the 95 % prediction uncertainty intervals associated with the posterior parameter uncertainty, whereas the light grey region depicts the total predictive uncertainty. Details on how to compute these uncertainty intervals can be found in Schoups and Vrugt (2010) and so will not be repeated herein. For completeness, the top panel plots the observed rainfall hyetograph (blue bars) and potential evaporation (red bars). The HYDRUS-1D model matched the in situ observations (red circles) very well, with a corresponding root mean square error of 0.009 cm 3 cm −3 , and model efficiency (Nash and Sutcliffe, 1970) of 0.87. These two measures of goodness of fit were calculated for the HYDRUS-1D simulation that best described the observational data, that is, using the parameter values that had maximum posterior density (θ r = 0.066, θ s = 0.445, α = 0.0048, n = 1.68, K s = 0.074, L = 0.63, and  Table 1. h LB = −128). Note that two observations (day of year 190 and 238) fall outside the 95 % prediction uncertainty bounds. This, however, is in good agreement with statistical expectation (2 out of 29 observations correspond to approximately 7 %). To convey an impression of the spatial variability of soil water content within the 50 m × 50 m plot, we added bars to the mean values that span the 2.5 % and 97.5 % quantiles of the spatially distributed observations. Figure 9 plots the 95 % uncertainty bounds of the soil hydraulic functions corresponding to the prior and posterior distributions. We also plotted the observed soil water content data. Note that the in situ data exhibited relative small variability and covered only a limited range of soil water states. This explains, at least in part, why in situ observations of soil water dynamics contain insufficient information to estimate accurately and precisely all VGM parameters. Note also that some of the data fall outside the posterior uncertainty bounds. The reason for this is that these bounds show the uncertainty in the soil hydraulic functions due to uncertainty in the soil hydraulic parameters only. They do not include uncertainty due to observational and model errors.
The likelihood function used in this study was based on assumptions of uncorrelated, homoscedastic, and normally distributed residuals. If any of these assumptions is violated then the posterior distribution and corresponding prediction uncertainty intervals are subject to error and should be revisited. Good statistical practice therefore constitutes checking whether the underlying assumptions of the likelihood model have been met (e.g. Schoups and Vrugt, 2010). To assess the validity of these assumptions we conducted three diagnostic tests, and the results of this are plotted in Fig. 10. The top panel measures the correlation among the residuals by plotting the autocorrelation function. The autocorrelation at given lag (red circles) remained within the theoretical 95 % significance interval of a time series of uncorrelated residuals (blue lines), indicating that the residuals are uncorrelated. The middle panel tests whether the magnitude of the residuals depends on the magnitude of the soil water content observation. The residuals appear homoscedastic, that is, independent of the magnitude of the observational data. Finally, the bottom panel presents a quantile-quantile plot and explores whether the residuals follow a normal distribution. Except 95% parameter uncertainty bounds 95% total model uncertainty bounds Fig. 8. 95 % predictive uncertainty bounds using real soil water content data: (a) meteorological forcing and (b) observed and predicted soil water content. The dark grey region represents the predictive uncertainty due to uncertainty in the estimated parameters. The light grey region denotes the predictive uncertainty due to the combined effects of parameter, model, boundary condition and calibration data uncertainty. The red circles mark the mean values of the spatially distributed soil water content observations, while the red bars span the 2.5 % and 97.5 % quantiles of these observations. Note that we depict daily values of potential evaporation and precipitation, but the HYDRUS-1D model was actually run with hourly values of these forcing variables.
for the two large residuals discussed above, the quantiles of the residuals were in good agreement with this assumption. We therefore concluded that the underlying assumptions of the likelihood model were met and that the posterior distribution and corresponding predictive uncertainty intervals were adequate.

Discussion
The results presented in this paper clearly indicate that the in situ observations of soil water content did not contain sufficient information to warrant an accurate and precise estimation of all parameters of interest. This finding is not new but has been reported previously (e.g. Jacques et al., 2002;Ritter et al., 2003). Naturally occurring boundary conditions display insufficient variability to result in a wide range of soil water states, which is a prerequisite to successfully estimating the VGM parameters (Vrugt et al., 2001(Vrugt et al., , 2002. The necessary information for some of the soil hydraulic parameters simply appears beyond the range of the actual soil water content observations, and these parameters are therefore difficult to constrain. There is different ways in which the information content of our data could have been increased. One approach is to increase measurement frequency and to measure directly after rainfall events (Fig. 8). This would have resulted in a larger range of observed soil water states. Another possibility is to consider the presence of vegetation, and consequently, root water uptake. This would have extended the range of observed soil water states to the dry end. The gain of information associated with this, however, would come at a cost. The simulation of root water uptake as a function of time and depth requires specification of additional parameters, which would need to be estimated simultaneously with the other parameters, introducing additional parameter and model uncertainty (e.g. Ines and Mohanty, 2008;Wollschläger et al., 2009).
The use of prior information about the soil hydraulic parameters substantially improved parameter identifiability. To achieve this improvement, however, it was necessary to include information about parameter correlation in the prior distribution. Using this additional information, our results with synthetic data demonstrated that the Bayesian approach is effective and robust, even in case of biased prior information. It is noteworthy that the bias we used in this study to test the robustness of the approach was rather moderate. Clearly, another situation occurs if the magnitude of the bias is actually so large that the prior information becomes incompatible with the information contained in the observational data, as illustrated and discussed in Hou and Rubin (2005). In this situation, any approach of using prior information is likely to fail. In practice, however, this failure becomes evident from inspection of the posterior distribution, with the maximum posterior density of at least some of the parameters located at the bounds of the feasible parameter space (Hou and Rubin, 2005).
Another option to alleviate problems with parameter nonidentifiability is to reduce the dimensionality of the model calibration problem and fix some of the soil hydraulic parameters at some a priori defined value (e.g. Jacques et al., 2002;Ritter et al., 2003). This approach is practical but may impair the ability of the soil hydraulic model to accurately describe the experimental data, in particular if parameters are fixed at nonoptimal values. Due to parameter correlation, fixing single parameters at nonoptimal values will likely corrupt the estimates of the remaining parameters. Moreover, the decision which parameter to fix is often rather arbitrary, without consideration of the actual information content of the data (but see, e.g. Ritter et al., 2003;Mertens et al., 2005, for some exceptions). This might results in fixing a parameter whose value is actually well defined by the calibration data. An example of this is the L parameter in the VGM model that is often fixed at some value taken from the literature because it is deemed unimportant or insensitive (e.g. Abbaspour et al., 2000;Wollschläger et al., 2009). In contrast, the results presented in this study demonstrated that the in situ observations contained valuable information that helped to substantially constrain this parameter (Fig. 6). If we still decide to fix some of the soil hydraulic parameters, it remains typically difficult to choose appropriate values. This is particularly true for θ r , θ s , and K s . These parameters are generally poorly defined by direct measurements and should therefore be considered as fitting parameters (e.g. van Genuchten and Nielsen, 1985). The Bayesian approach presented in this study avoids many of the problems associated with fixing parameters.
A recent review by Vereecken et al. (2010) discusses the strengths and limitations of existing pedotransfer functions to make accurate and reliable predictions of the parameters in the VGM model. They also make suggestions on how to improve the predictive capabilities of future pedotransfer functions. We like to add to their suggestions to include information on parameter correlation in the pedotransfer function output. The Monte Carlo approach we presented might prove useful in this regard. It is generally applicable and easy to implement. As shown in this study, information on correlation among the soil hydraulic parameters is highly beneficial in Bayesian inverse modelling of in situ data, but equally useful in many other applications such as stochastic modelling of soil hydraulic properties and soil water flow (e.g. Mishra and Parker, 1989;Mallants et al., 1996).

Summary and conclusions
Many contributions to the soil hydrological literature have demonstrated the limited information content of in situ measurements of soil water state variables under natural boundary conditions to estimate the soil hydraulic properties. A general approach, which has yet received very little attention in the soil hydrological literature, is to use an informative prior distribution of the soil hydraulic parameters and to combine this distribution with the in situ observations using Bayes' theorem. In this paper, we investigated to which degree prior information about the soil hydraulic parameters can help improve parameter identifiability in inverse modelling of in situ soil water dynamics under natural boundary conditions. We used percentages of sand, silt, and clay as input variables to the ROSETTA pedotransfer function that predicts the soil hydraulic parameters. Textural data constitute basic soil information that is readily available in most vadose zone studies. In addition to the standard ROSETTA prediction that provides the mean values and standard deviations of the predicted parameters, we tested a Monte Carlo approach to derive the correlation structure of the predicted parameters. It was one objective of this study to explore the value of this additional information on parameter correlation in Bayesian inverse modelling. Another objective was to test the robustness of the Bayesian approach in case of biased prior information. We formulated three different prior distributions that incorporate to different extents the prior information derived with ROSETTA. We illustrated our approach using synthetic and real-world observations of in situ soil water dynamics under natural boundary conditions.
The results of this study demonstrated that prior information about the soil hydraulic parameters significantly improved parameter identifiability in Bayesian inverse modelling of in situ soil water dynamics. The results also indicated that the Bayesian approach was effective and robust under the conditions tested in this study, even in case of biased prior information. For the Bayesian approach to be effective and robust, however, it was essential to incorporate information about the correlation structure of the soil hydraulic parameters in the prior distribution. The proposed algorithm to derive this additional information proved useful, yielding reasonable estimates of the correlation coefficients of the ROSETTA predicted parameters. Using the so derived full covariance matrix in Bayesian inverse modelling enabled us to successfully calibrate a one-dimensional effective vadose zone model using real-world data with limited information content.