Journal of Hydrology 624 (2023) 129822

Contents lists available at ScienceDirect

Journal of Hydrology
journal homepage: www.elsevier.com/locate/jhydrol

Research papers

An in-depth analysis of Markov-Chain Monte Carlo ensemble samplers forinverse vadose zone modeling
Giuseppe Brunetti a,e ,∗, Jiri Šimunek b , Thomas Wöhling c,d , Christine Stumpp aa Department of Water, Atmosphere and Environment, Institute of Soil Physics and Rural Water Management, University of Natural Resources and Life Sciences,

Vienna, Muthgasse, 18, Vienna, 1180, Austriab
Department of Environmental Sciences, University of California, Riverside, 900 University Ave, Riverside, 92521, CA, USAc
Chair of Hydrology, Technische Universität Dresden, Bergstrasse 66, Dresden, 01069, GermanydLincoln Agritech Ltd, Environmental Research, Hamilton, 3214, New Zealande
Department of Civil Engineering, University of Calabria, Rende, 87019, Italy

ARTICLEINFOThis manuscript was handled by Yuefei Huang,
Editor-in-Chief, with the assistance of Shuo Zhang,
Associate Editor.
Keywords:
Vadose zone
Markov-chain Monte Carlo
Ensemble samplers
Differential evolution
Affine invariant

ABSTRACTThis study elucidates the behavior of Markov-Chains Monte Carlo ensemble samplers for vadose zone inverse
modeling by performing an in-depth comparison of four algorithms that use Affine-Invariant (AI) moves orDifferential Evolution (DE) strategies to approximate the target density. Two Rosenbrock toy distributions, andone synthetic and one actual case study focusing on the inverse estimation of soil hydraulic parameters using
HYDRUS-1D, are used to compare samplers in different dimensions d. The analysis reveals that an ensemble
with 𝑁 = 𝑑 + 1 chains evolved using DE-based strategies converges to the wrong stationary posterior, while
AI does not suffer from this issue but exhibits delayed convergence. DE-based samplers regain their ergodic
properties when using 𝑁 ≥ 2𝑑 chains. Increasing the number of chains above this threshold has only minor
effects on the samplers’ performance, while initializing the ensemble in a high-likelihood region facilitates
its convergence. AI strategies exhibit shorter autocorrelation times in the 7d synthetic vadose zone scenario,
while DE-based samplers outperform them when the number of soil parameters increases to 16 in the actual
scenario. All evaluation metrics degrade as d increases, thus suggesting that sampling strategies based only oninterpolation between chains tend to become inefficient when the bulk of the posterior lays in increasingly
small portions of the parameters’ space.

1. Introduction
From simple linear regressions to complex numerical systems, models use mathematics, physics, and chemistry to translate real experience
into formalized thinking. The resulting array of concepts and formulations (i.e., model structure) is accepted if it can reproduce or predict
the behavior of the real system (i.e., measurements) with satisfactory
accuracy under multiple conditions. The model adherence to reality
depends not only on its structural adequacy, but also on a proper (andphysically realistic) definition of its parameters. The process by which
the model structure is selected and parameters are adjusted to minimize
the distance between model predictions and observations is broadly
defined as model calibration.
Model calibration plays a prominent role in vadose zone hydrology,
as relevant physical quantities (e.g., recharge and leaching fluxes) cannot be directly measured at the field scale. Hence, measurements of soil
hydraulic states (e.g., matric potential, volumetric water content) arefrequently used to calibrate vadose zone hydrological models, whichshould finally offer a comprehensive and experimentally validated

framework to estimate water and solute fluxes. Such indirect estimation
is feasible only if the model structure realistically ties together different
water and solute transport processes in the soil, thus leading to a wellconstrained theoretical framework. This is the case for mechanistic
models, such as HYDRUS (Šimůnek et al., 2016), SWAP (Van Dam et al.,2008), or MACRO (Jarvis and Larsbo, 2012).
Generally founded on the Richards and advection–dispersion equations, mechanistic models require the parametric definition of the soil
hydraulic (e.g. Durner, 1994; van Genuchten, 1980) and transport
properties. These parameters are often calibrated on transient fieldobservations, which are preferable to laboratory estimates in describing
effective properties at the field level due to better adherence to local
operating conditions (Ritter et al., 2003).
Deterministic calibration procedures focus on obtaining a single setof parameters that leads to the best attainable match between model
predictions and real measurements. However, measurement errors aswell as conceptual model approximations introduce uncertainty in∗ Corresponding author at: Department of Civil Engineering, University of Calabria, Rende, 87019, Italy.

https://doi.org/10.1016/j.jhydrol.2023.129822
Received 21 December 2022; Received in revised form 9 May 2023; Accepted 10 June 2023
Available online 26 June 2023
0022-1694/© 2023 Elsevier B.V. All rights reserved.

Journal of Hydrology 624 (2023) 129822

G. Brunetti et al.

densities, as the Gaussian transition kernel is not suited for thesekinds of distributions. The problem exacerbates if the covariance matrix
specified at the beginning of the analysis is not properly designed, thus
leading to low acceptance rates and slow convergence.
MCMC ensemble samplers overcome these limitations by simultaneously evolving multiple chains using different strategies (e.g., differential evolution, affine moves). Each chain is updated using the current
positions of all the other chains in the ensemble, thus assuming that the
distribution of the chains in the complementary ensemble carries useful
information about the target density, 𝑝(𝒙). In this way, it is possible toautomatically adapt the trial move to the target distribution without
any assumption on its shape. Different studies have shown that ensemble samplers outperform RWMH schemes in correlated anisotropic
and moderate-dimensional problems (Goodman and Weare, 2010; ter
Braak, 2006; Karamanis and Beutler, 2020). Their use for vadose zone
model calibration is growing, but it is still limited and not variegated, as
most of the existing research coupled mechanistic vadose zone models
with ensemble samplers based on the Differential Evolution strategy(DE-MC) (Huisman et al., 2010; Laloy et al., 2013; Wöhling and Vrugt,
2011; Vrugt and Beven, 2018; Lehmann et al., 2022; Kuhlman et al.,2022). Despite good results, Goodman and Weare (2010) developed
another ensemble sampler that uses affine invariant moves (AI) to deal
with skewed posteriors. The algorithm was successfully used by Li
et al. (2022) to estimate soil hydraulic parameters of a dual porous
system, and by Brunetti et al. (2019b) to calibrate the HYDRUS model.
However, a rigorous benchmark of different MCMC algorithms to provide guidelines for their application in vadose zone hydrological model
calibration is still missing.
This work intends to address this scientific gap and elucidate
the performance of different MCMC ensemble samplers in vadosezone model calibration. Only those ensemble samplers that have been
proven to be reversible and ergodic are considered (Appendix A).
These include: the Affine Invariant Ensemble Sampler (AI) as presented
in Goodman and Weare (2010), the Differential Evolution MCMC
(DE-MC) developed by ter Braak (2006), and the Differential Evolution with snooker update (DE-SNK) proposed by ter Braak and
Vrugt (2008). Multiple scenarios are used to evaluate algorithms interms of convergence, integrated autocorrelation time, effective sample
size, and acceptance rate. The problem is addressed in the following
way. First, n-dimensional toy likelihood functions that approximate
typical anisotropic target distributions encountered in vadose zone
modeling are used to test the MCMC samplers. Samplers’ convergenceis analyzed for different posterior shapes and dimensions, increasing
number of chains, and multiple initializations. Then, the HYDRUS
model is utilized to generate synthetic data mimicking measurementsfrom a weighable lysimeter. Data are combined with ensemble samplers
to calibrate HYDRUS and approximate the posterior distribution ofsoil hydraulic parameters. Finally, actual soil moisture observations
from a layered soil profile are used to compare the performance ofdifferent samplers in a real vadose zone model calibration scenario.
The resulting array of testing scenarios is used to clarify pitfalls andopportunities in the use of ensemble samplers for vadose zone model
calibration. Multimodal posterior distributions, typical of geostatistical
inverse problems (Reuschen et al., 2020, 2021), are not included in theanalysis, being less frequent in vadose zone hydrology.

the calibration, which question the assumption of a single true/bestparameter set. Instead, it turns out to be more scientifically realistic toassume the existence of a distribution of probable parameters’ values
that acknowledges the existence of an imperfect model, as well as errors
in the data. The estimation of this distribution is one of the objective of
the uncertainty analysis, which has become an integral part of vadosezone hydrological model calibration (Vereecken et al., 2016).
The Bayes theorem offers an appealing probabilistic framework
to uncertainty assessment in hydrological modeling. It constitutes astatistical approach for the estimation of a set of parameters 𝒙 in amodel 𝐻 for the data 𝐷:𝑝(𝐷|𝒙, 𝐻)𝑝(𝒙|𝐻)
𝑝(𝒙|𝐷, 𝐻) =(1)
𝑝(𝐷|𝐻)
where 𝑝(𝒙|𝐷, 𝐻) ≡ 𝑝(𝒙) is the posterior probability of the parameters,𝑝(𝐷|𝒙, 𝐻) ≡ 𝐿(𝒙) is the likelihood, 𝑝(𝒙|𝐻) ≡ 𝜋(𝒙) is the prior, and𝑝(𝐷|𝐻) ≡ 𝛺 is the marginal likelihood. By combining prior information with observations and model predictions, it became popularamong hydrologists as it enables them to infer parameter posterior
distributions, verify model adequacy, and assess the model’s predictive
uncertainty (Brunetti et al., 2020a,b; Huisman et al., 2010; Vrugt et al.,2008; Scharnagl et al., 2011; Wöhling and Vrugt, 2011; Schübl et al.,2022b). In particular, the posterior distribution is frequently the variable of interest for modelers, as it describes the epistemic uncertainty
of model parameters conditioned on measurements. However, its shape
can only be approximated using Monte Carlo sampling techniques,
since Eq. (1) has no closed-form solution for most real applications.
Sampling from the posterior distribution is statistically and numerically challenging, as it involves approximating a multidimensional
distribution whose shape is unknown a priori, but which is frequently
anisotropic and correlated (e.g. Brunetti et al., 2019b; Šimůnek and van
Genuchten, 1996). Markov-Chain Monte Carlo (MCMC) can solve thisissue and approximate the parameter posterior distribution, 𝑝(𝒙), which
can also be used for the calculation of sensitivity measures (Pérez et al.,2006; Xiao and Nowak, 2022).
MCMCs evolve step-by-step one or multiple chains to simulate aseries of values so that their density follows the underlying posterior
after a given amount of time. Briefly, starting from one location in
the parameter space, a proposal distribution is used to propose a newposition. The proposed position is accepted or rejected with sometransition probability (Hastings, 1970). If accepted, the chain moves tothe new proposed state. By repeating this procedure a large numberof times, the algorithm generates a chain of states where the next
proposed position depends only on the current position (Markov process). Assuming that both the so-called detailed balance and ergodicity
are preserved (Neal, 1993), the resulting scheme will sample from aunique and invariant distribution that is equal to the posterior 𝑝(𝒙)(Appendix A), and its sampling efficiency will mainly depend on how
the proposal distribution conforms to the target density.
MCMCs’ applications for uncertainty assessment in vadose zone
modeling are relatively new but highly promising. Some studies haveshown how these methods can be successfully combined with mechanistic vadose zone hydrological models to inversely estimate soil
hydraulic and transport parameters and assess their uncertainty (e.g.
Arora et al., 2012; Brunetti et al., 2019b; Eck et al., 2016; Jana andMohanty, 2012; Laloy et al., 2013; Vrugt et al., 2008; Vrugt andBeven, 2018; Wöhling and Vrugt, 2011; Wöhling et al., 2012; Laloyet al., 2012; Das et al., 2008). Altogether, these and other studies (e.g.
Brunetti et al., 2019b; Šimůnek and van Genuchten, 1996) highlighted
the peculiarity of inverse problems in vadose zone hydrology, whichis the presence of skewed, anisotropic, and correlated posterior parameters distribution. This poses several problems to classic MCMC
approaches, such as the single chain Random Walk Metropolis–Hastings
schemes (RWMH) (Hastings, 1970), which use a multivariate Gaussian
proposal to generate new states of the chain. These new states are
frequently rejected when sampling from anisotropic and correlated

2. Methods
2.1. Ensemble Monte Carlo Markov chains
A pseudocode of ensemble samplers is shown in Fig. 1. An ensemble
constituted by 𝑁 chains is first initialized in the parameters space,
and then evolved by moving one chain 𝒙𝑘 at a time using the current
positions of all the other remaining chains, which form the complementary ensemble. The ensemble size 𝑁 must be larger than the numberof dimensions, 𝑑, to guarantee that no lower-dimensional hyperplane2

Journal of Hydrology 624 (2023) 129822

