Next Article in Journal
Removal of Diesel Oil in Soil Microcosms and Implication for Geophysical Monitoring
Next Article in Special Issue
A Modeling Platform for Landslide Stability: A Hydrological Approach
Previous Article in Journal
Hydraulic Capacity and Efficiency of a Low-Speed Nonpressurized Coil Pump
Previous Article in Special Issue
Smoothing of Slug Tests for Laboratory Scale Aquifer Assessment—A Comparison among Different Porous Media
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
LHyGES, University de Strasbourg/EOST/ENGEES, CNRS, 67084 Strasbourg, France
2
LISAH, University Montpellier, INRA, IRD, Montpellier SupAgro, 34060 Montpellier, France
3
LMHE, Ecole Nationale d’Ingénieurs de Tunis, Tunis 1002, Tunisia
4
INRGREF, Laboratory of Rural Engineering, University of Carthage, Carthage 1054, Tunisia
*
Author to whom correspondence should be addressed.
Water 2019, 11(8), 1660; https://doi.org/10.3390/w11081660
Submission received: 28 July 2019 / Revised: 30 July 2019 / Accepted: 7 August 2019 / Published: 11 August 2019

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.

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,2,3,4,5,6,7,8,9,10,11,12,13]. Small scale laboratory column experiments were often used for the characterization 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,19,20,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 cm3 cm−3. 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 (Z1 = 2.5 cm), middle (Z2 = 17.5 cm), and bottom of the column (Z3 = 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 cm3/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.

2.2. Numerical Model

The 1D unsaturated flow through the column can be described using Richards equation (RE):
{ θ t = ( C s ( h ) + S s θ θ S ) h t = q d z q d = K ( h ) ( h z 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 ) = θ ( h ) θ r θ s θ r = { 1 ( 1 + | α h | n ) m h < 0 1 h 0 K ( S e ) = K s S e 1 / 2 [ 1 ( 1 S e 1 / m ) m ] 2
where θ r is the residual water content, Se is the effective saturation, K s is the hydraulic conductivity at saturation, m = 1 1 / n , and α 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 i n j = 0.08   cm / min ) during the injection period (24 min < t < 586 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:
( θ C ) t + ( q d C ) z z ( θ D C z ) = 0
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 cm2/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 θ 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]:
( θ m C m ) t + θ i m C i m t + ( q C m ) z z ( θ m D C m z ) = 0
θ i m C i m t = ω ( C m C i m )
θ = θ m + θ i m
where subscripts (m) and (im) represent the mobile and immobile phases, respectively, ω is the first-order exchange coefficient between phases, and θ 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 = θ m / θ ) is assumed to be constant. As a consequence, the fraction of immobile water θ i m 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 θ i m 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,31,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 m e s 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 , θ S , θ r , a l , α , and n for the soil and the saturated conductivity for the porous plate ( K S p ). Hence, the vector of unknown parameters is ξ = ( K S , K S p , θ S , θ r , a l , α , 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 ( ω ). Hence, in this case, it writes ξ = ( K S , K S p , θ S , θ r , a l , α , n , f , ω ) .
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,41,42,43,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 poor 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 σ = ( σ h , σ c ) the vector of hyperparameters.
Setting the calibration problem in a Bayesian framework yields the following posterior pdf,
p ( ξ , σ | y m e s ) σ h N h σ C N C exp ( S S h ( ξ ) 2 σ h 2 S S C ( ξ ) 2 σ C 2 )
where S S h ( ξ ) and S S C ( ξ ) 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, S S h ( ξ ) = k = 1 N h ( h m e s ( k ) h mod ( k ) ( ξ ) ) 2 includes the observed ( h m e s ( 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 = ( ξ i , σ 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 α H [46],
α H = min ( 1 , p ( ξ i , σ i | y m e s ) q ( x i | x i 1 ) p ( ξ i 1 , σ i | y m e s ) q ( x i 1 | x i ) ) .
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 s t a t 1.2 ) and if the chains are not autocorrelated.

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 Z1 = 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 Z2 = 17.5 cm and Z3 = 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, Z1 = 2.5 cm, Z2 = 17.5 cm, and Z3 = 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.
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).
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 θ 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 α 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).
Figure 2 and Figure 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 Figure 2 and Figure 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 θ r , are much narrower than their prior intervals. The histogram of marginal distribution of θ 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 ( α and n ). As previously, the parameters K S and K S p are moderately correlated ( r = 0.87 ), whereas the parameters α 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, θ 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 term 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.
Because the two transport models yield almost the same mean estimated values and uncertainty ranges for all hydraulic parameters ( K S , K S p , θ S , θ r , α , 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.

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 θ S , the residual water content θ r , the dispersivity coefficient a l , and the van Genuchten parameters ( α 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 ω . 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 θ 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. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar]
  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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [CrossRef]
  16. Mishra, S.; Parker, J.C. Parameter estimation for coupled unsaturated flow and transport. Water Resour. Res. 1989, 25, 385–396. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [CrossRef] [Green Version]
  22. Kaipio, J.; Somersalo, E. Statistical and computational inverse problems; Springer: New York, NY, USA, 2005; ISBN 978-0-387-22073-4. [Google Scholar]
  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. [Google Scholar] [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. [Google Scholar]
  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. [Google Scholar] [CrossRef]
  26. Mualem, Y. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 1976, 12, 513–522. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [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. [Google Scholar]
  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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [CrossRef] [Green Version]
  37. Hindmarsh, A.C. LSODE and LSODI, two new initial value ordinary differnetial equation solvers. ACM Signum Newsl. 1980, 15, 10–11. [Google Scholar] [CrossRef]
  38. Radhakrishnan, K. Description and Use of LSODE, the Livermore Solver for Ordinary Differential Equations; NASA report; NASA: Washington, DC, USA, 1993. [Google Scholar]
  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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [CrossRef]
  43. Hayek, M.; Younes, A.; Zouali, 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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar]
  46. Hastings, W.K. Monte Carlo sampling methods using Markov chains and their applications. Biometrika 1970, 57, 97–109. [Google Scholar] [CrossRef]
  47. Gelman, A.; Rubin, D.B. Inference from Iterative Simulation Using Multiple Sequences. Stat. Sci. 1992, 7, 457–472. [Google Scholar] [CrossRef]
Figure 1. Schematic of the laboratory infiltration experimental setup.
Figure 1. Schematic of the laboratory infiltration experimental setup.
Water 11 01660 g001
Figure 2. Observed (symbols), calibrated (lines), and predictive uncertainty (grey band) for the pressure head at Z1 = 2.5 cm, Z2 = 17.5 cm, and Z3 = 32.5 cm from the soil surface.
Figure 2. Observed (symbols), calibrated (lines), and predictive uncertainty (grey band) for the pressure head at Z1 = 2.5 cm, Z2 = 17.5 cm, and Z3 = 32.5 cm from the soil surface.
Water 11 01660 g002
Figure 3. Observed (symbols), calibrated (lines), and predictive uncertainty (grey band) for the effluent bromide concentration using the linear transport model.
Figure 3. Observed (symbols), calibrated (lines), and predictive uncertainty (grey band) for the effluent bromide concentration using the linear transport model.
Water 11 01660 g003
Figure 4. Convergence of the unsaturated soil parameters for linear transport: Variation of the Gelman–Rubin (Rstat) with the number of iterations. The convergence is achieved if the chains reach the threshold in broken-line (Rstat ≤ 1.2).
Figure 4. Convergence of the unsaturated soil parameters for linear transport: Variation of the Gelman–Rubin (Rstat) with the number of iterations. The convergence is achieved if the chains reach the threshold in broken-line (Rstat ≤ 1.2).
Water 11 01660 g004
Figure 5. Histograms of marginal distribution of parameters with the linear transport model.
Figure 5. Histograms of marginal distribution of parameters with the linear transport model.
Water 11 01660 g005
Figure 6. Convergence of the unsaturated soil parameters for the mobile–immobile transport model: Variation of the Gelman–Rubin (Rstat) with the number of iterations. The convergence is achieved if the chains reach the threshold in broken-line (Rstat ≤ 1.2).
Figure 6. Convergence of the unsaturated soil parameters for the mobile–immobile transport model: Variation of the Gelman–Rubin (Rstat) with the number of iterations. The convergence is achieved if the chains reach the threshold in broken-line (Rstat ≤ 1.2).
Water 11 01660 g006
Figure 7. Histograms of marginal distribution of parameters with the mobile–immobile transport model.
Figure 7. Histograms of marginal distribution of parameters with the mobile–immobile transport model.
Water 11 01660 g007
Figure 8. Observed (symbols), calibrated (lines), and predictive uncertainty (grey band) for the effluent bromide concentration using the mobile-immobile transport model.
Figure 8. Observed (symbols), calibrated (lines), and predictive uncertainty (grey band) for the effluent bromide concentration using the mobile-immobile transport model.
Water 11 01660 g008
Table 1. Ranges of the soil parameters.
Table 1. Ranges of the soil parameters.
ParameterLower BoundUpper Bound
K S (cm/min)0.010.3
K S p (cm/min)5 10−410−2
θ S (-)0.20.45
θ r (-)0.0.15
α (cm−1)0.0010.1
n (-)1.57.0
a l (cm)0.13.0
f (-)0.51.
ω (min−1)10−510−2
Table 2. Estimated mean values, confidence intervals (CIs), and size of the posterior CIs for the linear transport model.
Table 2. Estimated mean values, confidence intervals (CIs), and size of the posterior CIs for the linear transport model.
UnitMean Estimated ValueConfidence IntervalSize of the CI
K S (cm/min)0.2[0.19–0.21]0.013
K S p (cm/min)0.0022[0.0021–0.0023]10−4
θ S (-)0.33[0.32–0.34]0.016
θ r (-)0.1[0.02–0.15]0.13
α (cm−1)0.015[0.012–0.018]0.006
n (-)2.57 [1.96–3.4]1.43
a l (cm)2.36[1.98–2.91]0.93
Table 3. Estimated mean values, confidence intervals (CIs), and size of the posterior CIs for the mobile–immobile model.
Table 3. Estimated mean values, confidence intervals (CIs), and size of the posterior CIs for the mobile–immobile model.
UnitMean Estimated ValueConfidence IntervalSize of the CI
K S [cm/min]0.2[0.19–0.21]0.012
K S p [cm/min]0.0022[0.0021–0.0023]1 10−4
θ S [-]0.32[0.32–0.33]0.005
θ r [-]0.11[0.02–0.15]0.13
α [cm−1]0.015[0.012–0.018]0.007
n [-]2.43 [1.95–3.24]1.28
a l [cm]0.7[0.58–0.86]0.27
f [-]0.915[0.907–0.922]0.015
ω [min−1]5 10−4[4 10−4–5.9 10−4]1.9 10−4

Share and Cite

MDPI and ACS Style

Younes, A.; Zaouali, J.; Kanzari, S.; Lehmann, F.; Fahs, M. Bayesian Simultaneous Estimation of Unsaturated Flow and Solute Transport Parameters from a Laboratory Infiltration Experiment. Water 2019, 11, 1660. https://doi.org/10.3390/w11081660

AMA Style

Younes A, Zaouali J, Kanzari S, Lehmann F, Fahs M. Bayesian Simultaneous Estimation of Unsaturated Flow and Solute Transport Parameters from a Laboratory Infiltration Experiment. Water. 2019; 11(8):1660. https://doi.org/10.3390/w11081660

Chicago/Turabian Style

Younes, Anis, Jabran Zaouali, Sabri Kanzari, Francois Lehmann, and Marwan Fahs. 2019. "Bayesian Simultaneous Estimation of Unsaturated Flow and Solute Transport Parameters from a Laboratory Infiltration Experiment" Water 11, no. 8: 1660. https://doi.org/10.3390/w11081660

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop