Hydraulic and transport parameter assessment using column infiltration experiments

Abstract

The quality of statistical calibration of hydraulic and transport soil properties is studied for infiltration experiments in which, over a given period, tracer-contaminated water is injected into an hypothetical column filled with a homogeneous soil. The saturated hydraulic conductivity, the saturated and residual water contents, the Mualem–van Genuchten shape parameters and the longitudinal dispersivity are estimated in a Bayesian framework using the Markov chain Monte Carlo (MCMC) sampler. The impact of the kind of measurement sets (water content, pressure inside the column, cumulative outflow and outlet solute concentration) and that of the solute injection duration is investigated by analyzing the calibrated model parameters and their confidence intervals for different scenarios. The results show that the injection period has a significant effect on the quality of the estimation, in particular, on the posterior uncertainty range of the parameters. All hydraulic and transport parameters of the investigated soil can be well estimated from the experiment using only the outlet concentration and cumulative outflow, which are measured non-intrusively. An improvement of the identifiability of the hydraulic parameters is observed when the pressure data from measurements taken inside the column are also considered in the inversion.

Introduction

The soil parameters that influence water flow and contaminant transport in unsaturated zones are not generally known a priori and have to be estimated by fitting model responses to observed data. The unsaturated soil hydraulic parameters can be (more or less accurately) estimated from dynamic flow experiments (e.g., Hopmans et al., 2002; Vrugt et al., 2003a; Durner and Iden, 2011; Younes et al., 2013). Several authors have investigated different types of transient experiments and boundary conditions suited for a reliable estimation of soil hydraulic properties (e.g., van Dam et al., 1994; Šimůnek and van Genuchten, 1997; Inoue et al., 1998; Durner et al., 1999). Soil hydraulic properties are often estimated using inversion of one-step (Kool et al., 1985; van Dam et al., 1992) or multistep (Eching et al., 1994; van Dam et al., 1994) outflow experiments or controlled infiltration experiments (Hudson et al., 1996).

Kool et al. (1985) and Kool and Parker (1988) suggested that the transient experiments should cover a wide range of water contents to obtain a reliable estimation of the parameters. Van Dam et al. (1994) have shown that more reliable parameter estimates are obtained by increasing the pneumatic pressure in several steps instead of a single step. The multistep outflow experiments are the most popular laboratory methods (e.g., Eching and Hopmans, 1993; Eching et al., 1994; van Dam et al., 1994; Hopmans et al., 2002). However, their application is limited by expensive measurement equipment (Nasta et al., 2011).

Infiltration experiments have been investigated by Mishra and Parker (1989) to study the reliability of hydraulic- and transport-estimated parameters for a soil column of 200 cm using measurements of water content, concentration and water pressure inside the column. They showed that the simultaneous estimation of hydraulic and transport properties yields smaller estimation errors for model parameters than the sequential inversion of hydraulic properties from the water content and/or pressure head followed by the inversion of transport properties from concentration data (Mishra and Parker, 1989).

Inoue et al. (2000) performed infiltration experiments using a soil column of 30 cm. Pressure head and solute concentration were measured at different locations. A constant infiltration rate was applied to the soil surface and a balance was used to measure the cumulative outflow. They showed that both hydraulic and transport parameters can be assessed by the combination of flow and transport experiments.

Furthermore, infiltration experiments were often conducted in lysimeters for pesticide leaching studies. Indeed, lysimeter experiments are generally used to assess the leaching risks of pesticides using soil columns of around 1.2 m depth which is the standard scale for these types of experiments (Mertens et al., 2009; Kahl et al., 2015). Before performing the column leaching experiment, several infiltration–outflow experiments are often realized to estimate the soil hydraulic parameters (Kahl et al., 2015; Dusek et al., 2015).

The key objective of the present study is to evaluate the reliability of different experimental protocols for estimating hydraulic and transport parameters and their associated uncertainties for column experiments. We consider the flow and the transport of an inert solute injected into a hypothetical column filled with a homogeneous sandy clay loam soil. We assume that flow can be modeled by the Richards' equation (RE) and that the solute transport can be simulated by the classical advection–dispersion model. Furthermore, the Mualem and van Genuchten (MvG) models (Mualem, 1976; van Genuchten, 1980) are chosen to describe the retention curve and to relate the hydraulic conductivity of the unsaturated soil to the water content. The estimation of the flow and transport parameters through flow–transport model inversion is investigated for two injection periods of the solute and different data measurement scenarios.

Inverse modeling is often performed using local search algorithms such as the Levenberg–Marquardt algorithm (Marquardt, 1963). The latter is computationally efficient to evaluate the optimal parameter set (Gallagher and Doherty, 2007). Besides, the degree of uncertainty in the estimated parameters, expressed by their confidence intervals, is often calculated using a first-order approximation of the model near its minimum (Carrera and Neuman, 1986; Kool and Parker, 1988). However, as stated by Vrugt and Bouten (2002), parameter interdependence and model nonlinearity occurring in hydrologic models may violate the use of this first approximation to obtain accurate confidence intervals of each parameter. Therefore, in this work, the estimation of hydraulic and transport parameters is performed in a Bayesian framework using the Markov chain Monte Carlo (MCMC) sampler (Vrugt and Bouten, 2002; Vrugt et al., 2008). Unlike classical parameter optimization algorithms, the MCMC approach generates sets of parameter values randomly sampled from the posterior joint probability distributions, which are useful to assess the quality of the estimation. The MCMC samples can be used to summarize parameter uncertainties and to perform predictive uncertainty (Ades and Lu, 2003).