G. Brunetti et al.

Table 1
Probability of different moves in the four ensemble samplers analyzed.
Algorithm ID

AI
DE
AI-SK
DE-SK

SNK

0
1
0
0.9

0
0
0.1
0.1

2.1.3. Affine invariant ensemble sampler
Goodman and Weare (2010) developed an approach that uses only
one additional chain 𝒙𝑗 from the complementary ensemble to generatea proposal that is affine invariant:contains all the chains (ter Braak, 2006; Goodman and Weare, 2010).
A loop through 𝑁 chains constitutes one step of the ensemble. Thisalgorithmic structure is common for all ensemble samplers, whichinstead differ in proposing and accepting new states. The algorithm
will sample from a unique and invariant distribution that is equal to theposterior if, at each step, the move and acceptance scheme is reversible
and ergodic. It follows that it is possible to evolve the ensemble using
different moves, provided that all chains are evolved using the same
move scheme in the current step.

𝒚 = 𝒙𝑗 (𝑡) + 𝑍(𝒙𝑘 (𝑡) − 𝒙𝑗 (𝑡))(6)where 𝑍 is a scaling variable, which is random sampled from a symmetric distribution 𝑔(𝑍) of the form:{
𝑍 −1∕2 𝑎−1 ≤ 𝑍 < 𝑎𝑔(𝑍|𝑎) =(7)0
𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒

2.1.1. Differential evolution
The Differential Evolution Monte Carlo Markov Chain algorithm(DE-MC) was developed by ter Braak (2006). It combines the differential evolution strategy used for optimization with MCMC to address
the issue of choosing an appropriate scale and orientation for the
proposal distribution. The proposed move for the chain 𝒙𝑘 is generated
as follows:

where 𝑎 is typically equal to 2. The proposal is accepted with probability:()
𝑝(𝒚)𝑞 = 𝑚𝑖𝑛 1, 𝑍 𝑑−1(8)𝑝(𝒙𝑘 (𝑡))The main advantages of AI are its simplicity and its affine invarianceproperty, which makes it particularly suited for anisotropic and correlated target densities (Appendix B). The proof of its affine invarianceis detailed in Christen and Foxy (2010) and not reported here for
conciseness. Nevertheless, the scheme tends to become inefficient when
the dimensionality of the inverse problem increases (Huijser et al.,2022).

(2)

Two random chains, 𝒙𝑅1 and 𝒙𝑅2 , are selected from the complementaryensemble, and their difference is multiplied by a scale factor 𝛾, andadded to a random vector 𝝐 drawn from a symmetric distribution withsmall variance b compared to that of the target density. The resulting
vector is added to 𝑋𝑘 to generate the proposed move 𝒚. The proposal
is accepted with probability:
)
(
𝑝(𝒚)
(3)
𝑞 = 𝑚𝑖𝑛 1,
𝑝(𝒙𝑘 (𝑡))2.2. Ensemble samplers comparison
Previously described ensemble moves can be mixed to obtain hybrid
ensemble samplers. For example, ter Braak and Vrugt (2008) proposeda sampling scheme that, at each algorithm step, selects DE-MC and SNKwith 90% and 10% of probability, respectively. On this line, in thepresent study, the performance of four different ensemble samplers areinvestigated (Table 1). Algorithms are compared by assuming constantscaling parameters (i.e., 𝛾, 𝑎) during sampling, as their ergodicity andreversibility was proven in these settings. Their value was set according
to developers’ recommendations and other studies (Hogg and ForemanMackey, 2018; Goodman and Weare, 2010; ter Braak and Vrugt, 2008;
ter Braak, 2006).
Two 𝑡𝑜𝑦 distributions, and one synthetic and one actual case study
focusing on the inverse estimation of soil hydraulic parameters, areutilized to test and compare algorithms. Due to their analytical form
and limited computational cost, toy distributions are investigated fora broad array of conditions. These include: different posterior shapes
and dimensions, increasing number of chains, and multiple chains

The resulting scheme is ergodic and preserves the detailed balance,
as proven by ter Braakand Vrugt (2008). The same authors suggest√
setting 𝛾 = 2.38∕ 2𝑑 and 𝝐 ∼ 𝑁(0, 𝑏)𝑑 , where d is the number ofparameters. If the posterior is suspected to be multimodal, 𝛾 can beoccasionally set equal to 1 to allow jumps between disconnected modes.
2.1.2. Differential evolution with snooker update
ter Braak and Vrugt (2008) developed a variant of DE-MC that uses
a snooker update to evolve chains (SNK), and that is particularly suitedfor moderate-dimensional problems. The proposed move is generated asfollows:𝒚 = 𝒙𝑘 (𝑡) + 𝛾(𝒙𝑍1 (𝑡) − 𝒙𝑍2 (𝑡)) + 𝝐

DE-MC

1
0
0.9
0

projected orthogonally on to the line 𝒙𝑘 – 𝒙𝑍 , thus yielding 𝒙𝑍1 and𝒙𝑍2 . Their difference is multiplied by 𝛾 and added to 𝒙𝑘 to obtain the
new proposed position 𝒚, which is accepted with probability:()‖𝑑−1𝑝(𝒚) ‖‖𝒚 − 𝒙𝑍 (𝑡)‖𝑞 = 𝑚𝑖𝑛 1,(5)‖𝑑−1
𝑝(𝒙𝑘 (𝑡)) ‖‖𝒙𝑘 (𝑡) − 𝒙𝑍 (𝑡)‖ter Braak and Vrugt (2008) suggest randomly sampling 𝛾 froma bounded uniform distribution 𝛾 ∼ 𝑈 [1.2, 2.2]. SNK is generally
mixed with other moves to diversify the jumping possibilities, and it
is expected to be more efficient in moderate (10 ≤ 𝑑 < 30) to highdimensions (𝑑 ≥ 30) (ter Braak and Vrugt, 2008).

Fig. 1. Pseudocode of ensemble samplers.

𝒚 = 𝒙𝑘 (𝑡) + 𝛾(𝒙𝑅1 (𝑡) − 𝒙𝑅2 (𝑡)) + 𝝐

Moves’ probability (–)AI(4)

To update the position of chain 𝒙𝑘 , three chains from the complementary ensemble are selected, namely 𝒙𝑍 , 𝒙𝑅1 and 𝒙𝑅2 . 𝒙𝑅1 and 𝒙𝑅2 are
3

Journal of Hydrology 624 (2023) 129822

G. Brunetti et al.

Table 2
Benchmark scenarios, settings, and evaluation metrics, used to compare ensemble samplers listed in Table 1. The number of chains 𝑁 is set to respect the constraint 𝑁 > 𝑑.
Scenario

Type

𝑑

𝑁

Classic

10
30(11, 20, 50, 100, 200)
(31, 50, 100, 200)

Shifted-rotated1030(11, 20, 50, 100, 200)
(31, 50, 100, 200)

Synthetic
Actual7
16(50, 100, 200)
(50, 100, 200)

Rosenbrock

Vadose zone

initialization. Results are used to guide the following computationallyintensive analysis of ensemble samplers in vadose zone scenarios. The
full set of the investigated scenarios is listed in Table 2, and describedin detail in the following sections.

Initialization

Evaluation metric

High-likelihood
RandomLow likelihood

EMD, 𝜏𝑎𝑣𝑔 ,AR𝑎𝑣𝑔 , ESS𝑎𝑣𝑔

High-likelihood
High-likelihood

𝜏𝑎𝑣𝑔 , AR𝑎𝑣𝑔 , ESS𝑎𝑣𝑔
𝜏𝑎𝑣𝑔 , AR𝑎𝑣𝑔 , ESS𝑎𝑣𝑔

samplers exhibit ARs that decrease with time as chains approach thebulk of the target distribution. This is particularly evident in highdimensions, as the typical set is contained in a very small volume ofthe parameters’ space, which is difficult to sample efficiently.
Effective sample size. The ESS, which measures the informativeness ofan ensemble of chains, is calculated for each parameter, and is defined
as follows (Geyer, 1992):

2.2.1. Evaluation metrics
Four evaluation metrics are used to compare ensemble samplers
and assess their convergence: the Integrated Autocorrelation Time (𝜏),the Earth Mover’s Distance (EMD), the Acceptance Rate (AR), and the
Effective Sample Size (ESS). Albeit used in other studies, the Gelman–
Rubin statistic (Gelman et al., 2013) is not applied here since chains are
not formally independent in ensemble samplers, as described in Fig. 1.
An alternative approach is to use this metric to check the convergenceof multiple ensembles evolved simultaneously. However, this would
substantially increase the analysis’s computational burden, especiallyfor vadose zone scenarios. Therefore, EMD and 𝜏 are here considered as
the key indicators of convergence (Hogg and Foreman-Mackey, 2018).

𝑀 ⋅𝑁
(9)
𝜏
where M is the chains’ length, N is the number of chains, and 𝜏 is theintegrated autocorrelation time.
The integrated autocorrelation time and the effective sample sizeare calculated for each parameter of the distribution and then averaged
to have a single evaluation metric. Similarly, an average acceptance
rate is obtained by averaging the acceptance rates calculated for eachchain in the ensemble.
ESS =Integrated autocorrelation time. 𝜏 is the number of steps required bya chain to produce an independent sample, and represents the key
diagnostic for MCMC algorithms. Shorter 𝜏 are desirable as this implies
that the chain requires a lower number of steps to start sampling fromthe target density. In this study, 𝜏 is estimated for each parameter by
combining the procedure of Sokal (1997) and a mean normalized autocorrelation function, computed by averaging the normalized autocorrelation functions calculated for each chain in the ensemble (Brunetti
et al., 2019b). The transitory period required by the ensemble toconverge to the true posterior is called burn-in, and it is here assumed
equal to two times the maximum integrated autocorrelation time.

2.2.2. Toy distributions
The n-dimensional Rosenbrock distribution. 𝑇 𝑜𝑦 functions and distributions are frequently used to test and compare optimization andsampling strategies, respectively. Being computationally inexpensive,
they allow extensive testing and can be designed to mimic response
surfaces or posterior densities that most commonly arise from inverse
problems in real applications. In vadose zone hydrology, this would
mean a moderate-dimensional distribution (e.g., 𝑑 < 30), which isfrequently banana-shaped (Brunetti et al., 2019b; Šimůnek and van
Genuchten, 1996).
The Rosenbrock distribution conforms to these criteria, being characterized by a bulk lying in an elongated, narrow parabolic valley (Goodman and Weare, 2010). However, the banana-shaped valley
tends to disappear when extending the original 2-dimensional formulation to higher dimensions (Pagani et al., 2022), thus facilitating
sampling from MCMC algorithms. The complexity can still be retainedby shifting and rotating the distribution, thus making variables nonseparable and the bulk characterized by different scales in different
dimensions.
Therefore, in the present work, classic and shifted-rotated variantsof a d-dimensional Rosenbrock-type distribution are used to preliminarily compare ensemble samplers (Table 2). Sampling efficiency andconvergence is investigated for both variants in 10 and 30 dimensions.
This choice is made intentionally to analyze the performance of thealgorithms for a number of parameters that is potentially representativeof the average and worst-case dimensionalities encountered in vadosezone model calibration. The shifted-rotated distribution is obtained bymultiplying a modified classic Rosenbrock likelihood by a square dxd
rotation matrix 𝑹 and adding a d-dimensional translation vector 𝒐:Earth mover’s distance. EMD measures the distance between two probability distributions. More precisely, it is defined as the minimal costrequired to transform one histogram into the other, where there isa ‘‘ground distance’’ between the basic features that are aggregatedinto the histogram (Pele and Werman, 2009). In the present study,
EMD is used to measure the distance between a reference distribution
and the posterior approximated by ensemble samplers in a relatively
small number of steps. The latter is obtained by excluding the burn-in
period, while the reference target distribution is produced by running
the widely used implementation of the AI algorithm (Foreman-Mackey
et al., 2013) for ten million steps using 500 chains initialized aroundthe maximum likelihood. The use of a well-tested external package isintended to guarantee an unbiased estimation of the posterior. Low
EMD values are considered an indicator of the convergence of theensemble samplers to the target distribution. It must be noted that
this metric is used only for the toy distributions, as the calculationof a reference posterior for vadose zone problems would require an
impractically high computational cost.

{ 𝑑−1
]}
∑[
𝐿(𝒙) ∝ exp −100(𝒙2𝑖 − 𝒙𝑖+1 )2 + (𝒙𝑖 − 1)2 + (𝒙𝑖 − 1)2

Acceptance rate. AR is a qualitative indicator of the algorithm sampling
efficiency, as it describes the percentage by which a proposed stepis accepted. A moderately high acceptance rate (e.g., 20%–50%) isdesirable for ensemble samplers, as it indicates that proposals are
accepted fairly frequently. Higher ARs generally indicate that chains
movements are too small, thus leading to long autocorrelation times,
while very low ARs (e.g., <10%) suggest that the generated proposals
are unlikely to belong to the posterior distribution. In general, ensemble

𝑖=1

(10)

( ( 2.048(𝒙 − 𝒐) ))𝐿(𝒙) = 𝐿 𝑹
+1
(11)
100
The shift vectors and the rotation matrices are taken from Lianget al. (2014). 𝑹 is obtained from standard normally distributed entries
by Gram–Schmidt ortho-normalization, while 𝒐 is randomly distributedin the interval (−100,100). 𝑹 and 𝒐 can be downloaded at this link4

Journal of Hydrology 624 (2023) 129822

G. Brunetti et al.

(https://github.com/giubrune/Paper_EnsembleMCMC). The statisticalframework to analyze toy distributions is implemented in the Python
programming language. The resulting reference marginal distributions
for the 10-dimensional Classic and Shifted-Rotated Rosenbrock problems are shown in Figs. C.2 and C.3, respectively. Their 30-dimensionalcounterparts are not reported for visualization clarity, as they would
result in tight 2D plots.
The transformation is applied only to the likelihood, while theprior is kept unchanged, assumed uniform, and spanning an interval(−100, 100) for all parameters. Ensemble samplers are compared for
an increasing number of chains, by taking into account the requirement
𝑁 > 𝑑 (Table 2). Simultaneously, three different chains initializations
are investigated:

atmospheric forcing conditions in HYDRUS for both calibration scenarios (Fig. D.7). The first three months serve to generate a spin-up period
to relax initial conditions during model calibration, while the remaining
six months constitute the proper calibration window. Measured volumetric water contents and synthetic data (simulated for a homogeneous
soil profile) from 2014-04-01 to 2014-09-30 are thus used to inversely
estimate soil hydraulic parameters in the actual (Fig. D.9) and synthetic
scenario (Fig. D.8), respectively.
Water flow and root water uptake. The hydrological model HYDRUS1D (Šimůnek et al., 2016) is used to numerically solve the onedimensional Richards equation, which describes the variably-saturated
water flow in the soil:
[
(
)]
𝜕𝜃(ℎ)
𝜕
𝜕ℎ
=
𝐾(ℎ)
+ 1 − 𝑆(ℎ)
(12)
𝜕𝑡
𝜕𝑧
𝜕𝑧

1. High-Likelihood Region: Chains are randomly sampled from amultivariate normal distribution with mean equal to the known
maximum likelihood and diagonal covariance matrix (𝜎 2 = 0.1
for all dimensions)
2. Random: Chains are randomly sampled from the uniform prior.
3. Low Likelihood Region: Chains are randomly sampled in a regionof the parameters’ space characterized by a similar low value of𝐿(𝒙) for both the Classic and Shifted-Rotated likelihood.

where 𝜃 is the volumetric water content (L3 L−3 ), 𝑡 is time (T), 𝑧 is
the vertical coordinate (L), 𝐾 is the unsaturated hydraulic conductivity(LT−1 ), ℎ is the pressure head (L), and 𝑆(ℎ) is a sink term (L3 L−3 T−1 ),defined as a volume of water removed from a unit volume of soil perunit of time due to root water uptake. The unimodal van Genuchten
model (van Genuchten, 1980) parameterizes the soil water retention
curve 𝜃(ℎ):{
𝜃𝑠 −𝜃𝑟ℎ<0
𝜃𝑟 + (1+(𝛼|ℎ|)
𝑛 )𝑚
(13)
𝜃(ℎ) =
𝜃𝑠ℎ. ≥ 0The use of toy functions is intended as a computationally inexpensive but informative first step in comparing ensemble samplers.
Results are used to narrow the number of vadose zone scenarios by
excluding settings (e.g., initialization, number of chains) that lead to
MCMC convergence problems.

where 𝜃𝑟 is the residual water content (L3 L−3 ), 𝜃𝑠 is the saturated
water content (L3 L−3 ), and 𝛼 (L−1 ) and 𝑛 (–) are shape parameters with
𝑚 = 1 − 1∕𝑛. The unsaturated hydraulic conductivity 𝐾(ℎ) is describedusing the van Genuchten–Mualem(VGM) relationship:

2.2.3. Vadose zone hydrological model calibration
Two case studies involving the inverse estimation of soil hydraulicparameters using the hydrological model HYDRUS-1D (Šimůnek et al.,2016) are considered (Table 2): (1) a synthetic scenario simulatinga homogeneous weighable lysimeter, (2) an actual scenario involving
transient soil moisture observations at different depths in a layered
vadose zone profile in Austria. For simplicity, both scenarios share thesame set of atmospheric conditions, but differ in terms of the bottom
boundary condition, soil layering, and measurements used in the model
calibration.
The probabilistic Bayesian framework as well as algorithms listed inTable 1 have been implemented in FORTRAN and embedded into theoriginal HYDRUS source code to facilitate the workflow, and lay thefoundations for a future software release that would allow a widespreaduse of ensemble samplers in vadose zone modeling. In the next sections,the field site and common characteristics are defined first, followed bya description of the modeling theory and inversion scenarios.

[]2⎧⎪𝐾 𝛩𝐿 1 − (1 − 𝛩1∕𝑚 )𝑚𝐾(ℎ) = ⎨ 𝑠
⎪𝐾𝑠
⎩

ℎ<0

(14)

ℎ. ≥ 0where 𝛩 is the effective saturation (–), 𝐾𝑠 is the saturated hydraulicconductivity (LT−1 ), and L is a pore connectivity parameter (–), which
is usually set to 0.5 (van Genuchten, 1980).
The simulated vegetation water uptake (i.e., 𝑆(ℎ) in Eq. (12)) fromthe root zone 𝑅 includes the effect of root density 𝛿(𝑧) and root water
stress 𝛽(ℎ, 𝑧) on the actual transpiration:
𝑆(ℎ) =∫𝑅

𝑆(ℎ, 𝑧)𝑑𝑧 = 𝑆𝑝∫𝑅𝛽(ℎ, 𝑧)𝛿(𝑧)𝑑𝑧(15)

where 𝛽(ℎ, 𝑧) is a dimensionless water stress function that depends
on the soil pressure head ℎ, 𝛿(𝑧) is a spatial distribution functionthat accounts for root density, and 𝑆𝑝 is the potential transpiration
demand, which is calculated from the potential evapotranspiration rate
𝐸𝑇𝑝 using the Leaf Area Index (LAI). 𝛽(ℎ, 𝑧) is parameterized by thepiecewise linear stress function proposed by Feddes et al. (1978), whoseparameters are set according to Wesseling et al. (1991), while 𝛿(𝑧) is
assumed to decrease linearly with depth. For a thorough descriptionof evapotranspiration partitioning and root water uptake in HYDRUS,
please refer to Brunetti et al. (2019a). Governing equations are solved
numerically using Galerkin type linear finite element schemes coupledwith an implicit integration in time for both saturated and unsaturated
conditions.

Case study description. The data originates from a soil water monitoring
station in Elsbethen, Austria (47◦ 45′ 49′′ N, 13◦ 04′ 54′′ E), belonging to
the Austrian Soil Water Monitoring Program of the Federal Ministryof Agriculture, Regions and Tourism. Observations from this program
were used in Schübl et al. (2022a) to inversely estimate soil hydraulicparameters at different locations in Austria. Measured particle size
distribution revealed a 150 cm-deep layered soil profile, with two silty
soil horizons of different bulk densities 𝜌𝑏 overlying a sandy layer(Table D.3). The vegetation consists of grass, whose roots’ density is
assumed to decrease linearly with depth from the high-organic topsoil
to 𝑧 = −20 cm. Long-term measurements of precipitation, climate
variables, and volumetric water content at multiple depths (i.e., 𝑧 (cm)
= [−10, −30, −60, −90, −150]) have been recorded to characterizewater fluxes in the vadose zone. Soil moisture was measured using time
domain reflectometers with an acquisition frequency of 1 h. Daily solar
radiation, air humidity and temperature, and wind speed are used to
calculate the daily reference evapotranspiration demand 𝐸𝑇0 using thePenman-Monteith equation (Allen et al., 1998).
Nine-month long time series of daily precipitation and evapotranspiration (i.e., 2014-01-01 to 2014-09-30) are extracted and set as

Synthetic scenario: A weighable lysimeter. A weighable lysimeter is anexperimental device frequently used in soil physics and hydrology tounderstand vadose zone functioning in field conditions (Reth et al.,2021). By allowing precise measurement of water fluxes across boundaries, weighable lysimeters limit uncertainty in model calibration. Theirfeatures have been replicated in synthetic modeling scenarios to understand the informativeness of different types of measurements forvadose zone modeling (e.g., Schübl et al., 2022b; Schelle et al., 2012).
In their most basic configuration, only the mass of water leaving thelysimeter across the lower boundary and the total mass of the lysimeter5

Journal of Hydrology 624 (2023) 129822

G. Brunetti et al.

Likelihood function. Error residuals are assumed to be uncorrelated andnormally distributed with constant variance 𝜎 2 , thus leading to thefollowing likelihood function 𝐿(𝒙):are measured. As shown by Durner et al. (2008), these observations
can be used only to infer effective soil hydraulic parameters, but not toestimate them for each soil layer.
Building on these considerations, in the present study, a synthetic
experiment is used to compare the performance of ensemble samplers for the inverse estimation of effective soil hydraulic parametersfrom a weighable lysimeter. The lysimeter mimics the soil characteristics of the actual soil profile in Austria, but is modeled in HYDRUS
as one-dimensional and homogeneous. Soil hydraulic parameters are
generated by combining the weighted average of the particle size distributions and bulk densities reported in Table D.3 with the pedotransfer
function ROSETTA (Schaap et al., 2001), thus leading to the following
effective VGM parameters: 𝜃𝑟 = 0.032, 𝜃𝑠 = 0.39, 𝛼 (1/cm) = 0.033, 𝑛 =1.45, 𝐾𝑠 (cm/day) = 98.71. The numerical domain is discretized into
151 Finite Elements refined at the top, where an atmospheric boundary
condition is specified to simulate precipitation and evaporation fluxes.
A seepage face is set at the bottom of the domain. A constant and equalvalue of 2 of the LAI is used to partition evapotranspiration fluxes. The
roots’ density is assumed to decrease linearly from the soil surface down
to a maximum rooting depth of 20 cm.
A forward simulation is performed with the resulting model to generate 183 days-long time series of simulated average volumetric watercontent, 𝜃𝑎𝑣𝑔 , and cumulative water flux across the bottom boundary of
the lysimeter, Q. The execution time of a single simulation is on average
0.3 s on the available computer (Intel i7-8750H CPU 2.20 GHz, RAM
16 GB). Synthetic data are generated by perturbing these time series
with Gaussian noise characterized by null mean and standard deviation
𝜎 equal to 0.001 (Nolz et al., 2013) and 0.01 cm (Schelle et al., 2012)for 𝜃𝑎𝑣𝑔 and Q, respectively. Finally, synthetic measurements (Fig. D.8)
are combined with ensemble samplers (Table 1) and HYDRUS to estimate the posterior distribution of soil hydraulic parameters and infer
𝜎 for both measurement types, thus leading to a total of 7 calibrated
parameters. Uniform priors are assumed for all parameters, and theirbounds are reported in Table D.4.

𝐿(𝒙) =

𝑘
∏
𝑗=1

√

[
]2
1
exp − 𝜎 2 (𝐻𝑗 (𝒙) − 𝐷𝑗 )
2
2𝜋𝜎 2
1

(16)

For algebraic simplicity, the logarithm of 𝐿(𝒙) is used:
𝓁(𝒙) = −𝑘𝑘
1 ∑𝑘ln(2𝜋) − ln(𝜎 2 ) −(𝐻𝑗 (𝒙) − 𝐷𝑗 )22
2
2
2𝜎 𝑗=1(17)

The aggregated sum of log-likelihood is used in the analysis whendifferent types of measurements are assimilated in the Bayesian inference (e.g., synthetic scenario). As the HYDRUS execution time for boththe synthetic and actual scenario is relatively low, the inference does
not require the use of a surrogate model to emulate the response of thehydrological model.
3. Results
3.1. Toy distributions
3.1.1. The 10-dimensional Rosenbrock distributions
Figs. 2 and 3 show the calculated EMD and 𝜏𝑎𝑣𝑔 , respectively,
as a function of algorithm step for the 10-dimensional Classic andShifted-Rotated Rosenbrock distributions. EMD and 𝜏𝑎𝑣𝑔 are computedfor each ensemble sampler listed in Table 1 by initializing 𝑁 =(11, 20, 50, 100, 200) chains in high and low likelihood regions, and in
random positions. Results indicate that all algorithms in Table 1 canapproximate the stationary target distribution (i.e., EMD approachingzero) only under specific conditions. In particular, the convergenceis affected by both the number of chains in the ensemble and theirinitializations.
The most striking result is that using an ensemble of 𝑁 = 𝑑+1 chains(i.e., 11) leads constantly to convergence issues for DE-based samplers
independently from their initializations. The calculated integrated autocorrelation time is low for DE and DE-SNK samplers (Fig. 3), whichquickly start sampling from a stationary distribution. However, EMD
values indicate that this distribution is not the true one, thus hintingthat the scheme is reversible but not ergodic (Roberts and Rosenthaly,2004) (Appendix A). This behavior is detailed in Figs. C.2 and C.3,
which compares the marginal reference distribution (black contours)with that obtained using DE with 𝑁 = 𝑑 + 1 (red contours). The
ensemble remains trapped near the maximum likelihood region andfails to explore the true posterior landscape. Conversely, AI and AI-SNKwith 𝑁 = 𝑑 + 1 chains are still able to approximate the target posterior
when chains start in a high-likelihood region (blue contours in Figs. C.2
and C.3). Other types of initializations delay their convergence (Fig. 2),
as also proven by the increased 𝜏𝑎𝑣𝑔 .
The situation improves for all samplers when the ensemble size(i.e., number of chains) is increased, although the convergence rate
remains influenced by the initialization (Fig. 2). In particular, 𝑁 =20 leads to a relatively quick convergence only if chains start in ahigh-likelihood region. Vice-versa 𝜏𝑎𝑣𝑔 doubles and triples if chains are
initialized in random positions and low likelihood regions, respectively(Fig. 3). The effect of the initialization gradually mitigates when 𝑁 ≥50, although starting chains in a high-likelihood region results constantly in lower autocorrelation times, independently from the shapeof the posterior.
The average acceptance rate AR𝑎𝑣𝑔 (Fig. C.4) consistently decreases
with the number of steps as chains approach the typical set (i.e., areaof the parameter space that covers most of the mass of the distribution)
for all algorithms. However, the rate of decline appears more evident
for DE-based samplers (i.e., DE and DE-SNK), which generally exhibit
lower acceptance rates compared to AI-based samplers (i.e., AI and AISNK). Increasing the number of chains has only minor effects on the

Actual scenario: Vadose zone profile. The second scenario compares
ensemble samplers performance on actual observations of volumetricwater content at multiple depths. After a preliminary pre-processing
phase, measurements at z = −150 cm have been excluded from theanalysis due to accuracy issues, thus reducing the number of observation points to four (i.e., 𝑧 (cm) = [−10, −30, −60, −90]) (Fig. D.9).
The effect of the spatial heterogeneity of soil properties is not includedin the analysis, as only point measurements of soil water states areavailable. A 150 cm-deep soil profile consisting of three soil materials(i.e., Table D.3) is modeled in HYDRUS-1D using 151 Finite Elements
refined at the top, where an atmospheric boundary condition is assigned. Since the water table lies far below the domain, a unit pressure
head gradient is set at the bottom of the domain to simulate freedrainage. The execution time of a single simulation is on average 0.3 s
on the available computer (Intel i7-8750H CPU 2.20 GHz, RAM 16 GB).
As the textural composition of the superficial soil layers includes
a substantial amount of fine particles (Table D.3), Eqs. (13) and (14)are corrected with a fictitious air-entry value of −2cm to force theslope of the retention function to become zero when approachingsaturation (Vogel et al., 2000). The correction has negligible effects
on the macroscopic description of the water retention function of finetextured soils, but it substantially alters the shape of the hydraulicconductivity near saturation, making it less nonlinear and more realistic
for this type of soils. Finally, root water uptake is simulated using thesame settings as the synthetic scenario.
Ensemble samplers are combined with actual measurements and
HYDRUS to estimate the posterior distribution of VGM soil hydraulicparameters for each soil layer (𝜃𝑟 , 𝜃𝑠 , 𝛼, 𝑛, 𝐾𝑠 ), and infer the measurement error 𝜎𝜃 , thus leading to a total of 16 calibrated parameters. Thisapproach implies that the uncertainty in the applied inputs (e.g., precipitation, evapotranspiration) is not treated independently in the analysis,
but lumped into the posterior of soil hydraulic parameters. Again,
uniform priors are assumed for all parameters, and their bounds arereported in Table D.5.
6

Journal of Hydrology 624 (2023) 129822

G. Brunetti et al.

Fig. 2. Earth Mover’s Distance EMD as a function of algorithm step for the 10-dimensional Classic and Shifted-Rotated Rosenbrock distributions. EMD is calculated for eachensemble sampler listed in Table 1 by initializing 𝑁 = (11, 20, 50, 100, 200) chains in high and low likelihood regions, and in random positions.

Fig. 3. Average integrated autocorrelation time 𝜏𝑎𝑣𝑔 as a function of algorithm step for the 10-dimensional Classic and Shifted-Rotated Rosenbrock distributions. 𝜏𝑎𝑣𝑔 is calculated
for each ensemble sampler listed in Table 1 by initializing 𝑁 = (11, 20, 50, 100, 200) chains in high and low likelihood regions, and in random positions.
7

Journal of Hydrology 624 (2023) 129822

G. Brunetti et al.

Fig. 4. Earth Mover’s Distance EMD as a function of algorithm step for the 30-dimensional Classic and Shifted-Rotated Rosenbrock distributions. EMD is calculated for eachensemble sampler listed in Table 1 by initializing 𝑁 = (11, 20, 50, 100, 200) chains in high and low likelihood regions, and in random positions.

Results show that the random and low-likelihood initializations
are less problematic for the 30-dimensional shifted-rotated Rosenbrock
distribution. This is mainly caused by the shapes of the target densities (Figs. C.2 and C.3). In particular, the shifted-rotated posteriorcrosses a larger portion of the parameter’s space (Fig. C.3), while
the classic Rosenbrock distribution has many components lying in anarrow part of the domain (Fig. C.2). In these conditions, the random and low-likelihood initializations spread the chains far from thebulk of the classic Rosenbrock posterior, but relatively close to it
for its shifted-rotated variant, thus facilitating the convergence of thesamplers.
The integrated autocorrelation time shows the desired asymptotic
trend mainly when chains are initialized in a high-likelihood region,
while other starting positions lead to delayed convergence (Fig. 5) especially for DE and DE-SNK. This is confirmed by the average acceptancefraction (Fig. C.5), which drops significantly when DE chains start far
from the bulk of the posterior. The calculated 𝜏𝑎𝑣𝑔 strongly increases for
all algorithms in 30 dimensions, showing that 10,000 steps are often insufficient to obtain a precise estimate of the integrated autocorrelation
time. This is a well-known problem for ensemble samplers since theirsampling strategy becomes inefficient in high dimensions as the typical
set is contained in a tiny area of the parameter space. The comparison of
𝜏𝑎𝑣𝑔 indicates that AI and AI-SNK are less efficient than DE and DE-SNK
in high dimensions, when properly initialized (Huijser et al., 2022).

AR𝑎𝑣𝑔 , which remains in an acceptable range also when 𝑁 = 𝑑 + 1. This
indicates that AR is a qualitative indicator of samplers’ performance,
which cannot be used to detect MCMC convergence problems.
Finally, the average effective sample size (Fig. C.6) generally exhibits a linear behavior in the range 𝑁 ≥ 20, while its value below
20 should be considered as unreliable for DE and DE-SNK due to theconvergence problems previously explained. ESS values are constantly
higher for all samplers when chains are initialized in a high-likelihood
region, with DE-based samplers exhibiting the best performance. Conversely, the random initialization appears more detrimental for thesealgorithms than AI-based strategies.
3.1.2. The 30-dimensional Rosenbrock distributions
The calculated EMD and 𝜏𝑎𝑣𝑔 for the 30-dimensional Rosenbrock
distributions are shown in Figs. 4 and 5, respectively. In broad terms,
results replicate what was already observed for the 10-dimensional
distributions, though with some meaningful differences. On one hand,DE-based samplers with 𝑁 = 𝑑 + 1 (i.e., 31) chains exhibit low andasymptotic autocorrelation times (Fig. 5), but the analysis of EMD
again reveals their convergence to the wrong stationary distribution.
Relatively large EMD and 𝜏𝑎𝑣𝑔 values are instead observed for AI-based
strategies, thus suggesting a delayed convergence to the target density.
This suggests that 𝜏𝑎𝑣𝑔 is a good convergence indicator for AI also for asmall number of chains, but cannot detect the ergodicity issues of DE.
Convergence problems are exacerbated when chains are initialized
in random positions or low likelihood regions. The former and thelatter are particularly detrimental for DE-based and AI-based samplingstrategies, respectively. These behaviors were already encountered in
the 10-dimensional problems. At the same time, increasing the numberof chains (𝑁 ≥ 100) leads only to minor improvements.

3.1.3. Key findings• Ensembles constituted by 𝑁 = 𝑑+1 chains evolved using DE-based
strategies do not converge to the true stationary distribution,
independently from the posterior shape and dimensions, andinitialization.
8

Journal of Hydrology 624 (2023) 129822

G. Brunetti et al.

Fig. 5. Average integrated autocorrelation time 𝜏𝑎𝑣𝑔 as a function of algorithm step for the 30-dimensional Classic and Shifted-Rotated Rosenbrock distributions. 𝜏𝑎𝑣𝑔 is calculated
for each ensemble sampler listed in Table 1 by initializing 𝑁 = (11, 20, 50, 100, 200) chains in high and low likelihood regions, and in random positions.

• AI-based samplers maintain their ergodic properties also for 𝑁 =𝑑 + 1, but exhibit delayed convergence if initialized far from thebulk of the posterior.
• The convergence can be retained and improved for all ensemble
samplers by using a higher number of chains initialized in ahigh-likelihood region.
• In these conditions, DE-based samplers outperform AI-based
strategies, especially in high dimensions.

by algorithm AI with an ensemble of 100 chains is shown. The resulting
distribution is characterized by narrow correlated densities for soil hydraulic parameters. Conversely, as expected, 𝜎𝑄 and 𝜎𝜃 do not interact
visibly with other parameters.
Evaluation metrics for this scenario are summarized in Fig. 6, whichshows the average acceptance rate, integrated autocorrelation time,and effective sample size as calculated by selected samplers using a
different number of chains. AR𝑎𝑣𝑔 decreases for all samplers as thechains explore the bulk of the posterior, but remains consistently higher
for algorithms AI and AI-SNK. The very low acceptance rates for DEbased samplers indicates that proposals described by Eqs. (2) and (4)are not particularly suited for this low-dimensional inverse problem,
leading to autocorrelation times that are appreciably higher than theircounterparts for the affine invariant samplers. As a result, these latterare able to draw more samples from the target distribution [(G) in
Fig. 6]. Increasing the number of chains has only marginal effects on
the samplers’ efficiency for this modeling scenario, but increases theeffective sample size.
Notwithstanding the differences, all samplers satisfactorily retrieve
the original soil hydraulic parameters used to generate synthetic data
(Table D.4). Estimated percentiles differ only slightly among algorithms(Fig. D.10), with DE-based samplers exhibiting marginally wider posteriors (i.e., 𝐾𝑠 in Table D.4). This difference can be partially explainedby their higher 𝜏 [(D), (E), (F) in Fig. 6], and thus delayed convergence,
which makes their final posterior estimates less precise. It is worth
noting that 5000 steps were barely sufficient to observe the desired
asymptotic behavior of 𝜏 despite having started all samplers in the
known high-likelihood region.
Summarizing, the comparison between samplers for the synthetic
lysimeter scenarios suggests that:

3.2. Vadose zone hydrological model calibration
Based on the previous findings, the array of testing conditions forthe vadose zone modeling scenarios is reduced. In particular:• The number of chains is varied between 𝑁 = (50, 100, 200).
• In the synthetic scenario, chains are initialized around the position 𝒙0 = (𝜃𝑟 = 0.032, 𝜃𝑠 = 0.39, 𝛼 (1/cm) = 0.033, 𝑛 =1.45, 𝐾𝑠 (cm/day) = 98.71, 𝜎𝑄 (cm) = 0.01, 𝜎𝜃 = 0.001). Thisparameters’ set was used to generate synthetic data. In the actual scenario, chains are initialized around the global optimum
𝒙0 identified with the optimization strategy recently proposedby Brunetti et al. (2022).
• Chains are initialized by sampling from a multivariate normal distribution characterized by a diagonal covariance matrix, and withmean and standard deviation equal to 𝒙0 and 0.01𝒙0 , respectively.
3.2.1. Synthetic scenario
The resulting marginal distribution of soil hydraulic parameters andmeasurement errors 𝜎 for the synthetic lysimeter scenario is shown in
Fig. D.11. For visualization simplicity, only the distribution calculated
9

Journal of Hydrology 624 (2023) 129822

G. Brunetti et al.

Fig. 6. Synthetic scenario: average acceptance rate AR𝑎𝑣𝑔 [(A), (B), (C)] and integrated autocorrelation time 𝜏𝑎𝑣𝑔 [(D), (E), (F)] as a function of algorithm step calculated by theselected ensemble samplers (Table 1) using a different number of chains 𝑁; (G): Average effective sample size ESS𝑎𝑣𝑔 as a function of the chains number 𝑁 calculated by thedifferent ensemble samplers.

• Even for a low-dimensional calibration problem (i.e., 7 parameters), ensemble samplers require thousands of steps to achieve aformal convergence and correctly quantify the parameters’ uncertainty (i.e., a sample from the unique invariant distribution thatis equal to the posterior).
• Affine invariant ensemble samplers outperformed DE-based samplers, which exhibited delayed convergence.
• Increasing the number of chains had only marginal effects on thesamplers’ convergence, but increased the effective sample size.

Estimated parameters well reflect the textural composition of thevadose zone profile in Elsbethen (Table D.3). The fine-textured topsoil
layer is characterized by very low values of the shape parameters 𝛼 and𝑛, which increase with depth as the sand content becomes dominant
in soil layer 3. A similar trend is observed for the saturated hydraulicconductivity 𝐾𝑠 , whose low value for the topsoil layer is partiallyexplained by the high organic content and a low fraction of coarse
particles, which are also correlated with saturated water content values.
In particular, 𝜃𝑠3 is relatively low and should be considered more as
an effective value as it reflects the low variability of the measured
volumetric water content for this soil horizon [(D) in Fig. D.9]. Finally,
the estimated soil moisture measurement error 𝜎𝜃 is around 0.01, which
is a realistic value for the time domain reflectometer. Overall, estimatedparameter values appear plausible, thus confirming the reliability of thecalibration procedure and the adequacy of the model structure.
The model fitting performance is summarized in Fig. 7 which shows
a comparison between measured (black circles) and modeled (graylines) volumetric water contents at four different soil depths obtainedby a random sampling of 200 solutions from the posterior parameters distribution calculated by algorithm AI using 100 chains. Model
predictive checks highlight the good predictive accuracy and limited
predictive uncertainty of the calibrated model, thus further confirming
the reliability of the estimated parameters.

3.2.2. Actual scenario
As this scenario includes actual data, the discussion about thestatistical comparison between ensemble samplers is preceded by an assessment of the estimated parameters and a verification of the model’s
predictive accuracy to confirm that the chosen model structure isappropriate and that parameters are physically plausible.
Parameters assessment and model predictive performance. Fig. D.13
shows the marginal distribution of the measurement error 𝜎𝜃 and soil
hydraulic parameters of different layers for the actual scenario. As forthe synthetic case study, for visualization simplicity, only the marginal
distribution calculated by algorithm AI using 100 chains is shown. As
is often the case in vadose zone modeling, the estimated marginal
distributions are characterized by elongated narrow high-likelihood
regions that reflect the high correlation between some parameters(e.g., 𝜃𝑟1 − 𝑛1 , 𝛼2 − 𝑛1 , 𝜃𝑠3 − 𝐾𝑠3 ). Posterior parameters distributions
are generally updated from the uniform prior, thus indicating that soil
moisture measurements were sufficiently informative. Nevertheless, theuncertainty remains high for some parameters, such as 𝜃𝑟1 and 𝐾𝑠3 forwhich data were mainly informative of the wet and dry parts of theretention and conductivity curves, respectively.
A more quantitative assessment of the parameters’ uncertainty isreported in Table D.5, which lists the percentiles of the posterior
distribution for each parameter calculated by different samplers using
100 chains. Percentile values are generally similar among different
algorithms (Fig. D.12), thus confirming their convergence to the sametarget density. The highest differences are observed for the saturated
conductivities 𝐾𝑠 of soil layers 2 and 3, and for the residual watercontent 𝜃𝑟 of the topsoil layer.

Ensemble sampler comparison. Evaluation metrics for the actual scenario are summarized in Fig. 8, which shows the average acceptance
rate, integrated autocorrelation time, and effective sample size as calculated by selected samplers using a different number of chains. Essentially, the increased dimensionality of the actual scenario exacerbatessome of the convergence traits observed in the synthetic modeling. The
average acceptance rate AR𝑎𝑣𝑔 decreases exponentially for all samplers,but more prominently for DE-based samplers, which exhibit a very highrejection rate towards the end of the MCMC run. Algorithms AI and AISNK have slightly better ARs, but it is evident that the dimensionality
increase (from 7 to 16) strongly impairs their performance, as reflectedby the integrated autocorrelation time, which is now higher for AI and
AI-SNK.
As already noted in the other scenarios, a higher number of chains
increases the effective sample size [(G) in Fig. 8] but has only minor effects on the algorithm convergence, which remains difficult to formally
10

Journal of Hydrology 624 (2023) 129822

G. Brunetti et al.

Fig. 7. Actual scenario: a comparison between measured (black circles) and modeled (gray lines) volumetric water contents at four different soil depths obtained by a random
sampling without replacement of 200 solutions from the posterior parameters distribution calculated by algorithm AI using 100 chains.

Fig. 8. Actual scenario: average acceptance rate AR𝑎𝑣𝑔 [(A), (B), (C)] and integrated autocorrelation time 𝜏𝑎𝑣𝑔 [(D), (E), (F)] as a function of algorithm step calculated by theselected Ensemble samplers (Table 1) using a different number of chains 𝑁; (G): Average effective sample size ESS𝑎𝑣𝑔 as a function of the chains number 𝑁 calculated by thedifferent Ensemble samplers.

11

Journal of Hydrology 624 (2023) 129822

G. Brunetti et al.

who analyzed in-depth the behavior of the stretch move described in
Eq. (6). As the number of dimension increases, the volume exterior
to the typical set overwhelms the volume interior to it, and almost
every ensemble samplers proposal generates points lying outside ornear the tails of the typical set (Betancourt, 2017; Cover and Thomas,
2005; Robert and Casella, 2004). However, the density of these points
is so small that almost all proposals are rejected and the Markov chain
rarely moves (Fig. 8). This explains the performance degradation for all
samplers when transitioning from 7 to 16 parameters, and suggests thatensemble samplers should be used cautiously when the dimensions of
the target distributions exceeds 10 (rule of thumb).
DE-based samplers exhibit shorter autocorrelation times and lower
acceptance rate compared to affine-invariant samplers in the actual
scenario. This happens mainly because affine moves are less audacious,
but this timidity leads to delayed convergence (Hogg and ForemanMackey, 2018). This dual behavior can be altered by modifying thescale parameter 𝛾 and a in Eqs. (2) and (7), respectively. This hints
at the possibility that dynamically adapting scale parameters 𝛾 anda could improve samplers’ performance. Intuitively, these parameters
should be set to generate longer jumps in the parameter space atthe beginning of the Monte Carlo run to quickly move the ensemble
towards the bulk of the posterior (lowering 𝜏), but then they should
be gradually reduced to increase the acceptance rate, especially inhigh dimensions. However, the adaptation phase should have a limited
extension to preserve the ergodicity and the detailed balance (Andrieu
and Thoms, 2008).
Alternatively, other MCMCs, such as the Metropolis Adjusted
Langevin algorithm (Roberts and Tweedie, 1996), Hamiltonian Monte
Carlo (Neal, 2010), and the Barker proposal (Livingstone and Zanella,2022) can generate very efficient proposals by eliciting gradient information. However, a precise gradient estimation for nonlinear partial
differential equations such as the Richards equation (Eq. (12)) is
challenging, as it requires the use of full or adjoint sensitivity methods
that are difficult to implement (Bitterlich and Knabner, 2002). An approximation can be obtained through numerical differentiation, whichis, however, computationally expensive and prone to numerical errors.
Albeit these challenges, gradient-based samplers are supposed to outperform ensemble samplers for moderate and high-dimensional inverse
problems, though there have been no studies proving this for vadosezone model calibration. Recent advances in automatic differentiation
provided by machine learning allow reimplementing an existing model
to obtain a differentiable model version. This opens new opportunities
in using gradient-based samplers for model calibration (Shen et al.,2023) and sensitivity analysis (Pérez et al., 2006; Xiao and Nowak,
2022; Asheghi et al., 2020).
Finally, Tables D.4 and D.5 show that the largest differences inthe estimated parameters’ percentiles between samplers are observed
for those factors that exhibit the largest uncertainty (e.g., 𝐾𝑠 for thesynthetic scenario; 𝜃𝑟1 , 𝐾𝑠2 , 𝐾𝑠3 for the actual scenario). This behavior
is mainly caused by the initial position of the chains, which is distantfrom the target posterior for these parameters. Consequently, the burnin period is extended (larger 𝜏) and different among ensemble samplers,
which deliver slightly different estimates of the posterior’s percentiles.
Albeit irrelevant for this study’s purposes, these discrepancies might
become meaningful for real vadose zone modeling scenarios where anunder or overestimation of the uncertainty of some parameters caninfluence the decision-making process (e.g., overestimation of 𝐾𝑠 forsurface runoff risk assessment). Practically, it is recommended to use
a conservative approach when data are supposed to be uninformative
for some parameters, preferably by performing additional ad-hoc measurements (e.g., field disk infiltrometer to characterize 𝐾𝑠 ) when using
the model for predictive purposes that are far from the calibration
conditions. If this is not possible for some operative reasons, thenit is advisable to rely on a practical worst-case scenario and alter
the posterior (e.g., skew the 𝐾𝑠 distribution towards lower values forsurface runoff risk assessment). Ensemble samplers are very useful
to calibrate vadose zone models and quantify the uncertainty, however, their outcomes should be cautiously evaluated and then used for
decision-making by keeping in mind the purpose of the modeling,achieve. Indeed, 10,000 steps, which amount to 2 million HYDRUS
executions when using 200 chains, are not sufficient to have a precise
estimation of 𝜏𝑎𝑣𝑔 [(D), (E), (F) in Fig. 8]. This has some practical
implications for vadose zone modeling, which are discussed later.
The results of the comparison between ensemble samplers for theactual scenario can be summarized as follows:• The higher dimensionality of the field soil case has strong negative effects on the performance of all samplers, though the performance decline is more evident for the affine invariant ensemble
samplers, which are outperformed by DE-based strategies.
• Increasing the number of chains leads only to minor improvements.
• Despite the delayed convergence, all samplers retrieve similar
posterior distributions, and lead to successful model calibration.
4. Discussion
The analysis of the performance of ensemble samplers on multiple
scenarios opens some points of discussion on their applicability forvadose zone modeling. Their performance on toy distributions hasrevealed that differential evolution strategies tend to lose their ergodic
properties and converge to the wrong posterior if the ensemble iscomposed by 𝑁 = 𝑑 +1 chains, with 𝑑 being the number of dimensions.
Conversely, affine invariant samplers remain ergodic and reversible,
but their convergence is delayed. The number of chains has to be
increased to values around 𝑁 ≥ 2𝑑 to restore the ergodic propertiesof DE-based samplers.
Increasing the number of chains above this value has negligible
influence on samplers’ performance, thus indicating that the multichain
effect does not play an important role. This behavior is also encountered in vadose zone modeling scenarios. It follows that a high numberof chains should be preferred when parallel computing infrastructures
are available (Foreman-Mackey et al., 2013). Conversely, if the computational resources are limited, an ensemble constituted by 𝑁 ≈ 2𝑑
chains can be used, but has to be evolved for a longer number ofsteps to obtain comparable values of the effective sample size ESS.
Our experiences suggest that sampling from the invariant distribution(e.g., after burn-in) for few autocorrelation times is sufficient to obtaina reliable evaluation of the parameters’ uncertainty for vadose zone
problems.
Initializing chains in a high-likelihood region has proven to facilitate
the convergence of ensemble samplers to the true posterior. In these
conditions, the scheme simultaneously satisfies the detailed balance
and is geometrically ergodic, thus leading to fast chains’ mixing (Jones,
2004). It is therefore advisable to use first a global search algorithm
to identify the bulk of the posterior, and then initialize chains in thisregion. It can be speculated that this is particularly effective when theprior is uninformative (e.g., uniform), as it avoids starting chains in
low likelihood regions. In these conditions, chains’ trajectories cannot
be corrected by the prior, and will require multiple steps to traverse
the parameter space and reach the posterior. This delayed convergence
becomes particularly problematic when the model execution time increases. This happens frequently when the spatial domain becomes
multidimensional, or multiple processes are simulated (e.g., solutetransport). Surrogate models (Brunetti et al., 2017, 2020b) can partiallymitigate this problem, however their accuracy degrades when thenumber of dimensions increases (Razavi et al., 2012).
Another important result is that ensemble samplers are able tosuccessfully calibrate the HYDRUS model on both the synthetic and actual scenarios, although they exhibit slightly different behaviors. Affine
invariant samplers outperformed DE-based samplers in the synthetic
scenario (Fig. 6) (i.e., 7 parameters), but exhibited slower convergence
for the moderate-dimensional actual scenario (Fig. 8) (i.e., 16 parameters), despite the very low acceptance rates observed for DE-based
samplers. This problem was already reported by Huijser et al. (2022),
12

Journal of Hydrology 624 (2023) 129822

G. Brunetti et al.

Fig. B.1. (A) Contours of the Gaussian function as defined in Eq. (B.1); (B) Contours of the transformed Gaussian function as defined in Eq. (B.3).

Fig. C.2. A comparison between the reference Classic Rosenbrock distribution (Eq. (10)) (black contours), and the marginal distributions approximated by DE (red contours) andAI (blue contours) when using an ensemble of 11 chains. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of thisarticle.)13

Journal of Hydrology 624 (2023) 129822

G. Brunetti et al.

Fig. C.3. A comparison between the reference Shifted-Rotated Rosenbrock distribution (Eq. (11)) (black contours), and the marginal distributions approximated by DE (red contours)
and AI (blue contours) when using an ensemble of 11 chains. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of thisarticle.)• The use of ensemble samplers is discouraged for problems involving more than 20 calibrated parameters. If used, a thoroughassessment of autocorrelation time is recommended to obtainreliable estimates of the parameters’ uncertainty.
• High 𝑁 are preferable if parallel computing infrastructures areavailable. Vice versa, low values of 𝑁 can still be used, but chains
have to be evolved longer to obtain comparable values of theeffective sample size.

5. Recommendations
The discussion of the results allows us to propose some recommendations on the use of ensemble samplers for vadose zone modeling.
However, the following should be intended as general guidelines to be
tailored to the specific calibration task:• The Rosenbrock distribution is a valid surrogate benchmark totest the performance of MCMC methods for vadose zone problems. The Python code, as well as the shift vector and the rotation
matrices, can be downloaded at this link (https://github.com/
giubrune/Paper_EnsembleMCMC).
• The number of chains must be 𝑁 ≥ 2𝑑, where 𝑑 is the number ofdimensions.
• It is recommended to initialize chains in a high-likelihood region.
This can be achieved by using first a global search algorithm.
• The affine-invariant ensemble sampler should be preferred forvadose zone calibration problems whose dimensionality does notexceed 10 parameters, while DE-based strategies might be preferable when the number of the parameters is between 10 and
20.

6. Conclusions
The main objective of this study was to provide a comprehensive
analysis of the performance of ergodic and reversible Monte Carlo
Markov Chain ensemble samplers for vadose zone model calibration. A
variegated array of testing scenarios, encompassing toy functions, andsynthetic and actual experiments, was used to statistically compare four
different algorithms, and understand limitations and opportunities intheir use for vadose zone modeling.
The analysis revealed that all samplers could be successfully used
to calibrate soil hydraulic parameters, though their performance isinfluenced by the dimensionality of the inverse problem, the numberof chains, and their initializations. These issues can exacerbate the14

Journal of Hydrology 624 (2023) 129822

G. Brunetti et al.

Fig. C.4. Average acceptance rate AR𝑎𝑣𝑔 as a function of algorithm step for the 10-dimensional Classic and Shifted-Rotated Rosenbrock distributions. AR𝑎𝑣𝑔 is calculated for eachensemble sampler listed in Table 1 by initializing 𝑁 = (11, 20, 50, 100, 200) chains in high and low likelihood regions, and in random positions.

Writing – review & editing. Thomas Wöhling: Validation, Writing –review & editing. Christine Stumpp: Resources, Writing – review &editing.

computational burden of the model calibration, thus hindering thewidespread application of these algorithms to more complex problems(e.g., reactive solute transport, 2D/3D numerical domains). Results
suggest that more research is needed to develop solutions to overcome
these limitations (e.g., parallel computing) and extend the applicabilityof ensemble samplers to vadose zone model calibration. In addition,other algorithms, such as the Hamiltonian Monte Carlo or the Nested
sampling, should be tested and compared to ensemble samplers to
identify under which conditions one is preferable to the other. Thisis a limitation of the present study that should be addressed in future
research. In this view, the set of inverse problems should be broadened
to include more pronounced multimodal posteriors, which are difficult
to sample for MCMC algorithms.
Finally, we want to emphasize that ensemble samplers should be
intended as a statistically rigorous tool to calibrate the model, understand its adequacy, and assess its predictive capability. However,even the most comprehensive and powerful numerical tool cannot do
without accurate measurements, in-depth knowledge of processes, andabove all, a humble approach to modeling a complex system such as
the vadose zone.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared toinfluence the work reported in this paper.
Data availability
The manuscript includes a link to download data used for the
analysis of the Rosenbrock function and synthetic scenario
Acknowledgments
Meteorological data were provided by the Central Institution for
Meteorology and Geodynamics (ZAMG). The shift vector, rotation matrices, and initial position of the chains for the Rosenbrock function, as well as the HYDRUS-1D model used in the synthetic scenario
can be downloaded at this link (https://github.com/giubrune/Paper_
EnsembleMCMC).

CRediT authorship contribution statement
Giuseppe Brunetti: Conceptualization, Methodology, Writing –original draft, Formal analysis. Jiri Šimunek: Supervision, Software,
15

Journal of Hydrology 624 (2023) 129822

G. Brunetti et al.

Fig. C.5. Average acceptance rate AR𝑎𝑣𝑔 as a function of algorithm step for the 30-dimensional Classic and Shifted-Rotated Rosenbrock distributions. AR𝑎𝑣𝑔 is calculated for eachensemble sampler listed in Table 1 by initializing 𝑁 = (31, 50, 100, 200) chains in high and low likelihood regions, and in random positions.

Appendix A. Detailed balance and ergodicity
Table D.3
Measured particle size distribution and bulk density of different soil horizons at
Elsbethen monitoring station.

We are interested in building numerical schemes that are able tosample from a posterior distribution that is invariant and unique.
Invariant distribution. A distribution 𝜋(𝑥) is invariant (or stationary)over the states of a Markov chain with transition probabilities 𝑇𝑛 (𝑥, 𝑥′ )if:∑
𝜋(𝑥) =𝜋(𝑥)𝑇
̃ 𝑛 (𝑥,
̃ 𝑥)
(A.1)

Soil
horizon

Depth
(cm)

Sand
(%)

Silt
(%)

Clay
(%)

𝜌𝑏
(g/cm3 )Organic
content (%)

1
2
3

0–20
20–50
50–150

33
26.5
84

61
68.5
14

6
5
2

1.2
1.3
1.3

6.3
1.45
0.3

𝑥̃

This implies that 𝜋(𝑥) persists forever once it is reached. A Markov
chain that satisfies the more restrictive detailed balance condition, will
sample from an invariant distribution. Detailed balance requires that theprobabilities of transitioning from a state 𝑥 to state 𝑥′ and vice versa
are the same:
𝜋(𝑥)𝑇 (𝑥, 𝑥′ ) = 𝜋(𝑥′ )𝑇 (𝑥′ , 𝑥)𝑥̃MCMC methods with affine invariance properties exhibit performance that are independent of the aspect ratio when sampling highly
anisotropic distributions (Christen and Foxy, 2010). A classic example
is the following probability density:()−(𝑥1 − 𝑥2 )2 (𝑥1 + 𝑥2 )2𝜋(𝑥) ∝ exp
−
(B.1)
2𝜈
2

(A.2)

which implies that 𝜋 is invariant, since:
∑
∑
∑
𝜋(𝑥)𝑇
̃ (𝑥,
̃ 𝑥) =
𝜋(𝑥)𝑇 (𝑥, 𝑥)̃ = 𝜋(𝑥)𝑇 (𝑥, 𝑥)
̃ = 𝜋(𝑥)𝑥̃Appendix B. Affine invariant property

(A.3)

𝑥̃

Ergodic Markov chains. It is not sufficient to build a scheme thatconverges to stationary distribution. It is also required that the Markov
chain converges to the invariant distribution as 𝑛 → ∞ regardless of
the initial state of the chain. In this case, the chain is ergodic and will
sample from a unique distribution. The proof is reported in Neal (1993)
and not repeated here for conciseness.

Single-variable MCMC methods would be forced to make movements√of order 𝜈 to sample the distribution [(A) in Fig. B.1], while more
efficient MCMC schemes would use these perturbations only in one
16

Journal of Hydrology 624 (2023) 129822

G. Brunetti et al.

Fig. C.6. Average effective sample size ESS𝑎𝑣𝑔 as a function of the chains number 𝑁. ESS𝑎𝑣𝑔 is calculated for each ensemble samplers by initializing chains in high and low
likelihood regions, and in random positions.

Fig. D.7. Measured potential evapotranspiration (green line) and precipitation (blue line) used in HYDRUS for both synthetic and actual scenarios. The hatched area covers themodel spin-up period. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)17

