# Bayesian Simultaneous Estimation of Unsaturated Flow and Solute Transport Parameters from a Laboratory Infiltration Experiment

Anis Younes, Jabran Zaouali, Sabri Kanzari, Francois Lehmann, Marwan Fahs

► To cite this version:

Anis Younes, Jabran Zaouali, Sabri Kanzari, Francois Lehmann, Marwan Fahs. Bayesian Simultaneous Estimation of Unsaturated Flow and Solute Transport Parameters from a Laboratory Infiltration Experiment. Water, 2019, 11 (8), pp.1660. 10.3390/w11081660. [https://hal.science/hal-02390243](https://hal.science/hal-02390243)

HAL Id: hal-02390243

https://hal.science/hal-02390243

Submitted on 2 Dec 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Bayesian Simultaneous Estimation of Unsaturated Flow and Solute Transport Parameters from a Laboratory Infiltration Experiment

Anis Younes <fnref n="1" />,<fnref n="2" />,<fnref n="3" />, Jabran Zaouali <fnref n="1" />, Sabri Kanzari <fnref n="4" />, Francois Lehmann <fnref n="1" /> and Marwan Fahs <fnref n="1" />,*

Received: 28 July 2019; Accepted: 7 August 2019; Published: 11 August 2019

<figure><img src="image_2.png" /></figure>

check for
updates

**Abstract:** Numerical modeling has become an irreplaceable tool for the investigation of water flow and solute transport in the unsaturated zone. The use of this tool for real situations is often faced with lack of knowledge of hydraulic and soil transport parameters. In this study, advanced experimental and numerical techniques are developed for an accurate estimation of the soil parameters. A laboratory unsaturated flow and solute transport experiment is conducted on a large undisturbed soil column of around 40 cm length. Bromide, used as a nonreactive contaminant, is injected at the surface of the undisturbed soil, followed by a leaching phase. The pressure measurements at different locations along the soil column as well as the outflow bromide concentration are collected during the experiment and used for the statistical calibration of flow and solute transport. The Richards equation, combined with constitutive relations for water content and permeability, is used to describe unsaturated flow. Both linear and non-equilibrium mobile-immobile transport models are investigated for the solute transport. All hydraulic and mass transport parameters are inferred using a one-step Bayesian estimation with the Markov chain Monte Carlo sampler. The results prove that the pressure and concentration measurements are able to identify almost all hydraulic and mass transport parameters. The mobile-immobile transport model better reproduces the infiltration experiment. It produces narrower uncertainty intervals for soil parameters and predictive output concentrations.

**Keywords:** infiltration laboratory experiment; Richards equation; advection-dispersion equation; mobile-immobile model; predictive uncertainty; Markov chain Monte Carlo

## 1. Introduction

Numerical modeling of water flow and solute transport through the vadoze zone is essential to estimate the amount of recharge and/or to prevent soil and subsurface water contamination. It is becoming increasingly important for several applications as water quality management, water resources planning in urban zones, and for mitigating groundwater pollution. The processes of flow and mass transport in the unsaturated zone depend on the driving forces (matric and gravitational potential) and on the soil properties. Thus, to perform real field simulations, the numerical models require an accurate prior knowledge of the soil parameters. These parameters are usually determined using laboratory unsaturated flow and solute transport experiments performed using soil samples collected from the field. Inversion methods can then be applied to fit the model outputs to the laboratory measurements and to assess the soil parameters.

Different laboratory experiments have been conducted for the estimation of the soil parameters e.g., [1–13]. Small scale laboratory column experiments were often used for the characterization

<footnote n="1">¹ LHyGES, University de Strasbourg/EOST/ENGEES, CNRS, 67084 Strasbourg, France</footnote>

<footnote n="2">² LISAH, University Montpellier, INRA, IRD, Montpellier SupAgro, 34060 Montpellier, France</footnote>

<footnote n="3">³ LMHE, Ecole Nationale d'Ingénieurs de Tunis, Tunis 1002, Tunisia</footnote>

<footnote n="4">⁴ INRGREF, Laboratory of Rural Engineering, University of Carthage, Carthage 1054, Tunisia</footnote>

<footnote>* Correspondence: fahs@unistra.fr; Tel.: +33-3-6885-0448</footnote>
of hydraulic and mass transport properties [1,2,11,14,15]. Mishra and Parker [16] showed that the simultaneous estimation of hydraulic and mass transport parameters yields more accurate results than the two-steps estimation in which hydraulic and transport parameters are determined sequentially. Inoue et al. [17] used electrical conductivity and matric pressure head measurements at different depths for the flow-transport inversion with the local search Levenberg-Marquardt algorithm. Recently, several Bayesian approaches, where measurements are combined with prior parameter information to provide posterior parameter distributions, have been investigated for the estimation of the unsaturated hydraulic soil parameters, among others [10,12,18–21]. The term Bayesian is used to describe statistical inversion by considering [22]: (i) That model variables are random, (ii) that randomness describes the degree of information for their realization, and (iii) that the solution of the estimation problem is the posterior probability distribution from which several statistics can be obtained. Bayesian hydraulic soil parameter estimation was investigated for a laboratory drainage experiment in [10]. In [12], it is shown that the infiltration laboratory experiment yields more accurate hydraulic parameters than the percolation-drainage experiment. Bayesian inference of both hydraulic and transport unsaturated soil parameters has been investigated for hypothetical experiments in the case of linear transport [21] and in the case of pesticide transport [23] in unsaturated soils. Moreira et al. [11] considered the estimation of hydraulic and solute transport parameters based on a Bayesian inversion of a laboratory infiltration experiment performed on a small column of 7.3 cm length. One tensiometric sensor was used to monitor the pressure head in the middle of the column. The pressure head and outflow concentrations measurements were employed for parameters estimation. Because of the insufficient collected data, some hydraulic soil parameters were not included in the inversion procedure to prevent identification difficulties. These parameters were considered as known (measured) and they were assigned informative Gaussian priors. Furthermore, Moreira et al. [11] used a sequential procedure in which the hydraulic parameters were firstly estimated using the pressure head measurements. Then, in a second step, these estimates were employed as prior information with Gaussian distributions to determine the mass transport parameters from the concentration measurements.

In the present study, a laboratory infiltration experiment is conducted in order to estimate soil hydraulic and mass transport parameters. Bromide is used as a nonreactive contaminant and injected at the surface of the soil column. The current research differs from the literature by (i) considering a large column of around 40 cm length of an undisturbed sandy soil, (ii) performing a simultaneous estimation of all hydraulic and mass transport soil parameters using the outflow breakthrough curve of bromide and the pressure head measurements at three different locations, (iii) assuming no prior information for all hydraulic and solute transport parameters for which we attribute uniform prior distributions with large intervals to reflect our poor prior knowledge of their values, and (iv) estimating all parameters using a one-step Bayesian procedure via the Markov chain Monte Carlo (MCMC) method [24].

The direct problem is based on the Richards equation coupled with the Darcy's law and the mass conservation equation [25]. The Mualem and van Genuchten models are used as constitutive relations [26,27]. The linear advection-dispersion and the non-equilibrium mobile-immobile models are investigated for bromide transport. Bayesian inversion results are discussed in terms of fitting between simulations and measurements and in terms of uncertainty ranges associated with the estimated parameter values and to the model outputs.

## 2. Materials and Methods

### 2.1. *Laboratory Experiment*

A cylindrical PVC tube of 60 cm length, beveled at its end to form a cutting edge, was sunk into the soil in the region of Bekalta in the east center of Tunisia. A homogeneous soil column of 38.4 cm length and 15 cm diameter has been carefully extracted. The soil, essentially formed by sand, was placed over a porous plate of 6 mm high with a saturated water content of 0.42 cm³ cm⁻³. The soil column (sand + porous plate) is set on a container filled with gravel (Figure 1). An orifice,
located face to the soil bottom, allows effluent recovery and fixing the pressure head h [L] at the lower boundary to zero. The vertical axis Z is downward oriented with the origin located at the soil surface. Three tensiometric sensors, connected to a data logger and a computer, allow monitoring the pressure head near the top (Z₁ = 2.5 cm), middle (Z₂ = 17.5 cm), and bottom of the column (Z₃ = 32.5 cm). The outflow bromide concentration is measured using an ion selective electrode (HI4102 from HANNA instruments). Initially, hydrostatic equilibrium is obtained by verifying that each tensiometer indicates a suction equal to its distance to the bottom of the column. A constant water flow rate of 14 cm³/min was applied at the top surface of the soil using a high-accuracy Masterflex® L/S peristaltic pump. Cotton fiber wicks were deposited at the surface of the soil to obtain a uniform distribution of the flow rate over the soil area. The tensiometric data and the cumulative outflow were monitored until t = 660 min. Infiltration occurs during a period of 562 min, it starts at t = 24 min and was stopped at t = 586 min. During a first period (24 min < t < 220 min), the pump injects a potassium bromide (KBr) solution from a reservoir with a concentration of 0.08 g/L. The period of injection ensures around one pore volume of injected KBr solution. Then, clean water was injected at the surface of the soil with the same flow rate for the period 220 min < t < 586 min. This ensures around two pore volumes of injected clean water. To limit evaporation, the upper end of the column was sealed with the exception of a small orifice, which ensures equilibrium with atmospheric pressure.

<figure><img src="image_3.png" /><figcaption>Figure 1. Schematic of the laboratory infiltration experimental setup.</figcaption></figure>

## 2.2. Numerical Model

The 1D unsaturated flow through the column can be described using Richards equation (RE):

$$ \left\{ \begin{array}{l} \frac{\partial \theta}{\partial t} = (C_s(h) + S_s \frac{\theta}{\partial s}) \frac{\partial h}{\partial t} = - \frac{\partial q_d}{\partial z} \\ q_d = -K(h) \left( \frac{\partial h}{\partial z} - 1 \right) \end{array} \right. \qquad (1) $$

where θ is the water content, t is time, C_s(h) represents the specific moisture capacity, S_s is the specific storage, θ_s is the saturated water content, q_d is the water velocity, and K(h) is the hydraulic conductivity.
Equation (1) is closed using the constitutive relations of Mualem [26] and van Genuchten [27] expressing the hydraulic conductivity in terms of saturation and the saturation in terms of pressure head, respectively.

$$ S_e(h) = \frac{\theta(h) - \theta_r}{\theta_s - \theta_r} = \begin{cases} \frac{1}{(1 + |\alpha h|^m)^{\alpha}} & h < 0 \\ 1 & h \ge 0 \end{cases} \qquad (2) $$ $$ K(S_e) = K_s S_e^{1/2} \left[ 1 - \left( 1 - S_e^{1/m} \right)^m \right]^2 $$

where $\theta_r$ is the residual water content, $S_e$ is the effective saturation, $K_s$ is the hydraulic conductivity at saturation, $m = 1 - 1/n$, and $\alpha$ and $n$ are the shape parameters specific to the soil.

The flow system of Equations (1) and (2) is subject to the following boundary conditions. At the upper boundary, we impose a constant flux ($q_{inj} = 0.08$ cm/min) during the injection period ($24 \text{ min} < t < 586 \text{ min}$); otherwise, the upper flux is zero. A zero-pressure head (Dirichlet boundary condition) is imposed at the column lower boundary. Initial conditions correspond to a hydrostatic pressure distribution ($h = z - L$).

Two transport models are used to investigate bromide transport through the studied soil column. The first model is the well-known linear transport model, based on the advection-dispersion equation:

$$ \frac{\partial(\theta C)}{\partial t} + \frac{\partial(q_d C)}{\partial z} - \frac{\partial}{\partial z}\left(\theta D \frac{\partial C}{\partial z}\right) = 0 \qquad (3) $$

where C is the bromide concentration, $D = a_l q_d + d_m$ is the dispersion coefficient, $d_m$ is the molecular diffusion coefficient fixed to a small value of $3.10^{-4}$ cm$^2$/min as common in literature, and $a_l$ is the soil dispersion coefficient. The transport and flow equations are coupled together through the water content $\theta$ and the Darcy velocity $q_d$.

The second transport model is based on the non-equilibrium approach. It is usually called mobile-immobile model. This model proceeds by considering that the liquid phase in the porous medium can exist as mobile and immobile phases. No advective solute transport is assumed in the immobile zone and a first-order kinetic diffusion transport process occurs between the mobile and immobile zones. The mobile-immobile model in variably saturated porous media writes [28,29]:

$$ \frac{\partial(\theta_m C_m)}{\partial t} + \theta_{im} \frac{\partial C_{im}}{\partial t} + \frac{\partial(q_m C_m)}{\partial z} - \frac{\partial}{\partial z} \left( \theta_m D \frac{\partial C_m}{\partial z} \right) = 0 \qquad (4) $$

$$ \theta_{im} \frac{\partial C_{im}}{\partial t} = \omega(C_m - C_{im}) \qquad (5) $$

$$ \theta = \theta_m + \theta_{im} \qquad (6) $$

where subscripts (m) and (im) represent the mobile and immobile phases, respectively, $\omega$ is the first-order exchange coefficient between phases, and $\theta$ is the total water content. We assume that the mobile and immobile water contents can vary in space and time but the ratio of mobile water content to the total water content ($f = \theta_m / \theta$) is assumed to be constant. As a consequence, the fraction of immobile water $\theta_{im}$ is also constant during the transient flow (see Equation (6)). This assumption is less restrictive and more representative of reality than assuming a constant immobile water content $\theta_{im}$ during the unsaturated flow [11].

For both transport models, initial conditions correspond to a zero-solute concentration inside the whole column. Boundary upper condition ($z = 0$) corresponds to a Dirichlet condition with the concentration fixed to that of the bromide solution during pollutant injection and to zero during clean water injection. The lower boundary condition ($z = L$) corresponds to a zero concentration gradient.

Numerical solution of flow and transport equations in unsaturated porous media is a challenging problem [30–32]. The systems of equations in both linear transport and mobile-immobile transport cases are solved with the standard finite volume method. A uniform mesh of 192 cells of 0.2 cm is employed for the discretization of the soil column and porous plate. This spatial discretization
yields mesh independent solutions. Temporal discretization is performed with a variable order ODE (ordinary differential equation) solver via the method of lines (MOL), as in Fahs et al. [33,34]. MOL has been successfully used for solving flow and/or solute transport in unsaturated porous media in many studies e.g., [10,12,21,35,36]. The ODE solver used in this work (DLSODIS: Double-precision Livermore Solver for ordinary differential equations-implicit form and sparse matrix) allows variable time-step and variable order integration with time error control [37,38].

## 2.3. Bayesian Parameter Inference

The vector of observations $y_{mes}$ is formed by the measurements of the pressure head at three different locations inside the column and the breakthrough curve of bromide at the outflow. The specific storage is fixed to a very small value ($S_s = 10^{-8}cm^{-1}$). All hydraulic and solute transport soil parameters are assumed unknowns and sought by inverse modeling. In the case of linear transport model, the estimated parameters are: $K_S$, $\theta_S$, $\theta_r$, $a_l$, $\alpha$, and $n$ for the soil and the saturated conductivity for the porous plate ($K_S^p$). Hence, the vector of unknown parameters is $\xi = (K_S, K_S^p, \theta_S, \theta_r, a_l, \alpha, n)$. In the case of mobile-immobile model, the vector of unknown parameters contains, in addition, the fraction of mobile water content ($f$) and the exchange coefficient ($\omega$). Hence, in this case, it writes $\xi = (K_S, K_S^p, \theta_S, \theta_r, a_l, \alpha, n, f, \omega)$.

The soil parameters for both linear and mobile-immobile situations are estimated using a Bayesian approach where the measurements are combined with prior parameter information to determine the posterior probability distribution function (pdf). In the following, the pdfs are inferred using DREAM(ZS) software [39] based on the Markov chain Monte Carlo (MCMC) method. The MCMC method has been used by several authors in hydrogeology, e.g., [11,40–44]. With MCMC, random sequences of parameter sets are generated and converge asymptotically toward the target distribution [45]. The latter can then be used to evaluate mean estimated parameter values and confidence intervals to characterize parameter uncertainty. Due to the absence of prior information on the parameter values, independent uniform prior distributions are assumed for all soil parameters. With these non-informative distributions, all possible parameter values are equally likely. Furthermore, the parameter intervals are chosen sufficiently large (Table 1) to reflect our usual prior knowledge of the model parameter values.

Error measurements of the tensiometric pressure sensors and of the ion selective electrode are assumed to be normally and independently distributed with unknown standard deviations. The latter are considered as hyperparameters and are simultaneously estimated with the physical parameters. We denote by $\sigma = (\sigma_h, \sigma_c)$ the vector of hyperparameters.

Setting the calibration problem in a Bayesian framework yields the following posterior pdf,

$$p(\xi, \sigma | y_{mes}) \propto \sigma_h^{-N_h} \sigma_c^{-N_c} \exp\left(-\frac{SS_h(\xi)}{2\sigma_h^2} - \frac{SS_c(\xi)}{2\sigma_c^2}\right) \quad (7)$$

where $SS_h(\xi)$ and $SS_c(\xi)$ are defined as the sums of the square of differences between the observed and predicted data of the pressure head and concentration, respectively.

For instance, $SS_h(\xi) = \sum_{k=1}^{N_h} (h_{mes}^{(k)} - h_{mod}^{(k)}(\xi))^2$ includes the observed ($h_{mes}^{(k)}$) and predicted ($h_{mod}^{(k)}$) pressure heads at time $t_k$ for the number of pressure head observations ($N_h$).

With the MCMC sampler, a new candidate $x^i = (\xi^i, \sigma^i)$ is generated from a proposal distribution $q(x^i|x^{i-1})$, which only depends on the previous accepted candidate. The new candidate can be accepted or rejected depending on the Hasting ratio $\alpha_H$ [46],

$$\alpha_H = \min\left(1, \frac{p(\xi^i, \sigma^i | y_{mes}) q(x^i|x^{i-1})}{p(\xi^{i-1}, \sigma^i | y_{mes}) q(x^{i-1}|x^i)}\right). \quad (8)$$
The last parameter sets that adequately fit the model onto observations are used to appraise the quality of parameter estimate (such as the mean and the standard deviation of the estimated parameters). The MCMC sampler is used hereafter with three parallel chains. Results are considered stationary if the criterion, defined by Gelman and Rubin [47], is verified ($R_{stat} \le 1.2$) and if the chains are not autocorrelated.

<table><caption>**Table 1.** Ranges of the soil parameters.</caption><thead><tr><th>Parameter</th><th>Lower Bound</th><th>Upper Bound</th></tr></thead><tbody><tr><td>K<sub>S</sub> (cm/min)</td><td>0.01</td><td>0.3</td></tr><tr><td>K<sub>S</sub><sup>p</sup> (cm/min)</td><td>5 10<sup>-4</sup></td><td>10<sup>-2</sup></td></tr><tr><td>θ<sub>S</sub> (-)</td><td>0.2</td><td>0.45</td></tr><tr><td>θ<sub>r</sub> (-)</td><td>0.</td><td>0.15</td></tr><tr><td>α (cm<sup>-1</sup>)</td><td>0.001</td><td>0.1</td></tr><tr><td>n (-)</td><td>1.5</td><td>7.0</td></tr><tr><td>a<sub>l</sub> (cm)</td><td>0.1</td><td>3.0</td></tr><tr><td>f (-)</td><td>0.5</td><td>1.</td></tr><tr><td>ω (min<sup>-1</sup>)</td><td>10<sup>-5</sup></td><td>10<sup>-2</sup></td></tr></tbody></table>

### 3. Results and Discussion

The surface injection induces an increase of the pressure head at the tensiometer locations. When the surface injection is stopped, the measured pressure heads decrease as the water front moves below the tensiometers (Figure 2). For each tensiometer, the pressure head is almost constant during the injection period. Note that since the tensiometers are separated by the same distance (15 cm between each other), the jump in the pressure levels during the injection period between the first and second tensiometers is similar to that between the second and third ones (Figure 2). The results of Figure 2 show that only the tensiometer located at $Z_1 = 2.5$ cm from the soil surface continuously measures negative pressure heads. Hence, it is the only tensiometer that remains under unsaturated conditions after the arrival of the water front. The other tensiometers (located at $Z_2 = 17.5$ cm and $Z_3 = 32.5$ cm from the soil surface) measure positive pressure heads indicating that the latter are under the water table. Note that after the end of the infiltration period (i.e., for $t > 586$ min), desaturation occurs and the pressure head at each tensiometer decreases until reaching its initial hydrostatic negative value corresponding to -35.9 cm, -20.9 cm, and -5.8 cm at, respectively, $Z_1 = 2.5$ cm, $Z_2 = 17.5$ cm, and $Z_3 = 32.5$ cm from the soil surface. The pressure measurements during this desaturation period (Figure 2) are prone to contain useful information for the characterization of the unsaturated soil parameters. The measured concentration at the column outflow is given in Figure 3. This figure depicts the logarithm of the bromide concentration in the effluent since the potential measured by the ion-selective electrode is directly proportional to the logarithm of the ion concentration.

Parameter identifiability is assessed for the linear transport model using MCMC with 30,000 model runs. Figure 4 shows that the chains have converged after 14,000 model runs. The values of the mean estimated parameters and the corresponding confidence intervals, for the linear transport model, are depicted in Table 2.
<table><caption>**Table 2.** Estimated mean values, confidence intervals (CIs), and size of the posterior CIs for the linear transport model.</caption><thead><tr><th></th><th>Unit</th><th>Mean Estimated Value</th><th>Confidence Interval</th><th>Size of the CI</th></tr></thead><tbody><tr><td>K<sub>S</sub></td><td>(cm/min)</td><td>0.2</td><td>[0.19–0.21]</td><td>0.013</td></tr><tr><td>K<sub>S</sub><sup>p</sup></td><td>(cm/min)</td><td>0.0022</td><td>[0.0021–0.0023]</td><td>10<sup>-4</sup></td></tr><tr><td>θ<sub>S</sub></td><td>(-)</td><td>0.33</td><td>[0.32–0.34]</td><td>0.016</td></tr><tr><td>θ<sub>r</sub></td><td>(-)</td><td>0.1</td><td>[0.02–0.15]</td><td>0.13</td></tr><tr><td>α</td><td>(cm<sup>-1</sup>)</td><td>0.015</td><td>[0.012–0.018]</td><td>0.006</td></tr><tr><td>n</td><td>(-)</td><td>2.57</td><td>[1.96–3.4]</td><td>1.43</td></tr><tr><td>a<sub>l</sub></td><td>(cm)</td><td>2.36</td><td>[1.98–2.91]</td><td>0.93</td></tr></tbody></table>

The results of this table show that, except the residual water content, all the parameters are well estimated since the posterior intervals of the estimated parameters (Table 2) are strongly narrowed, as compared to their prior intervals (given in Table 1). The significant difference between the prior and posterior distributions of a parameter testifies the high model sensitivity of this parameter. Furthermore, the marginal distributions of all soil parameters (except $θ_r$) exhibit almost bell-shaped posterior distributions (Figure 5).

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Observed (symbols), calibrated (lines), and predictive uncertainty (grey band) for the pressure head at $Z_1$ = 2.5 cm, $Z_2$ = 17.5 cm, and $Z_3$ = 32.5 cm from the soil surface.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Observed (symbols), calibrated (lines), and predictive uncertainty (grey band) for the effluent bromide concentration using the linear transport model.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Convergence of the unsaturated soil parameters for linear transport: Variation of the Gelman–Rubin ($R_{stat}$) with the number of iterations. The convergence is achieved if the chains reach the threshold in broken-line ($R_{stat} \le 1.2$).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Histograms of marginal distribution of parameters with the linear transport model.</figcaption></figure>

The saturated conductivities of the soil $K_S$ and that of the porous plate $K_S^p$ are accurately estimated and show very narrow uncertainty intervals. The sand soil is around 100 times more permeable than the porous plate. Note that $K_S$ and $K_S^p$ are moderately correlated ($r = -0.85$). The saturated water content of the soil is also well estimated with an associated uncertainty less than 5%. The residual water content has no influence on the pressure head nor on the output concentration. Indeed, the posterior interval of $\theta_r$ is similar to its prior interval, which indicate that this parameter is not influent and, hence, cannot be well estimated. This non identifiability is due to the fact that dry conditions were not attained during the experiment. The results of the Figure 5 show that the parameters $\alpha$ and $n$ are quite well estimated, although they are strongly correlated ($r = 0.99$). Finally, the dispersivity transport parameter $a_l$ is quite well estimated with a high mean estimated value of 2.33 cm and a large confidence interval (Figure 5).

Figures 2 and 3 depict the pressure head and the outflow concentration evaluated at the measurement points using results of the MCMC inversion. A good agreement is observed between the predicted and simulated pressure heads using the mean estimated parameter values. However, a less satisfactory matching is observed between the predicted and measured concentrations. The disagreement between concentration profiles is especially notable at the leaching phase. Indeed, the measured concentrations show an asymmetric and long-tailed breakthrough curve, which cannot be reproduced by the symmetrical simulated concentration profile. The 95% uncertainty bands corresponding to the 97.5 and 2.5 percentiles of the distributions are shown in Figures 2 and 3. In these figures, the green regions represent the total predictive uncertainty, which takes into account both parametric uncertainty and measurement errors. All pressure and concentration observations fit
within the 95% total uncertainty range. Small uncertainties were obtained for pressure head responses, while very large uncertainties were assigned to the concentration outputs. This phenomenon is related to the bad matching between measured and simulated concentrations. This bad matching is responsible of the large dispersivity coefficient found by the MCMC sampler in order to reproduce the long-tailed breakthrough curve of the observed concentrations.

The Bayesian inversion of the mobile–immobile model, which includes nine unknown parameters, was performed with the MCMC sampler using 40,000 model runs. The convergence of the MCMC chains requires more simulations to reach the convergence than with the linear transport model. As reported in Figure 6, the overall chains have converged after 24,000 model runs. The mean values of the estimated parameters and the associated confidence intervals are summarized in Table 3. As previously, the posterior confidence intervals of all parameters, except $\theta_r$, are much narrower than their prior intervals. The histogram of marginal distribution of $\theta_r$ shows almost a uniform distribution similar to its prior one. The histograms of marginal distributions of the rest of parameters show almost bell-shaped distributions (Figure 7) significantly different from their prior uniform distributions which indicate the high sensitivity of the parameters.

The results of Table 3 and Figure 7 show that the new two parameters *f* and *ω* are well estimated. Note that the fraction of the mobile water content to the total water content *f* is very accurately estimated with an associated uncertainty less than 2%. A moderate correlation (*r* = −0.73) is observed between the parameters *f* and *ω*. The comparison between the results of the linear transport model (Table 2) and the mobile–immobile one (Table 3) show that the two models yield almost similar mean parameter estimated values and similar confidence intervals for the saturated hydraulic conductivities ($K_S$ and $K_S^p$) and the van-Genuchten parameters ($\alpha$ and $n$). As previously, the parameters $K_S$ and $K_S^p$ are moderately correlated (*r* = −0.87), whereas the parameters $\alpha$ and $n$ are very strongly correlated (*r* = 0.99). The mean estimated values of the saturated water content of the two models are almost identical. However, the associate uncertainty intervals are significantly different. Indeed, the inversion of the mobile–immobile model yields an uncertainty interval three times narrower than with the linear transport model. Hence, $\theta_S$ is much more accurately estimated with the mobile–immobile transport model, although both models yield the same mean estimated parameter value. As with the linear transport model, the residual water content is irrelevant and, hence, cannot be identified. Among the parameters estimated by the inversion of the two transport models, the dispersivity coefficient shows the most significant differences both in terms of mean estimated value and in terms of associate uncertainty range. Indeed, the inversion of the mobile–immobile transport model yields a mean estimated $a_l$ value that is 3.3 times less than the one obtained with the linear transport model. The confidence interval is also more than 3 times narrower with the mobile–immobile transport model.

<table><caption>**Table 3.** Estimated mean values, confidence intervals (CIs), and size of the posterior CIs for the mobile–immobile model.</caption><thead><tr><th></th><th>Unit</th><th>Mean Estimated Value</th><th>Confidence Interval</th><th>Size of the CI</th></tr></thead><tbody><tr><td>K<sub>S</sub></td><td>[cm/min]</td><td>0.2</td><td>[0.19–0.21]</td><td>0.012</td></tr><tr><td>K<sub>S</sub><sup>p</sup></td><td>[cm/min]</td><td>0.0022</td><td>[0.0021–0.0023]</td><td>1 10<sup>-4</sup></td></tr><tr><td>θ<sub>S</sub></td><td>[-]</td><td>0.32</td><td>[0.32–0.33]</td><td>0.005</td></tr><tr><td>θ<sub>r</sub></td><td>[-]</td><td>0.11</td><td>[0.02–0.15]</td><td>0.13</td></tr><tr><td>α</td><td>[cm<sup>-1</sup>]</td><td>0.015</td><td>[0.012–0.018]</td><td>0.007</td></tr><tr><td>n</td><td>[-]</td><td>2.43</td><td>[1.95–3.24]</td><td>1.28</td></tr><tr><td>a<sub>l</sub></td><td>[cm]</td><td>0.7</td><td>[0.58–0.86]</td><td>0.27</td></tr><tr><td>f</td><td>[-]</td><td>0.915</td><td>[0.907–0.922]</td><td>0.015</td></tr><tr><td>ω</td><td>[min<sup>-1</sup>]</td><td>5 10<sup>-4</sup></td><td>[4 10<sup>-4</sup>–5.9 10<sup>-4</sup>]</td><td>1.9 10<sup>-4</sup></td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Convergence of the unsaturated soil parameters for the mobile–immobile transport model: Variation of the Gelman–Rubin ($R_{stat}$) with the number of iterations. The convergence is achieved if the chains reach the threshold in broken-line ($R_{stat} \le 1.2$).</figcaption></figure>

Because the two transport models yield almost the same mean estimated values and uncertainty ranges for all hydraulic parameters ($K_S, K_S^p, \theta_S, \theta_r, \alpha, n$), the matching between measured and calibrated pressure heads by inverting the mobile–immobile transport model (not shown) is similar to the matching depicted in Figure 2 obtained by inverting the linear transport model. Likewise, the total predictive uncertainty associated with the pressure head output with the mobile–immobile model (not shown) is similar to the one observed with the linear transport model (grey bands in Figure 2).

Measured and simulated output concentrations using the mean estimated parameter values obtained from the inversion of the mobile–immobile transport model are plotted in Figure 8. Contrary to the results of the linear transport model (Figure 3), a very good matching is observed between measured and simulated concentrations during the whole experiment (including the leaching phase). Furthermore, the total uncertainty (grey bands) characterizing the concentration output is much narrower with the mobile–immobile transport model (Figure 8) than with the linear one (Figure 3). Hence, for the investigated undisturbed soil, the mobile–immobile transport model better reproduces bromide infiltration and yields much more accurate output concentrations than the linear one.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Histograms of marginal distribution of parameters with the mobile–immobile transport model.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Observed (symbols), calibrated (lines), and predictive uncertainty (grey band) for the effluent bromide concentration using the mobile-immobile transport model.</figcaption></figure>

**4. Conclusions**

In this study, a laboratory infiltration experiment was conducted by injecting bromide solution at the surface of an undisturbed soil column. The local pressure head measurements inside the column and the outflow concentration measurements were used to simultaneously estimate all hydraulic and transport soil parameters. Unsaturated water flow is simulated using the Richards equation coupled with the Mualem/van Genuchten models. The linear advection-dispersion and the non-equilibrium mobile-immobile models were investigated for modeling the bromide propagation. A Bayesian inversion was performed using the pressure head and concentration measurements. In the case of linear transport model, seven soil parameters were estimated: The saturated conductivity $K_S$, the saturated water content $\theta_S$, the residual water content $\theta_r$, the dispersivity coefficient $a_l$, and the van Genuchten parameters ($\alpha$ and $n$) of the soil and the saturated conductivity $K_S^p$ of the porous plate. In the case of mobile-immobile model, the vector of unknown parameters contains, in addition, the fraction of mobile water content $f$ and the exchange coefficient $\omega$. All hydraulic and transport soil parameters were estimated in a one-step procedure using the MCMC sampler starting from non-informative uniform distributions for all parameters. The obtained results are summarized hereafter:

- The flow through the investigated soil can be well reproduced by the Richards equation combined with the Mualem/van-Genuchten models.

- The pressure head and outflow concentration measurements are able to identify all the unsaturated hydraulic and solute transport soil parameters, except the residual water content $\theta_r$ which was not identified with either of the models since dry conditions were not attained during the laboratory experiment.

- significant difference was observed between the two models for the estimation of the dispersivity coefficient. Indeed, with the mobile-immobile transport model, the estimated value of $a_l$ was three times lower than the one obtained with the linear transport model. The associated uncertainty was also three times smaller with the mobile-immobile transport model.

- The linear transport model yielded a good agreement between the measured and simulated pressure heads, but a less satisfactory matching was observed between measured and simulated
concentrations. Furthermore, small uncertainties were obtained for pressure head responses, while very large uncertainties were assigned to the output concentrations.

- The mobile-immobile transport model better reproduced the infiltration experiment. Indeed, a very good matching was obtained between measured and simulated concentrations. Furthermore, the concentration uncertainty region was much narrower than with the linear transport model.

Finally, the results of this work point out the necessity of testing different models for interpreting real experiments and to employ accurate Bayesian methods for model calibration. These methods provide accurate information concerning the estimated parameters, their confidence intervals, and model output uncertainties.

**Author Contributions:** J.Z. and S.K. worked on the laboratory experiment; F.L. worked on the numerical simulations; A.Y. and M.F. supervised the research and finalized the manuscript.

**Funding:** This research was partially supported by the French National Research Agency (ANR Hydrocrizsto) and the Tunisian-French joint international laboratory NAILA (http://www.lmi-naila.com/).

**Acknowledgments:** The authors express their thanks to Radhouane Hamdi for his careful assistance during the experimental work.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Kool, J.B.; Parker, J.C.; van Genuchten, M.T. Determining Soil Hydraulic Properties from One-step Outflow Experiments by Parameter Estimation: I. Theory and Numerical Studies. *Soil Sci. Soc. Am. J.* **1985**, *49*, 1348-1354. [CrossRef]

2. Parker, J.C.; Kool, J.B.; van Genuchten, M.T. Determining Soil Hydraulic Properties from One-step Outflow Experiments by Parameter Estimation: II. Experimental Studies. *Soil Sci. Soc. Am. J.* **1985**, *49*, 1354-1359. [CrossRef]

3. Van Dam, J.C.; Stricker, J.N.M.; Droogers, P. Inverse Method to Determine Soil Hydraulic Functions from Multistep Outflow Experiments. *Soil Sci. Soc. Am. J.* **1994**, *58*, 647-652. [CrossRef]

4. Durner, W.; Schultze, B.; Zurmühl, T. *State-of-the-Art in Inverse Modeling of Inflow/Outflow Experiments*; University of California: Riverside, CA, USA, 1999.

5. Beydoun, H.; Lehmann, F. Expériences de drainage et estimation de paramètres en milieu poreux non saturé. *C. R. Geosci.* **2006**, *338*, 180-187. [CrossRef]

6. Puhlmann, H.; Von Wilpert, K.; Lukes, M.; Dröge, W. Multistep outflow experiments to derive a soil hydraulic database for forest soils. *Eur. J. Soil Sci.* **2009**, *60*, 792-806. [CrossRef]

7. Kumar, S.; Sekhar, M.; Reddy, D.V.; Mohan Kumar, M.S. Estimation of soil hydraulic properties and their uncertainty: Comparison between laboratory and field experiment. *Hydrol. Process.* **2010**, *24*, 3426-3435. [CrossRef]

8. Schelle, H.; Iden, S.C.; Peters, A.; Durner, W. Analysis of the Agreement of Soil Hydraulic Properties Obtained from Multistep-Outflow and Evaporation Methods. *Vadose Zone J.* **2010**, *9*, 1080-1091. [CrossRef]

9. Durner, W.; Iden, S.C. Extended multistep outflow method for the accurate determination of soil hydraulic properties near water saturation. *Water Resour. Res.* **2011**, *47*, 13. [CrossRef]

10. Younes, A.; Mara, T.A.; Fajraoui, N.; Lehmann, F.; Belfort, B.; Beydoun, H. Use of Global Sensitivity Analysis to Help Assess Unsaturated Soil Hydraulic Parameters. *Vadose Zone J.* **2013**, *12*, 1-12. [CrossRef]

11. Moreira, P.H.S.; van Genuchten, M.T.; Orlande, H.R.B.; Cotta, R.M. Bayesian estimation of the hydraulic and solute transport properties of a small-scale unsaturated soil column. *J. Hydrol. Hydromech.* **2016**, *64*, 30-44. [CrossRef]

12. Younes, A.; Zaouali, J.; Fahs, M.; Slama, F.; Grunberger, O.; Mara, T.A. Bayesian soil parameter estimation: Results of percolation-drainage vs infiltration laboratory experiments. *J. Hydrol.* **2018**, *565*, 770-778. [CrossRef]
13. Szabó, N.P. A genetic meta-algorithm-assisted inversion approach: Hydrogeological study for the determination of volumetric rock properties and matrix and fluid parameters in unsaturated formations. *Hydrogeol. J.* **2018**, *26*, 1935–1946. [CrossRef]

14. Abbaspour, K.C.; Johnson, C.A.; van Genuchten, M.T. Estimating Uncertain Flow and Transport Parameters Using a Sequential Uncertainty Fitting Procedure. *Vadose Zone J.* **2004**, *3*, 1340–1352. [CrossRef]

15. Abbaspour, K.C.; Schulin, R.; van Genuchten, M.T. Estimating unsaturated soil hydraulic parameters using ant colony optimization. *Adv. Water Resour.* **2001**, *24*, 827–841. [CrossRef]

16. Mishra, S.; Parker, J.C. Parameter estimation for coupled unsaturated flow and transport. *Water Resour. Res.* **1989**, *25*, 385–396. [CrossRef]

17. Inoue, M.; Šimůnek, J.; Shiozawa, S.; Hopmans, J.W. Simultaneous estimation of soil hydraulic and solute transport parameters from transient infiltration experiments. *Adv. Water Resour.* **2000**, *23*, 677–688. [CrossRef]

18. Laloy, E.; Weynants, M.; Bielders, C.L.; Vanclooster, M.; Javaux, M. How efficient are one-dimensional models to reproduce the hydrodynamic behavior of structured soils subjected to multi-step outflow experiments? *J. Hydrol.* **2010**, *393*, 37–52. [CrossRef]

19. Diamantopoulos, E.; Iden, S.C.; Durner, W. Inverse modeling of dynamic nonequilibrium in water flow with an effective approach. *Water Resour. Res.* **2012**, *48*, 16. [CrossRef]

20. Mara, T.A.; Delay, F.; Lehmann, F.; Younes, A. A comparison of two Bayesian approaches for uncertainty quantification. *Environ. Model. Softw.* **2016**, *82*, 21–30. [CrossRef]

21. Younes, A.; Mara, T.; Fahs, M.; Grunberger, O.; Ackerer, P. Hydraulic and transport parameter assessment using column infiltration experiments. *Hydrol. Earth Syst. Sci.* **2017**, *21*, 2263–2275. [CrossRef]

22. Kaipio, J.; Somersalo, E. *Statistical and computational inverse problems*; Springer: New York, NY, USA, 2005; ISBN 978-0-387-22073-4.

23. Younes, A.; Mara, T.A.; Voltz, M.; Guellouz, L.; Musa Baalousha, H.; Fahs, M. A new efficient Bayesian parameter inference strategy: Application to flow and pesticide transport through unsaturated porous media. *J. Hydrol.* **2018**, *563*, 887–899. [CrossRef]

24. Lee, P.M. *Bayesian Statistics: An Introduction*, 4th ed.; Wiley: Chichester, West Sussex, UK; Hoboken, NJ, USA, 2012; ISBN 978-1-118-33257-3.

25. Zhu, H.; Liu, T.; Xue, B.; Wang, G. Modified Richards' Equation to Improve Estimates of Soil Moisture in Two-Layered Soils after Infiltration. *Water* **2018**, *10*, 1174. [CrossRef]

26. Mualem, Y. A new model for predicting the hydraulic conductivity of unsaturated porous media. *Water Resour. Res.* **1976**, *12*, 513–522. [CrossRef]

27. Van Genuchten, M.T. A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils. *Soil Sci. Soc. Am. J.* **1980**, *44*, 892–898. [CrossRef]

28. Van Genuchten, M.T.; Wagenet, R.J. Two-Site/Two-Region Models for Pesticide Transport and Degradation: Theoretical Development and Analytical Solutions. *Soil Sci. Soc. Am. J.* **1989**, *53*, 1303–1310. [CrossRef]

29. Clothier, B.E.; Vogeler, I.; Green, S.R.; Scotter, D.R. Transport in unsaturated soil: Aggregates, macropores, and exchange. In *Physical Nonequilibrium in Soils: Modeling and Application*; Ann Arbor Press: Chelsea, MI, USA, 1998; pp. 273–295, ISBN 9781575040493.

30. Farthing, M.W.; Ogden, F.L. Numerical Solution of Richards' Equation: A Review of Advances and Challenges. *Soil Sci. Soc. Am. J.* **2017**, *81*, 1257–1296. [CrossRef]

31. Ku, C.-Y.; Liu, C.-Y.; Xiao, J.-E.; Yeih, W. Transient Modeling of Flow in Unsaturated Soils Using a Novel Collocation Meshless Method. *Water* **2017**, *9*, 954. [CrossRef]

32. Zha, Y.; Yang, J.; Zeng, J.; Tso, C.-H.M.; Zeng, W.; Shi, L. Review of numerical solution of Richardson-Richards equation for variably saturated flow in soils. *WIREs Water* **2019**, *6*, e1364. [CrossRef]

33. Fahs, M.; Younes, A.; Lehmann, F. An easy and efficient combination of the Mixed Finite Element Method and the Method of Lines for the resolution of Richards' Equation. *Environ. Model. Softw.* **2009**, *24*, 1122–1126. [CrossRef]

34. Fahs, M.; Koohbor, B.; Belfort, B.; Ataie-Ashtiani, B.; Simmons, C.; Younes, A.; Ackerer, P. A Generalized Semi-Analytical Solution for the Dispersive Henry Problem: Effect of Stratification and Anisotropy on Seawater Intrusion. *Water* **2018**, *10*, 230. [CrossRef]

35. Miller, C.T.; Abhishek, C.; Farthing, M.W. A spatially and temporally adaptive solution of Richards' equation. *Adv. Water Resour.* **2006**, *29*, 525–545. [CrossRef]
36. Younes, A.; Zaouali, J.; Lehmann, F.; Fahs, M. Sensitivity and identifiability of hydraulic and geophysical parameters from streaming potential signals in unsaturated porous media. *Hydrol. Earth Syst. Sci.* **2018**, *22*, 3561–3574. [CrossRef]

37. Hindmarsh, A.C. LSODE and LSODI, two new initial value ordinary differential equation solvers. *ACM Signum Newsl.* **1980**, *15*, 10–11. [CrossRef]

38. Radhakrishnan, K. *Description and Use of LSODE, the Livermore Solver for Ordinary Differential Equations*; NASA report; NASA: Washington, DC, USA, 1993.

39. Laloy, E.; Vrugt, J.A. High-dimensional posterior exploration of hydrologic models using multiple-try DREAM (ZS) and high-performance computing. *Water Resour. Res.* **2012**, *48*, W01526. [CrossRef]

40. Vrugt, J.A.; ter Braak, C.J.F.; Clark, M.P.; Hyman, J.M.; Robinson, B.A. Treatment of input uncertainty in hydrologic modeling: Doing hydrology backward with Markov chain Monte Carlo simulation: FORCING DATA ERROR USING MCMC SAMPLING. *Water Resour. Res.* **2008**, *44*, W00B09. [CrossRef]

41. Linde, N.; Ginsbourger, D.; Irving, J.; Nobile, F.; Doucet, A. On uncertainty quantification in hydrogeology and hydrogeophysics. *Adv. Water Resour.* **2017**, *110*, 166–181. [CrossRef]

42. Younes, A.; Delay, F.; Fajraoui, N.; Fahs, M.; Mara, T.A. Global sensitivity analysis and Bayesian parameter inference for solute transport in porous media colonized by biofilms. *J. Contam. Hydrol.* **2016**, *191*, 1–18. [CrossRef]

43. Hayek, M.; Younes, A.; Zaouali, J.; Fajraoui, N.; Fahs, M. Analytical solution and Bayesian inference for interference pumping tests in fractal dual-porosity media. *Comput. Geosci.* **2018**, *22*, 413–421. [CrossRef]

44. Rajabi, M.M.; Ataie-Ashtiani, B.; Simmons, C.T. Model-data interaction in groundwater studies: Review of methods, applications and future directions. *J. Hydrol.* **2018**, *567*, 457–477. [CrossRef]

45. Gelman, A. *Bayesian Data Analysis*, 3rd ed.; Chapman & Hall/CRC Texts in Statistical Science; CRC Press: Boca Raton, FL, USA, 2014; ISBN 978-1-4398-4095-5.

46. Hastings, W.K. Monte Carlo sampling methods using Markov chains and their applications. *Biometrika* **1970**, *57*, 97–109. [CrossRef]

47. Gelman, A.; Rubin, D.B. Inference from Iterative Simulation Using Multiple Sequences. *Stat. Sci.* **1992**, *7*, 457–472. [CrossRef]

<figure><img src="image_11.png" /></figure>

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).