J. Hydrol. Hydromech., 64, 2016, 1, 30–44
DOI: 10.1515/johh-2016-0002

Bayesian estimation of the hydraulic and solute transport propertiesof a small-scale unsaturated soil column
Paulo H. S. Moreira, Martinus Th. van Genuchten, Helcio R. B. Orlande*, Renato M. Cotta
DEM/PEM – Politécnica/COPPE, Federal University of Rio de Janeiro, UFRJ, Caixa Postal 68503, Cidade Universitária, Rio de Janeiro,
21941-972, Brazil.
*
Corresponding author. E-mail: helcio@mecanica.coppe.ufrj.br

Abstract: In this study the hydraulic and solute transport properties of an unsaturated soil were estimated simultaneously
from a relatively simple small-scale laboratory column infiltration/outflow experiment. As governing equations we used
the Richards equation for variably saturated flow and a physical non-equilibrium dual-porosity type formulation forsolute transport. A Bayesian parameter estimation approach was used in which the unknown parameters were estimatedwith the Markov Chain Monte Carlo (MCMC) method through implementation of the Metropolis-Hastings algorithm.
Sensitivity coefficients were examined in order to determine the most meaningful measurements for identifying the
unknown hydraulic and transport parameters. Results obtained using the measured pressure head and soluteconcentration data collected during the unsaturated soil column experiment revealed the robustness of the proposed
approach.
Keywords: Markov Chain Monte Carlo method; Metropolis-Hastings algorithm; Hydraulic parameters; Transportparameters; Soil column experiment.
INTRODUCTION
Efficient methods for accurately estimating physical andgeochemical parameters of mathematical models for the analysis of soil, hydrologic and contaminant transport problems are
critically needed. A large number of increasingly sophisticatedpublic domain and commercial software packages are now
available for a broad range of applications, including coupledvadose zone-groundwater systems (Miller et al., 2013; Šimůnek
and Bradford, 2008). While a range of experimental techniques
are now available for measuring the physical properties of soils(Dane and Topp, 2002), many of the standard procedures are
based on assumptions that simplify the mathematical analysis,
such as one-dimensional flow, steady-state flow conditions
and/or homogeneous media properties. At the same time, more
complicated and often nonlinearly coupled physical and biogeochemical processes are being included in the models for relevant applications. Therefore, proper formulations and solutionsof the associated parameter estimation problems (Beck and
Arnold, 1977) are needed to provide estimates for the required
subsurface media and contaminant transport properties.
Many methods and techniques are available for estimating
flow and transport soil properties. In this study a Bayesianestimator is used to identify the conditional probability distribution of the unknown properties given the measurements (i.e.,
the posterior probability density). The measurement model thatincorporates the measurement errors and related uncertainties iscalled the likelihood (i.e., the conditional probability of themeasurements given the parameters). The model for the unknowns that reflects all of the uncertainty in the parameters
without the information conveyed by the measurements iscalled the prior model (Beck and Arnold, 1977; Kaipio andSomersalo, 2004; Lee, 2004; Tan et al., 2006). The formal
mechanism to combine new information (measurements) withpreviously available information (prior) is known as Bayes’
theorem (Beck and Arnold, 1977; Kaipio and Somersalo, 2004;
Lee, 2004; Tan et al., 2006). Therefore, the term Bayesian isoften used to describe statistical inversion approaches based on

30

the principles that (Kaipio and Somersalo, 2004): (i) all variables in the model are treated as random variables, (ii) the randomness describes the degree of information concerning theirrealization, which is coded in probability distributions, and (iii)the solution of the inverse problem is the posterior probability
distribution from which point estimates and other statistics arecomputed.
This study aims at characterizing the unsaturated hydraulicand solute transport properties of variably saturated soils. Similar studies have been carried out in the past using both smallscale laboratory soil column experiments (Abbaspour et al.,2001, 2004; Kool et al., 1985; Parker et al., 1985) as well as formore elaborate field scale settings (Hosseini et al., 2011; Kohne
et al., 2006; Si and Kachanoski, 2000; Vrugt et al., 2006), in
which minimization of the maximum likelihood objective function was used for estimating the unknown parameters. Various
studies have used stochastic methods for studying fluid flow
and contaminant transport processes using synthetic experimental data (Franssen et al., 2003a, Fu and Gómez-Fernandez,2009), including particle filters for real-time estimations(Franssen et al., 2003b; Li et al., 2012; Vrugt et al., 2013; Xu etal., 2013; Zhou et al., 2011). Additionally, Bayesian estimators
have been used for estimating unsaturated soil hydraulic parameters (Diamantopoulos et al., 2012; Laloy et al., 2010).
The present work differs from other studies in the literatureby considering the estimation of soil hydraulic and transport
properties by means of a Bayesian estimator using observed
data of the pressure head and outflow solute concentration
obtained from a small-scale unsaturated soil column experiment. Selected unsaturated soil hydraulic and transport parameters affecting fluid flow and contaminant transport are estimated with a Markov Chain Monte Carlo (MCMC) method
through implementation of the Metropolis-Hastings algorithm(Hastings, 1970; Kaipio and Somersalo, 2004; Lee, 2004;
Metropolis et al., 1953; Tan et al., 2006). The Bayesian approach permits one to include a priori information available forthe parameters, such as from previous experimental runs orother experimental setups, or even from theoretical predictions.

Bayesian estimation of the hydraulic and solute transport properties of a small-scale unsaturated soil column

The flow problem in our study is described in terms of thestandard Richards equation assuming uniform variably-saturated
flow model in a single-porosity system, and a two-region (dualporosity) type physical non-equilibrium model for transport(Šimůnek et al., 2003; van Genuchten and Wierenga, 1976).
The present study extends earlier work by Orlande et al.
(2009) who had used only simulated data in the inverse analysis. For the results presented here, actual pressure head andsolute concentration experimental measurements were collected
in a small-scale laboratory soil column experiment and used toestimate the unsaturated soil hydraulic and transport parameters. We analyzed the sensitivity coefficients in order to examine which parameters can be estimated simultaneously. The
accuracy of the Bayesian approach was also evaluated by reproducing the estimated unsaturated soil hydraulic parameters
obtained in the earlier studies by Kool et al. (1985) and Parker
et al. (1985) involving the use of cumulative flux measurements. Below we present the mathematical formulation for thephysical problem considered here, the solution of the inverse
problem, a description of the experimental procedure and theparameter estimation results.
PHYSICAL PROBLEM AND MATHEMATICALFORMULATION
The physical problem considered here concerns the flow ofwater and a dissolved tracer through a vertical unsaturated laboratory soil column of length L, as shown in Figure 1. At the initialtime (t = 0), the pressure head, hi(z), inside the soil column isassumed to vary linearly as a function of the position, z, in thecolumn (but not necessarily following equilibrium conditions).
For t > 0, a tracer solution with known concentration, c0, is injected at the upper end of the column at a constant flow rate, q0. The
column is connected at its bottom to a vacuum chamber withpressure control, thus making it possible to keep a constant pressure head, hb, at the outflow boundary. Values of the pressure
head and solute concentration are assumed to be uniform in eachcross section.