Journal of Hydrology 624 (2023) 129822

G. Brunetti et al.

Fig. D.8. Synthetic scenario: synthetic average volumetric water content 𝜃𝑎𝑣𝑔 (A) and cumulative bottom outflow 𝑄 (B) used to inversely estimate soil hydraulic parameters andmeasurements’ errors 𝜎𝜃 and 𝜎𝑄 .

Table D.4
Synthetic scenario: bounds of the uniform prior distributions, and percentiles of the posterior distribution calculated by different samplers using 100 chains. Statistics are calculated
by excluding the burn-in phase.
Synthetic scenario
Parameter

Prior bounds𝜃𝑟 (cm3 ∕cm3 )𝜃𝑠 (cm3 ∕cm3 )𝛼 (1∕cm)𝑛 (−)
𝐾𝑠 (cm∕day)𝜎𝑄 (cm)𝜎𝜃 (cm3 ∕cm3 )(0, 0.1)(0.3, 0.5)
(0.001, 0.1)
(1.01, 2.5)
(10, 500)(0.001, 0.1)
(1E−05, 1E−02)

AI

DE

AI-SNK

DE-SNK

5%

50%

95%

5%

50%

95%

5%

50%

95%

5%

50%

95%

0.032
0.388
0.032
1.447
89.2
0.009
0.001

