Towards a New Calibration Module for HYDRUS:Pitfalls and Opportunities in the Use of Markov-Chain
Monte Carlo Ensemble Samplers
Giuseppe Brunettia, Jirka Šimůnekb , Thomas Wohlingc, and Christine Stumppd
a Department of Civil Engineering, University of Calabria, Rendeb Department of Environmental Sciences, University of California, Riversidec Chair of Hydrology, Technical University of Dresden, Dresdend Institute of Soil Physics and Rural Water Management, University of Natural Resources and Life Sciences, Vienna

7th International Conference
HYDRUS Software Applications to Subsurface Flow and Contaminant Transport Problems
20 April 2023, Prague

WHYWHATHOWProbabilistic
calibration of vadosezone models isfrequent butchallenging

Performance ofMarkov Chain Monte
Carlo ensemble
samplers

Analytical functions,
synthetic and actual
vadose zone scenarios,
using different metrics

THE SCIENTIFIC WORKFLOW

DETERMINISTICCALIBRATIONSELECT A FITNESSMETRICITERATIVE COMPARISON
OF MODEL PREDICTIONSAND OBSERVATIONSBEST FITNESSONE SET OFPARAMETERSUNCERTAINTY

Measurement uncertainty

PROBABILISTIC CALIBRATION

Numerical errors

DISTRIBUTION OF PROBABLE PARAMETERS

LikelihoodPrior

BAYESIAN INFERENCE

Pr ( D | Ω ) Pr ( Ω )Pr ( Ω | D ) =Pr ( D )posteriorMarginal Likelihood

We have to approximate an unknown
distribution in d dimensions

MARKOV CHAIN MONTE CARLO"Markov Chain Monte Carlo schemes that are reversible and ergodic approximate theinvariant and unique posterior distribution in a finite number of steps"

https://blog.stata.com/2016/11/15/introduction-to-bayesianstatistics-part-2-mcmc-and-the-metropolis-hastings-algorithm/

CHALLENGES IN VADOSE ZONE MODELINGCHALLENGES INVADOSE ZONEMODELING
Posteriors are frequently
banana-shaped
Classic MCMC schemes
become inefficient
Gradient is difficult tocalculate
Brunetti, G., Sˇim˚unek, J., Bogena, H., Baatz, R., Huisman, J.A., Dahlke, H., Vereecken, H.,
2019. On the Information Content of Cosmic-Ray Neutron Data in the Inverse Estimation ofSoil Hydraulic Properties. Vadose Zone Journal 18, 0. doi:10.2136/vzj2018.06.0123.

ComputationallyintensiveENSEMBLE OF
MONTE CARLOMARKOV CHAINS
Differential evolution (DE)(Ter Braak, 2006)

Evolve N chains
simultaneously
Each chains is evolved
using the position ofother chains

Differential evolution Snooker (SNK)(Ter Braak and Vrugt, 2008)

Gradient-freeAffine Invariant move (AI)(Goodman and Weare, 2010)

Can be parallelizedARE THEY SUITED FOR VADOSE ZONE MODELING?

WHICH ONE IS BETTER?
WHAT ARE THE KEY FACTORS?
HOW TO USE THEM?

BENCHMARKSCENARIOSENSEMBLE SAMPLERS

Differential evolution (DE)(Ter Braak, 2006)

Differential evolution Snooker (SNK)(Ter Braak and Vrugt, 2008)

Affine Invariant move (AI)(Goodman and Weare, 2010)

BENCHMARKSCENARIOSANALYTICAL FUNCTIONS
d-dimensional classic Rosenbrock(d=10, 30 and Uniform priors)d-dimensional shifted-rotated Rosenbrock(d=10, 30 and Uniform priors)

R: Rotation matrix
o: Shift vector*Reference solution IS available  Emcee (Foreman-Mackey, 2013) with 500 chains for 10 million steps

BENCHMARKSCENARIOSVADOSE ZONE MODELING
Synthetic: Weighing Lysimeter*Uniform priors

https://israelagri.com/optimizingirrigation-and-fertilization-usinglysimeters/*Reference solution IS NOT available

Actual: Vadose zone profile in Austria*Uniform priors

BENCHMARKSCENARIOSSUMMARYEVALUATION METRICS:• EMD: Earth Movers’ Distance (Convergence to the reference distribution)• τavg: Average Integrated Autocorrelation Time (Steps required to start sampling from the truetarget distribution)• ARavg: Average Acceptance Rate (percentage by which a proposed step is accepted)• ESSavg: Average Effective Sample Size (informativeness of an ensemble of chains)

RESULTS
ANALYTICAL FUNCTIONS

d=10

•

Algorithms are sampling from astationary distribution (asymptotic τavg)•

DE-based samplers decorrelates fasterthan AI•

Increasing the number of chains leads tominor improvements•

Initializing chains in a high-likelihood
region reduces autocorrelation

RESULTSANALYTICAL FUNCTIONS

d=10•

Convergence issues for DE-based samplers
if N = d+1 chains (i.e., I1)•These issues are not detected by τavg

d=30

RESULTSANALYTICAL FUNCTIONSCORE FINDINGSANALYTICAL FUNCTIONS•Ensembles constituted by N = d + 1 chains evolved using DE-based
strategies do not converge to the true stationary distribution,
independently from the posterior shape and dimensions, andinitialization.

•The convergence can be retained and improved for all ensemble
samplers by using a higher number of chains initialized in a highlikelihood region.

•In these conditions, DE-based samplers outperform AI-based
strategies, especially in high dimensions.

•Performance of Ensemble samplers degrade in high-dimensions.

RESULTSVADOSE ZONE MODELING
Synthetic: Weighing Lysimeterhttps://israelagri.com/optimizingirrigation-and-fertilization-usinglysimeters/

Actual: Vadose zone profile in Austria

RESULTSSYNTHETIC LYSIMETER•High computational cost (up to 1
million HYDRUS executions)•AI samplers > DE samplers.

•Increasing the number of chains marginal effects on the samplers’convergence, but increases ESS.

RESULTSSYNTHETIC LYSIMETERRESULTSACTUAL VADOSE ZONE•

More parameters lead to delayedconvergence•

AI is more affected andoutperformed by DE•

Increasing the number of chainsleads only to minor
improvements.

RESULTSSYNTHETIC LYSIMETER•

Largest differences are observed for Ks and θs•

All samplers retrieve similar posterior distributions,and lead to successful model calibration.

Recommendations• The Rosenbrock distribution is a valid surrogate benchmark• The number of chains must be N≥2d, where d is the number of dimensions.
• Initialize chains in a high-likelihood regions.
• if d≤10:choose AIelse if 10<d≤20:
choose DEelse:
be careful• High N is preferable if parallel computing infrastructures are available.

THANK YOU FOR YOUR ATTENTION