The problem is simulated by assuming uniform flow as
described by the standard Richards equation as follows:∂θ∂  ∂h=+KK∂t ∂ z  ∂ z
in 0 ≤ z ≤ L and t > 0(1)where θ is the volumetric water content (L3L–3), h is the pressure head (L), K is the hydraulic conductivity (L T – 1 ), z is the
distance from the bottom boundary upwards (L), and L is thelength of the column (L).
Equation (1) was solved subjected to the following boundary
conditions:h = hb at z = 0, for t > 0

(2a)

 ∂h

− K
+ K  = q0 at z = L, for t > 0
 ∂z


(2b)

and the initial condition:h = hi(z) in 0 < z < L, at t = 0(3)which varies linearly between the values
hi(0) = h0 – Δh at z = 0(4a)hi(L) = h0 +Δh at z = L(4b)where h0 is the initial pressure head in the middle of the column
and Δh some positive value of the pressure head difference
between the middle of the column and its boundaries.
The relation between θ and h, as well as the hydraulic conductivity function, K, are described here by (van Genuchten,1980):

θ ( h) = θ r +θs − θr1 + α h n 



(

(5)

m

)

m
K ( Se ) = K s Sel 1 − 1 − Se1/ m 



2

(6)

where θr and θs are the residual and saturated water contents(L3L–3), respectively, Ks is the hydraulic conductivity at saturation (L T – 1 ), α (L–1) and n (–) are empirical shape parameters,
m = 1–1/n, l = 0.5 is a pore connectivity parameter, and Se iseffective saturation given bySe =θ ( h) − θ r
θs − θr(7)

A dual-porosity (two-region) type physical non-equilibrium
model was used for the solute transport problem. By assuming
linear equilibrium sorption, the equations governing transport in avariably-unsaturated soil are given by (van Genuchten and
Wierenga, 1976):
Fig. 1. Schematic of the one-dimensional flow and solute transport
problem considered in this study.

∂c∂[(θ m + ρb fK d )cm ] + [θim + ρb (1 − f ) K d ] im =
∂t
∂t
∂cm
∂ 

− qcm θ m Dm
∂z 
∂z


(8a)

31

Paulo H. S. Moreira, Martinus Th. van Genuchten, Helcio R. B. Orlande, Renato M. Cotta[θ im + ρb (1 − f ) K d ]
t>0

∂cim
= α m [cm − cim ]∂t

in 0 ≤ z ≤ L and(8b)where the subscripts m and im refer to mobile and immobile
liquid regions, respectively, c is the solute concentration (ML–3),
f defines the fraction of sorption sites in contact with mobile
water (–), ρb is the soil bulk density (ML–3), Kd is the
distribution coefficient (L–3M), Dm is the dispersion coefficient(L2T – 1) and αm is a mass transfer coefficient determining thesolute transfer rate between the mobile and immobile regions(T – 1 ). The dual-porosity formulation was used here since many
tropical soils (often oxisols) of the type used in this study
exhibit multi-modal pore- or particle-size distributions (e.g.,
Carducci et al., 2011; Sommer et al., 2003; Spohrer et al., 2006)
favoring the existence of relatively immobile water pockets(Nkedi-Kizza et al., 1984; Melamed et al., 1994).
A third-type condition was used at the inflow boundary
and a zero diffusive flux at the outflow boundary, that is,∂cm= 0 at z = 0 for t > 0∂z−θ m Dm∂ cm+ q0 cm = q0 c0 at z = L for t > 0
∂z

(9a)

(9b)

The initial condition is of the form:m = im =0 in

0 < z < L for t = 0(9c)

The dispersion coefficient Dm was defined following Bear’sclassical formulation with the tortuosity given by Millington
and Quirk's model (Kohne et al., 2006) as follows:θ m Dm = Ddθ m10/3
+λ q
θ s2

(10)

where Dd is the molecular diffusion coefficient (L 2 T – 1 ), λ isthe dispersivity (L), and q is the volumetric fluid flux (L T – 1 )given by the Darcy-Buckingham law as follows: ∂h 
q = − K ( h)  + 1  ∂z 

(11)

Following Nkedi-Kizza et al. (1984), we assumed that thefraction of sorption sites in contact with mobile water (f) is thesame as the ratio of the mobile water content (θm) to the total
water content (θ), that is:

θ
f = m
θ

(12a)

with

θ = θ m + θim(12b)

By using Eqs. (12), Eqs. (8a, b) can be rewritten as:f∂∂  ∂c( Rθ cm ) =  D* m − qcm  − α m [cm − cim ]
∂t
∂z 
∂z


32

(13a)

(1 − f ) Rθ∂cim
= α m [cm − cim ]
∂t

(13b)

where

D* =

Dd f 10/3

θ s2

and R = 1 +θ 10/3 + λ q

ρb K d
θ

(13c)

(13d)

A unique feature of Eqs. (13a, b) is the fact that the fractionof immobile water (and hence also f) is assumed to be constantduring transient variably-saturated flow conditions, in which
the mobile and immobile water contents can vary in time andspace. We believe that this assumption is more realistic, andcertainly less restrictive, than having a constant immobile watercontent (θim) in time, as assumed for example in the HYDRUS
software package (Šimůnek et al., 2013), unless a dual-porosity
type water flow exchange term is used between the mobile andimmobile liquid regions. Limited experimental data (e.g., De
Smedt and Wierenga, 1984; Maraqa et al., 1997) also support
the assumption of having a constant immobile water fraction, f,
independent of the degree of saturation. An additional advantage of the dual-porosity formulation used here is that themodel simplifies immediately to the standard advectiondispersion equation when f = 1. Hence, no a priori assumption
is needed as to which formulation (the standard advection dispersion equation or the dual-porosity model) is most optimalfor the transport problem being considered.
The forward problem, given by Equations (1) to (7), (9) and(13), with all model parameters, initial and boundary conditions
known, was solved using an implicit finite volume discretization approach with the WUDS (Raithby and Torrance, 1974)
interpolation scheme. The numerical solution of the forward
problem was implemented in a Fortran code, which was verified against several reference solutions for different test-cases.
Results presented below were obtained with the code developedin this work, by using converged grids.
PARAMETER ESTIMATIONThe combined unsaturated flow and solute transport problem
defined above involves the following vector of parameters:
T
T
PT = [PHydraulic, PTransport
]

(14a)

which contains separate vectors of hydraulic and solute transportparameters as considered in our analysis:TPHydraulic = [ K s ,θ s ,θ r ,α , n , h0 , Δh ]

(14b)

T
PTransport= [ f , ( ρb K d ), Dd , λ , α m ](14c)Note that the parameters ρb and Kd appear only in the formof a product in the transport problem, and hence cannot betreated independently. However, ρb can be measured with relatively standard techniques, which would leave only Kd as one of
the unknowns in the vector of parameters.
Parameter estimation within the Bayesian framework is performed by treating all parameters in the mathematical formulation of the physical problem as random variables, with prior

Bayesian estimation of the hydraulic and solute transport properties of a small-scale unsaturated soil column

information coded in the form of probability distributions(Kaipio and Somersalo, 2004). Consider the vector of parameters appearing in the physical model formulation asPT ≡ [P1, P2, ..., PN], where N is the number of parameters.
Consider also that a vector of measurements Y is available for
estimation of the parameters P. Bayes’ theorem can then be
stated as (Kaipio and Somersalo, 2004; Lee, 2004; Tan et al.,2006):

π posterior (P) = π(P Y) =π(P)π(Y P)

(15)

π(Y)

where πposterior(P) is the posterior probability density, π(P) is theprior density, π(Y|P) is the likelihood function, and π(Y) is themarginal probability density of the measurements, which is anormalizing constant.
In our study we assumed that the measurement errors are
Gaussian random variables, with zero means and constant variance, and that the measurement errors are additive and independent of the unknowns. With these assumptions, the likelihood function can be expressed as (Kaipio and Somersalo,
2004; Lee, 2004; Tan et al., 2006):
|

= 2

⁄

| |

⁄

−

−

−
(16)

where M is the number of measurements, Σ is the covariancematrix of the measurement errors and T(P) is the solution of thedirect problem obtained with the vector of parameters P. We
note that parameters estimated through maximization of Equation (16) are referred to as maximum likelihood estimates (Beckand Arnold, 1977).
The unknown parameters were estimated by using the Metropolis-Hastings algorithm for the Markov Chain Monte Carlo(MCMC) method (Hastings, 1970; Kaipio and Somersalo,
2004; Lee, 2004; Metropolis et al., 1953; Tan et al., 2006).
Implementation of the Metropolis-Hastings algorithm starts
with the selection of a proposal distribution p(P*,P(t–1)), which
is used to draw a new candidate state P*, given the current stateP(t–1) of the Markov chain. Once this moving distribution has
been selected, the Metropolis-Hastings sampling algorithm canbe implemented by repeating the following steps:
1. Sample a Candidate Point P* from the proposal distributionp(P*,P(t-1)).
2. Calculate the acceptance factor:π(P* | Y) p(P (t −1) , P* ) (t −1)| Y) p (P* , P (t −1) )  π(Pα = min 1,(17)

3. Generate a random value U which is uniformly distributedon (0, 1).
tt
4. If U ≤ α, set P = P*. Otherwise, set P = P(t-1).
5. Return to step 1 in order to generate the sequence12n{P , P , …, P }.
For more details on the theoretical aspects of the MetropolisHastings algorithm and MCMC methods, readers are referred to
Kaipio and Somersalo (2004), Lee (2004) and Tan et al. (2006).
Proposal densities (Kaipio and Somersalo, 2004; Lee, 2004;
Tan et al., 2006) for implementation of the Metropolis-Hastings
algorithm were taken in this work as random walks, with
Gaussian perturbations for the parameters with Gaussian priors,
and with uniform perturbations for the parameters having uniform priors. Gaussian perturbations had zero means andstandard deviations of 0.01% of the prior means. Uniform
perturbations involved maximum steps of 0.5% of the parameter range established for the uniform prior. These proposal
densities were selected based on numerical experiments.
Estimation of the unknown parameters must be preceded byan analysis of the sensitivity coefficients (Beck and Arnold,
1977; Özisik and Orlande, 2000). The sensitivity matrix isdefined as:
 ∂T1
 ∂P
 1
 ∂T2
T
 ∂ TT (P ) J (P) = 
 =  ∂P1
P
∂


 

 ∂TM
 ∂P
 1

∂T1
∂P2

∂T1
∂P3

∂T2
∂P2

∂T2
∂P3


∂TM
∂P2


∂TM
∂P3

where the sensitivity coefficients J ij =

∂T1 
∂PN 

∂T2 


∂PN  (18)
 

∂TM 

∂PN 
∂Tigive the sensitivity of∂Pj

the response of problem Ti (i.e., the dependent variable being thesolution of the direct problem) with respect to changes in theparameter Pj. A small value of Jij indicates that large changes in Pj
yield small changes in Ti. Estimation of parameter Pj is typicallyvery difficult in that case since very similar values of Ti would beobtained for a wide range of Pj values. When the sensitivity
coefficients are small, J T J ≈ 0 and the inverse problem is illconditioned. Furthermore,

JT J

is null when a column of J

becomes a linear combination of the other columns (Beck and
Arnold, 1977; Özisik and Orlande, 2000). It is therefore
desirable to have linearly independent sensitivity coefficients Jijwith large magnitudes so that the parameter estimation problem is
not very sensitive to measurement errors, and hence accurate
parameter estimates can be obtained. A comparison of themagnitude of the sensitivity coefficients, as well as an analysis ofpossible linear dependencies, is more easily performed by using
reduced sensitivity coefficients instead of the original ones. The
reduced sensitivity coefficients are obtained by multiplying theoriginal sensitivity coefficients, Jij, by the parameters that theyrefer to. They have the same units as the measured variables, thus
providing a basis for the comparisons. Since the sensitivity
coefficients are functions of the model parameters in thepresent case, the parameter estimation problem is nonlinear andonly a local analysis of the sensitivity coefficients is thus
possible (Beck and Arnold, 1977; Özisik and Orlande, 2000).
Sensitivity coefficients in our study were calculated using
forward finite differences; convergence of this first orderapproximation was verified for all cases analyzed.
EXPERIMENTAL PROCEDURE
Laboratory soil column inflow/outflow experiments werecarried out in order to obtain the data needed to estimate theunknown soil hydraulic and solute transport parameters ofconcern in our study, given by Eqs. (14). The apparatus (Figure 2)
used for the experiment consisted mostly of equipment obtained
from Soil Measurement Systems (Tucson, AZ), designed forunsaturated flow experiments. We used an acrylic column with aninternal diameter of 5.1 cm and length of 7.3 cm, packed with an
undisturbed Brazilian Latosol soil sample, taken from a study area
near Caetité, BA, Brazil. The bottom of the column was connected to a sealed vacuum chamber, which in turn was connected33

Paulo H. S. Moreira, Martinus Th. van Genuchten, Helcio R. B. Orlande, Renato M. Cotta

Fig. 2. Schematic of the experimental apparatus.

to a vacuum pump. The pressure inside the vacuum chamber wasmeasured with a precision mercury U-tube manometer attached tothe chamber. The pressure was controlled by a vacuum regulator,
which allowed us to maintain a constant negative pressure head at
the outflow boundary of the soil column.
Using a high-accuracy syringe pump, a potassium bromide(KBr) solution with a concentration of 2.5 g/L was injected from areservoir to the top of the column, at a constant rate of
3.25 x 10–4 cm/s. Variations in the pressure head in the center ofthe vertical column, at z = 3.5 cm, were measured with atensiometer attached to a pressure transducer and connected to adata logger and computer. A fraction collector, placed inside thevacuum chamber and containing a large number of samplingtubes, was used to collect the outflow fluid, thus allowing us toobtain a series of outflow rates and associated fluid
concentrations.
The tensiometer halfway down the column consisted of a
10-mm diameter, 60-mm long acrylic tube. One of its ends was
attached to a 6-mm diameter ceramic bar in direct contact withsoil inside the column, while the other end was connected to apressure transducer. The tensiometer was carefully installed suchthat outside air could not enter the column. Since transientpressure heads inside the column were recorded in terms of anelectrical signal (mV), the pressure transducer needed to becalibrated first. For this purpose, both the pressure transducer andthe mercury manometer were connected to the vacuum chamberso that it was possible to record the voltages corresponding to aseries of pressure head values.
KBr solute concentrations were assumed to be linearly related
to the electric conductivity (μS/cm) of the outflow solution. This
assumption was verified by measuring the electric conductivity of
KBr solutions within the expected range of concentrations. It wasthen possible to immediately convert electric conductivities
measured for each outflow sample, as well as of the inflow
solution, into reduced solute concentrations, c(t)/c0, with c0 being
the concentration of the injected solution. Based on values
obtained during calibration, the standard deviations were 2.0
cm for the pressure head and 0.01 (1%) for the normalized
solute concentrations.
Relatively standard experiments (Dane and Topp, 2002) werecarried out to obtain prior information about the soil bulk density,
ρb, the saturated hydraulic conductivity, Ks, and the saturated
water content, θs. Mean values obtained with these measurements
are listed in Table 1. Standard deviations of the measured data
were approximately 1% of the mean values. The priordistributions for ρb, Ks, and θs were for these reasons assumed tohave Gaussian distributions, centered at the mean values
presented in Table 1 and with standard deviations of 1.0% of themean values. Values of the residual water content θr, and the34water retention parameters α and n, were estimated from themeasured soil textural distribution and the bulk density using the
Rosetta pedotransfer functions (Schaap et al., 2001). Estimatedvalues for these parameters are also shown in Table 1.
Table 1. Properties of the soil used in this study.
Parameterθr (cm3/cm3)
θs (cm3/cm3)

ρb (g/cm3)
Ks (cm/s)
α (cm–1)
n (–)

Value
0.042
0.408
1.637
9.37 x 10–4
0.040
2.632

RESULTS AND DISCUSSIONS
Before analysis of the infiltration/outflow parameter estimation problem considered in this study, we evaluated the implemented MCMC method by estimating unsaturated soil hydraulic parameters obtained with the one-step column outflow
experiment described by Kool et al. (1985) and Parker et al.
(1985). That experiment considered an undisturbed soil core of
5.40 cm diameter and 4.00 cm length (Parker et al., 1985). For
times t > 0, a constant pressure head of 10 m was applied at the
top of the column, causing water stored in the soil column tomove out of the medium, to be collected at the bottom of thecolumn, which was maintained at atmospheric pressure (Parker
et al., 1985). For the simulations and analysis of the sensitivity
coefficients, the spatial domain was defined as a 3.95 cm longsoil core placed over a 0.55 cm thick ceramic plate, resulting in
a L = 4.50 cm long column consisting of two materials having
different physical properties (the soil and the ceramic plate).
The inverse problem dealt with cumulative outflow flux simulated measurements for the estimation of hydraulic parametersof the soil. In order to avoid an “inverse crime”, as explained by
Kaipio and Somersalo (2007), the one-step outflow experiment
was simulated with the HYDRUS-1D code (Šimůnek et al.,2013) to generate the simulated data using a different meshfrom the finite volume numerical solution we used for the forward problem. The finite volume solution for the forward
problem was used for analysis of the sensitivity coefficients
based on the parameter values shown in Table 2 (Kool et al.,
1985; Parker et al., 1985), as well as for the parameter estimation problem using the simulated measurements. The simulated
measurements contained additive, uncorrelated Gaussian errors,with zero mean and a constant standard deviation of 2% of the
maximum cumulative flux.

Bayesian estimation of the hydraulic and solute transport properties of a small-scale unsaturated soil column
Table 2. Parameter values used for analysis of the sensitivity coefficients in the experiment of Parker et al. (1985).
Parameter
Residual water content of the soil core
Saturated water content of the soil core
Saturated hydraulic conductivity of the soil corevan Genuchten α parameter of the soil corevan Genuchten n parameter of the soil core
Residual water content of the ceramic plate
Saturated water content of the ceramic plate
Saturated hydraulic conductivity of the ceramic plate
van Genuchten α parameter of the ceramic plate
van Genuchten n parameter of the ceramic plate

Variations in the cumulative outflow, as well as in the
reduced sensitivity coefficients, are presented in Figure 3.
Figure 3a shows that the sensitivity coefficients had verysmall values compared to the variations in the cumulative
outflow. Figure 3b provides a close-up of the reduced
sensitivity coefficients versus time. The variations areapproximately three orders of magnitude smaller than the
cumulative outflow rates, meaning that variations in thehydraulic parameters have very little effect on the cumulative
outflow. The results also indicate a linear dependence between
the parameters θr and θs, as well as among the parameters Ks, αand n. Since the values of θs and Ks were directly measured byParker et al. (1985), we defined their a priori distributions to be
Gaussian, centered on the values obtained for each parameter.
Nevertheless, the parameters α and n were found to be linearly
dependent and their curves for the sensitivity coefficients werecoincident. This linear dependency made it necessary to fix aninformative a priori distribution for either α and n in order tomake the estimation possible. Van
Genuchten(1978)
previously showed a graphical procedure to obtain thecoefficient n from observed water content data, independently
from α. Other studies (e.g., Leij et al., 1997; van Genuchten andNielsen, 1985; Yates et al., 1992) showed that for a given soil
the values of α and n generally can be obtained independently.
For these reasons, we assumed a Gaussian a priori distribution
for n, centered on the value obtained by Parker et al. (1985).
Since θs, Ks, and n could be considered as known parameters,
their prior distributions were defined with a small standard
deviation. For the parameters θr and α, the a priori conditions
were defined as uniform distributions defined from the initialvalues used in Parker et al. (1985). All of the prior distributions
are listed in Table 3. The parameter values estimated with the
MCMC method are presented in Table 4, including theirrespective confidence intervals, and compared to the values
obtained by Parker et al. (1985). Notice that the values obtained
by Parker et al. (1985) are inside the confidence intervals
obtained in our study. Figure 4 shows a comparison between
the experimental cumulative outflow curve and the outflow
curve obtained with the estimated MCMC parameter values.
Results indicate excellent agreement between the experimental
and estimated curves, with small and practically uncorrelated
residuals (difference between experimental and estimatedvalues).
We next focus on the inverse analysis using the experimental
setup described earlier (Fig. 2). The following measurements
were used to estimate the sensitivity coefficients: (i) transientvariations in the pressure head halfway down the column;(ii) volumetric fluid fluxes at the outlet of the column; andValueθr,soil = 0.161
θs,soil = 0.388
Ks,soil = 1.5 x 10–3 cm/s

αsoil= 0.033 cm–1
nsoil = 1.387

θr,plate = 0.0
θs,plate = 1.0
Ks,plate = 8.7 x 10–7 cm/s

αplate= 1.0 x 10–20 cm–1
nplate = 1.1(iii) solute concentrations of the outflowing fluid. For the
analysis of the sensitivity coefficients with respect to theunknown hydraulic and solute transport parameters, we usedthe values presented in Table 5. These values were selectedbased on direct measurements (see Table 1) and on expectedvalues of the parameters.

(a)

(b)
Fig 3. Variations in the cumulative outflow flux and its sensitivity
coefficients (a), and close-up of the variations in the sensitivity
coefficients with respect to the hydraulic parameters (b) in theexperiment of Parker et al. (1985).

35

Paulo H. S. Moreira, Martinus Th. van Genuchten, Helcio R. B. Orlande, Renato M. Cotta(a)Fig. 4. Variations in time of the simulated measurements, estimatedcumulative flux and the residuals of the experiment by Parker et al.
(1985).
Table 3. Prior distributions considered for the MCMC evaluationof the experiment of Parker et al. (1985).
Parameterθr

A priori
Distributions(MCMC)
UniformθsGaussian

Distribution data
Interval: (0.6 ̅ r , 1.4 ̅ r)
̅ r = 0.200
Mean: θ s = 0.388
Standard Deviation: 0.01 θ s

Ks

Gaussianα

Uniform

Interval: (0.6α , 1.4 ) cm–1
α = 2.5 x 10–2cm–1nGaussian

Mean: n =1.387
Standard Deviation: 0.01 n

Mean: K s = 1.5 x 10–3 cm/s
Standard Deviation: 0.01 K s

Figure 5a shows transient variations in the pressure head in
the middle of the soil column and the reduced sensitivity
coefficients with respect to all the hydraulic and solutetransport parameters, using the values presented in Table 5.
Note that the reduced sensitivity coefficients are allpractically zero relative to the pressure head values. The
reduced sensitivity coefficients with respect to the hydraulicparameters are about four orders of magnitude smaller than thepressure head, as illustrated in Figure 5b, which is amagnification of Figure 5a within the region of interest foranalysis of the sensitivity coefficients. An analysis of Figure 5b
reveals a linear dependence among the sensitivity coefficientswith respect to Ks, α and n, as well as between θr and θs. The
sensitivity coefficients with respect to h0 and h werepractically zero. Figure 5b also shows that the sensitivity
coefficients become constant when steady state conditions areapproached. For the transport parameters, the sensitivity
coefficients were zero, as shown in Figure 5c. This was expected
since the transport parameters do not appear in the formulationof the flow problem given by Equations (1)–(4). The reverse
situation, however, is the case in that solute transport
information can improve the estimation of soil hydraulicparameters.

36

(b)

(c)
Fig. 5. Variations in the pressure head in the middle of the column
and in its sensitivity coefficients (a), and close-ups of variations in thesensitivity coefficients with respect to the hydraulic (b) and transport(c) parameters.

Bayesian estimation of the hydraulic and solute transport properties of a small-scale unsaturated soil column
Table 4. Evaluation of the MCMC analysis of the experiment of Parker et al. (1985).
Parameter

Values estimated by Parker et al. (1985)

Estimated values using MCMC

θrθs
Ks(cm/s)α (cm–1)
n

0.161
0.388
1.500 x 10–3
3.350 x 10–2
1.387

0.163
0.388
1.503 x 10–3
2.990 x 10–2
1.388

99% confidence intervals(MCMC)(0.162 , 0.164)(0.385 , 0.391)(1.500 x 10–3 , 1.506 x 10–3)(2.90 x 10–2 , 3.10 x 10–2)(1.386 , 1.390)

(a)

(a)

(b)

(b)

(c)
Fig. 6. Variations in the fluid flux exiting the column and its
sensitivity coefficients (a), and close-ups of variations in thesensitivity coefficients with respect to the hydraulic (b) and transport(c) parameters.

(c)Fig. 7. Variations in the solute concentration of the fluid exiting thesoil columns and its sensitivity coefficients (a), and magnificationsof variations in the sensitivity coefficients with respect to thehydraulic (b) and transport (c) parameters.

37

Paulo H. S. Moreira, Martinus Th. van Genuchten, Helcio R. B. Orlande, Renato M. Cotta

Variations in the volumetric fluid flux exiting the column, as
well as its sensitivity coefficients, are presented in Figure 6a.
The parameter values listed in Table 5 were also used in thisanalysis. Similarly to the pressure head measurements
examined in Figures 5, the sensitivity coefficients of the fluid
flux are practically zero. Figure 6b shows that the sensitivity
coefficients of the hydraulic parameters are again four orders ofmagnitude smaller than the measured values for our test case.
Table 5. Parameter values used for the analysis of the sensitivity
coefficients.
Parameter

Value

Column height

L = 7.3 cm

Infiltration flux

q0 = 3.050 x 10–4 cm/s

Head at the bottom of the column

hb = –60 cm

Constant initial head in the column

hi = –60 cm

Residual water content

θr = 0.042

Saturated water content

θs = 0.408

Saturated hydraulic conductivity
van Genuchten α parameter

Ks = 9.730 x 10–4 cm/s

α = 0.040 cm–1

van Genuchten n parameter
Fraction of sorption sites in contact
with mobile water
Soil bulk density

ρb = 1.637 g/cm3

Distribution coefficient

Kd = 1.0 cm3/g

Dispersivity

n = 2.632
f = 0.500

λ = 1.0 cm

Distribution coefficient

Dd = 1.0 x 10–5 cm2/s

Mass transfer coefficient

αm = 1.0 x 10–5 s–1

Normalized initial concentration
Normalized concentration of theinjected tracer

ci = 0.0
c0 = 1.0

Additionally, the sensitivity coefficients with respect to thehydraulic parameters are all linearly dependent, while thosewith respect to the transport parameters are all zero, as can beobserved in Figure 6c. This indicates that measurements of thefluid flux did not contribute much useful information for estimation of the unknown hydraulic and transport parameters.
Variations in the solute concentrations of the outflow, as
well as its sensitivity coefficients, are presented in Figure 7a forthe parameter values specified in Table 5. Once again, thereduced sensitivity coefficients are four orders of magnitude
smaller than the solute concentrations. The sensitivity
coefficients with respect to the hydraulic parameters in
Figure 7b show that α and n are linearly dependent, as well as
θs and Ks. The sensitivity coefficients are zero for h0 and Δh.
Figure 7c further shows that the transport parameters f, (ρbKd), λ1–0.1622 –0.1241 –0.0395–0.16221
0.5484 0.6076–0.1241 0.5484
1
0.9169cor θr ,θs ,ks ,α,n,h0 ,h = –0.0395 0.6076 0.9169
1
0.2029 –0.3273 –0.8413 –0.5649–0.1136 0.3412 0.1271 0.4137–0.0018 –0.1045 –0.3321 –0.2272

38

and n are linearly independent, while the sensitivity coefficientis zero for the dispersion coefficient, Dd. The latter is due to thedominant effects of hydrodynamic dispersion in the dispersion
coefficient as compared to molecular diffusion (Equation 10).
The above analysis revealed that the parameter estimation
problem considered in this study is difficult due to small values
and linear dependencies of the sensitivity coefficients. Only
measurements of the pressure head in the middle of the column
and of the outflow solute concentration provided meaningful
information for estimation of several parameters.
The Bayesian estimation approach is capable of dealing withthese difficulties by using informative prior distributions for theparameters that could be independently measured or estimated.
The values of θs and Ks were experimentally measured, as discussed above. The initial pressure head at the observation position point, h0, was read by the tensiometer situated in the soilcolumn, while the pressure head difference, Δh, was obtained
from the initial pressures at the surfaces of the column. It
was thus possible to estimate the vector of hydraulic parameters assuming uniform or Gaussian prior distributions, according to the quality of the a priori information available for eachparameter. Since θs, Ks, h0 and Δh were directly measured,
Gaussian prior distributions were adopted for these parameters.
For θr, α and n, prior uniform distributions were used. For Dd,
with practically null sensitivity coefficient, a Gaussian priordistribution was defined based on the known range of effective
diffusion coefficients for unsaturated soils (e.g., Shackelford,1991).
The estimation process was carried out in two steps. In the
first
step,
the
hydraulicparameters
PTHydraulic = θr , θs , Ks , α, n, h0 , h were estimated by using
the measured data of the pressure head in the middle of thecolumn. In the second step we used estimates for the flow parameters as prior information for the coupled flow-transport
problem
in
which
the
transport
parameters
T
PTransport = [ f , ( ρb K d ), Dd , λ , α m ] were estimated from themeasured concentration data of the outflow and the pressure
heads in the middle of the column. Since the hydraulic parameters were treated as random in the second step, they were also
allowed to vary in the Markov chain, in accordance with theirprior distributions.
The prior distributions used for the hydraulic parameters in
the first step are presented in Table 6. In order to estimate thehydraulic parameters, pressure head measurements in the middleof the column were considered until a time equal to 3300 s for thesolution of the inverse problem, with experimental errors assumed
to have additive, uncorrelated Gaussian distributions with zeromean and standard deviation corresponding to 5% of the difference between the minimum and maximum measured values(2.0 cm, for our example). Table 7 lists the values estimated foreach hydraulic parameter, while Equation (19) presents the obtained correlation matrix, where strong correlations among someparameters can be observed, as expected.
0.2029 –0.1136–0.3273 0.3412–0.8413 0.1271–0.5649 0.4137
1
0.33310.33311
0.3979 0.1959–0.0018–0.1045–0.3321–0.2272
0.39790.19591(19)

Bayesian estimation of the hydraulic and solute transport properties of a small-scale unsaturated soil column
Table 6. Prior distributions of the hydraulic parameters used in thefirst step.
Parameter

Distribution

Distribution information

θr

Uniform

Interval: (0 , θs )

θsGaussian

Mean: θs = 0.408 (cm3/cm3)
Standard deviation: 0.01 θs (cm3/cm3)

Ks

Gaussianα

Uniform

Interval: (0.8α , 1.2α) cm–1
α = 0.018 cm–1n

Uniform

Interval: (0.8 , 1.2 )n = 1.29

h0GaussianΔhGaussian

Mean: = –66.0cm
Standard deviation: 0.05cm
Mean: = 3.5cm
Standard deviation: 0.05cm

Mean: K s = 7.0 x 10-4 cm/s
Standard deviation: 0.2 K s(a)Table 7. Estimated values of the unsaturated soil hydraulicparameters.
Parameter
3

Mean
3

99% confidence intervalθ r (cm /cm )0.118 x 10

θ s (cm3/cm3)

0.409

K s (cm/s)

α (cm )
–1

1.07 x 10–2(0 , 0.366 x10–2)(0.408 , 0.410)–3(0.956 x 10–3 , 1.196 x10–3)–2(2.009 x 10–2 , 2.204 x10–2)

2.107 x 10

n (–)

1.236

(1.228 , 1.243)

h0 (cm)–59.308(–59.536 , –59.081)

Δh(cm)

3.523(3.516 , 3.530)

With the mean values of the parameters as presented inTable 7, it was possible to simulate the variations in the pressure
head halfway down the column using Equations (1)–(4). Figure 8ashows a comparison of the experimental and simulated curves forthe pressure head inside the soil column. The simulated curve
generally falls within the range of the experimental errors, except
at the very beginning of the experiment when the effects of theestimation of the initial pressure head and its gradient are moresignificant. Figure 8b further shows that the pressure headresiduals (i.e. the difference between experimental andestimated values), although small, are correlated, which indicatesthat the mathematical model used to describe the fluid flow
problem did not fully represent the physics of the problem.
Figure 9 presents the number of accepted states and values of
the hydraulic parameters at each state of the Markov chain. The
Markov chain was generated with 50,000 states, considering thelast 30,000 states for computation of the statistics. The
convergence of the chain was tested by calculating the meanvalues of the parameters in different ranges inside the interval
considered for the calculations (Geweke, 1992). The chain
converged for all parameters, despite the linear dependence ofθs and Ks, with respect to θr and n, respectively. The acceptanceratio of the Markov chain was approximately 14%.
In the second step, the measured pressure head and soluteconcentration data were used to estimate the hydraulic andtransport parameters of the combined vector of parameters
TT
PT = [PHydraulic, PTransport] . The estimated values of the(b)Fig. 8. Comparison of experimental and simulated curves for thepressure head halfway down the column (a) and the residuals (b).

hydraulic parameters of step 1 were then used as prior information with Gaussian distributions centered at the mean values
presented in Table 7 and with the covariance resulting from theestimation in step 1. The prior distributions used for thetransport parameters are defined in Table 8. The value of the
product ρbKd of the distribution coefficient was assumed to bevery small, that is, KBr was taken approximately as aconservative tracer (Equation 13d), which is characteristic ofsoils such as the one used in this work at neutral pH (Goldberg
and Kabengi, 2010).
The estimated means of the marginal posterior distributions
for the transport parameters are presented in Table 9, together
with their 99% confidence intervals. With the mean values forthe parameters it was possible to compare the experimental andtheoretical solute concentration curves of the column outflow,
as seen in Figure 10a. Excellent agreement was obtained for thetwo curves. Although correlated, the concentration residuals areof the same order as the measurement uncertainties, as shownby Figure 10b.

39

Paulo H. S. Moreira, Martinus Th. van Genuchten, Helcio R. B. Orlande, Renato M. Cotta

(a)

(e)

(b)

(f)

(c)

(g)

(d)
(h)
Fig. 9. States of the Markov chain for the unsaturated soil hydraulic parameters: (a) accepted states, (b) θr, (c) θs, (d) Ks, (e) α , (f) n, (g) h0,
and (h) Δh.

40

Bayesian estimation of the hydraulic and solute transport properties of a small-scale unsaturated soil column

(a)

(b)

Fig. 10. Comparison of experimental and simulated concentrations of the outflow flux (a) and the residuals (b).

(a)

(d)

(b)

(e)

(c)

(f)

Fig. 11. States of the Markov chain for the solute transport parameters: (a) accepted states, (b) f, (c) ρbKd, (d) λ, (e) Dd , and (f) αm.

41

Paulo H. S. Moreira, Martinus Th. van Genuchten, Helcio R. B. Orlande, Renato M. Cotta

The evolution of accepted states and the Markov chains forthe transport parameters are presented in Figure 11. A total of
25,000 states were considered for the Markov chains, with thelast 10,000 states used to obtain the mean values and standard
deviations of the marginal posterior distributions. The
convergence of the chains was verified by using the method
proposed by Geweke (1992), with the acceptance ratio
corresponding to about 15% of the states of the Markov chain.
Computational times for estimation of the hydraulic andtransport parameters were of the order of 5 hours each, on acomputer with an I5 2.27 GHz processor and 4.0 GByte of
RAM memory.
Table 8. Prior distributions of the transport parameters used in
the second step.
Parameter

Distribution

Distribution informationf

Uniform

Interval: (0 , 1)ρb Kd

Uniform

Interval: (0.8 ρb Kd , 1.2 ρb Kd )λ (cm)

Uniform

Dd (cm²/s)

Gaussianαm (s–1)

Uniformρb Kd = 1.637 x 10–11

Interval: (0.2 λ , 1.8 λ)λ = 7.0 x 10–1
Mean: Dd = 1.00 x 10–5
Standard Deviation: 0.01 Dd
Interval: (0.01 α m , 1)

α m = 1.00 x 10–2
Table 9. Estimated values of the transport parameters.
Parameter

Mean

f

0.535

ρb Kd

1.655 x 10

99% Confidence Interval(0.514 , 0.557)–11(1.639 , 1.671) x 10–11

λ (cm)

0.224(0.222 , 0.226)

Dd (cm2/s)

1.001 x 10–5

(9.998 x 10–6 , 1.003 x 10–5)

α m (s–1)

0.276(0.273 , 0.279)

Due to correlation of the residuals observed in Figures 8 and
10, a final analysis was performed in order to investigate thepossibility of the measurements themselves being correlated.
The parameters were again estimated by successively using
only half of the measurements, that is, by halving the frequencyof the original measurements. However, the estimatedparameters and the residuals obtained with 50% and 25% of theoriginal measurements were practically the same as thoseobserved when 100% of the measurements were used.
Therefore, the correlated residuals can be attributed to thestrong linear dependency of the physical parameters shown in
Figures 5–7, as well as to some mismatch between the adopted
mathematical model and the physical processes of the problem.
On the other hand, the small residuals observed for themeasured concentration validate the use of the dual-porosity
physical non-equilibrium model, in particular the hypotheses of
having a constant f (see Equations 12a,b) and a non-reactive
tracer (KBr) (see Equation 13d) for our experiment.
CONCLUSIONSIn this study we used a Bayesian approach for solving aparameter estimation problem involving the simultaneous

42

characterization of soil hydraulic and solute transport parametersby means of a small-scale laboratory soil column experiment. The
Metropolis-Hastings algorithm of the Markov Chain Monte Carlo
method was implemented for estimation of the unknown
parameters. The sensitivity coefficients with respect to the model
parameters were analyzed, revealing that pressure head (at the
middle of the column) and outflow solute concentration
measurements provided meaningful information for estimation ofsome of the unknowns. Even for such measurements, thesensitivity coefficients were generally very small, thus making theestimation procedure difficult and sensitive to measurement
errors. For those parameters with small or linearly dependent
sensitivity coefficients, Gaussian priors were utilized, based on
direct measurements with other techniques. For the otherparameters, non-informative uniform distributions were used aspriors.
The estimation approach proposed in this work consisted oftwo steps. In the first step, the hydraulic parameters were
estimated by using the measured data of the pressure head in the
middle of the column. In the second step, these estimates of theflow parameters served as prior information for the coupled flowtransport problem, in which the transport parameters were
estimated from the measured concentration data of the column
discharge and of the pressure head in the middle of the column.
Since the flow parameters were treated as random variables in thesecond step, they were also allowed to vary in the Markov chainin accordance with their prior distributions. The convergence ofthe Markov chain was verified by calculating the difference
between the means in different ranges of the chain. Results
revealed the accuracy and robustness of the estimation approachby producing small, although correlated, residuals.

Acknowledgements. The authors would like to thank Capes and
CNPq for supporting this research.
REFERENCES
Abbaspour, K., Schulin, R., van Genuchten, M.Th., 2001. Estimating unsaturated soil hydraulic parameters using ant
colony optimization. Adv Water Resour., 24, 827–841.
Abbaspour, K., Johnson, C., van Genuchten, M.Th., 2004.
Estimating uncertain flow and transport parameters using asequential uncertainty fitting procedure. Vadose Zone J., 3,
1340–1352.
Beck, J., Arnold, K., 1997. Parameter Estimation in Engineering and Science. Wiley Interscience, New York.
Carducci, C.E., de Oliviera, G.C., da Costa, S., Zeviani, W.M.,
2011. Modeling the water retention curve in Oxisols using
the double van Genuchten equation. R. Bras. Ci. Solo., 35,
77–86. (In Portuguese.)
Dane, J.H., Topp, G.C., 2002. Methods of Soil Analysis, part 4,
Physical Methods. Soil Sci. Soc. Am., Madison, WI.
De Smedt, F., Wierenga, P.J., 1984. Solute transfer throughcolumns of glass beads. Water Resour. Res., 20, 225–232.
Diamantopoulos, E., Iden, S.C., Durner, W., 2012. Inverse
modeling of dynamic nonequilibrium in water flow with aneffective approach. Water Resour. Res., 48, W03503,
doi: 10.1029/2011WR010717.
Franssen, H.J.H., Gómez-Hernández, J.J., Sahuquillo, A.,
2003a. Coupled inverse modeling of groundwater flow andmass transport and the worth of concentration data. J.
Hydrol., 281, 281–295.
Franssen, H.J.H., Stauffer, F., Kinzelbach, W., 2003b. Joint
estimation of transmissivity and recharges-application:

Bayesian estimation of the hydraulic and solute transport properties of a small-scale unsaturated soil column

Stochastic characterization of well capture zones. J. Hydrol.,294, 1–3, 87–102.
Fu, J., Gómez-Hernández, J.J., 2009. Uncertainty assessment
and data worth in groundwater flow and mass transportmodeling using a blocking Markov chain Monte Carlo
method. J. Hydrol., 364, 3–4, 328–341.
Geweke, J., 1992. Evaluating the accuracy of sampling-based
approaches to the calculation of posterior moments. In: Bernardo, J., Berger, J., Dawid, A., Smith, A. (Eds): Bayesian
Statistics. Oxford University Press, London.
Goldberg, S., Kabengi, N.J., 2010. Bromide adsorption byreference minerals and soils. Vadose Zone J., 9, 780–786.
Hastings, W. K., 1970. Monte Carlo Sampling Methods Using
Markov Chains and Their Applications, Biometrika, 57, 97–109.
Hosseini, A.H., Deutsch, C.V., Mendoza, C.A., Biggar, K.W.,
2011. Inverse modeling for characterization of uncertainty intransport parameters under uncertainty of source geometry inheterogeneous aquifers. J Hydrol., 405, 3–4, 402–416.
Kaipio, J., Somersalo, E., 2004. Statistical and Computational
Inverse Problems. Applied Mathematical Sciences 160,
Springer-Verlag.
Kaipio, J., Somersalo, E., 2007. Statistical inverse problems:
discretization, model reduction and inverse crimes. J. Comp.
Appl. Math, 198, 493–504.
Kohne, J., Mohanty, B., Simunek, J., 2006. Inverse dualpermeability modeling of preferential water flow in a soilcolumn and implications for field-scale solute transport.
Vadose Zone J., 5, 59–76.
Kool, J.B., Parker J.C., van Genuchten, M.T., 1985. Determining soil hydraulic properties from one-step outflow experiments by parameter estimation, I. Theory and numericalstudies. Soil Sci. Soc. Am. J., 49, 1348–1354.
Laloy, E., Weynants, M., Bielders, C.L., Vanclooster, M., Javaux, M., 2010. How efficient are one-dimensional models
to reproduce the hydrodynamic behavior of structured soils
subjected to multi-step outflow experiments. J. Hydrol., 393,
1–2, 37–52, doi: 10.1016/ j.jhydrol.2010.02.017.
Lee, P., 2004. Bayesian Statistics. Oxford University Press,
London.
Leij, F.J., Russell, W.B., Lesch, S.M., 1997. Closed-form expressions for water retention and conductivity data.
Groundwater, 35, 5, 848–858.
Li, L., Zhou, H., Gómez-Hernández, J.J., Franssen, H.J.H.,
2012. Jointly mapping hydraulic conductivity and porosity
by assimilating concentration data via ensemble Kalman filter. J. Hydrol., 428–429, 152–169.
Maraqa, M.A., Wallace, R.B., Voice, T.C., 1997. Effect ofdegree of saturation on dispersivity and immobile water insandy columns. J. Contam. Hydrol., 25, 199–218.
Melamed, R., Jurinak, J.J., Dudley, L.M., 1994. Anion exclusion – pore water velocity interaction affecting transport ofbromine through an Oxisol. Soil Sci. Soc. Am. J., 58, 1405–
1410.
Metropolis, N., Rosenbluth, A., Rosenbluth, M., Teller, A.,
Teller, E., 1953. Equation of state calculation by fast computing machines. J. Chemical Phys., 21, 1087–1092.
Miller, C.T., Dawson, C.N., Farthing, M.W., Hou, T.Y., Huang,J., Kees, C.E., Kelley, C.T., Langtangen, H.P., 2013. Numerical simulation of water resources problems: Models, methods, and trends. Adv. Water Resour., 51, 405–437.
Nkedi-Kizza, P., Biggar, J.W., Selim, H.M., van Genuchten,
M.Th., Wierenga, P.J., Davidson, J.M., Nielsen, D.R., 1984.
On the equivalence of two conceptual models for describing
ion exchange during transport through an aggregated oxisol.
Water Resour. Res., 20, 8, 1123–1130.

Orlande, H.R.B., van Genuchten, M.Th., Cotta, R.M., Moreira,P.H., 2009. Bayesian estimation of hydraulic and solutetransport parameters from laboratory soil column experiments. In: Proc. Int. Symp. Convective heat and mass transfer in sustainable energy, Hammamet, Tunisia, 20 p.
Özisik, M.N., Orlande, H.R.B., 2000. Inverse Heat Transfer:
Fundamentals and Applications. Taylor and Francis, NewYork.
Parker, J.C., Kool, J.B., van Genuchten, M.Th., 1985. Determining soil hydraulic properties from one-step outflow experiments by parameter estimation. II. Experimental studies.
Soil Sci Soc. Am. J., 49, 1354–1359.
Raithby, G.D., Torrance, K.E., 1974. Upstream-weighted differencing scheme and their application to elliptic problems
involving fluid flow. Computers & Fluids, 2, 191–206.
Schaap, M.G., Leij, F.J., van Genuchten, M.Th., 2001. Rosetta:
A computer program for estimating hydraulic parameterswith hierarchical pedotransfer functions. J Hydrol., 251,
163–176.
Si, B.C., Kachanoski, R.G., 2000. Estimating soil hydraulicproperties during constant flux infiltration: Inverse procedures. Soil Sci. Soc. Am. J., 64, 2, 439–449.
Šimůnek, J., Jarvis, N., van Genuchten, M.Th., Gardenas, A.,
2003. Review and comparison of models for describing nonequilibrium and preferential flow and transport in the vadosezone. J. Hydrol., 272, 14–35.
Šimůnek J., Bradford S.A., 2008. Vadose zone modeling: Introduction and importance. Vadose Zone J., 7, 2, 581–586.
Šimůnek, J., Šejna, M., Saito, H., Sakai, M., van Genuchten, M.
Th., 2013. The HYDRUS-1D Software Package for Simulating the Movement of Water, Heat, and Multiple Solutes in
Variably Saturated Media, Version 4.16, HYDRUS Software
Series 3, Department of Environmental Sciences, Universityof California Riverside, Riverside, California, USA, 340 p.
(http://www.pcprogress.com/Downloads/Pgm_hydrus1D/HYDRUS1D4.16.pdf).
Shackelford, C.D., 1991. Laboratory diffusion testing for waste
disposal – A review. J. Contam. Hydrol., 7, 177–217.
Sommer, R., Fölster, H., Vielhauer, K., Maklouf, E.J., Vlek,
P.J.G., 2003. Deep soil water dynamics and depletion bysecondary vegetation in the Eastern Amazon. Soil Sci. Soc.
Am. J., 67, 1672–1686.
Spohrer, K., Herrmann, L., Ingwersen, J., Stahr, K., 2006.
Applicability of uni- and bimodal retention functions for water flow modeling in a tropical Acrisol. Vadose Zone J., 5,
48–58.
Tan, S., Fox, C., Nicholls, G., 2006. Inverse Problems. Course
Notes for Physics 707. University of Auckland, Auckland.
van Genuchten, M.Th., Wierenga, P.J., 1976. Mass transfer
studies in sorbing porous media, I. Analytical solutions. Soil
Sci. Soc. Am. J., 40, 4, 473–480.
van Genuchten, M.Th., 1978. Calculating the unsaturated hydraulic conductivity with a new closed-form analytical model. Hydrology Document Number 412. Department of Civil
Engeneering, Princeton University, Princeton, New Jersey,USA.
van Genuchten, M.Th., 1980. A closed form equation for predicting the hydraulic conductivity of unsaturated soils. Soil
Sci. Soc. Am. J., 44, 5, 892–898.
van Genuchten, M.Th., Nielsen, D.R., 1985. On describing andpredicting the hydraulic properties of unsaturated soils. Annales Geophysicae, 3, 5, 615-628.
Vrugt, J.A., Gupta, H.V., Dekker, S.C., Sorooshian, S.,
Wagener, T., Bouten, W., 2006. Application of stochastic

43Paulo H. S. Moreira, Martinus Th. van Genuchten, Helcio R. B. Orlande, Renato M. Cotta

parameter optimization to the Sacramento soil moisture accounting model. J. Hydrol., 325, 288–307.
Vrugt, J.A., ter Braak, C.J.F., Diks, C.D.H., Schoups, G., 2013.
Hydraulic data assimilation using particle Markov chain
Monte Carlo simulation: Theory, concepts and applications.
Adv. Water Res., 51, 457–478.
Xu, T., Gómez-Hernández, J.J., Zhou, H., Li, L., 2013. The
power of transient piezometric head data in inverse modeling: An application of the localized normal-score EnKF withcovariance inflation in a heterogenous bimodal hydraulicconductivity field. Adv. Water Res., 53, 100–118.

44

Yates, S.R., van Genuchten, M.Th., Warrick, A.W., Leij, F.J.,
1992. Analysis of measured, predicted, and estimatedhydraulic conductivity using the RETC computer program.
Soil Sci. Soc. Am. J., 56, 2, 347–354.
Zhou, H., Gómez-Hernández, J.J., Franssen, H.J.H., Li, L.,2011. An approach to handling non-Gaussianity of parameters and state variables in ensemble Kalman filtering. Adv.
Water Res., 34, 7, 844–864.
Received 4 November 2014
Accepted 15 September 2015