0.032
0.39
0.033
1.451
97.5
0.01
0.001

0.033
0.391
0.033
1.46
102.6
0.011
0.001

0.027
0.388
0.032
1.425
89
0.009
0.001

0.032
0.389
0.033
1.453
95.9
0.01
0.001

0.034
0.398
0.038
1.46
140.9
0.013
0.001

0.032
0.388
0.032
1.447
89.5
0.009
0.001

0.032
0.39
0.033
1.452
96.9
0.01
0.001

0.033
0.391
0.034
1.46
102.7
0.011
0.001

0.03
0.388
0.032
1.433
89.8
0.009
0.001

0.032
0.39
0.033
1.45
98.4
0.01
0.001

0.033
0.394
0.036
1.459
120.2
0.011
0.001

Table D.5
Actual scenario: bounds of the uniform prior distributions, and percentiles of the posterior distribution calculated by different samplers using 100 chains. Statistics are calculated
by excluding the burn-in phase. Subscripts 1, 2, and 3 refer to the soil horizons as indicated in Table D.3.
Actual scenario
Parameter

Prior bounds𝜃𝑟1 (cm3 ∕cm3 )𝜃𝑠1 (cm3 ∕cm3 )𝛼1 (1∕cm)
𝑛1 (−)
𝐾𝑠1 (cm∕day)𝜃𝑟2 (cm3 ∕cm3 )𝜃𝑠2 (cm3 ∕cm3 )𝛼2 (1∕cm)
𝑛2 (−)
𝐾𝑠2 (cm∕day)𝜃𝑟3 (cm3 ∕cm3 )𝜃𝑠3 (cm3 ∕cm3 )𝛼3 (1∕cm)
𝑛3 (−)
𝐾𝑠3 (cm∕day)𝜎𝜃 (cm3 ∕cm3 )(0, 0.2)(0.2, 0.6)
(0.001, 1)(1.01, 3)(0.1, 2000)(0, 0.2)(0.2, 0.6)
(0.001, 1)(1.01, 3)(0.1, 2000)
0, 0.2)(0.2, 0.6)
(0.001, 1)(1.01, 3)(0.1, 10000)(0.001, 0.05)