Hypothetical infiltration experiments are considered for a column of 120 cm depth, initially under hydrostatic conditions, free of solute and filled with a homogeneous sandy clay loam soil. Continuous flow and solute injection are performed during a time period $T_{\mathrm{inj}}$ at the top of the column and with a zero-pressure head at the bottom. The unknown parameters for the water flow are the hydraulic parameters: $k_{\mathrm{s}}$ (L T${}^{-\mathrm{1}}$), the saturated hydraulic conductivity; ${\mathit{\theta }}_{\mathrm{s}}$ (L${}^{\mathrm{3}}$ L${}^{-\mathrm{3}}$), the saturated water content; ${\mathit{\theta }}_{\mathrm{r}}$ (L${}^{\mathrm{3}}$ L${}^{-\mathrm{3}}$), the residual water content; and $\mathit{\alpha }$ (L${}^{-\mathrm{1}}$) and $n$ ($-$), the MvG shape parameters. The only unknown parameter of the tracer transport is the longitudinal dispersivity, $a_{\mathrm{L}}$(L).

Several scenarios corresponding to different sets of measurements are investigated to address the following questions:

• Can we obtain an appropriate estimation of all flow and transport parameters from tracer-infiltration experiments, even though a limited range of water contents is covered (only moderately dry conditions are obtained because of gravity drainage conditions prescribed at the bottom of the soil column)?

• What is the optimal set of measurements for the estimation of all the parameters? Can we use only non-intrusive measurements (cumulative outflow and concentration breakthrough curve) or are intrusive measurements of pressure heads and/or water contents inside the column unavoidable?

• Is there an optimal design for the tracer injection?

For this purpose, synthetic scenarios are considered in the sequel in which data from numerical simulations are used to avoid the uncontrolled noise of experiments that could bias the conclusions.

The paper is organized as follows. The mathematical models describing flow and transport in the unsaturated zone are detailed in Sect. 2. Section 3 describes the MCMC Bayesian parameter estimation procedure used in the DREAM${}_{\left(\mathrm{ZS}\right)}$ sampler. Section 4 presents the different investigated scenarios and discusses the results of the calibration in terms of mean parameter values and uncertainty ranges for each scenario. Conclusions are given in Sect. 5.

Unsaturated flow–transport model