AI

DE

AI-SNK

50%

95%

5%

50%

95%

5%

50%

95%

5%

50%

95%

0.09
0.424
0.003
1.11
0.27
0.004
0.388
0.004
1.42
16.4
0.048
0.21
0.19
1.51
2716
0.01

0.14
0.426
0.004
1.13
0.32
0.053
0.395
0.005
1.58
43.7
0.055
0.24
0.39
1.65
5197
0.011

0.19
0.428
0.005
1.17
0.35
0.135
0.401
0.007
1.77
101.8
0.06
0.29
0.72
1.87
8871
0.011

0.03
0.425
0.003
1.08
0.29
0.004
0.389
0.004
1.43
25.4
0.048
0.21
0.15
1.51
2922
0.01

0.12
0.426
0.004
1.12
0.32
0.039
0.394
0.005
1.57
76.5
0.055
0.24
0.43
1.65
6226
0.011

0.19
0.429
0.005
1.16
0.36
0.1
0.403
0.007
1.73
197.6
0.059
0.29
0.89
1.89
9490
0.011

0.05
0.425
0.003
1.11
0.21
0.005
0.387
0.004
1.42
18.6
0.049
0.21
0.16
1.52
2768
0.01

0.14
0.426
0.004
1.13
0.31
0.051
0.393
0.005
1.58
38.5
0.054
0.23
0.32
1.63
4781
0.011

0.19
0.428
0.004
1.18
0.33
0.114
0.4
0.006
1.83
89.2
0.059
0.27
0.56
1.81
7449
0.011

0.014
0.425
0.003
1.08
0.22
0.007
0.387
0.004
1.45
26.7
0.049
0.21
0.19
1.52
3240
0.01

0.11
0.427
0.004
1.11
0.32
0.082
0.394
0.005
1.59
88.9
0.055
0.24
0.57
1.61
6256
0.011

0.19
0.428
0.005
1.17
0.37
0.146
0.402
0.007
1.79
290.5
0.06
0.28
0.95
1.81
9292
0.011

direction, while performing larger steps in the other. The affine trans-

which can be solved without any particular customization by classic
MCMC methods (e.g., Metropolis–Hastings). MCMC with affine invariance properties view the two distributions as equally difficult, thus
avoiding the need to tune the proposal strategy.

formation:𝑦1 =𝑥2 − 𝑥2√ ,𝜈

𝑦2 = 𝑥1 + 𝑥2

(B.2)

Appendix C. The Rosenbrock distribution
See Figs. C.2 and C.6.

turns the distribution in Eq. (B.1) into [(B) in Fig. B.1]:(
𝜋𝐴 (𝑦) ∝ exp−(𝑦21 + 𝑦22 )2DE-SNK

5%

Appendix D. Vadose zone model calibration

)
(B.3)