We consider a uniform soil profile in the column and an injection of a solute tracer such as bromide, as described in Mertens et al. (2009). The unsaturated water flow in the vertical soil column is modeled with the one-dimensional pressure head form of the RE: $$\left\{\begin{array}{l}\left(c,\left(h\right),+,S_{\mathrm{s}},{\displaystyle \frac{\mathit{\theta }}{{\mathit{\theta }}_S}}\right){\displaystyle \frac{\partial h}{\partial t}}={\displaystyle \frac{\partial q}{\partial z}}\\ q=K\left(h\right)\left({\displaystyle \frac{\partial h}{\partial z}},-,\mathrm{1}\right)\end{array}\right.,$$ where $h$ (L) is the pressure head; $q$ (L T${}^{-\mathrm{1}}$) is the Darcy velocity; $z$ (L) is the depth, measured as positive in the downward direction; $S_{\mathrm{s}}$ ($-$) is the specific storage; $\mathit{\theta }$ and ${\mathit{\theta }}_{\mathrm{s}}$ (L${}^{\mathrm{3}}$ L${}^{-\mathrm{3}}$) are the actual and saturated water contents, respectively; $c\left(h\right)$ (L${}^{-\mathrm{1}}$) is the specific moisture capacity; and $K\left(h\right)$ (L T${}^{-\mathrm{1}}$) is the hydraulic conductivity. The latter two parameters are both functions of the pressure head. In this study, the relations between the pressure head, conductivity and water content are described by the following standard models of Mualem (1976) and van Genuchten (1980): $$\begin{array}{l}{S}_{\mathrm{e}}\left(h\right)={\displaystyle \frac{\mathit{\theta }\left(h\right)-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}}=\left\{\begin{array}{l}{\displaystyle \frac{\mathrm{1}}{{\left(\mathrm{1},+,{\left|\mathit{\alpha },h\right|}^n\right)}^m}}h<\mathrm{0}\\ \mathrm{1}h\ge \mathrm{0}\end{array}\right.\\ K\left(S_{\mathrm{e}}\right)=K_{\mathrm{s}}{S}_{\mathrm{e}}^{\mathrm{1}/\mathrm{2}}{\left[\mathrm{1},-,{\left(\mathrm{1},-,S_{\mathrm{e}}^{\mathrm{1}/m}\right)}^m\right]}^{\mathrm{2}}\end{array},$$ where $S_{\mathrm{e}}$ ($-$) is the effective saturation, ${\mathit{\theta }}_{\mathrm{r}}$ (L${}^{\mathrm{3}}$ L${}^{-\mathrm{3}}$) is the residual water content, $K_{\mathrm{s}}$ (L T${}^{-\mathrm{1}}$) is the saturated hydraulic conductivity, and $m=\mathrm{1}-\mathrm{1}/n$, $\mathit{\alpha }$ (L${}^{-\mathrm{1}}$) and $n$ ($-$) are the MvG shape parameters.

The tracer transport is governed by the following convection–dispersion equation: $${\displaystyle \frac{\partial \left(\mathit{\theta },C\right)}{\partial t}}+{\displaystyle \frac{\partial \left(q,C\right)}{\partial z}}-{\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\theta },D,{\displaystyle \frac{\partial C}{\partial z}}\right)=\mathrm{0},$$ where $C$ (M L${}^{-\mathrm{3}}$) is the concentration of the tracer, $D$ (L${}^{\mathrm{2}}$ T${}^{-\mathrm{1}}$) is the dispersion coefficient in which $D=a_{l}q+d_{\mathrm{m}}$ and $a_l$ (L) is the dispersivity coefficient of the soil and $d_{\mathrm{m}}$ (L${}^{\mathrm{2}}$ T${}^{-\mathrm{1}}$) is the molecular diffusion coefficient, which is set as 1.04 10${}^{-\mathrm{4}}$ cm${}^{\mathrm{2}}$ min${}^{-\mathrm{1}}$.

The transport Eq. (3) is coupled with the flow Eq. (1) by the water content $\mathit{\theta }$ and Darcy's velocity $q$. The initial conditions are as follows: a hydrostatic pressure distribution with zero-pressure head at the bottom of the column $\left(z,=,L\right)$ and a solute concentration of zero inside the whole column. An infiltration with a flux $q_{\mathrm{inj}}$ of contaminated water with a concentration $C_{\mathrm{inj}}$ is then applied at the upper boundary condition ($z=$ 0) during a period $T_{\mathrm{inj}}$. Hence, the boundary conditions at the top of the column can be expressed as $$\begin{array}{cc}& {\displaystyle }& {\displaystyle }\mathrm{for}\mathrm{0}<t\le T_{\mathrm{inj}}\left\{\begin{array}{l}K\left({\displaystyle \frac{\partial h}{\partial z}},-,\mathrm{1}\right)=q_{\mathrm{inj}}\\ \mathit{\theta }D{\displaystyle \frac{\partial C}{\partial z}}+qC=q_{\mathrm{inj}}{C}_{\mathrm{inj}}\end{array}\right.{\displaystyle }& {\displaystyle }\mathrm{for}t>T_{\mathrm{inj}}\left\{\begin{array}{l}K\left({\displaystyle \frac{\partial h}{\partial z}},-,\mathrm{1}\right)=\mathrm{0}\\ C_{\mathrm{inj}}=\mathrm{0}\end{array}\right..\end{array}$$ A zero-pressure head is maintained at the lower boundary $\left(z,=,L\right)$ of the column and a zero-concentration gradient is used as the lower boundary condition for the solute transport, namely $${\left(h\right)}_{z=l}=\mathrm{0}{\left({\displaystyle \frac{\partial C}{\partial z}}\right)}_{z=l}=\mathrm{0}.$$ In the sequel, the infiltration rate and the injected solute concentration are $q_{\mathrm{inj}}=\mathrm{0.015}$ cm min${}^{-\mathrm{1}}$ and $C_{\mathrm{inj}}=\mathrm{1}$ g cm${}^{-\mathrm{3}}$, respectively. The system (Eqs. 1–5) is solved using the standard finite difference method for both flow and transport spatial discretization. A uniform mesh of 600 cells is employed. Temporal discretization is performed with the high-order method of lines (MOL) (e.g., Miller et al., 1998; Tocci et al., 1997; Younes et al., 2009; Fahs et al., 2011). Error checking, robustness, order selection and adaptive time step features, available in sophisticated solvers, are applied to the time integration of partial differential equations (Tocci et al., 1997). The MOL has been successfully used to solve RE in many studies (e.g., Farthing et al., 2003; Miller et al., 2006; Li et al., 2007; Fahs et al., 2009). Details on the use of the MOL for solving RE are described in Fahs et al. (2009).

Bayesian parameter estimation

The vector of unknown parameters that has to be identified by model calibration is $\mathit{\xi }=\left(k_{\mathrm{s}},,,{\mathit{\theta }}_{\mathrm{s}},,,{\mathit{\theta }}_{\mathrm{r}},,,\mathit{\alpha },,,n,,,a_{\mathrm{L}}\right)$. To analyze the performance of the model calibration procedures, a reference solution is generated by simulating the flow–transport problem (Eqs. 1–5) using the following parameter values (corresponding to a sandy clay loam soil): $k_{\mathrm{s}}=$ 50 cm day${}^{-\mathrm{1}}$, ${\mathit{\theta }}_{\mathrm{s}}=\mathrm{0.43}$, ${\mathit{\theta }}_{\mathrm{r}}=\mathrm{0.09}$, $\mathit{\alpha }=$ 0.04 cm${}^{-\mathrm{1}}$, $n=\mathrm{1.4}$ and $a_l=\mathrm{0.2}$ cm. Four types of variables are extracted from the results of the simulation: the pressure head and water content 5 cm below the top of the column, the cumulative outflow and the solute breakthrough concentration at the outflow of the column. These four data series are modified by adding a normally distributed white noise using the following standard deviations: ${\mathit{\sigma }}_{\mathrm{h}}=$ 1 cm for the pressure head, ${\mathit{\sigma }}_{\mathit{\theta }}=\mathrm{0.02}$ for the water content, ${\mathit{\sigma }}_Q=\mathrm{0.1}$ cm for the cumulative outflow and ${\mathit{\sigma }}_{\mathrm{C}}=\mathrm{0.01}$ g cm${}^{-\mathrm{3}}$ for the exit concentration. These perturbations mimic measurement errors and the resulting values of water pressure, water content, cumulative outflow and solute breakthrough concentration are considered as measurements in the following.

The flow–transport model is used to analyze the effects of different measurement sets on parameter identification. For this purpose, we adopt a Bayesian approach that involves the parameter joint posterior distribution (Vrugt et al., 2008). The latter is assessed with the DREAM${}_{\left(\mathrm{ZS}\right)}$ MCMC sampler (Laloy and Vrugt, 2012). This software generates random sequences of parameter sets that asymptotically converge toward the target joint posterior distribution (Gelman et al., 1997). Thus, if the number of runs is sufficiently high, the generated samples can be used to estimate the statistical measures of the posterior distribution, such as the mean and variance, among other measures.

The Bayes theorem states that the probability density function of the model parameters conditioned onto data can be expressed as $$p\left(\mathit{\xi },\mathrm{|},y_{\mathrm{mes}}\right)\propto p\left(y_{\mathrm{mes}},\mathrm{|},\mathit{\xi }\right)p\left(\mathit{\xi }\right),$$ where $p\left(\mathit{\xi },\mathrm{|},y_{\mathrm{mes}}\right)$ is the likelihood function measuring how well the model fits the observations $y_{\mathrm{mes}}$, and $p\left(\mathit{\xi }\right)$ is the prior information about the parameter before the observations are made. Independent uniform priors within the ranges reported in Table 1 are chosen. In this work, a Gaussian distribution defines the likelihood function because the observations are simulated and corrupted with Gaussian errors. Hence, the parameter posterior distribution is expressed as $$p\left(\mathit{\xi },\mathrm{|},y_{\mathrm{mes}}\right)\propto \exp \left(-,{\displaystyle \frac{{\mathrm{SS}}_{\mathrm{h}}\left(\mathit{\xi }\right)}{\mathrm{2}{\mathit{\sigma }}_{\mathrm{h}}^{\mathrm{2}}}},-,{\displaystyle \frac{{\mathrm{SS}}_{\mathit{\theta }}\left(\mathit{\xi }\right)}{\mathrm{2}{\mathit{\sigma }}_{\mathit{\theta }}^{\mathrm{2}}}},-,{\displaystyle \frac{{\mathrm{SS}}_Q\left(\mathit{\xi }\right)}{\mathrm{2}{\mathit{\sigma }}_Q^{\mathrm{2}}}},-,{\displaystyle \frac{{\mathrm{SS}}_{\mathrm{C}}\left(\mathit{\xi }\right)}{\mathrm{2}{\mathit{\sigma }}_{\mathrm{C}}^{\mathrm{2}}}}\right),$$ where SS${}_{\mathrm{h}}\left(\mathit{\xi }\right)$, SS${}_{\mathit{\theta }}\left(\mathit{\xi }\right)$, SS${}_Q\left(\mathit{\xi }\right)$ and SS${}_{\mathrm{C}}\left(\mathit{\xi }\right)$ are the sums of the squared differences between the observed and modeled data of the pressure head, water content, cumulative outflow and output concentration, respectively. For instance, SS${}_{\mathrm{h}}\left(\mathit{\xi }\right)={\sum }_{k=\mathrm{1}}^{N_{\mathrm{h}}}{\left(h_{\mathrm{mes}}^{\left(k\right)},-,h_{\mathrm{mod}}^{\left(k\right)},\left(\mathit{\xi }\right)\right)}^{\mathrm{2}}$, which includes the observed $h_{\mathrm{mes}}^{\left(k\right)}$ and predicted $h_{\mathrm{mod}}^{\left(k\right)}$ pressure heads at time $t_k$ for the number of pressure head observations $N_{\mathrm{h}}$.

Prior lower and upper bounds of the uncertainty parameters and reference values.

Parameters	Lower	Upper	Reference
	bounds	bounds	values
$k_{\mathrm{s}}$ (cm min${}^{-\mathrm{1}}$)	0.025	0.1	0.0347
${\mathit{\theta }}_{\mathrm{s}}$ ($-$)	0.3	0.5	0.43
${\mathit{\theta }}_{\mathrm{r}}$ ($-$)	0.05	0.2	0.09
$\mathit{\alpha }$ (cm${}^{-\mathrm{1}}$)	0.01	0.3	0.04
$n$ ($-$)	1.2	5	1.4
$a_l$ (cm)	0.05	0.6	0.2

Bayesian parameter estimation is performed hereafter with the DREAM${}_{\left(\mathrm{ZS}\right)}$ software (Laloy and Vrugt, 2012), which is an efficient MCMC sampler. DREAM${}_{\left(\mathrm{ZS}\right)}$ computes multiple sub-chains in parallel to thoroughly explore the parameter space. Archives of the states of the sub-chains are stored and used to allow a strong reduction of the “burn-in” period in which the sampler generates individuals with poor performances. Taking the last 25 % of individuals of the MCMC (when the chains have converged) yields multiple sets of parameters, $\mathit{\xi }$, that adequately fit the model onto observations. These sets are then used to estimate the updated parameter distributions, the pairwise parameter correlations and the uncertainty of the model predictions. As suggested in Vrugt et al. (2003b), we consider that the posterior distribution is stationary if the Gelman and Ruban (1992) criterion is $\le \mathrm{1.2}$.

Results and discussion

In this section, the identifiability of the parameters is investigated for seven different scenarios of measurement sets (Table 1). In the first scenario, only measured pressure heads and cumulative outflow are used for the calibration. Scenarios 2 to 5 investigate the benefit of adding measured water contents and/or solute outlet concentrations to pressure heads and outflow. The last scenarios (6, 7) investigate the use of measured cumulative outflow and concentration breakthrough at the column outflow because these measurements do not require intrusive techniques. Scenarios 5 to 7 investigate the effects of solute injection duration on the identifiability of the parameters as well.

In all cases, the MCMC sampler was run with three simultaneous chains for a total number of 50 000 runs. Depending on the scenario, the MCMC required between 5000 and 20 000 model runs to reach convergence and was terminated after 30 000 runs. The last 25 % of the runs that adequately fit the model onto observations are used to estimate the updated probability density function (pdf).

The data sets for parameter estimation

The data sets obtained from solving the flow–transport problems (Eqs. 1–5) using the parameters given in Sect. 2 are shown in Fig. 1. The pressure head at 5 cm from the top of the column (Fig. 1a) increases from its initial hydrostatic negative value ($-$115 cm) and reaches a plateau ($-$1.75 cm) in less than 100 min during the injection period. After the injection is finished, it progressively decreases due to the drainage caused by the gravity effect. A similar behavior is observed for the water content at the same location (Fig. 1b), where the value of the plateau is close to the saturation value. The cumulative outflow (Fig. 1c) starts to increase at approximately 1000 min after the beginning of the injection. It shows an almost linear behavior until 5500 min. It then slowly increases with an asymptotic behavior due to the natural drainage after the end of the injection period. Fig. 1d displays the water saturation as a function of the pressure head. It is worth noting that only a few parts of this curve are described during the infiltration experiment. Indeed, only moderate dry conditions are established because the minimum pressure head reached in the column is $-$120 cm, which corresponds to the initial pressure head at the top of the column.

(a) Pressure head at 5 cm below the soil surface, (b) water content at 5 cm below the soil surface, (c) cumulative outflow, (d) retention curve, (e) output concentration for $T_{\mathrm{inj}}=$ 5000 and (f) for $T_{\mathrm{inj}}=$ 3000 min. Solid lines represent model outputs and dots represent the sets of perturbed data serving as conditioning information for model calibration.

[Figure omitted. See PDF]

The breakthrough concentration curve (Fig. 1e) shows a sharp front, which starts shortly after 3000 min. Note that if the injection of both water and contaminant are stopped once the solute reaches the output. For an injection period of 3000 min, the breakthrough curve exhibits a smoother progression (Fig. 1f).

The data considered as measurements, which are used as conditioning information for model calibration, are also shown in Fig. 1. In Fig. 1b, the water content seems to be more affected by the perturbation of data than the pressure head and cumulative outflow. This phenomenon is due to the relative importance of the measurement errors of the water content often observed with time-domain reflectometry probes and to the weak variations of the water content during the infiltration experiment. The perturbation of the breakthrough curve is relatively small because of the low added noise since output concentrations can be accurately measured. The perturbations of the pressure head and cumulative outflow seem weak because of the large variation of these variables during the experiment.

MCMC solutions for the transport scenario 1. The diagonal plots represent the inferred posterior probability distribution of the model parameters. The off-diagonal scatterplots represent the pairwise correlations in the MCMC drawing.

[Figure omitted. See PDF]

MCMC solutions for transport scenario 2 (see Fig. 2 caption).

[Figure omitted. See PDF]

Results of the parameter estimation

The uncertainty model parameters are assumed to be distributed uniformly over the ranges reported in Table 1. This table also lists the reference values used to generate data observations before perturbation. Seven scenarios are considered, corresponding to different sets of measurements for the estimation of the hydraulic and transport soil parameters (Table 2).

Measurement sets and injection periods for the different scenarios. The pressure head $h$ and the water content $\mathit{\theta }$ are measured at 5 cm from the top of the column. The cumulative outflow $Q$ and the concentration $C$ are measured at the exit of the column.

Scenario	Measured variables				Injection period
	$h$	$\mathit{\theta }$	$Q$	$C$	$T_{\mathrm{inj}}=\mathrm{5000}$ min	$T_{\mathrm{inj}}=\mathrm{3000}$ min
1	$\mathit{\nu }$		$\mathit{\nu }$		$\mathit{\nu }$
2	$\mathit{\nu }$	$\mathit{\nu }$	$\mathit{\nu }$		$\mathit{\nu }$
3	$\mathit{\nu }$	$\mathit{\nu }$	$\mathit{\nu }$	$\mathit{\nu }$	$\mathit{\nu }$
4	$\mathit{\nu }$		$\mathit{\nu }$	$\mathit{\nu }$	$\mathit{\nu }$
5	$\mathit{\nu }$		$\mathit{\nu }$	$\mathit{\nu }$		$\mathit{\nu }$
6			$\mathit{\nu }$	$\mathit{\nu }$	$\mathit{\nu }$
7			$\mathit{\nu }$	$\mathit{\nu }$		$\mathit{\nu }$

The MCMC results of the seven studied scenarios are given in Figs. 2–8. The “on-diagonal” plots in these figures display the inferred parameter distributions, whereas the “off-diagonal” plots represent the pairwise correlations in the MCMC sample. If the draws are independent, non-sloping scatterplots should be observed. However, if a good value of a given parameter is conditioned by the value of another parameter, then their pairwise scatterplot should show a narrow sloping stripe. The sensitivity of parameters is obtained by comparing prior to posterior parameter distribution. A significant difference between the two distributions for a parameter indicates high model sensitivity to that parameter (Dusek et al., 2015).

MCMC solutions for transport scenario 3 (see Fig. 2 caption).

[Figure omitted. See PDF]

MCMC solutions for transport scenario 4 (see Fig. 2 caption).

[Figure omitted. See PDF]

MCMC solutions for transport scenario 5 (see Fig. 2 caption).

[Figure omitted. See PDF]

MCMC solutions for transport scenario 6 (see Fig. 2 caption).

[Figure omitted. See PDF]

MCMC solutions for transport scenario 7 (see Fig. 2 caption).

[Figure omitted. See PDF]

To facilitate the comparison between the different scenarios, Figs. 9–14 show the mean and the 95 % confidence intervals of the final MCMC sample that adequately fit the model onto observations for each scenario, and Table 3 summarizes the pairwise parameter correlations.

Summary of the pairwise parameter correlations.

Scenario
1	$\left(k_{\mathrm{s}},,,n\right)=-\mathrm{0.97}$	$\left(k_{\mathrm{s}},,,\mathit{\alpha }\right)=\mathrm{0.94}$			$\left({\mathit{\theta }}_{\mathrm{r}},,,{\mathit{\theta }}_{\mathrm{s}}\right)=\mathrm{0.96}$
2	$\left(k_{\mathrm{s}},,,n\right)=-\mathrm{0.98}$	$\left(k_{\mathrm{s}},,,\mathit{\alpha }\right)=\mathrm{0.94}$	$\left(k_{\mathrm{s}},,,{\mathit{\theta }}_{\mathrm{r}}\right)=-\mathrm{0.94}$	$\left({\mathit{\theta }}_{\mathrm{r}},,,n\right)=\mathrm{0.98}$
3	$\left(k_{\mathrm{s}},,,n\right)=-\mathrm{0.97}$	$\left(k_{\mathrm{s}},,,\mathit{\alpha }\right)=\mathrm{0.91}$	$\left(k_{\mathrm{s}},,,{\mathit{\theta }}_{\mathrm{r}}\right)=-\mathrm{0.94}$	$\left({\mathit{\theta }}_{\mathrm{r}},,,n\right)=\mathrm{0.99}$
4	$\left(k_{\mathrm{s}},,,n\right)=-\mathrm{0.98}$	$\left(k_{\mathrm{s}},,,\mathit{\alpha }\right)=\mathrm{0.95}$	$\left(k_{\mathrm{s}},,,{\mathit{\theta }}_{\mathrm{r}}\right)=-\mathrm{0.96}$	$\left({\mathit{\theta }}_{\mathrm{r}},,,n\right)=\mathrm{0.99}$
5	$\left(k_{\mathrm{s}},,,n\right)=-\mathrm{0.96}$	$\left(k_{\mathrm{s}},,,\mathit{\alpha }\right)=\mathrm{0.93}$	$\left(k_{\mathrm{s}},,,{\mathit{\theta }}_{\mathrm{r}}\right)=-\mathrm{0.91}$	$\left({\mathit{\theta }}_{\mathrm{r}},,,n\right)=\mathrm{0.98}$
6	$\left(k_{\mathrm{s}},,,n\right)=-\mathrm{0.95}$			$\left({\mathit{\theta }}_{\mathrm{r}},,,n\right)=\mathrm{0.95}$
7	$\left(k_{\mathrm{s}},,,n\right)=-\mathrm{0.95}$			$\left({\mathit{\theta }}_{\mathrm{r}},,,n\right)=\mathrm{0.94}$

Figure 2 shows the inferred distributions of the parameters identified with the MCMC sampler using only the pressure and cumulative outflow measurements (scenario 1). The parameters $k_{\mathrm{s}}$, $\mathit{\alpha }$ and $n$ are well estimated; their prior intervals of variation are strongly narrowed and they essentially show bell-shaped posterior distributions. This shows the high sensitivity of the model responses to these parameters.

The parameter $k_{\mathrm{s}}$ is strongly correlated to $\mathit{\alpha }$ (0.94) and $n$ ($-$0.97). These results confirmed the results of Eching et al. (1994) from multistep outflow experiments where it was found that the inverse solution technique is greatly improved when both cumulative outflow and pressure head data from some positions inside the column are used. The two water-content-related parameters are strongly correlated (0.96) and cannot be identified accurately because the water retention relationship depends on the difference between ${\mathit{\theta }}_{\mathrm{s}}$ and ${\mathit{\theta }}_{\mathrm{r}}$, and only this difference is identifiable. Note that the prior intervals of ${\mathit{\theta }}_{\mathrm{r}}$ and ${\mathit{\theta }}_{\mathrm{s}}$, which are, respectively, $\left[\mathrm{0.05},,,\mathrm{0.2}\right]$ and $\left[\mathrm{0.3},,,\mathrm{0.5}\right]$, have changed to the posterior intervals $\left[\mathrm{0.05},,,\mathrm{0.16}\right]$ and $\left[\mathrm{0.39},,,\mathrm{0.5}\right]$ because the target difference should be ${\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}=\mathrm{0.34}.$ In the literature, van Genuchten and Nielsen (1985), Eching and Hopmans (1993) and Zurmühl (1996) considered that saturated water content is determined independently and considered only ${\mathit{\theta }}_{\mathrm{r}}$ to be an empirical parameter that should be fitted to the data.

The dispersivity coefficient $a_l$ has not been identified in this first scenario.

The MCMC results in Fig. 3 show that water content measurements throughout the experiment (scenario 2) allow the estimation of both the residual and saturated water contents. The parameter ${\mathit{\theta }}_{\mathrm{r}}$ strongly correlates to $k_{\mathrm{s}}$ ($-$0.94) and $n$ (0.98) and the parameter $k_{\mathrm{s}}$ remains strongly related to $\mathit{\alpha }$ (0.94) and $n$ ($-$0.98). Although the water content data are subject to relatively high measurement errors, a good estimation is obtained for ${\mathit{\theta }}_{\mathrm{s}}$ and ${\mathit{\theta }}_{\mathrm{r}}$. The parameters $k_{\mathrm{s}}$, $\mathit{\alpha }$ and $n$ are estimated with the same accuracy as for the first scenario. All parameters (except the dispersivity coefficient) are highly sensitive since their posterior intervals of variations are strongly reduced compared to the prior intervals. Moreover, the prior uniform distributions give place to almost Gaussian posterior distributions. These results show that, although Kool et al. (1985) and Kool and Parker (1988) suggested that the transient experiments should cover a wide range of water contents, an appropriate estimation of all parameters can be obtained with the infiltration experiment even though a limited range of water contents is covered.

When the concentration measurements are also considered in the inversion (scenario 3), the results depicted in Fig. 4 show very significant correlations between $k_{\mathrm{s}}$ and ${\mathit{\theta }}_{\mathrm{r}}$ ($-$0.94), $k_{\mathrm{s}}$ and $\mathit{\alpha }$ (0.91), $k_{\mathrm{s}}$ and $n$ ($-$0.97) and $n$ and ${\mathit{\theta }}_{\mathrm{r}}$ (0.99). The posterior uncertainty ranges of $k_{\mathrm{s}}$, $\mathit{\alpha }$, $n$ and ${\mathit{\theta }}_{\mathrm{r}}$ are similar to the previous scenarios. Those of ${\mathit{\theta }}_{\mathrm{s}}$ and $a_l$ are strongly reduced, leading to a good identification of these parameters when using $C$ measurements (Figs. 10 and 14). A better estimate of the saturated water content is obtained because advective transport is a function of this variable.

Posterior mean values and 95 % confidence intervals of the saturated hydraulic conductivity for the different scenarios.

[Figure omitted. See PDF]

Posterior mean values and 95 % confidence intervals of the saturated water content for the different scenarios.

[Figure omitted. See PDF]

Posterior mean values and 95 % confidence intervals of the residual water content for the different scenarios.

[Figure omitted. See PDF]

In the inversion procedure of scenario 4, the measurements of the water content are not considered. This scenario leads to the same quality of the estimation for the parameters $k_{\mathrm{s}}$, ${\mathit{\theta }}_{\mathrm{r}}$, $\mathit{\alpha }$ and $n$ (Figs. 9, 11, 12, 13) and similar correlations between the parameters as in the previous scenario. This result shows that the intrusive water content measurements, which are subject to more significant measurement errors than the output concentration, are not required if the output concentration is measured. Compared with the results of scenario 2, it can be concluded that better parameter estimations are obtained using $h$, $Q$ and $C$ data than using $h$, $Q$ and $\mathit{\theta }$ data, especially for ${\mathit{\theta }}_{\mathrm{s}}$. Therefore, using $C$ instead of $\mathit{\theta }$ measurements in combination with $h$ and $Q$ measurements allows the estimation of $a_l$ and yields a better estimate of ${\mathit{\theta }}_{\mathrm{s}}$.

The pressure head, cumulative outflow and concentration measurements are used in the estimation procedure of scenario 5, but the injection period is now reduced to $T_{\mathrm{inj}}=$ 3000 min. The obtained results (Fig. 6) show the same correlations between the parameters as for $T_{\mathrm{inj}}=$ 5000 min. For the parameters $k_{\mathrm{s}}$, ${\mathit{\theta }}_{\mathrm{s}}$, ${\mathit{\theta }}_{\mathrm{r}}$, $\mathit{\alpha }$ and $n$, almost the same mean estimates are obtained as for scenario 4. However, the parameters are better identified (Figs. 9–13). Indeed, the uncertainty of these parameters is smaller because the credible interval is reduced by a factor of 25 % for $k_{\mathrm{s}}$, 8 % for ${\mathit{\theta }}_{\mathrm{s}}$, 26 % for ${\mathit{\theta }}_{\mathrm{r}}$, 10 % for $\mathit{\alpha }$ and 25 % for $n$ when compared to the results obtained using $T_{\mathrm{inj}}=\mathrm{5000}$ min. The parameter $a_l$ is also much better estimated than in the previous scenario. Its mean value approaches the reference solution and the posterior uncertainty range is reduced by approximately 75 % (Fig. 14).

In scenario 6, the pressure head measurements are removed and only non-intrusive measurements ($Q$ and $C$ data) are used for the calibration with an injection period of $T_{\mathrm{inj}}=$ 5000 min. These kinds of non-intrusive measures have been used by Mertens et al. (2009) to estimate some of the hydraulic and pesticide leaching parameters. The results depicted in Fig. 7 show high correlations only between $k_{\mathrm{s}}$ and $n$ ($-$0.95) and ${\mathit{\theta }}_{\mathrm{r}}$ and $n$ (0.95). On the one hand, these results show that all the parameters are well estimated since, as compared to the prior intervals (given in Table 1), the confidence intervals of the estimated parameters (plotted in Figs. 9–14) are strongly reduced, especially for the parameters $\mathit{\alpha }$, $n$ and ${\mathit{\theta }}_{\mathrm{s}}$. On the other hand, compared to the results of scenario 4, which also considers pressure data, $k_{\mathrm{s}}$ is not as well estimated (the mean value is not as close to the reference value and the confidence interval is 27 % larger). The mean estimated values for ${\mathit{\theta }}_{\mathrm{r}}$ and $n$ also degraded (not as close to the reference solution), although their confidence intervals are similar to those of scenario 4 (Figs. 11, 13). The estimated mean value of the parameter $\mathit{\alpha }$ is similar to that in scenario 4. However, its uncertainty is much larger because the credible interval is 77 % larger (Fig. 14). The parameters ${\mathit{\theta }}_{\mathrm{s}}$ and $a_l$ are estimated as well as in scenario 4 (in terms of mean estimated value and credible interval).

The last scenario (scenario 7) is similar to the previous one, but the injection period is reduced to $T_{\mathrm{inj}}=$ 3000 min. The results depicted in Fig. 8 show similar correlations between the parameters as for $T_{\mathrm{inj}}=$ 5000 min. However, a significant improvement is observed for the mean estimated values, which approach the reference solution for $k_{\mathrm{s}}$, ${\mathit{\theta }}_{\mathrm{r}}$, $n$ and $a_l$ (Figs. 9, 11, 13, 14). The uncertainties of $k_{\mathrm{s}}$, $\mathit{\alpha }$ and $a_l$ are also reduced by approximately 40, 15 and 70 %, respectively. The parameter ${\mathit{\theta }}_{\mathrm{s}}$ is estimated as well as in scenario 6. The improvement of the parameter estimation in this last scenario compared to the previous one can be explained by the fact that the injection of water and solute contaminant is stopped once the concentration reaches the column outlet. Hence, the injected volume (0.015 $\times$ 3000 $=$ 45 cm${}^{\mathrm{3}}$ cm${}^{-\mathrm{2}})$ is slightly less than the pore volume (120 $\times$ 0.43 $=$ 51 cm${}^{\mathrm{3}}$ cm${}^{-\mathrm{2}})$. Thus, when the injection is stopped, the column is not fully saturated and the outlet flux strongly reduces (see the asymptotic behavior of the cumulative outflow when the injection is stopped in Fig. 1c). As a consequence, the concentration profile increases smoothly (see Fig. 1f) until reaching its maximum value, in contrast to the sharp front observed for $T_{\mathrm{inj}}=\mathrm{5000}min$ in scenario 6 (see Fig. 1e). Hence, the breakthrough curve obtained with $T_{\mathrm{inj}}=$ 3000 min is more affected by the hydraulic parameters than the breakthrough curve obtained with $T_{\mathrm{inj}}=$ 5000 min. This explains why a better estimation of the parameters is observed for the last scenario compared to scenario 6.

Posterior mean values and 95 % confidence intervals of the shape parameter $\mathit{\alpha }$ for the different scenarios.

[Figure omitted. See PDF]

Posterior mean values and 95 % confidence intervals of the shape parameter $n$ for the different scenarios.

[Figure omitted. See PDF]

Posterior mean values and 95 % confidence intervals of dispersivity for the different scenarios.

[Figure omitted. See PDF]

Conclusions

In this work, estimation of hydraulic and transport soil parameters have been investigated using synthetic infiltration experiments performed in a column filled with a sandy clay loam soil, which was subjected to continuous flow and solute injection over a period $T_{\mathrm{inj}}$.

The saturated hydraulic conductivity, the saturated and residual water contents, the Mualem–van Genuchten shape parameters and the longitudinal dispersivity are estimated in a Bayesian framework using the MCMC sampler. Parameter estimation is performed for different scenarios of data measurements.

The results reveal the following conclusions:

• All hydraulic and transport parameters can be appropriately estimated from the described infiltration experiment. However, the accuracy differs and depends on the type of measurement and the duration of the injection $T_{\mathrm{inj}}$, even if the water content remains close to saturated conditions.

• The use of concentration measurements at the column outflow, in addition to traditional measured variables (water content, pressure head and cumulative outflow), reduces the hydraulic parameter uncertainties, especially those of the saturated water content (comparison between scenarios 2 and 3).

• The saturated hydraulic conductivity is estimated with the same order of accuracy, independent of the observed variables.

• The estimation of the dispersivity is sensitive to the injection duration. Scenarios 5 and 7 with $T_{\mathrm{inj}}=$ 3000 min yield much more accurate dispersivity estimations than scenarios 4 and 6 with $T_{\mathrm{inj}}=$ 5000 min due to the extended spreading of the solute observed for $T_{\mathrm{inj}}=$ 3000 min.

• A better identifiability of the soil parameters is obtained using $C$ instead of $\mathit{\theta }$ measurements, in combination with $h$ and $Q$ data (comparison between scenarios 2 and 4).

• Using only non-intrusive measurements (cumulative outflow and output concentration) yields satisfactory estimation of all parameters (scenario 7). The uncertainty of the parameters significantly decreases when the injection of water and solute is maintained for a limited period (comparison between scenarios 6 and 7).

This last point has practical applications for designing simple experimental setups dedicated to the estimation of hydrodynamic and transport parameters for unsaturated flow in soils. The setup has to be appropriately equipped to measure the cumulative water outflow (e.g., weighing machine) and the solute breakthrough at the column outflow (e.g., flow through electrical conductivity). The injection should be stopped as soon as the solute concentration reaches the outflow. The accuracy of the estimation of ${\mathit{\theta }}_{\mathrm{r}}$, $\mathit{\alpha }$ and $n$ improves by adding pressure measurements inside the column, close to the injection.

These results are of course related to the models and experimental conditions we used. This work can be extended to different types of soils, water retention and/or relative permeability functions to evaluate the interest of coupling flow and transport for parameter identification. This work can also be extended to reactive solutes.

No data sets were used in this article.

The authors declare that they have no conflict of interest.

Acknowledgements

The authors are grateful to the French National Research Agency, which funded this work through the program AAP Blanc – SIMI 6 project RESAIN (no. ANR-12-BS06-0010-02). Edited by: H. Cloke Reviewed by: four anonymous referees