See Figs. D.7–D.13 and Tables D.3 and D.5.
18

Journal of Hydrology 624 (2023) 129822

G. Brunetti et al.

Fig. D.9. Actual scenario: measured volumetric water contents at four different depths used to inversely estimate the soil hydraulic parameters of different soil horizons and themeasurement error 𝜎𝜃 .

Fig. D.10. Synthetic scenario: Box plots of the soil hydraulic parameters and measurements’ error calculated by ensemble samplers in Table 1 using 100 chains. Statistics arecalculated by excluding the burn-in phase.

19

Journal of Hydrology 624 (2023) 129822

G. Brunetti et al.

Fig. D.11. Synthetic scenario: marginal distribution of the soil hydraulic parameters and measurements’ error calculated by algorithm AI using 100 chains. Statistics are calculated
by excluding the burn-in phase.

20

Journal of Hydrology 624 (2023) 129822

G. Brunetti et al.

Fig. D.12. Actual scenario: Box plots of the soil hydraulic parameters and measurements’ error calculated by ensemble samplers in Table 1 using 100 chains. Statistics are calculated
by excluding the burn-in phase.

21

Journal of Hydrology 624 (2023) 129822

G. Brunetti et al.

Fig. D.13. Actual scenario: marginal distribution of the soil hydraulic parameters and measurement error calculated by algorithm AI using 100 chains. Statistics are calculated by
excluding the burn-in phase.

References

Brunetti, G., Šimůnek, J., Glöckler, D., Stumpp, C., 2020b. Handling model complexity
with parsimony: Numerical analysis of the nitrogen turnover in a controlled aquifer
model setup. J. Hydrol. 584, 124681. http://dx.doi.org/10.1016/j.jhydrol.2020.
124681.

Allen, R.G., Pereira, L.S., Raes, D., Smith, M., 1998. FAO Irrigation and Drainage Paper
No. 56: Crop Evapotranspiration. FAO, Rome, p. 300. http://dx.doi.org/10.1016/j.
eja.2010.12.001, arXiv:arXiv:1011.1669v3.

Brunetti, G., Šimůnek, J., Turco, M., Piro, P., 2017. On the use of surrogate-basedmodeling for the numerical analysis of Low Impact Development techniques. J.
Hydrol. 548, 263–277. http://dx.doi.org/10.1016/j.jhydrol.2017.03.013, URL: http://linkinghub.elsevier.com/retrieve/pii/S0022169417301506.

Andrieu, C., Thoms, J., 2008. A tutorial on adaptive MCMC. Stat. Comput. 18 (4),343–373. http://dx.doi.org/10.1007/s11222-008-9110-y.
Arora, B., Mohanty, B.P., McGuire, J.T., 2012. Uncertainty in dual permeability model
parameters for structured soils. Water Resour. Res. 48 (1), http://dx.doi.org/10.
1029/2011WR010500.

Brunetti, G., Stumpp, C., Šimůnek, J., 2022. Balancing exploitation and exploration: A
novel hybrid global-local optimization strategy for hydrological model calibration.
Environ. Model. Softw. 150, 105341. http://dx.doi.org/10.1016/j.envsoft.2022.
105341.

Asheghi, R., Hosseini, S.A., Saneie, M., Shahri, A.A., 2020. Updating the neural network
sediment load models using different sensitivity analysis methods: A regional
application. J. Hydroinform. 22 (3), 562–577. http://dx.doi.org/10.2166/hydro.
2020.098.

Christen, J.A., Foxy, C., 2010. A general purpose sampling algorithm forcontinuous distributions (the t-walk). Bayesian Anal. 5 (2), 263–282. http://dx.doi.org/10.1214/10-BA603, URL: https://projecteuclid.org/journals/bayesiananalysis/volume-5/issue-2/A-general-purpose-sampling-algorithm-for-continuousdistributions-the-t/10.1214/10-BA603.full.

Betancourt, M., 2017. A conceptual introduction to Hamiltonian Monte Carlo. URL:
http://arxiv.org/abs/1701.02434, arXiv:1701.02434.
Bitterlich, S., Knabner, P., 2002. An efficient method for solving an inverse problem
for the Richards equation. J. Comput. Appl. Math. 147 (1), 153–173. http://dx.
doi.org/10.1016/S0377-0427(02)00430-2.

Cover, T.M., Thomas, J.A., 2005. Elements of Information Theory. John Wiley and Sons,
pp. 1–748. http://dx.doi.org/10.1002/047174882X, URL: https://onlinelibrary.
wiley.com/doi/book/10.1002/047174882X.

Brunetti, G., Kodešová, R., Šimůnek, J., 2019a. Modeling the translocation andtransformation of chemicals in the soil-plant continuum: A dynamic plant uptake
module for the HYDRUS model. Water Resour. Res. http://dx.doi.org/10.1029/
2019WR025432.

Das, N.N., Mohanty, B.P., Njoku, E.G., 2008. A Markov chain Monte Carlo algorithm for
upscaled soil-vegetation-atmosphere-transfer modeling to evaluate satellite-based
soil moisture measurements. Water Resour. Res. 44 (5), http://dx.doi.org/10.1029/
2007WR006472, URL: http://doi.wiley.com/10.1029/2007WR006472.

Brunetti, G., Papagrigoriou, I.-A., Stumpp, C., 2020a. Disentangling model complexity
in green roof hydrological analysis: A Bayesian perspective. Water Res. 115973.
http://dx.doi.org/10.1016/j.watres.2020.115973.

Durner, W., 1994. Hydraulic conductivity estimation for soils with heterogeneous
pore structure. Water Resour. Res. 30 (2), 211–223. http://dx.doi.org/10.1029/
93WR02676.

Brunetti, G., Šimůnek, J., Bogena, H., Baatz, R., Huisman, J.A., Dahlke, H.,
Vereecken, H., 2019b. On the information content of cosmic-ray neutron datain the inverse estimation of soil hydraulic properties. Vadose Zone J. 18 (1),
http://dx.doi.org/10.2136/vzj2018.06.0123.

Durner, W., Jansen, U., Iden, S.C., 2008. Effective hydraulic properties of layeredsoils at the lysimeter scale determined by inverse modelling. Eur. J. Soil Sci.
http://dx.doi.org/10.1111/j.1365-2389.2007.00972.x.
22

Journal of Hydrology 624 (2023) 129822

G. Brunetti et al.
Eck, D.V., Qin, M., Hirmas, D.R., Giménez, D., Brunsell, N.A., 2016. Relating quantitative soil structure metrics to saturated hydraulic conductivity. Vadose Zone J. 15(1), 1–11. http://dx.doi.org/10.2136/vzj2015.05.0083.
Feddes, R.A., Kowalik, P.J., Zaradny, H., 1978. Simulation of Field Water Use and Crop
Yield. PUDOC, Wageningen, p. 189.
Foreman-Mackey, D., Hogg, D.W., Lang, D., Goodman, J., 2013. emcee : The MCMC
Hammer. Publ. Astron. Soc. Pac. 125 (925), 306–312. http://dx.doi.org/10.1086/
670067, arXiv:1202.3665.
Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A., Rubin, D.B., 2013.
Bayesian Data Analysis, third ed..
van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulicconductivity of unsaturated soils1. Soil Sci. Am. J. 44 (5), 892–898. http://dx.
doi.org/10.2136/sssaj1980.03615995004400050002x.
Geyer, C.J., 1992. Practical markov chain monte carlo. Statist. Sci. 7 (4), 473–483.
http://dx.doi.org/10.1214/ss/1177011137.
Goodman, J., Weare, J., 2010. Ensemble samplers with affine invariance. Commun.
Appl. Math. Comput. Sci. 5 (1), 65–80. http://dx.doi.org/10.2140/camcos.2010.5.
65.
Hastings, W.K., 1970. Monte carlo sampling methods using Markov chains and theirapplications. Biometrika 57 (1), 97–109. http://dx.doi.org/10.1093/biomet/57.1.
97.
Hogg, D.W., Foreman-Mackey, D., 2018. Data analysis recipes: Using Markov Chain
Monte Carlo. Astrophys. J. Suppl. Ser. 236 (1), 11. http://dx.doi.org/10.3847/
1538-4365/aab76e, arXiv:1710.06068.
Huijser, D., Goodman, J., Brewer, B.J., 2022. Properties of the affine-invariant ensemble
sampler’s ‘stretch move’ in high dimensions. Aust. N. Z. J. Stat. 64 (1), 1–26.
http://dx.doi.org/10.1111/anzs.12358, arXiv:1509.02230.
Huisman, J., Rings, J., Vrugt, J., Sorg, J., Vereecken, H., 2010. Hydraulic properties ofa model dike from coupled Bayesian and multi-criteria hydrogeophysical inversion.
J. Hydrol. 380 (1–2), 62–73. http://dx.doi.org/10.1016/J.JHYDROL.2009.10.023.
Jana, R.B., Mohanty, B.P., 2012. A comparative study of multiple approaches to soil
hydraulic parameter scaling applied at the hillslope scale. Water Resour. Res. 48(2), http://dx.doi.org/10.1029/2010WR010185.
Jarvis, N., Larsbo, M., 2012. MACRO (v5.2): Model use, calibration, and validation.
Trans. ASABE 55 (4), 1413–1423. http://dx.doi.org/10.13031/2013.42251.
Jones, G.L., 2004. On the markov chain central limit theorem. Probab. Surv.
1(1),
299–320.
http://dx.doi.org/10.1214/154957804100000051,
URL:
https://projecteuclid.org/journals/probability-surveys/volume-1/issue-none/Onthe-Markov-chain-central-limit-theorem/10.1214/154957804100000051.full,arXiv:0409112.
Karamanis, M., Beutler, F., 2020. Ensemble Slice Sampling: Parallel, black-box andgradient-free inference for correlated and multimodal distributions. Stat. Comput.
31 (5), http://dx.doi.org/10.1007/s11222-021-10038-2, arXiv:2002.06212.
Kuhlman, K.L., Mills, M.M., Heath, J.E., Paul, M.J., Wilson, J.E., Bower, J.E., 2022. Parameter estimation from spontaneous imbibition into volcanic tuff. Vadose Zone J.
21 (2), e20188. http://dx.doi.org/10.1002/vzj2.20188, URL: https://onlinelibrary.
wiley.com/doi/10.1002/vzj2.20188.
Laloy, E., Linde, N., Vrugt, J.A., 2012. Mass conservative three-dimensional water
tracer distribution from Markov chain Monte Carlo inversion of time-lapse groundpenetrating radar data. Water Resour. Res. 48 (7), http://dx.doi.org/10.1029/
2011WR011238, URL: http://doi.wiley.com/10.1029/2011WR011238.
Laloy, E., Rogiers, B., Vrugt, J.A., Mallants, D., Jacques, D., 2013. Efficient posterior
exploration of a high-dimensional groundwater model from two-stage Markov chain
Monte Carlo simulation and polynomial chaos expansion. Water Resour. Res. 49(5), 2664–2682. http://dx.doi.org/10.1002/wrcr.20226.
Lehmann, F., Rajabi, M.M., Belfort, B., Delay, F., Fahs, M., Ackerer, P., Younes, A.,
2022. Bayesian inversion of laboratory experiments of transport through limestone
fractures. J. Contam. Hydrol. 249, 104045. http://dx.doi.org/10.1016/j.jconhyd.
2022.104045.
Li, X., Wen, Z., Zhan, H., Zhu, Q., Jakada, H., 2022. On the bimodal radialsolute transport in dual-permeability porous media. Water Resour. Res. 58(11), e2022WR032580. http://dx.doi.org/10.1029/2022WR032580, URL: https://
onlinelibrary.wiley.com/doi/10.1029/2022WR032580.
Liang, J.J., Qu, B.Y., Suganthan, P.N., 2014. Problem Definitions and Evaluation
Criteria for the CEC 2014 Special Session and Competition on Single Objective
Real-Parameter Numerical Optimization.
Livingstone, S., Zanella, G., 2022. The Barker proposal: Combining robustness andefficiency in gradient-based MCMC. J. R. Stat. Soc. Ser. B Stat. Methodol. 84 (2),496–523. http://dx.doi.org/10.1111/rssb.12482, URL: https://academic.oup.com/
jrsssb/article/84/2/496/7056132, arXiv:1908.11812.
Neal, R.M., 1993. Probabilistic Inference using Markov Chain Monte Carlo Methods.
Technical Report, Dept. of Computer Science, University of Toronto, Toronto.
Neal, R.M., 2010. MCMC using Hamiltonian dynamics. In: Handbook of Markov Chain
Monte Carlo, Vol. 54. pp. 113–162.
Nolz, R., Kammerer, G., Cepuder, P., 2013. Interpretation of lysimeter weighing data
affected by wind. J. Plant Nutr. Soil Sci. 176 (2), 200–208. http://dx.doi.org/10.
1002/jpln.201200342.
Pagani, F., Wiegand, M., Nadarajah, S., 2022. An <i>n</i> -dimensional Rosenbrock
distribution for Markov chain Monte Carlo testing. Scand. J. Stat. 49 (2), 657–680.
http://dx.doi.org/10.1111/sjos.12532.

Pele, O., Werman, M., 2009. Fast and robust earth mover’s distances. In: Proceedingsof the IEEE International Conference on Computer Vision. pp. 460–467. http://dx.doi.org/10.1109/ICCV.2009.5459199.
Pérez, C.J., Martín, J., Rufo, M.J., 2006. MCMC-based local parametric sensitivity
estimations. Comput. Statist. Data Anal. 51 (2), 823–835. http://dx.doi.org/10.
1016/j.csda.2005.09.005.
Razavi, S., Tolson, B.A., Burn, D.H., 2012. Review of surrogate modeling inwater resources. Water Resour. Res. 48 (7), n/a. http://dx.doi.org/10.1029/
2011WR011527, URL: http://doi.wiley.com/10.1029/2011WR011527.
Reth, S., Perez-Priego, O., Coners, H., Nolz, R., 2021. Lysimeter. In: Springer Handbooks. Springer Science and Business Media Deutschland GmbH, pp. 1583–1593.
http://dx.doi.org/10.1007/978-3-030-52171-4_58.
Reuschen, S., Jobst, F., Nowak, W., 2021. Efficient discretization-independent Bayesian
inversion of high-dimensional multi-Gaussian priors using a hybrid MCMC. Water
Resour. Res. 57 (8), e2021WR030051. http://dx.doi.org/10.1029/2021WR030051,
URL: https://onlinelibrary.wiley.com/doi/10.1029/2021WR030051.
Reuschen, S., Xu, T., Nowak, W., 2020. Bayesian inversion of hierarchical geostatisticalmodels using a parallel-tempering sequential Gibbs MCMC. Adv. Water Resour. 141,
103614. http://dx.doi.org/10.1016/j.advwatres.2020.103614.
Ritter, A., Hupet, F., Muñoz-Carpena, R., Lambot, S., Vanclooster, M., 2003. Using
inverse methods for estimating soil hydraulic properties from field data as analternative to direct methods. Agricult. Water Manag. http://dx.doi.org/10.1016/
S0378-3774(02)00160-9.
Robert, C.P., Casella, G., 2004. Monte Carlo Statistical Methods. In: Springer Texts
in Statistics, Springer New York, New York, NY, http://dx.doi.org/10.1007/978-14757-4145-2, URL: http://link.springer.com/10.1007/978-1-4757-4145-2.
Roberts, G.O., Rosenthaly, J.S., 2004. General state space markov chains and
MCMC algorithms. Probab. Surv. 1 (1), 20–71. http://dx.doi.org/10.1214/
154957804100000024,
URL:
https://projecteuclid.org/journals/probabilitysurveys/volume-1/issue-none/General-state-space-Markov-chains-and-MCMCalgorithms/10.1214/154957804100000024.full.
Roberts, G.O., Tweedie, R.L., 1996. Exponential convergence of Langevin distributions
and their discrete approximations. Bernoulli 2 (4), 341–363. http://dx.doi.org/10.
2307/3318418.
Schaap, M.G., Leij, F.J., van Genuchten, M.T., 2001. Rosetta: A computer program forestimating soil hydraulic parameters with hierarchical pedotransfer functions. J.
Hydrol. 251 (3–4), 163–176. http://dx.doi.org/10.1016/S0022-1694(01)00466-8.
Scharnagl, B., Vrugt, J.A., Vereecken, H., Herbst, M., 2011. Inverse modelling of insitu soil water dynamics: investigating the effect of different prior distributionsof the soil hydraulic parameters. Hydrol. Earth Syst. Sci. 15 (10), 3043–3059.
http://dx.doi.org/10.5194/hess-15-3043-2011.
Schelle, H., Iden, S.C., Fank, J., Durner, W., 2012. Inverse estimation of soil hydraulicand root distribution parameters from lysimeter data. Vadose Zone J. http://dx.
doi.org/10.2136/vzj2011.0169.
Schübl, M.A., Brunetti, G., Fuchs, G., Stumpp, C., 2022a. From soil water monitoringdata to vadose zone water fluxes: a comprehensive example of reverse hydrology.
Hydrol. Earth Syst. Sci. Discuss. 2022, 1–30. http://dx.doi.org/10.5194/hess-2022261, URL: https://hess.copernicus.org/preprints/hess-2022-261/.
Schübl, M., Stumpp, C., Brunetti, G., 2022b. A Bayesian perspective on the information
content of soil water measurements for the hydrological characterization of thevadose zone. J. Hydrol. 613, 128429. http://dx.doi.org/10.1016/j.jhydrol.2022.
128429.
Shen, C., Appling, A.P., Gentine, P., Bandai, T., Gupta, H., Tartakovsky, A., BaityJesi, M., Fenicia, F., Kifer, D., Li, L., Liu, X., Ren, W., Zheng, Y., Harman, C.J.,
Clark, M., Farthing, M., Feng, D., Kumar, P., Aboelyazeed, D., Rahmani, F.,
Beck, H.E., Bindas, T., Dwivedi, D., Fang, K., Höge, M., Rackauckas, C., Roy, T.,
Xu, C., Lawson, K., 2023. Differentiable modeling to unify machine learning andphysical models and advance geosciences. arXiv:2301.04027.
Šimůnek, J., Genuchten, M.T., Šejna, M., 2016. Recent developments and applicationsof the HYDRUS computer software packages. Vadose Zone J. 15 (7), 1–25. http://dx.doi.org/10.2136/vzj2016.04.0033.
Šimůnek, J., van Genuchten, M.T., 1996. Estimating unsaturated soil hydraulic properties from tension disc infiltrometer data by numerical inversion. Water Resour.
Res. http://dx.doi.org/10.1029/96WR01525.
Sokal, A., 1997. Monte Carlo Methods in Statistical Mechanics: Foundations and New
Algorithms. Springer, Boston, MA, pp. 131–192. http://dx.doi.org/10.1007/978-14899-0319-8-6.
ter Braak, C.J., 2006. A Markov Chain Monte Carlo version of the genetic algorithm
Differential Evolution: Easy Bayesian computing for real parameter spaces. Stat.
Comput. 16 (3), 239–249. http://dx.doi.org/10.1007/s11222-006-8769-1.
ter Braak, C.J., Vrugt, J.A., 2008. Differential Evolution Markov Chain with snooker
updater and fewer chains. Stat. Comput. 18 (4), 435–446. http://dx.doi.org/10.
1007/s11222-008-9104-9.
Van Dam, J.C., Groenendijk, P., Hendriks, R.F., Kroes, J.G., 2008. Advances of modelingwater flow in variably saturated soils with SWAP. Vadose Zone J. 7 (2), 640.
http://dx.doi.org/10.2136/vzj2007.0060.
Vereecken, H., Schnepf, A., Hopmans, J., Javaux, M., Or, D., Roose, T., Vanderborght, J., Young, M., Amelung, W., Aitkenhead, M., Allison, S., Assouline, S.,
Baveye, P., Berli, M., Brüggemann, N., Finke, P., Flury, M., Gaiser, T., Govers, G.,
Ghezzehei, T., Hallett, P., Hendricks Franssen, H., Heppell, J., Horn, R., Huisman, J., Jacques, D., Jonard, F., Kollet, S., Lafolie, F., Lamorski, K., Leitner, D.,23

Journal of Hydrology 624 (2023) 129822

G. Brunetti et al.

Wesseling, J., Elbers, J., Kabat, P., den Broek, B.V., 1991. SWATRE: instructions for
input.
Wöhling, T., Bidwell, V.J., Barkle, G.F., 2012. Dual-tracer, non-equilibrium mixing cell
modelling and uncertainty analysis for unsaturated bromide and chloride transport.
J. Contam. Hydrol. 140–141, 150–163. http://dx.doi.org/10.1016/j.jconhyd.2012.
08.001.
Wöhling, T., Vrugt, J.A., 2011. Multiresponse multilayer vadose zone model calibration
using Markov chain Monte Carlo simulation and field water retention data. Water
Resour. Res. 47 (4), http://dx.doi.org/10.1029/2010WR009265.
Xiao, S., Nowak, W., 2022. Reliability sensitivity analysis based on a two-stage Markov
chain Monte Carlo simulation. Aerosp. Sci. Technol. 130, 107938. http://dx.doi.
org/10.1016/j.ast.2022.107938.

McBratney, A., Minasny, B., Montzka, C., Nowak, W., Pachepsky, Y., Padarian, J.,
Romano, N., Roth, K., Rothfuss, Y., Rowe, E., Schwen, A., Šimůnek, J., Tiktak, A.,
Van Dam, J., van der Zee, S., Vogel, H., Vrugt, J., Wöhling, T., Young, I., 2016.
Modeling soil processes: Review, key challenges, and new perspectives. Vadose
Zone J. 15 (5), http://dx.doi.org/10.2136/vzj2015.09.0131, vzj2015.09.0131.
Vogel, T., van Genuchten, M.T., Cislerova, M., 2000. Effect of the shape of the soil
hydraulic functions near saturation on variably-saturated flow predictions. Adv.
Water Resour. 24 (2), 133–144. http://dx.doi.org/10.1016/S0309-1708(00)000373.
Vrugt, J.A., Beven, K.J., 2018. Embracing equifinality with efficiency: Limits ofAcceptability sampling using the DREAM(LOA) algorithm. J. Hydrol. 559, 954–971.
http://dx.doi.org/10.1016/j.jhydrol.2018.02.026.
Vrugt, J.A., Stauffer, P.H., Wöhling, T., Robinson, B.A., Vesselinov, V.V., 2008.
Inverse modeling of subsurface flow and transport properties: A review with newdevelopments. Vadose Zone J. 7 (2), 843–864. http://dx.doi.org/10.2136/vzj2007.
0078.

24