Journal of Hydrology 584 (2020) 124681

Contents lists available at ScienceDirect

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

Research papers

Handling model complexity with parsimony: Numerical analysis of thenitrogen turnover in a controlled aquifer model setup

T⁎

Giuseppe Brunettia, , Jiří Šimůnekb, David Glöcklerc, Christine StumppaaInstitute for Soil Physics and Rural Water Management, University of Natural Resources and Life Sciences, Vienna 1190, Austria
Department of Environmental Sciences, University of California, Riverside, CA 92521, USAc
Institute of Hydrochemistry, Technical University of Munich, Munich 80333, Germanyb

A R T I C LE I N FO

A B S T R A C T

This manuscript was handled by P. Kitanidis,
Editor-in-Chief, with the assistance of Ming Ye,
Associate Editor

Highly parameterized numerical models are frequently used in environmental sciences to interpret physicochemical processes. Despite their popularity, unjustiﬁed model complexity often undermines model generalizability, resulting in poor model predictive capabilities. A trade-oﬀ between model performance and complexity
can be achieved using model selection techniques, which are largely unexplored in biogeochemical modeling.
Our study aims to address this scientiﬁc gap by investigating the main advantages and ﬁndings of model selection in a typical biogeochemical modeling application. In particular, the study focuses on the numericaldescription, at ﬁve diﬀerent levels of complexity, of the reactive transport processes of nitrogen species in acontrolled aquifer model laboratory setup. To this end, the mechanistic model HYDRUS is coupled with themultimodal Nested Sampling algorithm for the estimation of the Bayesian model evidence and posterior parameter distributions. After successfully calibrating the HYDRUS-2D model against measurements from a tracer
test, a lower-ﬁdelity surrogate-based analysis is carried out to tackle the problem of computational cost, as wellas to compare ﬁve diﬀerent models which combine ﬁrst-order and Monod-based parameterization of the microbial nitrogen degradation. Results demonstrate that measured data contain enough information content tosupport an increase in the complexity of the Monod model, which can reproduce the transient accumulation of
nitrite at the beginning of the experiment. The Bayesian analysis reveals that such behavior is related to theparticular structure of measured data, which plays a fundamental role in the model selection. The analysis
concludes with a validation of the surrogate-based analysis, which conﬁrms the reliability of the proposed
numerical approach, opening up new perspectives on the use of multi-ﬁdelity surrogates for the Bayesian model
selection in environmental modeling.

Keywords:
Modeling
Model selection
Lower-ﬁdelity surrogate model
HYDRUS
Nitriﬁcation

1. Introduction
Water pollution is a global challenge that has increased in bothdeveloped and developing countries, undermining economic growth aswell as the physical and environmental health of billions of people. The
European Environmental Agency (EEA) has recently reported that only
38% of European water bodies met chemical pollution standards(European Environment Agency, 2018), with groundwater sites
showing high exposure to nutrient pollution from agricultural activities.
A similar trend was observed in China and the United States, where
agriculture was identiﬁed as the main source of pollution of surfacewater bodies (Mateo-Sagasta et al., 2017). Extensive use of organic andinorganic fertilizers, as well as waste management from livestock, isdirectly related to increased levels of nitrogen pollution in water bodies.
Preventative measures taken in the last few decades, such as optimized⁎

Corresponding author.

https://doi.org/10.1016/j.jhydrol.2020.124681
Received 15 December 2019; Accepted 11 February 2020
Available online 14 February 20200022-1694/ © 2020 Elsevier B.V. All rights reserved.

fertigation strategies and appropriate tillage and manure management,
have resulted in an appreciable decrease of nitrate concentrations insurface-water bodies, though their eﬀect on groundwater is still limited(EEA, 2018).
Groundwater vulnerability to nitrogen byproducts is mainly inﬂuenced by reactive transport processes in the vadose zone. Subsurface
redox conditions induced by the coupled water-solute transport regulate the microbial transformation of nitrogen species, such as nitriﬁcation and denitriﬁcation. An accurate numerical description ofthese multiple interacting physical and biological processes is crucial
for developing eﬀective and reliable strategies to mitigate groundwater
nitrogen pollution. In the last few decades, several models have beendeveloped and applied to describe and predict nitrogen turnover in soils
both at small and large scales. For example, the hydrological model,
HYDRUS (Šimůnek et al., 2016), has been successfully used to simulateJournal of Hydrology 584 (2020) 124681

G. Brunetti, et al.

contradictory model rankings, which complicate the model selection.
Lu et al. (2011) further conﬁrmed these ﬁndings by highlighting signiﬁcant diﬀerences in the geostatistical model selection between classical information criteria and the Monte Carlo Markov Chain (MCMC)
analysis. Similar results were reported by Schöniger et al. (2014), who
compared nine ways to calculate the Bayesian model evidence in ahydrological application, concluding that numerical integration is at
the basis of the Bayesian model selection.
To this aim, multiple studies have focused on developing or improving algorithms for the MC determination of the marginal likelihood. The most popular MC approximation is the Arithmetic Mean of
the joint likelihood evaluated for parameter samples obtained by various methods (e.g., MCMC, generalized likelihood uncertainty estimation methods). However, this approximation is aﬀected by slow convergence and overestimation of the marginal likelihood, especially if
the posterior distribution is much more concentrated than that of theprior (Kass and Raftery, 1995). To overcome this limitation, Pajor(2017) proposed a Corrected Arithmetic Mean Estimator, in which theprior sample is trimmed to eliminate regions of the parameter spacewith very low likelihood, similarly as in the Corrected Harmonic Mean
estimator (Lenk, 2009). Various other numerical methods have been
proposed in the literature to estimate the Bayesian model evidence
directly or to compute ratios of two marginal likelihoods (i.e., a Bayes
factor). Popular methods are based on the reversible jump algorithm(Green, 1995) and bridge sampling methods (Meng and Wong, 1996).
Recently, Volpi et al. (2017) proposed a Gaussian Mixture Importance(GAME) sampling, which estimates the evidence via multidimensional
numerical integration of the posterior parameter distribution using
bridge sampling. First, a large collection of samples is generated fromthe posterior parameter distribution using MCMC simulation with theDREAM algorithm (Vrugt, 2016) by conditioning separately each candidate model on the observed data. Then, a multivariate mixture distribution is ﬁtted to this MCMC collection of posterior samples using the
likelihood theory to estimate the marginal likelihood of each competing
hypothesis via bridge sampling. A popular alternative is to estimate themarginal likelihood from the Gibbs sampling output (Chib, 1995), later
extended to the Metropolis-Hastings output (Chib and Jeliazkov, 2001).
However, both variants are computationally ineﬃcient. Friel and Pettitt(2008) proposed a thermodynamic integration method, inspired by thepath sampling of Gelman and Meng (1998). Recently, Wei et al. (2019)
applied this method to calculate the marginal likelihood and discriminate between diﬀerent evapotranspiration models. In anotherstudy, Liu et al. (2016) showed that thermodynamic integration could
outperform semi-analytical solutions, Arithmetic Mean, Harmonic
Mean, and nested sampling, in terms of the solution accuracy andconsistency for repeated independent runs. A diﬃculty with thismethod is the choice of the optimal temperature schedule. To circumvent this issue, Weinberg (2012) proposed a novel approach based on
Lebesgue integration. However, the thermodynamic integration generally requires a high computational cost, which restraints its use inenvironmental modeling. Elshall and Ye (2019) proposed a variant of
the steppingstone sampling algorithm, which is less sensitive to sampling errors and more computationally eﬃcient than the thermodynamic integration, yet this method is mainly limited to moderately
complex problems.
From this literature background, it emerges that an accurate estimation of the marginal likelihood requires a high number of model
executions, which may not be possible for biogeochemical models dueto their high computational cost. A viable solution is to use highly efﬁcient MC algorithms, which provide a balanced trade-oﬀ between
marginal likelihood accuracy and computational eﬀort. This will allow
to have a statistical basis for biogeochemical model selection and uncertainty assessment, which is of great importance in biogeochemical
modeling. Skilling (2006) proposed a one-dimensional re-parametrization of the integral of the Bayesian model evidence, which signiﬁcantly
reduced the computational cost of the Bayesian analysis. This approach

nitrogen balance in several agricultural applications (e.g., Filipović
et al., 2016; Gärdenäs et al., 2005; Hanson et al., 2006; Li et al., 2015).
HYDRUS numerically solves the Richards and advection-dispersion-reaction equations describing water ﬂow and solute transport in soils,
respectively. Nitrogen reactions are generally simulated using a sequential chain with ﬁrst-order coeﬃcients, which reduces the numberof model parameters but simpliﬁes the description of the microbial
degradation. On the other hand, the model DAISY (Hansen et al., 1991)
reﬁnes the description of nitrogen turnover in the root zone by simulating the eﬀect of microbial growth in multiple conceptual pools. A
similar approach is used in the models RZWQM (Ma et al., 2012) and
CENTURY (Parton, 1996). Many studies have proven the accuracy ofthis approach (e.g., Cameira et al., 1998; Kröbel et al., 2010; Ma et al.,2007; Salazar et al., 2017).
However, the choice of an appropriate level of model complexity isa constant problem in biogeochemical modeling (e.g., Homann et al.,
2000; Kimmins et al., 2008). A large number of optimized parameters
can pose model calibration issues, such as signiﬁcant parameter uncertainty and limited model-predictive capabilities that undermine its
generalizability (Ajami and Gu, 2010; Elsheikh et al., 2013; Freni et al.,2011; Schoups et al., 2008). Since model uncertainty plays a crucialrole in optimal decision making for groundwater management (PeñaHaro et al., 2011; Wagner and Gorelick, 1987), the selected model must
be as simple as possible, as complex as necessary (Höge et al., 2018) toaccurately reproduce observations and reduce predictive uncertainty.
This approach should be ubiquitous in any modeling application. The
high conceptual parameterization of complex models often results inthe overﬁtting of measured data and parameter non-uniqueness, unlesslarge and informative observations are used for calibration (Pandeet al., 2015). In contrast, simpler models generally have lower predictive uncertainty since their few parameters are easily constrained byobservations, although the ﬁt quality deteriorates. Thus, a less accurate,
simple model ranks better than a more complex one that does notprovide a good enough ﬁt to justify its higher complexity (Myung,2000). This trade-oﬀ between quality of ﬁt and complexity embodies
the principle of Occam's razor, which states that “entities should not be
multiplied beyond necessity,” and is at the basis of model selection
techniques.
The Bayesian model selection is a formal approach that allows acomparison of alternative models on a statistically rigorous basis inlight of observations and any prior information available. This evaluation is based on the calculation of the Bayesian model evidence,
often named in the literature as the marginal likelihood. Computation of
the Bayesian model evidence is numerically challenging, as it involves a
multi-dimensional integration over the entire parameter space often
described by a multimodal peaked likelihood surface. Such problems
are exacerbated for environmental and biogeochemical models characterized by high nonlinearity and strong parameter interactions(Ajami and Gu, 2010; Brunetti et al., 2019, 2016; Dai et al., 2019). To
avoid the time-consuming Monte Carlo (MC) analysis, several authors
have suggested diﬀerent analytical approximations of the Bayesianmodel evidence. For example, the Akaike Information Criterion(Akaike, 1973) stems from the information theory and is frequently
applied for its ease of implementation (Guthery et al., 2003). In anotherstudy, Neuman (2003) proposed a Maximum Likelihood BayesianModel Averaging approach based on Kashyap's information criterion,
which reduces computational eﬀort. A further simpliﬁed index is theBayesian Information Criterion ﬁrst introduced by Schwarz (1978).
Recently, Watanabe (2013) proposed a Widely Applicable BayesianInformation Criterion (WBIC), which starts with the computed log
pointwise posterior predictive density and then adds a correction for aneﬀective number of parameters to adjust for overﬁtting. As discussed byGelman et al. (2014), this criterion has the advantage of averaging overthe posterior distribution rather than conditioning on a point estimate.
However, Friel et al., (2017) demonstrated that the WBIC can overestimate the Bayesian model evidence. Such inaccuracies can lead to2

Journal of Hydrology 584 (2020) 124681

G. Brunetti, et al.

2014; Laloy et al., 2013; Razavi et al., 2012). While most of the existing
literature focuses on data-driven surrogates, previous studies highlighted how lower-ﬁdelity surrogates overcome many of the limitationsof data-driven surrogates (Asher et al., 2015; Razavi et al., 2012).
Lower-ﬁdelity surrogates are computationally cheap models constructedfrom the complex model by reducing geometrical complexity, numerical resolution, or removing processes, and generally, yield eﬃcientcalibration providing that their approximation is limited and assessedagainst the original model. The combination of eﬃcient Bayesian model
selection techniques (i.e., Nested Sampling) with lower-ﬁdelity surrogates is still largely unexplored for biogeochemical modeling applications. More research is thus needed in this direction. Obtained results
could apply to other scientiﬁc ﬁelds as well.
Our study aims to cover these scientiﬁc gaps by exploring the main
advantages and ﬁndings of model selection in a typical biogeochemical
application. In particular, the study focuses on the numerical description, at diﬀerent levels of complexity, of the reactive transport processes of nitrogen species in a controlled aquifer model laboratory
setup. The aquifer model is ﬁrst subject to a tracer test to characterize
its hydraulic functioning and then to an application of cattle slurry toassess biogeochemical reactions involved in nitrogen transformations.
The problem is addressed in the following way. Initially, the mechanistic model HYDRUS (Šimůnek et al., 2016) is externally coupled with
the MULTINEST algorithm for the estimation of the Bayesian model
evidence and posterior parameter distributions. Next, the coupled twodimensional HYDRUS-2D model is employed to inversely estimate soil
hydraulic and solute transport parameters using measurements from a
tracer experiment. Estimated parameters are then used in the one-dimensional HYDRUS-1D model, which simulates transport and biogeochemical reactions in the system during the application of cattle slurry.
HYDRUS-1D is used as a lower-ﬁdelity surrogate to overcome the impractically high computational cost associated with the 2D model in thecomputationally intensive Bayesian model selection. Five levels ofcomplexity of the biogeochemical model are investigated. Either constant ﬁrst-order decay rates or Monod kinetics are alternatively usedand combined to interpret ammoniﬁcation and nitriﬁcation processes in
the aquifer. High spatial resolution measurements are used in theBayesian analysis to constrain the inverse problem. Finally, results fromthe surrogate-based analysis are projected in the 2D model for validation. It must be noted that the combined use of the MULTINEST algorithm and the lower ﬁdelity surrogate for the Bayesian model selection
represents a novelty in biogeochemical modeling, underlining that thedevelopment of eﬃcient model selection frameworks for such applications is of crucial importance.

is at the basis of the Nested Sampling (NS) algorithm, which has been
applied successfully in several research ﬁelds (Elsheikh et al., 2014,2013; Mukherjee et al., 2006; Paulson et al., 2019; Schöniger et al.,2014; Shaw et al., 2007). Liu et al. (2016) evaluated the convergence
and accuracy of NS using multiple numerical experiments and concluded that it is computationally more eﬃcient than thermodynamic
integration and similarly accurate. The low stability in the estimation of
the marginal likelihood stems from the necessity to draw eﬃciently
unbiased samples from the likelihood-constrained prior. Several variants of the original algorithm were developed to address this issue.
Brewer et al. (2011) proposed a Diﬀusive Nested Sampling algorithm
that uses MCMC to explore a mixture of progressively constrained
distributions. However, for highly correlated and multimodal target
distributions, multiple particles needs to be considered, which increases
the computational cost. Furthermore, multiple tuning parameters haveto be set, making the algorithm very sensitive. Handley et al. (2015)
developed the PolyChord NS algorithm, which uses the slice sampling
to sample within the nested isolikelihood contours, thus making itsuited for a high-dimensional parameter space. A recent study of Cao
et al. (2018) reported how PolyChord outperformed other NS-basedalgorithms in terms of accuracy and stability in high-dimensional inverse problems, even though it is computationally expensive. In anotherstudy, Zeng et al. (2018) used the DREAMZS algorithm to sample fromthe likelihood-constrained prior and compared it against the classic NSMetropolis-Hastings. Results of that study indicate that the NSDREAMZS is more eﬃcient and robust, which indicates that this algorithm is due to its higher eﬃciency preferable to the thermodynamic
integration when computational resources are inadequate. To improvethe computational eﬃciency of the NS, Higson et al. (2019) introduceda Dynamic Nested Sampling algorithm, in which the number of “livepoints” varies, allocating samples more eﬃciently. This results in a
speedup of computations by a factor of up to ~7 for marginal likelihood
calculations. Feroz et al. (2009) proposed an ellipsoidal rejection
sampling scheme (MULTINEST) that allows us to identify and isolate
near-Gaussian posterior modes in small overlapping ellipsoids. Once
identiﬁed, these modes can be evolved independently, thereby automatically accommodating posteriors degeneracies while maintaining
high eﬃciency for simpler problems. MULTINEST has proven to behighly eﬃcient and accurate for low dimensional problems (i.e., < 20parameters), while its accuracy degrades for highly parameterized
analyses. Later, Feroz et al. (2013) proposed an Importance Nested
Sampling, which overcomes such limitations and signiﬁcantly increases
the accuracy of MULTINEST.
Despite its computational eﬃciency, the NS still requires thousandsof model runs to calculate the Bayesian model evidence. In this perspective, recent studies on the surrogate-based analysis have demonstrated their eﬃciency for computationally expensive model calibration
applications (Asher et al., 2015; Brunetti et al., 2017; Elsheikh et al.,Fig. 1. A cross-section of the experimental facility with the spatial distribution of measurement points.
3

Journal of Hydrology 584 (2020) 124681

G. Brunetti, et al.

parameters. The tracer test was conducted under steady-state ﬂow
conditions by injecting deuterium (2H) enriched groundwater into thesystem. In particular, a tracer solution volume of 1.6 L with a concentration of 243 mg/L was applied at the soil surface for 173 min.
BTCs were obtained by regularly sampling pore water for ﬁve days at
16 selected sampling ports (Fig. 1). In total, approximately 40 mL ofpore water was extracted per sampling port over a period of up to fourdays, which is small compared to the total water volume (30 L), and not
signiﬁcant enough to inﬂuence water ﬂow within the system.

Table 1
Measured irrigation rates, nitrogen species concentrations in the inﬂow, andwater table depths during the cattle slurry injection.
Time
(day)

Irrigation
rate (L/day)

Water Table
depth (cm)

Norg
(mg/
L)

NH4+
(mg/L)

NO2−(mg/L)

NO3−(mg/L)

9
23
51
54
74
89

13.35
13.35
6.50
6.50
6.50
0.00

25.00
25.00
25.00
25.00
[16.00,34.00]
25.00

10.92
11.89
12.15
1.00
1.00
1.00

10.92
11.89
12.15
0.03
0.03
0.03

0.21
0.00
1.26
4.29
4.29
4.29

2.24
0.07
0.01
51.46
51.46
51.46

2.1.3. Experiments on the fate of ammonium and nitrate
Nitrogen mineralization was triggered by injecting cattle slurry into
the aquifer model for 89 days. The main aim of this experiment was toinvestigate nitrogen reactive processes in a combined saturated-unsaturated system under diﬀerent operating conditions.
The experiment required four diﬀerent reservoirs with a volume of
200 L. Two liters of preconditioned liquid manure were diluted in 198 Lof pure groundwater to reach a suﬃciently high concentration of nutrients. Biogeochemical measurements were carried out in both thereservoirs and at the injection needles when connecting a new reservoir
to the system. Because, as mentioned previously, it was not possible tomeasure organic nitrogen, its concentration was established using
measured NH4+ concentrations and compared with literature values. In
particular, previous studies indicated that ammonium generally represents about 50% of the nitrogen in the liquid cattle manure (e.g.,
Beauchamp et al., 2010; van Middelkoop and Holshof, 2017). Therefore, the organic nitrogen concentration in the inﬂow is assumed to be
the same as the ammonium during the ﬁrst 51 days of the experiments,
considering that nitrite and nitrate concentrations were considerably
smaller (Table 1). The last cattle manure reservoir was connected to theinjection system on day 51. Measured inﬂow concentrations indicatethat NO3− was the predominant N species (Table 1). This behavior isdue to strong mineralization processes likely triggered by preliminary
sieving and centrifugation of cattle slurry used to separate liquid andsolid fractions. Indeed, previous studies have demonstrated that suchtreatment results in a monodisperse particle size distribution of the liquid manure characterized by a majority of particles below 25 µm(Peters et al., 2011). Such particles are characterized by low carbon/
nitrogen ratios and a high nitrogen mineralization potential (Fangueiro
et al., 2008). Thus, a Norg concentration of 1 mg/L was assumed for thelast reservoir, considering that DOC was still detected in the inﬂow.
Nevertheless, preliminary numerical simulations showed that this value
plays a negligible role in the nitriﬁcation process due to its relatively
low magnitude.
The experiment was divided into three stages (Table 1). Steady-state
ﬂow conditions were induced for the ﬁrst 22 days by setting the recharge rate and the water table to 13.35 L/day and 25 cm, respectively.
On day 23, the recharge rate was decreased to 6.5 L/day to increasetransit times. Between days 54 and 73, transient ﬂow conditions weregenerated by simulating diurnal water table ﬂuctuations. The applied
amplitude and periodicity were consistent with water table ﬂuctuations
observed in natural environments (Gribovszki et al., 2010, 2008). In
particular, the automated piston device was programmed to raise andlower the outﬂow tube from 16 to 34 cm, providing a complete ﬂuctuation cycle over 24 h. On day 73, the ﬂuctuations were terminatedwhen the water table reached its starting position at a depth of 25 cm.
At that time, the peristaltic pumps were deactivated, resulting in stagnant water ﬂow. Measured dissolved oxygen in the aquifer model wasalways above 2 mg/L, which indicates suitable aerobic conditions fornitriﬁcation and negligible denitriﬁcation.

2. Materials and methods
2.1. Case study description
2.1.1. Experimental setup
A schematic of the experimental aquifer setup is shown in Fig. 1.
The aquifer model (106 cm long * 50 cm high * 20 cm wide) was divided into two chambers: a 100 cm wide chamber for the sediment
packing and a 6 cm wide chamber for drainage, separated by a waterpermeable grid. The soil was packed up to a height of 50 cm. To
guarantee uniform irrigation, water inﬂow (at the soil surface) was
pumped from a reservoir to the aquifer model using 4 peristaltic pumps
connected to 60 silicone tubes equipped with an injection needle. The
outﬂow tube from the drainage chamber was used to simulate watertable ﬂuctuations by adjusting its height via a self-built automatedpiston device. Outﬂow water was then collected in a waste containerand pumped to the sink.
The sediment used in the experiment originated from a gravel pit in
Bruckmühl near Munich, Germany. The particle size distribution analysis indicated a medium-coarse sandy sediment characterized by anegligible silt and clay content (i.e., 1.25%). Groundwater used for theexperiments was taken directly from a well at the Institute for
Groundwater Ecology in the Helmholtz Zentrum München, Germany,
where all experiments were conducted.
Soil water samplers (Rhizon SMS, Eijkelkamp) and optical dissolved
oxygen sensors (PSt3, Presenns) were installed in four vertical proﬁles(indicated by letters A, B, C, and D; vertical dashed lines in Fig. 1) witheight ports and spots each, respectively. Pore water samples extracted
during the tracer test were analyzed using an isotopic water analyzer(Picarro L2120-i) to reconstruct breakthrough curves (BTCs). Ammonium (NH4+), nitrite (NO2−), and nitrate (NO3−) concentrations insamples and groundwater were determined using an ion chromatography system (ICS 1100, Dionex), while a total organic carbon analyzer (TOC-5000A, Shimadzu) was used to measure the Dissolved Organic Carbon (DOC). No measurements of organic nitrogen (Norg) wereavailable due to operational issues. The soil pressure head was monitored by eleven tensiometers (black circles in Fig. 1) (T5-5, UMS)
distributed in three vertical proﬁles with an acquisition frequency of
30 min.
The experiment consisted of a tracer test, which was conducted
before the manure injection. The system was equilibrated by irrigatingwith pure groundwater for 15 weeks before starting with the manure
injection. In particular, during this period, the irrigation rate and thewater table were set to 13.35 L/day and 25 cm, respectively.

2.1.2. Tracer experiments
The tracer experiment took place between weeks 6 and 15 of theequilibration period in order to avoid disturbances inside the model
aquifer during experiments with cattle slurry (i.e., nitrogen). The main
objective of the tracer test was to accurately characterize the solutetransport processes in the aquifer and to use BTC data in the Bayesianframework to inversely estimate soil hydraulic and transport

2.2. Modeling theory
2.2.1. Water ﬂow and solute transport
The hydrological model HYDRUS (Šimůnek et al., 2016) is used tomodel the hydraulic behavior of the aquifer. HYDRUS is a ﬁnite4

Journal of Hydrology 584 (2020) 124681

G. Brunetti, et al.

bacteria, Ammonium Oxidizing Bacteria (AOB) and Nitrite Oxidizing
Bacteria (NOB). Such reactions can be numerically described using two
diﬀerent formulations: either as constant ﬁrst-order decay reactions or
as Monod Kinetics. The former neglects the bacteria growth and assumes constant reaction rates, while the latter adjusts the reaction depending on the simulated bacteria population. When considering constant ﬁrst-order decay, the reaction chain in the solid and liquid phases
is described as follow:

element model for simulating the movement of water, heat, and multiple solutes in variably-saturated porous media. It numerically solves
the Richards equation for multi-dimensional variably-saturated waterﬂow:∂θ= ∇ [K ·∇ (h − z )]∂t(1)

where t is time [T], z is the vertical coordinate [L], θ is the volumetric water content [L3L−3], K is the unsaturated hydraulic conductivity [LT−1], and h is the pressure head [L]. The unimodal van
Genuchten–Mualem (VGM) model (van Genuchten, 1980) is used todescribe the soil hydraulic properties:

ϕk = μL, k θCk + μS, k ρskwhere µL,k and µS,k are the ﬁrst-order degradation coeﬃcients for thekth solute, providing a connection between chain species in the liquidand solid phases, respectively [T−1]. This numerical approach has beenused extensively in the literature (e.g., Hanson et al., 2006; Li et al.,2015; Ramos et al., 2012). The use of this model is also justiﬁed bymeasured inﬂow concentrations (Table 1), which exhibit a complete
shift towards nitrate after day 54, suggesting a plausible negligible effect of reaction processes towards the end of the experiment.
On the other hand, biomass growth and decay are explicitly modeled using the Monod kinetics and the ﬁrst-order reaction, respectively:1Θ=⎧ (1 + (α | h |)n)m if h ⩽ 0⎨
1 if h > 0⎩θ−θΘ = θ − θr
s

(2)

r

( (

m 2

) ) ⎤⎦ if h<0⎧ Ks Θl ⎡ 1 − 1 − Θ m1
K=⎣⎨
Ks if h>0⎩1
m=1− ndΓCk=γ
Γ − bΓ
dt
S + Ck(3)
3 −3

where Θ is the eﬀective saturation [L L ], α is a shape parameter
related to the inverse of the air-entry pressure head [L−1], θs and θr are
the saturated and residual water contents, respectively [L3L−3], n andm are pore-size distribution indices [–], Ks is the saturated hydraulicconductivity [LT−1], and l is the tortuosity and pore-connectivity
parameter [–].
The reactive transport of nitrogen species in the aquifer is modeled
as a sequential reaction chain, in which organic nitrogen is ﬁrst converted into ammonium by heterotrophic organisms, then oxidized by
autotrophic bacteria into nitrite and nitrate. Denitriﬁcation processes
are neglected due to the persistent oxic conditions observed during theexperiment. The transport of the kth solute in the aquifer is describedusing the advection–dispersion-reaction equation, assuming that solutes
can exist only in the solid and liquid phases:

∂ρsk
∂θCk
= ∇ (θDkW ·∇Ck ) − ∇ (qCk ) − ϕ+∂t∂twhere Γ is the biomass concentration [ML ], γ is the maximum
aerobic growth rate [T−1], S is the half-saturation constant [MM−1],
and b is the ﬁrst-order coeﬃcient representing lysis [T−1]. This approach assumes that bacterial degradation only occurs in the liquid
phase (De Wilde et al., 2009) and that bacteria are immobile (De Wildeet al., 2009; Holden and Fierer, 2005; Wehrer et al., 2012), leading to:

ϕk =

Γ⎛
Ck ⎞
θ
γ
Y ⎝ S + Ck ⎠
⎜

⎟

(8)

where Y is the yield coeﬃcient [MM−1]. The Monod kinetics, whileproviding a more coherent description of nitrogen turnover, requires
more parameters and is computationally more expensive, and therefore
its use should be statistically justiﬁed in light of the information content
of measured data.
A preliminary analysis of measured eﬄuent concentrations showeda transient accumulation of nitrites and a delayed production of nitrates
in the aquifer at the beginning of the experiment, suggesting the bacterial growth had a signiﬁcant inﬂuence on the mineralization process(Fig. 2). Ammonium in injected water was readily available for oxidation into nitrites by AOBs. As the AOBs population grew, the nitrite
substrate became suﬃcient to trigger NOBs growth, which reduced

where Ck is the concentration of the kth solute in the liquid phase[ML−3], sk is the concentration of the kth solute in the solid phase[MM−1], ρ is the sediment density [ML−3], DkW is the dispersion tensor
for the kth solute in water [L2T−1], q is the water ﬂux [LT−1], and ϕkrepresents the reaction sink/source term [ML−3T−1]. No reactive processes are assumed for deuterium and nitrate.
Linear adsorption to the solid phase is considered for Norg and NH4+
(Gärdenäs et al., 2005; Li et al., 2015; Mariano et al., 2016), while 2H,NO2−, and NO3− exist only in the liquid phase:
(5)
3

−1

where Kd,k is the distribution coeﬃcient of the kth solute [L M
components of the dispersion tensor Dij,kW are expressed as:

θDijW, k = λT |q| δij + (λL − λT )

(7)
−3

(4)

Sk = K d, k Ck

(7)

qj qi
|q|

+ θDW , k τW δij]. The(6)where λT and λL are the transverse and longitudinal dispersivities [L],
δij is the Kronecker delta function [–], DW,k is the molecular diﬀusion
coeﬃcient of solute k in free water [L2T−1], and τW is the tortuosity
factor in the liquid phase [–] (e.g., Millington and Quirk, 1960).
2.2.2. Microbial degradation of nitrogen
Heterotrophic and autotrophic bacteria in the porous media are
responsible for nitrogen mineralization. Heterotrophic bacteria convert
organic nitrogen into ammonium (i.e., ammoniﬁcation), which is used
as a substrate in the next nitriﬁcation process operated by autotrophic

Fig. 2. Measured concentrations of NH4+, NO2−, and NO3− in outﬂow water
during the experiment showing the transient accumulation of NO2−.
5

Journal of Hydrology 584 (2020) 124681

G. Brunetti, et al.

neglected due to its limited eﬀect on the experiment (cover on top). A
special “Time-Variable Pressure Head” BC (blue line in Fig. 3) is used
for the right part of the domain exposed to air and water table ﬂuctuations. This BC assigns the speciﬁed pressure head value to the lowest
nodal point of the boundary, while pressure heads at other nodes areadjusted based on the z-coordinate. At the same time, a “Seepage Face”
BC is assigned to nodes with calculated negative pressure heads. The
concentration ﬂux along the top and right boundaries is simulated using
the Cauchy-type BC. Zero water and concentration ﬂux boundary
conditions (black lines in Fig. 3) are applied to all remaining boundariesof the domain.
Preliminary tests have shown that the execution time to simulate
both ammoniﬁcation and nitriﬁcation using the Monod kinetics in
HYDRUS-2D was around 10,000 s on a PC equipped with a CPU Intel
Core i7-8750H, 16 GB of RAM, and SSD. Since no high-performancecomputing cluster was readily available, such computational cost would
make the Bayesian model selection impractical. Therefore, HYDRUS-1D
is used as a lower-ﬁdelity surrogate model due to its limited run time(i.e., ≈8 s). This assumption appears acceptable for several reasons: 1)
HYDRUS-1D simpliﬁes only the geometry of the problem, while it
shares the same parameterization and simulates the same physicochemical processes as HYDRUS-2D, 2) the inﬁltration process in theunsaturated part of the domain is considered to be strictly vertical, andthus, the two models are expected to produce similar inﬁltration
pathways, 3) since the majority of reactive processes take place near the
sediment surface, the similarity between HYDRUS-1D and HYDRUS-2Din this zone guarantees a good description of the reactive transport also
in the surrogate, and 4) preliminary simulations showed that the velocities are higher in the saturated zone of the system than in the unsaturated zone. Due to the limited length of the aquifer model (i.e.,
100 cm), the delay induced by horizontal ﬂow in the saturated zone israther limited, and the assumption of instantaneous discharge in thesurrogate is thus tolerable. However, it must be emphasized that theseassumptions and the used approach are not general but restricted to thisparticular case study. There exist multiple numerical approaches toformally address the discrepancies between low and high-ﬁdelity
models, such as correction functions, space-mapping, and hybrid strategies (Razavi et al., 2012), which are, however, beyond the scope of
the present work. Nevertheless, ﬁnal model validation is carried out toassess the results of the surrogate-based analysis.
Similarly, the HYDRUS-1D vertical domain is discretized in 100 onedimensional elements which are reﬁned near the soil surface. An“Atmospheric” BC is speciﬁed on the top boundary, while a “TimeVariable Pressure Head” BC (the red line in Fig. 3) is used for thebottom boundary. The quality of the FE mesh is assessed by checking
the mass balance error reported by HYDRUS at the end of the simulation. Mass balance errors, which in this simulation were always below
1%, are generally considered acceptable at these low levels.
Initial pressure head conditions are assumed to be in hydrostatic
equilibrium from the lowest located nodal point (i.e., z = 0 and
h = 25 cm). Measured concentrations from diﬀerent sampling points

Table 2Models analyzed in the Bayesian model selection framework and their description of biogeochemical processes.
Biogeochemical
Process

Ammoniﬁcation
Ammonium Oxidation
Nitrite Oxidation

Model ID
M1

M2

M3

M4

M5

FirstorderFirstorderFirstorderFirst-orderFirst-orderFirst-order

Monod

Monod

First-order

Monod

Monod

First-order

Monod

Monod

Monod

nitrites by producing nitrates. After this transient period, both bacterial
populations were equilibrated, and the mineralization process wascharacterized by fast nitrate production and transport. This preliminary
data screening theoretically supports the use of Monod kinetics for theinterpretation of the nitrogen turnover in the aquifer. However, themost appropriate description of each process (e.g., ammoniﬁcation,ammonium oxidation, nitrite oxidation) must be identiﬁed to avoid
unjustiﬁed complexity. To this aim, the ﬁrst-order decay reactions andMonod kinetics are combined to obtain a set of ﬁve diﬀerent biogeochemical models of increasing complexity (Table 2). In particular, thesimplest M1 model describes all nitrogen reactions using only ﬁrstorder reactions, while the most complex M5 model simulates bacteria
growth in all reaction processes. Diﬀerent model compositions aim not
only to identify the most appropriate level of model complexity but also
to explain whether bacterial growth plays a crucial role in all reaction
steps.

2.2.3. Numerical domain and boundary conditions
In this section, model settings used in the following Bayesian analysis are described. The two-dimensional model HYDRUS-2D is used in
combination with tracer data for the Bayesian estimation of soil hydraulic and transport properties. Throughout this manuscript, HYDRUS2D is considered the high-ﬁdelity model due to its comprehensive description of Boundary Conditions (BCs). The two-dimensional domain
has a length and a depth of 100 and 50 cm, respectively. The domain isdiscretized into two-dimensional triangular elements reﬁned in theright part of the domain to accommodate velocity gradients at theoutﬂow, thus reducing mass balance errors. The generated FE mesh has
2142 nodes and 4100 two-dimensional elements.
An “Atmospheric” BC (green line in Fig. 3) is speciﬁed on the top
boundary of the model domain to simulate the irrigation of the aquifer.
This BC type can exist in three diﬀerent states: (a) irrigation and/or
potential evaporation ﬂuxes, (b) a zero pressure head (full saturation)
during ponding when both inﬁltration and surface runoﬀ occur, and (c)an equilibrium between the soil surface pressure and the atmosphericwater vapor pressure when the atmospheric evaporative demandcannot be met by evaporation from the soil proﬁle. Soil evaporation isFig. 3. The spatial distribution of applied boundary conditions in both HYDRUS-2D and HYDRUS-1D.
6

Journal of Hydrology 584 (2020) 124681

G. Brunetti, et al.

likelihood as the prior volume is shrunk. Once the evidence Z is estimated, posterior parameter distributions and associated weights can beinferred using the full sequence of discarded points from the nested
sampling process. In this study, the active set consisted of 100 points for
all Bayesian analyses. This value was selected based on a preliminary
sensitivity analysis on the eﬀect of active points on the calculated
marginal likelihood. Brieﬂy, the Coeﬃcient of Variation (CV) of themarginal likelihood was estimated by running the NS algorithm ﬁvetimes using a diﬀerent number of active points (i.e., 20, 50, 100, 200),
as described in Zeng et al. (2018). Only a negligible reduction in the CV
was obtained when using more than 100 points. This analysis was
carried out used to assess the stability of the estimator and its potential
eﬀect on models’ ranking.
The most challenging task in the NS algorithm is sampling from theprior with the hard constraint L > L0. To overcome this problem,
several techniques have been proposed in the past. In the present study,
the MULTINEST algorithm proposed by Feroz et al. (2009) is used.
MULTINEST is a combination of the ellipsoidal nested sampling(Mukherjee et al., 2006) with the recursive clustering (Shaw et al.,2007) that is particularly suited for degenerated multimodal posterior
distributions, which are often encountered when inversely estimating
soil hydraulic and transport parameters (e.g., Brunetti et al., 2019).
Indeed, in such circumstances, a single ellipsoid also includes lowlikelihood regions between high-likelihood clusters making the algorithm extremely ineﬃcient. The solution is to detect these clusters andbind them in separate ellipsoids. To achieve this, a recursive processwith 2-means clustering is adopted. If the resulting two ellipsoids havea signiﬁcantly lower total volume than the parent ellipsoid (less thanhalf), the split is accepted and the clustering and volume test are repeated on each of the two subsets of points. This process continues
recursively. Alternatively, if the total ellipse volume is signiﬁcantly
greater than expected, this indicates that there may be more than two
clusters and that 2-means clustering was not appropriate. Therefore, thealgorithm still tries to subdivide the clusters recursively. However, it
only accepts the ﬁnal clustering if the total volume decrease is signiﬁcant.
It can be assumed that the posterior bulk was explored suﬃciently
when the active points, in sum, incremented the evidence by only somesmall fraction of the total accumulated evidence. More speciﬁcally, thecurrently active point with the highest likelihood gives an estimate of
the highest likelihood in the region. The remaining evidence at iteration
i is estimated according to:were used to set initial solute concentrations in the solute domain. Since
no measured data were available, the microbial biomass initial concentration in model M2 was assumed to decrease with depth. In particular, a value of 10 mg/L is assumed for the top 10 cm, while 1 mg/L
is set in all remaining nodes of the numerical domain. This is compatible with what was already observed in other studies (e.g., Henrichs
et al., 2007).
The standard HYDRUS code was modiﬁed to simulate bacterial
degradation in model M2. To this end, a separate FORTRAN ﬁle wascreated to solve explicitly in time Eq. (7) and to calculate zero-order
reaction rates, which are then used to assemble the system of algebraic
equations resulting from the discretization of the governing advectiondispersion partial diﬀerential equations.
2.3. Model calibration and uncertainty assessment
2.3.1. MULTINEST algorithm
In this section, the MULTINEST algorithm used in the model selection and inverse parameter estimation is described. For simplicity,
only a general overview of the method is provided, while detailed descriptions about the theoretical background and numerical implementation can be found in Skilling (2006), Mukherjee et al. (2006),
and Feroz et al. (2009).
The Bayesian inference provides a statistically coherent approachfor the estimation of a set of parameters Ω in a model H for the data D:Pr(Ω|D,H) =Pr(D|Ω,H) Pr(Ω|H)
Pr(D|H)(9)

where Pr(Ω|D, H) ≡ P(Ω) is the posterior probability of the parameters, Pr(D|Ω, H) ≡ L(Ω) is the likelihood, Pr(Ω|H) ≡ π(Ω) is theprior, and Pr(D|H) ≡ Z is the Bayesian model evidence, which can be
expressed as:Z=∫ L (Ω) π (Ω) dDΩ

(10)where D is the dimensionality of the parameter space. The Bayesianevidence is the average of the likelihood over the prior, and it is expected to be smaller for a model with large areas in its parameter spacehaving low likelihood, even if the likelihood function is highly peaked,
and larger for a model with higher likelihood of model parameters in alarger fraction of its parameter space. The numerical estimation of Z is
at the basis of the model selection and is the core of the NS algorithm.
NS was ﬁrst presented by Skilling (2006) as an eﬃcient Monte Carlo
method to calculate the evidence and estimate posterior distributions as
by-products. The core of the method is the transformation of the multidimensional integral in Eq. (10) to a one-dimensional integral. To thisend, the prior mass X can be deﬁned as the cumulative prior mass
covering all likelihood values greater than the iso-likelihood contour ω:X (ω) =∫L (Ω) > ω π (Ω) dDΩ

where Xi is the remaining prior volume. The stopping criterion is thenbased on the ratio between the estimated total evidence and the current
evidence:log(Zi + Zest ) − log Zest < 0.5(11)2.3.2. Likelihood function and prior distributions
If we assume that error residuals are uncorrelated and normally
distributed with a constant variance, σ2, the likelihood function L(Ω)
can be written as:

1∫ L (X ) dX0(14)

To avoid the premature stopping of the algorithm, the maximum
number of algorithm iterations is set to 100,000. Such threshold value
has never been reached during the simulations.

For the case of normalized prior probability, the prior mass X has avalue 1 if ω = 0 and decreases to 0 as ω approaches 1. Deﬁning theinverse function L(X) as the likelihood which bounds a prior mass X ,one can simplify Eq. (10) from a multi-dimensional integral over theprior volume to a one-dimensional integral over the prior mass:

Z=

(13)

Zest = L max Xi(12)

The NS algorithm provides an iterative process to solve Eq. (12). At
each iteration i, N “active” samples are drawn from the prior. Samples
are then sorted in order of their likelihood, and the smallest (withlikelihood L0) is removed from the active set and replaced by a point
drawn from the prior subject to the constraint that the point has alikelihood L > L0. The iterative process is stopped when the stopping
criterion is satisﬁed. The algorithm thus travels through nested shells ofk

L (Ω) = ∏
i=1

1
2πσ 2

1
exp ⎡− σ 2 (Hi (Ω) − y~i )2⎤
⎦
⎣ 2

(15)

where Hi(Ω) and ỹi are the ith model realization and its corresponding
measured value, respectively. For algebraic simplicity, the log-likelihood ℓ(Ω) is used:7

Journal of Hydrology 584 (2020) 124681

G. Brunetti, et al.
kℓ(Ω) = −kk
1ln(2π ) − ln(σ 2) −∑ (Hi (Ω) − y~i )22
2
2σ 2 i = 1

Table 3The bounds of the soil hydraulic, transport, and reactive parameters used in theBayesian analyses.

(16)

In the inverse modeling, σ is treated as a nuisance parameter andintegrated out from the inference equation (e.g., Kavetski et al., 2006).
The new likelihood is then written as:2kℓ(Ω) = −k ⎡⎤ln ∑ (Hi (Ω) − y~i )2⎥2 ⎢ i=1⎣⎦(17)

Uniform prior distributions are used in the Bayesian analysis.
Parameter bounds are provided in the following sections.
2.3.3. Inverse estimation of soil hydraulic and solute transport parameters
Before the Bayesian model selection, the MULTINEST algorithm isused in conjunction with tracer data to inversely estimate soil hydraulicand transport parameters. This modeling scenario is denoted as TR. A
high-resolution dataset consisting of measured BTCs from 16 sampling
points is used to calibrate HYDRUS-2D. Since the tracer breakthroughcurves did not show any anomalous behavior, such as tailing (Knorret al., 2016; Pedretti et al., 2013), the eﬀects of soil heterogeneity havebeen neglected. To have a better agreement between HYDRUS-1D and
HYDRUS-2D, a weight of 0.2 has been assigned to measured deuterium
concentrations in the saturated zone where the eﬀect of the horizontal
ﬂow is appreciable. As such, the optimized soil hydraulic parameters
are mostly representative of the unsaturated zone, which is similarly
described by the high and lower ﬁdelity models. For such conditions,estimated soil hydraulic and transport parameters can be directlytransferred to HYDRUS-1D since they are representative of the soil
properties of the aquifer. The residual water content θr is ﬁxed to 0.0,
which is a common value for coarse-textured soils (e.g., Brunetti et al.,2017). The pore connectivity parameter l is set to 0.5, and the transverse dispersivity λT is assumed 0.1λL. Due to these simpliﬁcations, thenumber of unknown parameters is reduced to 5: the saturated watercontent θs, the VGM shape parameters α and n, the saturated hydraulicconductivity Ks, and the longitudinal dispersivity λL. The parameterbounds used in the analysis are reported in Table 2.

Parameter

Scenario ID

Lower Bound

Upper Boundθs (cm3/cm3)
α (1/cm)n (–)
Ks (cm/day)λL (cm)K d Norg (cm3/g)

TR
TR
TR
TR
TR
M1, M2, M3, M4, M5

0.20
0.05
1.50
100.00
0.10
0.00

0.60
0.30
4.00
1000.00
0.50
1000.00

μ L Norg (1/day)

M1, M2, M3, M4

0.01

3.00

μ S Norg (1/day)

M1, M2, M3, M4

0.01

3.00

K dNH 4 + (cm3/g)

M1, M2, M3, M4, M5

0.00

20.00
3.00

μ LNH 4 + (1/day)

M1, M30.01

μ SNH 4 + (1/day)

M1, M30.01

3.00

μ LNO2 − (1/day)

M1, M20.01

3.00

γNorg (1/day)

M5

0.20

10.0

YNorg (–)

M5

0.05

0.25

S Norg (mg/L)

M5

0.50

8.00

b Norg (1/day)

M5

0.01

0.15

γNH4+ (1/day)
YNH4+ (–)
SNH4+ (mg/L)
bNH4+ (1/day)
γNO2- (1/day)
YNO2- (–)
SNO2− (mg/L)
bNO2− (1/day)

M2, M4, M5
M2, M4, M5
M2, M4, M5
M2, M4, M5
M3, M4, M5
M3, M4, M5
M3, M4, M5
M3, M4, M5

0.20
0.05
0.50
0.01
0.20
0.05
0.50
0.01

0.90
0.25
8.00
0.15
0.90
0.25
8.00
0.15

model run. Modeled nitrite and nitrate concentrations at four diﬀerent
depths are compared with measured data to validate the results of thesurrogate-based analysis. This aim of this step is to identify plausible
signiﬁcant deviations between the higher and lower ﬁdelity models,
and to provide a general assessment of the surrogate-based analysis.
3. Results
3.1. Inverse estimation of soil hydraulic and solute transport parameters

2.3.4. Biogeochemical model selection
50th percentile values of the posterior distributions of soil hydraulicand transport parameters obtained in the scenario TR are then used in
HYDRUS-1D for the Bayesian model selection. Measured nitrite andnitrate concentrations from diﬀerent vertical proﬁles are averaged toobtain average concentrations for 8 observation points at diﬀerent
depths (Fig. 1). The number of unknown parameters in biogeochemical
models varies between 7 (M1) and 14 (M5) (Table 3). A preliminary
sensitivity analysis to identify unimportant factors is avoided due to therelatively low dimensionality of the inverse problem. In such circumstances, the Bayesian analysis normally reveals which parameters areinﬂuential and can be inferred from measured data. The soil bulk
density is set to 1.5 g/cm3, which is a common value for sandy sediments. Molecular diﬀusion coeﬃcients are assumed 1.92 cm2 day−1 for
all nitrogen species (Langergraber and Simunek, 2006). Parameterbounds reported in Table 2 are based on a literature review (e.g., Bitton,
2010; Hanson et al., 2006; Leslie Grady et al., 2011) and slightly expanded to detect particular patterns in the Bayesian optimization. Indeed, posterior parameter distributions located near the bounds of thesearch space generally indicate the theoretical inadequacy of the model,
which can only describe measured data by using unphysical parametervalues. Finally, the calculated Bayesian model evidence is used for themodel selection and to summarize model diﬀerences.

The joint and marginal posterior distributions of the soil hydraulicand transport parameters obtained from Scenario TR are shown in
Fig. 4 and have a markedly unimodal pattern characterized by a leptokurtic behavior. Consequently, the conﬁdence intervals (Table 3) arenarrow, thus indicating good identiﬁability of optimized parameters.
The Bayesian analysis reveals a well-posed inverse problem characterized by low uncertainty, mainly stemming from the simultaneous use of
informative measured concentrations at diﬀerent depths, which constrain the inverse estimation problem well. This is further conﬁrmed by
the limited correlation shown in the joint posterior distributions. Estimated parameter ranges are in agreement with what has been reported
in the literature for sandy soils (e.g., Thoma et al., 2014). Only thesaturated hydraulic conductivity exhibits a relatively lower value (e.g.,
Carsel and Parrish, 1988), indicating the compaction of the aquifer
model during the experiment. The longitudinal dispersivity λL ranges
between 0.11 and 0.19 cm, which corresponds with the values reported
in Sato et al. (2003) and Vanderborght and Vereecken (2007). Water
ﬂow is characterized by vertical inﬁltration in the unsaturated zone,
and by its gradual transition to horizontal ﬂow when approaching thesaturated zone. Such behavior is particularly evident near the outﬂow.
Posterior predictive checks (grey lines in Fig. 5) conﬁrm the accuracy of the numerical model in describing the coupled water-solutetransport in the aquifer. Fig. 5 shows a slight underestimation of the
tracer concentration at z = −20 cm in the vertical proﬁle D (Fig. 1).
This can be explained by lower weights assigned to the observation
points in the saturated zone (i.e., z < −25 cm) during the Bayesiananalysis, thus their less accurate ﬁt by the model. The overall accuracy

2.3.5. Model validation
50th percentile values of the posterior distributions obtained in theBayesian surrogate-based analysis (i.e., HYDRUS-1D) are used in thehigh-ﬁdelity model (i.e., HYDRUS-2D) to perform a single forward
8

Journal of Hydrology 584 (2020) 124681

G. Brunetti, et al.

Fig. 4. Joint (below diagonal) and marginal (diagonal) posterior distributions of the soil hydraulic and transport parameters for Scenario TR. The isolines and thedashed lines indicate the 12, 39, 68, and 84% probability levels and the 95% conﬁdence intervals, respectively. The axis ranges correspond to the parameter bounds
reported in Table 2.

3.2. Biogeochemical model selectionof the model predictions conﬁrms the assumption of homogeneous soil
and the negligible inﬂuence of heterogeneity on water and solute
movement. The estimated 50th percentile values of soil hydraulic andtransport parameters (Table 3) are used in the next Bayesian Model
Selection, where only adsorption and reaction parameters are investigated.

3.2.1. Marginal likelihoods
The MULTINEST algorithm required from 1211 (M1) to 2671 (M4)
iterations to satisfy the convergence criterion (Eq. (14)). In particular,
scenario M5 required, on average, approximately 3 days of computation
to perform a single Bayesian analysis. This behavior is related to anincrease in both the dimensionality of the inverse problem (i.e., 14
parameters) and the mathematical complexity of the numerical model.
Indeed, the addition of multiple Monod-type reactions (Eq. (8)) makes

Fig. 5. A comparison between the model predictions (grey lines) obtained by random sampling of 100 solutions from the posterior distributions (Scenario TR) andmeasured concentrations at four diﬀerent observation points.
9

Journal of Hydrology 584 (2020) 124681

G. Brunetti, et al.

Table 4
Median values and 95% conﬁdence intervals of the soil hydraulic, transport, and reactive parameters for six modeling scenarios (TR, M1, M2, M3, M4, and M5). The
estimated Bayesian model evidence and the coeﬃcient of variation are reported in the last two rows.
Parameter

Scenario ID
TR

M1

M2

M3

M4

M5

5%

50%

95%

5%

50%

95%

5%

50%

95%

5%

50%

95%

5%

50%

95%

5%

50%

95%

θs (cm /cm )α (1/cm)n (–)
Ks (cm/day)λL(cm)
K d Norg (cm3/g)

0.37
0.21
3.28
193
0.11
–

0.38
0.22
3.47
216
0.14
–

0.39
0.25
3.72
243
0.19
–

–
–
–
–
–
77

0.38
0.22
3.47
216.1
0.14
398

–
–
–
–
–
738

–
–
–
–
–
32

0.38
0.22
3.47
216.1
0.14
86

–
–
–
–
–
628

–
–
–
–
–
631

0.38
0.22
3.47
216
0.14
772

–
–
–
–
–
897

–
–
–
–
–
229

0.38
0.22
3.47
216
0.14
273

–
–
–
–
–
330

–
–
–
–
–
134

0.38
0.22
3.47
216
0.14
646

–
–
–
–
–
963

μ L Norg (1/day)

–

–

–

0.20

1.45

2.80

0.22

1.51

2.80

0.24

1.62

2.83

0.28

1.61

2.76

–

–

–

μ S Norg (1/day)

–

–

–

0.07

0.14

2.30

1.35

2.64

2.97

0.08

0.09

0.11

2.38

2.78

2.98

–

–

–

K dNH 4 + (cm3/g)

–

–

–

4.50

8.05

9.80

2.30

4.74

9.41

7.26

9.04

9.90

3.90

4.97

6.38

1.25

2.06

4.74

μ LNH 4 + (1/day)

–

–

–

0.22

1.47

2.82

–

–

–

0.21

1.25

2.69

–

–

–

–

–

–

3

3

μ SNH 4 + (1/day)

–

–

–

0.18

0.30

0.78

–

–

–

0.30

0.38

0.53

–

–

–

–

–

–

μ LNO2 − (1/day)

–

–

–

1.29

2.28

2.94

1.44

2.23

2.91

–

–

–

–

–

–

–

–

–
9.77

γNorg (1/day)

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

4.71

7.87

YNorg (–)

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

0.06

0.14

0.23

S Norg (mg/L)

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

0.71

2.55

4.01

–

–

–

–

–

–

–

–

0.01

0.06

0.13

0.82
0.25
5.07
0.11
–
–
–
–

–
–
–
–
–
–
–
–
0.29
0.39
0.06
0.14
2.25
4.80
0.09
0.13
248.6 ± 0.35
0.002
1569

–
–
–
–
0.48
0.23
7.50
0.15

0.64
0.70
0.16
0.22
0.51
0.65
0.08
0.12
0.32
0.36
0.05
0.06
5.39
7.05
0.12
0.14
294.09 ± 0.45
0.002
2671

0.75
0.25
0.98
0.14
0.40
0.09
7.92
0.15

0.34
0.37
0.05
0.05
0.97
2.94
0.02
0.05
0.50
0.53
0.14
0.19
7.40
7.80
0.14
0.15
307.8 ± 0.43
0.003
2532

b Norg (1/day)

–

–

–

–

–

–

–

γNH4+ (1/day)
YNH4+ (–)
SNH4+ (mg/L)
bNH4+ (1/day)
γNO2- (1/day)
YNO2- (–)
SNO2− (mg/L)
bNO2− (1/day)Z
CV
Iterations

–
–
–
–
–
–
–
–

–
–
–
–
–
–
–
–

–
–
–
–
–
–
–
–

–
–
–
–
–
–
–
–
–
–
–
–
–
–
–
–
131.19 ± 0.23
0.001
1211

–
–
–
–
–
–
–
–

0.36
0.69
0.10
0.21
0.65
1.26
0.01
0.05
–
–
–
–
–
–
–
–
233.8.09 ± 0.29
0.001
1299

the system of partial diﬀerential equations (Eq. (4)) numerically stiﬀ,
thus requiring smaller time steps to achieve convergence. The coeﬃcient of variation for multiple MULTINEST runs is very low for allmodeling scenarios (Table 4) and does not aﬀect the model ranking,thus indicating an overall satisfactory accuracy and stability of the algorithm. This is also conﬁrmed by the low error in the estimatedmarginal likelihoods (Table 4).
The biogeochemical models M1 and M5 exhibit the lowest (i.e.,
131.19 ± 0.23) and the highest (i.e., 307.8 ± 0.43) marginal likelihoods, respectively, indicating that the information content of themeasured data supports the use of the full Monod model to represent
nitrogen turnover in the aquifer. A closer look at the intercomparison
between diﬀerent models highlights interesting patterns. First, the useof the Monod kinetics instead of the ﬁrst-order reaction to represent the
AOBs reaction results in a signiﬁcant increase in the marginal likelihood, which almost doubles from 131.19 (M1) to 233.8 (M2). Similar
behavior is encountered when simulating the NOBs bacterial growth(M3). This provides the ﬁrst indication that microbial growth and degradation play a fundamental role in the analyzed case study.
Interestingly, Z is higher for model M3 compared to model M2, although they both share the same level of theoretical model complexity.
This behavior mainly indicates that the NOBs reaction has a larger
impact on nitrogen turnover in the aquifer than ammonia reduction,
though the diﬀerence in the estimated Bayesian evidence can also be
partially attributed to a higher number of parameters of model M3. The
use of the Monod kinetics for both nitrifying bacteria (M4) leads to afurther, but less substantial, increase of the marginal likelihood (i.e.,
294.09 ± 0.45). This increase is even less signiﬁcant when includingin the model ammonifying bacteria growth and degradation (M5), thus
indicating a saturation of the model complexity. Such behavior can beclearly inferred from Fig. 6, which shows how the marginal likelihood

0.40
0.06
6.57
0.09
0.56
0.24
8.00
0.15

Fig. 6. Estimated marginal likelihoods against the number of parameters. The
dashed grey line represents a second-order polynomial ﬁt line.

grows nonlinearly with the number of parameters, and how such dynamics seems to follow a second-order polynomial trend characterized
by a maximum around 14 parameters, which corresponds to the parameterization of model M5. From this perspective, it can be hypothesized that the use of a more complex model would result in overparameterization and worse model generalizability.
3.2.1.1. Posterior distributions. Marginal posterior distributions andmedian values and 95% conﬁdence intervals of adsorption and
10

Journal of Hydrology 584 (2020) 124681

G. Brunetti, et al.

μ LNH 4 + , are characterized by ﬂat posteriors and high uncertainty, suggesting a low inﬂuence of these parameters on the simulated process
and that their values cannot be inferred from measured data. This behavior is common in all models. Conversely, the ﬁrst-order decay
coeﬃcients in the solid phase, μ SNorg and μ SNH 4 + , are well estimated,thus indicating their preeminent role in the numerical description of thenitrogen turnover in the aquifer. The posterior distributions of the nitriﬁcation rate in the liquid phase, μ LNO2 − , are highest near the rightboundary for both models M1 and M2, conﬁrming an overall model
inadequacy again.
The Monod kinetics parameters for the ammoniﬁcation process inmodel M5 are generally characterized by high uncertainty, suggestingthat the information content of the data is not suﬃcient to well constrain their values. The uniform prior is reduced only for the maximum
growth rate and the half-saturation constant, while the yield and thelysis coeﬃcient are only partially inferred. This behavior can be anindication of potential overﬁtting.
Microbial reaction parameters for AOBs show mixed posterior patterns. The maximum aerobic growth rate, γNH4+, is characterized by
leptokurtic distributions in biogeochemical models M2, M4, and M5,
and this behavior exacerbates with an increase in model complexity.
The growth rate of AOBs, γNH4+, in model M4 ranges between 0.64 and0.75 1/day, which is in line with typical growth rates at 20 °C reported
in the literature (e.g., Leslie Grady et al., 2011), while its range appears
to be underestimated for model M5. Again, this behavior can be explained by a high correlation shown in Fig. 7, which indicates potential
overﬁtting. Posterior distributions of half-saturation constants for AOBs
are right-skewed and accumulated near the boundaries. The uncertainty is higher for models M2 and M5, while its range shrinks andvaries between 0.51 and 0.98 mg/L for model M4, which is in agreement with what has been reported in the literature (Bitton, 2010; Leslie
Grady et al., 2011) and conﬁrms that AOBs have a lower half-saturation
constant compared to NOBs (Sharma and Ahlert, 1977). Posterior distributions of the yield coeﬃcient, YNH4+, show a diﬀerent behavior
among models. In particular, models M2 and M4 have peaks near theleft boundary, while model M5 is characterized by very low yield values. The latter appears to contradict what was observed in other studies (e.g., Bitton, 2010), which reported a high value of YNH4+, thus
against suggesting potential overﬁtting of model M5. A behavior this,which is further conﬁrmed by the relatively low value of the lysis
coeﬃcient for AOBs observed in models M2 and M5.
On the other hand, microbial reaction parameters for NOBs show
similar posterior patterns. The maximum aerobic growth rate γNO2- formodels M3, M4, and M5 is well identiﬁed and characterized by limited
uncertainty. The high positive correlation observed in Fig. 7 for model

reaction parameters for all investigated biogeochemical models areshown in Fig. 8 and Table 4, respectively. At the ﬁrst inspection,
posteriors exhibit mixed patterns. Generally, parameters’ uncertainty
tends to decrease when the model complexity increases, indicating thetheoretical inadequacy of simple parameterizations in reproducing
nitrogen turnover in the aquifer. For model M1, the uniform prior is
reduced for four parameters, and among them, only two factors, μ SNorgand μ SNH 4 + , are well identiﬁed. On the other hand, the reduction ofalmost all priors for models M4 and M5 suggests that bacterial growth
dynamics is the mechanism responsible for reaction processes and thatthe ﬁrst-order reactions are not suﬃcient to explain thisbiogeochemical process.
The posterior distribution of the adsorption coeﬃcient for organic
nitrogen, K d Norg , varies signiﬁcantly between diﬀerent biogeochemical
models. On average (Table 4), its value is high and indicates strong
sorption of organic nitrogen to the solid phase, which is compatible
with observations in the transport studies on organic compounds insoils (Mariano et al., 2016). Interestingly, the posterior distributions of
K d Norg for models M3 and M4 are leptokurtic and unimodal, though their
median values (Table 4) diﬀer signiﬁcantly, while higher uncertainty is
encountered for models M1 and M5. While the behavior of the former
can be explained by an overall theoretical inadequacy of the parameterization, the diﬀerence between M4 and M5 is surprising, considering the similar level of complexity and estimated marginal likelihoods. However, correlation matrices (Fig. 7) can help in explaining
this aspect. In particular, Fig. 7 shows a strong negative correlation
between K d Norg and K dNH 4 + for model M5, thus implying that low sorption of organic nitrogen is compensated in the simulation by higher
ammonium sorption. Overall, there is an increase in parameter correlations in model M5, which can indicate overﬁtting, as already discussed for Fig. 6.
On the other hand, the posterior distribution of the adsorption
coeﬃcient for ammonium, K dNH 4 +, is skewed in all biogeochemical
models, thus suggesting its prominent role in nitrogen turnover. The
uncertainty range decreases with an increase in model complexity andshifts towards more realistic sorption values. In particular, models M4
and M5 present leptokurtic posterior distributions, which are in linewith values reported in the literature, generally ranging between 0.5
and 3.5 cm3/g (Hanson et al., 2006; Jellali et al., 2010; Li et al., 2015).
On the other hand, models M1 and M3 exhibit a dense probability region near the right boundary of the parameter space (i.e.,
K dNH 4 + ≈ 10 cm3/g), which is signiﬁcantly overestimated. As a consequence of these high sorption rates, the ﬁrst-order decay coeﬃcientsin the liquid phase for organic nitrogen and ammonium, μ L Norg and

Fig. 7. Heatmaps of the correlation matrices for models M4 (left) and M5 (right) calculated using the samples obtained by the MULTINEST algorithm.
11

Journal of Hydrology 584 (2020) 124681

G. Brunetti, et al.

Fig. 8. Marginal posterior distributions of the adsorption and reaction parameters for the biogeochemical modeling scenarios M1, M2, M3, M4, and M5. The dashedlines indicate 95% conﬁdence intervals. The axis ranges correspond to the parameter bounds reported in Table 2.

3.2.1.2. Model predictive checks. To better investigate the models’
behavior, the parameter uncertainty is propagated to generate
posterior predictive checks, which are used to “look for systematic
discrepancies between real and simulated data” (Gelman et al., 2004),
thus representing a fundamental tool for model criticism. Posterior
predictive checks of simulated nitrite concentrations for all
biogeochemical models are obtained by random sampling of 100
solutions from the posterior parameters’ distributions (Fig. 9). An
increase and a decrease in the quality of ﬁtting and predictive
uncertainty, respectively, is evident at the ﬁrst inspection. In
particular, the use of constant ﬁrst-order decay reactions in model M1
prevents the model from correctly simulating transient accumulation of
nitrites and delayed nitrate production at the beginning of theM5 can explain the slight overestimation of the growth rate. The posterior distribution of the half-saturation constant, SNO2−, is left-skewed,
and its uncertainty decreases with model complexity. All ranges are inagreement with what has been reported in the literature (Bitton, 2010;
Leslie Grady et al., 2011) and conﬁrm that NOBs have a higher halfsaturation constant (Sharma and Ahlert, 1977). The NOBs yield coefﬁcient, YNO2−, is signiﬁcantly higher for model M5 and suggest that
NOBs were present in higher number compared to AOBs, which contradicts what has been reported in the literature (Bitton, 2010). Finally,
lysis rates show a left-skewed posterior distribution, with bNO2− characterized by a lower uncertainty. This higher value of the death rate for
NOBs is compatible with what was reported in Copp and Murphy(1995), who reported a typical value of 0.14 1/day.
12

Journal of Hydrology 584 (2020) 124681

G. Brunetti, et al.

Fig. 9. Measured average NO2− (black circles) at four diﬀerent depths and posterior predictive checks (grey lines) for the biogeochemical modeling scenarios M1,
M2, M3, M4, and M5 (from top to bottom) obtained by random sampling of 100 solutions from the posterior distributions.

characterized by a slightly lower predictive uncertainty than model
M5, which again can indicate potential overﬁtting.
Similarly, posterior predictive checks of simulated nitrate concentrations for all biogeochemical models are obtained by random
sampling of 100 solutions from the posterior parameter distributions(Fig. 10). Predictive checks of nitrate concentrations conﬁrm andhighlight the inaccuracy of the model that uses ﬁrst-order reactions to
describe nitrogen reactions in the aquifer. The model bias is particularly
evident during the ﬁrst 50 days of the numerical simulation when theinjected manure is characterized by high organic nitrogen and ammonium concentrations that favor mineralization processes and microbial
growth (Table 1). After this period, nitrates are dominant in the inﬂow,
thus shifting the physical process towards a pure advection–dispersion
process with limited reactions. As expected under such circumstances,
the simple model M1 accurately reproduced NO3− concentrations, evencapturing the ﬁnal concentration rise at z = −10 cm. To summarize,
posterior predictive checks indicate that the model M1 can only partially reproduce measured data and that it cannot fully account forprocesses occurring early during the experiment.
Model M2 improves the ﬁtting but exhibits a very high uncertainty,

experiment. The situation slightly improves for model M2, which canreproduce the nitrites peak at the beginning of the experiment, but
signiﬁcantly overestimate their concentration in the remaining part ofthe experiment. This overestimation is particularly evident in theunsaturated zone, and it indicates high AOBs growth, which is notcompensated by the simultaneous growth of NOBs in model M2. The
situation improves for model M3, which slightly underestimates nitrite
accumulation, and which is characterized by low predictive
uncertainty. The anticipation in nitrite production is mainly due tothe assumption of a constant ﬁrst-order reaction for AOBs, which does
not account for the delay in the bacteria growth. The combination ofmultiple Monod kinetics in models M4 and M5 improves the ﬁtting andlimits the uncertainty, thus suggesting that bacterial growth andreaction were the mechanisms driving nitrogen turnover in theaquifer. This conﬁrms the ﬁndings reported in Cheyns et al. (2010)
and De Wilde et al. (2009), who suggested that Monod kinetics
performed better when describing pesticide biodegradation in alaboratory mineralization experiment. A more precise estimation of
the bacteria population at the beginning of the experiment would havefurther improved model accuracy. Interestingly, model M4 is
13

Journal of Hydrology 584 (2020) 124681

G. Brunetti, et al.

Fig. 10. Measured average NO3− (black circles) at four diﬀerent depths and posterior predictive checks (grey lines) for the biogeochemical modeling scenarios M1,
M2, M3, M4, and M5 obtained by random sampling of 100 solutions from the posterior distributions.

compared to model M4. Both these factors concur with the small difference in marginal likelihood between models M4 and M5.

while predictions of model M3 are much more generalizable. This isinteresting since it indicates that NOBs growth has a greater inﬂuenceon nitrogen turnover than AOBs growth. However, both models M2 andM3 are aﬀected by signiﬁcant underestimation of nitrate concentrations
at the beginning of the experiment. This behavior disappears in modelM4 and to a greater extent, in model M5, which accurately reproduces
measured nitrate concentrations. The uncertainty band encompasses
measured data and shows good model generalizability. Predictive uncertainty in model M5 tends to increase at deeper depths compared tomodel M4, which is instead more conservative. The eﬀect of predictive
uncertainty is more evident during the middle time period of the simulation. At that stage of the simulation, the bacteria population isalready established, and reaction rates become extremely sensitive tokinetic parameters, which cannot be precisely inferred from measured
data due to their limited temporal resolution during this period. Indeed,
no observations were available between days 37 and 54. Few measurements in that period would likely reduce the model predictive
uncertainty. The improvement in the quality of the ﬁtting is the main
reason for the highest estimated marginal likelihood of model M5,
which is, on the other hand, limited by the slightly higher uncertainty

3.3. Model validation
The estimated 50th percentile values of reaction parameters formodel M5 (Table 3) are used in the high-ﬁdelity model, HYDRUS-2D,
for validation. A comparison between measured and modeled nitrite
and nitrate concentrations at four diﬀerent vertical proﬁles and depths
is shown in Fig. 11. HYDRUS-2D can accurately reproduce measured
data at several locations in the aquifer, thus indicating good predictive
capabilities. The model correctly describes the transient accumulationof nitrites at the beginning of the experiment, although with someunderestimation and delay, which was already observed with HYDRUS1D. This can be related to potential uncertainties in the estimation of
the soil hydraulic and transport parameters, as well as to biased initialconditions. Indeed, such initial deviations between modeled and measured concentrations are also noticeable for nitrates, which are overestimated during the ﬁrst few days of the simulated period. Even so, thiseﬀect is small and limited in time. After this initial period, the model
14

Journal of Hydrology 584 (2020) 124681

G. Brunetti, et al.

Fig. 11. A comparison between measured (circles) and HYDRUS-2D modeled (solid lines) nitrite (top) and nitrate (bottom) concentrations at four diﬀerent vertical
proﬁles (A, B, C, D) and depths.

As mentioned in the introduction, this study presents the ﬁrst application of Nested Sampling for the Bayesian calibration of a biogeochemical mechanistic model. Results suggest that this algorithm can be
successfully used to simultaneously estimate the Bayesian model evidence and parameter posterior distributions in typical vadose zone
modeling applications, thus providing a statistical basis for the modelcomparison and uncertainty assessment in the light of measured data.
The possibility of having this combined statistical outcome at no extra
cost is one of the main advantages of Nested Sampling compared to
traditional MCMC techniques such as ensemble samplers (e.g.,Goodman and Weare, 2010; Vrugt, 2016), which are frequently used in
combination with mechanistic models (e.g., Brunetti et al., 2019;
Wöhling and Vrugt, 2011). In our analysis, it is ﬁrst used to infer theposterior distributions of soil hydraulic and transport parameters, thenin a classic model selection framework to discriminate between alternative parameterization of reactive processes. Furthermore, as thenumber of iterations demonstrates, it is computationally eﬃcient, requiring a lower number of model executions. This behavior has alreadybeen observed in past studies focused on the comparison of Nested
Sampling with other MCMC samplers (Allison and Dunkley, 2013). The
coeﬃcients of variation for multiple MULTINEST runs were generally
very low, indicating good stability of the algorithm. As such, future
research should focus on providing a rigorous assessment of the advantages and drawbacks of Nested Sampling compared to classical
MCMC samplers in high and low dimensions for typical inverse estimation problems in vadose zone modeling.

more suitable model parameterization for reproducing measured data.
In particular, the higher marginal likelihood for model M5 indicatesthat bacterial growth played a major role in the entire nitrogen degradation chain during the laboratory experiment, even though its effect is more signiﬁcant for nitriﬁcation, as indicated by the low difference in the estimated Bayesian evidence between models M4 and M5(Table 4). On the other hand, model M1, which uses the ﬁrst-order
reactions, underﬁtted measured data (Fig. 9 and Fig. 10), especially atthe beginning of the experiment. This mainly depends on the measured
data, which indicates that diﬀerent geochemical processes are dominant in diﬀerent parts of the experiment. During the ﬁrst part of theexperiment, the mineralization processes are dominant, thus tipping thebalance towards complex models, which can describe the transientaccumulation of nitrite. Model M1, which cannot account for theseprocesses, is instead able to satisfactorily describe the second part of theexperiment, characterized by negligible mineralization. Another factor
that reduces the ﬁtting quality of the complex model is the uncertainty
stemming from the hypothesized initial bacteria population. This further conﬁrms the point that the choice of a complex model must be fully
justiﬁed by the information content of measured data, which govern themodel selection. The higher the information content of observations,the higher the model complexity allowed to the modeler.
In summary, this study conﬁrms that Monod kinetics is the mostsuited parameterization for laboratory mineralization experiments,
which are generally characterized by a strong inﬂuence of microbial
growth, especially in the ﬁrst part of the experiment. This does not
imply that this model parameterization is a priori better than the one
with the traditional ﬁrst-order coeﬃcients, but it is simply more suited
to reproduce available measured data for this case study, and theBayesian analysis statistically proves this. Such a numerical approachshould be used ubiquitously in any model calibration to have a statistical basis for model selection. Claiming the appropriateness of a model
parameterization by relying only on the results obtained using a single
set of optimized parameters, as frequently done, can be highly misleading and lead to poor model calibration, especially for biogeochemical modeling applications.

4.2. Biogeochemical model parameterization

4.3. Importance of measured data

predicts nitrate concentrations at diﬀerent locations well. Some bias isevident at z = −40 cm, but the overall quality of the ﬁtting is satisfactory. Nevertheless, it must be noted that the ﬁt quality does notdegenerate but improves compared to HYDRUS-1D, suggesting theoverall reliability of the Bayesian calibration.
4. Discussion
4.1. Nested sampling performance

The Bayesian model selection suggests that Monod kinetics is theConsequently, measured data are crucial in the model calibration.
15

Journal of Hydrology 584 (2020) 124681

G. Brunetti, et al.

(Forrester et al., 2007).

They must be exhaustive and representative of the investigated process.
To maximize the information content of measured data, the model canbe used in a prognostic way by performing a model-based Bayesian
optimal experimental design, which gives the researcher indications ofpotential experimental designs that reduce the model predictive uncertainty. For example, in our case, few measurements about the initial
distribution of bacteria in the aquifer could have improved the predictions of the complex model. This is particularly important whennovel model parameterizations are investigated since overﬁtting canundermine their generalizability, as well as their applicability forpractical applications where the model predictive capability is fundamental for designing purposes.

4.6. Surrogate-based analysis
The combination of the multimodal Nested Sampling algorithmwith a surrogate model can signiﬁcantly increase the overall computational eﬃciency of the analysis, though a rigorous assessment of thecalculated Bayesian model evidence and posterior distributions is recommended. In this perspective, the use of synthetic data is suggested
to avoid the inﬂuence of unpredictable sources of uncertainty often
encountered in measured data.
5. Conclusions

4.4. Importance of priors and posteriors assessment
The aim of this study was to investigate the main advantages andﬁndings of model selection in a typical biogeochemical modeling application dealing with nitrogen turnover in a controlled aquifer model
setup. Results from a multilevel numerical analysis raised several interesting points of discussion and provided future research directions on
how to handle model complexity in biogeochemical model calibrationfrom both a statistical and computational point of view. Certainly, thisresearch did not pretend to be exhaustive nor conclusive for such a
dynamic ﬁeld as environmental modeling, but instead aimed to be ascientiﬁc contribution, improving the current understanding of model
complexity and placing emphasis on the role of parsimony in model
selection. Research advancements in environmental sciences and engineering allow the modeler to simulate multiple physical processes
with diﬀerent levels of complexity. This modeling ﬂexibility is undoubtedly a valuable tool but should always be counterbalanced by an
overall sense of parsimony, which is established by the aim and circumstances of the numerical analysis. We want to underline that thisconcept should be central in any modeling application, and especiallyin environmental sciences, which are, and will be even more in the next
future, heavily dependent on numerical models.

The Bayesian model selection framework included ﬁve parameterizations, which combine Monod kinetics and ﬁrst-order reactions.
The dynamics of the estimated marginal likelihood (Fig. 6) and thecorrelation analysis (Fig. 7) suggest a plausible saturation of theavailable model complexity and potential overﬁtting, respectively, thus
indicating that a further increase in the model complexity could not bejustiﬁable. However, it is not possible to reject a priori the possibility
that measured data could be better described in terms of ﬁtting andgeneralizability, by other models that were not considered in this study.
Nevertheless, the prior model set was selected based on the available
prior information obtained from a preliminary analysis of measured
data. A new level of model complexity (e.g., denitriﬁcation) should be
carefully justiﬁed based on available information, considering thecomputational eﬀort required to calculate the marginal likelihood.
Since a single case study cannot be entirely representative of thephysical process investigated, and since diﬀerent unpredictable uncertainties can aﬀect measured data, both prior and posterior parameter
distributions must be carefully inspected and compared with the existing body of literature to assess their physical meaning. Priors arefundamental to rule out unphysical solutions and constrain the inverse
problem and must be set based on both a thorough literature survey andpreliminary experiments. For example, in our case, batch experimentson the adsorption of ammonium to the solid phase could have beenused to obtain a more informative prior about the distribution coeﬃcient. In contrast, posteriors with peaks near the boundary of the search
space are generally correlated with model inadequacy, a behavior that
we encountered mainly in the simple model in our analysis. In thisregard, the Bayesian analysis can also be viewed as a diagnostic tool todetect modeling issues.

Declaration of Competing Interest
The authors declare that they have no known competing ﬁnancial
interests or personal relationships that could have appeared to inﬂuence the work reported in this paper.
Acknowledgments
We thank Benjamin Gralher and Günter Hinreiner for the custommade design of the experimental setup, Niklas Gassen for the assistanceof implementing the experiments, and Petra Seibel for the analysis oftracer concentrations.

4.5. Computational cost
Computational cost is the bottleneck in the Bayesian calibration ofenvironmental models that cannot be performed without carrying out
an uncertainty assessment during model calibration. In our study, theuse of HYDRUS-1D as a lower-ﬁdelity surrogate has proven to be extremely useful. It allowed us to perform multiple Bayesian analyses in areasonable amount of time. Furthermore, the HYDRUS-2D validation
against measured data suggests that the results of the surrogate-based
analysis are reliable and that the simpliﬁcation induced by spatial reduction is acceptable.
Interestingly, results in terms of ﬁtting are even better when surrogate-based optimized parameters are used in the high-ﬁdelity model(i.e., HYDRUS-2D). This indicates that the two models have similar
likelihood landscapes and that a proper multi-ﬁdelity analysis based on
the use of hierarchical models with correction mapping (Robinson
et al., 2008) could have improved the Bayesian model selection both interms of accuracy and computational cost. In this view, future researchshould explore the use of the multi-ﬁdelity surrogate-based analysis forenvironmental model calibration, in particular the co-kriging techniques to combine data from both the higher and lower-ﬁdelity models

References
Ajami, N.K., Gu, C., 2010. Complexity in microbial metabolic processes in soil nitrogenmodeling: a case for model averaging. Stoch. Environ. Res. Risk Assess. 24, 831–844.
https://doi.org/10.1007/s00477-010-0381-4.
Akaike, H., 1973. Information Theory and an Extension of the Maximum Likelihood
Principle BT – Selected Papers of Hirotugu Akaike. In: Second International
Symposium on Information Theory, https://doi.org/10.1007/978-1-4612-1694-0_15.
Allison, R., Dunkley, J., 2013. Comparison of sampling techniques for Bayesian parameter
estimation. https://doi.org/10.1093/mnras/stt2190.
Asher, M.J., Croke, B.F.W., Jakeman, A.J., Peeters, L.J.M., 2015. A review of surrogate
models and their application to groundwater modeling. Water Resour. Res. 51,
5957–5973. https://doi.org/10.1002/2015WR016967.
Beauchamp, E.G., Kidd, G.E., Thurtell, G., 2010. Ammonia volatilization from liquid dairy
cattle manure in the ﬁeld. Can. J. Soil Sci. https://doi.org/10.4141/cjss82-002.
Bitton, G., 2010. Wastewater Microbiology: Fourth Edition, Wastewater Microbiology:
fourth ed., https://doi.org/10.1002/9780470901243.
Brewer, B.J., Pártay, L.B., Csányi, G., 2011. Diﬀusive nested sampling. Comput. Stat.
https://doi.org/10.1007/s11222-010-9198-8.
Brunetti, G., Šimůnek, 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 of soil hydraulic properties. Vadose Zo. J. 18. https://doi.org/10.2136/vzj2018.

16

Journal of Hydrology 584 (2020) 124681

G. Brunetti, et al.

Gribovszki, Z., Kalicz, P., Szilágyi, J., Kucsara, M., 2008. Riparian zone evapotranspiration estimation from diurnal groundwater level ﬂuctuations. J. Hydrol. https://doi.
org/10.1016/j.jhydrol.2007.10.049.
Gribovszki, Z., Szilágyi, J., Kalicz, P., 2010. Diurnal ﬂuctuations in shallow groundwaterlevels and streamﬂow rates and their interpretation – a review. J. Hydrol. https://doi.
org/10.1016/j.jhydrol.2010.02.001.
Guthery, F.S., Burnham, K.P., Anderson, D.R., 2003. Model selection and multimodel
inference: a practical information-theoretic approach. J. Wildl. Manage. https://doi.
org/10.2307/3802723.
Handley, W.J., Hobson, M.P., Lasenby, A.N., 2015. POLYCHORD: next-generation nested
sampling. Mon. Not. R. Astron. Soc. https://doi.org/10.1093/mnras/stv1911.
Hansen, S., Jensen, H.E., Nielsen, N.E., Svendsen, H., 1991. Simulation of nitrogen dynamics and biomass production in winter wheat using the Danish simulation model
DAISY. Fertil. Res. https://doi.org/10.1007/BF01051131.
Hanson, B.R., Šimůnek, J., Hopmans, J.W., 2006. Evaluation of urea–ammonium–nitratefertigation with drip irrigation using numerical modeling. Agric. Water Manage. 86,
102–113. https://doi.org/10.1016/J.AGWAT.2006.06.013.
Henrichs, M., Langergraber, G., Uhl, M., 2007. Modelling of organic matter degradationin constructed wetlands for treatment of combined sewer overﬂow. Total Environ.
Sci. https://doi.org/10.1016/j.scitotenv.2006.11.044.
Higson, E., Handley, W., Hobson, M., Lasenby, A., 2019. Dynamic nested sampling: an
improved algorithm for parameter estimation and evidence calculation. Stat.
Comput. https://doi.org/10.1007/s11222-018-9844-0.
Höge, M., Wöhling, T., Nowak, W., 2018. A primer for model selection: the decisive roleof model complexity. Water Resour. Res. 54, 1688–1715. https://doi.org/10.1002/
2017WR021902.
Holden, P.A., Fierer, N., 2005. Microbial processes in the vadose zone. Vadose Zo. J.
https://doi.org/10.2136/vzj2005.0001.
Homann, P.S., McKane, R.B., Sollins, P., 2000. Belowground processes in forest-ecosystembiogeochemical simulation models. Ecol. Manage For. https://doi.org/10.1016/
S0378-1127(00)00408-4.
Jellali, S., Diamantopoulos, E., Kallali, H., Bennaceur, S., Anane, M., Jedidi, N., 2010.
Dynamic sorption of ammonium by sandy soil in ﬁxed bed columns: evaluation ofequilibrium and non-equilibrium transport processes. J. Environ. Manage. https://doi.org/10.1016/j.jenvman.2009.11.006.
Kass, R.E., Raftery, A.E., 1995. Bayes factors. J. Am. Stat. Assoc. https://doi.org/10.1080/
01621459.1995.10476572.
Kavetski, D., Kuczera, G., Franks, S.W., 2006. Bayesian analysis of input uncertainty inhydrological modeling: 1. Theory. Water Resour. Res. 42. https://doi.org/10.1029/
2005WR004368.
Kimmins, J.P., Blanco, J.A., Seely, B., Welham, C., Scoullar, K., 2008. Complexity inmodelling forest ecosystems: how much is enough? Ecol. Manage. For. https://doi.
org/10.1016/j.foreco.2008.03.011.
Knorr, B., Maloszewski, P., Stumpp, C., 2016. Quantifying the impact of immobile water
regions on the fate of nitroaromatic compounds in dual-porosity media. J. Contam.
Hydrol. https://doi.org/10.1016/j.jconhyd.2016.05.002.
Kröbel, R., Sun, Q., Ingwersen, J., Chen, X., Zhang, F., Müller, T., Römheld, V., 2010.
Modelling water dynamics with DNDC and DAISY in a soil of the North China Plain: acomparative study. Environ. Model. Softw. 25, 583–601. https://doi.org/10.1016/J.
ENVSOFT.2009.09.003.
Laloy, E., Rogiers, B., Vrugt, J.A., Mallants, D., Jacques, D., 2013. Eﬃcient posterior
exploration of a high-dimensional groundwater model from two-stage Markov chain
Monte Carlo simulation and polynomial chaos expansion. Water Resour. Res. 49,
2664–2682. https://doi.org/10.1002/wrcr.20226.
Langergraber, G., Simunek, J., 2006. The Multi-component Reactive Transport Module
CW2D for Constructed Wetlands for the HYDRUS Software Package, Manual –Version 1.0.
Lenk, P., 2009. Simulation pseudo-bias correction to the harmonic mean estimator of
integrated likelihoods. J. Comput. Graph. Stat. https://doi.org/10.1198/jcgs.2009.
08022.
Leslie Grady, C.P., Daigger, G., Love, N., Filipe, C., 2011. Biological Wastewater
Treatment, third ed. CRC Press.
Li, Y., Šimůnek, J., Zhang, Z., Jing, L., Ni, L., 2015. Evaluation of nitrogen balance in a
direct-seeded-rice ﬁeld experiment using Hydrus-1D. Agric. Water Manage. 148,
213–222. https://doi.org/10.1016/J.AGWAT.2014.10.010.
Liu, P., Elshall, A.S., Ye, M., Beerli, P., Zeng, X., Lu, D., Tao, Y., 2016. Evaluating marginal
likelihood with thermodynamic integration method and comparison with several
other numerical methods. Water Resour. Res. https://doi.org/10.1002/
2014WR016718.
Lu, D., Ye, M., Neuman, S.P., 2011. Dependence of bayesian model selection criteria andﬁsher information matrix on sample size. Math. Geosci. https://doi.org/10.1007/
s11004-011-9359-0.
Ma, L., Ahuja, L.R., Nolan, B.T., Malone, R.W., Trout, T.J., Qi, Z., Midwest, U.S., 2012.
Root zone water quality model (RZWQM2): model use, calibration, and validation.
Trans. ASABE.
Ma, L., Malone, R.W., Heilman, P., Jaynes, D.B., Saseendran, S.A., Kanwar, R.S., Ascough,J.C., 2007. RZWQM simulated eﬀects of crop rotation, tillage, and controlled drainage on crop yield and nitrate-N loss in drain ﬂow. Geoderma 140, 260–271. https://doi.org/10.1016/J.GEODERMA.2007.04.010.
Mariano, E., Jones, D.L., Hill, P.W., Trivelin, P.C.O., 2016. Mineralisation and sorption ofdissolved organic nitrogen compounds in litter and soil from sugarcane ﬁelds. Soil
Biol. Biochem. https://doi.org/10.1016/j.soilbio.2016.10.004.
Mateo-Sagasta, J., Marjani, S., Turral, H., Burke, J., 2017. Water pollution from agriculture: a global review. FAO and IWMI.
Meng, X.L., Wong, W.H., 1996. Simulating ratios of normalizing constants via a simple
identity: a theoretical exploration. Stat. Sin.

06.0123.
Brunetti, G., Šimůnek, J., Piro, P., 2016. A comprehensive numerical analysis of the hydraulic behavior of a permeable pavement. J. Hydrol. 540, 1146–1161. https://doi.
org/10.1016/j.jhydrol.2016.07.030.
Brunetti, G., Šimůnek, J., Turco, M., Piro, P., 2017. On the use of surrogate-based modeling for the numerical analysis of Low Impact Development techniques. J. Hydrol.
548, 263–277. https://doi.org/10.1016/j.jhydrol.2017.03.013.
Cameira, M.R., Sousa, P.L., Farahani, H.J., Ahuja, L.R., Pereira, L.S., 1998. Evaluation of
the RZWQM for the simulation of water and nitrate movement in level-basin, fertigated maize. J. Agric. Eng. Res. 69, 331–341. https://doi.org/10.1006/JAER.1997.
0253.
Cao, T., Zeng, X., Wu, J., Wang, D., Sun, Y., Zhu, X., Lin, J., Long, Y., 2018. Integrating
MT-DREAMzs and nested sampling algorithms to estimate marginal likelihood andcomparison with several other methods. J. Hydrol. https://doi.org/10.1016/j.
jhydrol.2018.06.055.
Carsel, R.F., Parrish, R.S., 1988. Developing joint probability distributions of soil water
retention characteristics. Water Resour. Res. 24, 755–769. https://doi.org/10.1029/
WR024i005p00755.
Cheyns, K., Mertens, J., Diels, J., Smolders, E., Springael, D., 2010. Monod kinetics rather
than a ﬁrst-order degradation model explains atrazine fate in soil mini-columns:
Implications for pesticide fate modelling. Pollut. Environ. https://doi.org/10.1016/j.
envpol.2009.12.041.
Chib, S., 1995. Marginal likelihood from the Gibbs output. J. Am. Stat. Assoc. https://doi.
org/10.1080/01621459.1995.10476635.
Chib, S., Jeliazkov, I., 2001. Marginal likelihood from the metropolis–hastings output. J.
Am. Stat. Assoc. https://doi.org/10.1198/016214501750332848.
Copp, J.B., Murphy, K.L., 1995. Estimation of the active nitrifying biomass in activatedsludge. Water Res. 29, 1855–1862. https://doi.org/10.1016/0043-1354(94)00347-A.
Dai, H., Chen, X., Ye, M., Song, X., Hammond, G., Hu, B., Zachara, J.M., 2019. Using
bayesian networks for sensitivity analysis of complex biogeochemical models. Water
Resour. Res. 55, 3541–3555. https://doi.org/10.1029/2018WR023589.
De Wilde, T., Spanoghe, P., Mertens, J., Sniegowksi, K., Ryckeboer, J., Jaeken, P.,
Springael, D., 2009. Characterizing pesticide sorption and degradation in macro scale
biopuriﬁcation systems using column displacement experiments. Pollut. Environ.
https://doi.org/10.1016/j.envpol.2008.11.032.
Elshall, Ye, 2019. Making stepping stones out of stumbling blocks: a bayesian model
evidence estimator with application to groundwater transport model selection.
Water. https://doi.org/10.3390/w11081579.
Elsheikh, A.H., Hoteit, I., Wheeler, M.F., 2014. Eﬃcient Bayesian inference of subsurface
ﬂow models using nested sampling and sparse polynomial chaos surrogates. Methods
Appl. Mech. Eng. Comput. https://doi.org/10.1016/j.cma.2013.11.001.
Elsheikh, A.H., Wheeler, M.F., Hoteit, I., 2013. Nested sampling algorithm for subsurface
ﬂow model selection, uncertainty quantiﬁcation, and nonlinear calibration. Water
Resour. Res. 49, 8383–8399. https://doi.org/10.1002/2012WR013406.
European Environment Agency (EEA), 2018. European waters Assessment of status andpressures 2018. Luxembourg. https://doi.org/10.2800/303664.
Fangueiro, D., Chadwick, D., Bol, R., 2008. Assessment of the potential N mineralizationof diﬀerent particle-size fractions in two dairy cattle slurries. J. Plant Nutr. Soil Sci.
171, 313–315. https://doi.org/10.1002/jpln.200700226.
Feroz, F., Hobson, M.P., Bridges, M., 2009. MultiNest: An eﬃcient and robust Bayesianinference tool for cosmology and particle physics. MNRAS 398, 1601–1614. https://doi.org/10.1111/j.1365-2966.2009.14548.x.
Feroz, F., Hobson, M.P., Cameron, E., Pettitt, A.N., 2013. Importance Nested Sampling
and the MultiNest Algorithm.
Filipović, V., Romić, D., Romić, M., Borošić, J., Filipović, L., Mallmann, F.J.K., Robinson,
D.A., 2016. Plastic mulch and nitrogen fertigation in growing vegetables modify soil
temperature, water and nitrate dynamics: experimental results and a modeling study.
Agric. Water Manage. 176, 100–110. https://doi.org/10.1016/J.AGWAT.2016.04.
020.
Forrester, A.I., Sóbester, A., Keane, A.J., 2007. Multi-ﬁdelity optimization via surrogate
modelling. Proc. R Soc. A Math. Phys. Eng. Sci. 463, 3251–3269. https://doi.org/10.
1098/rspa.2007.1900.
Freni, G., Mannina, G., Viviani, G., 2011. Assessment of the integrated urban water
quality model complexity through identiﬁability analysis. Water Res. 45, 37–50.
https://doi.org/10.1016/J.WATRES.2010.08.004.
Friel, N., McKeone, J.P., Oates, C.J., Pettitt, A.N., 2017. Investigation of the widely applicable Bayesian information criterion. Comput. Stat. https://doi.org/10.1007/
s11222-016-9657-y.
Friel, N., Pettitt, A.N., 2008. Marginal likelihood estimation via power posteriors. J. R.
Stat. Soc. Ser. B Stat. Methodol. https://doi.org/10.1111/j.1467-9868.2007.00650.x.
Gärdenäs, A.I., Hopmans, J.W., Hanson, B.R., Šimůnek, J., 2005. Two-dimensionalmodeling of nitrate leaching for various fertigation scenarios under micro-irrigation.
Agric. Water Manag. 74, 219–242. https://doi.org/10.1016/J.AGWAT.2004.11.011.
Gelman, A., Carlin, J.B., Stern, H.S., Rubin, D.B., 2004. Bayesian Data Analysis – Third
Edition, Chapman Texts in Statistical Science Series. https://doi.org/10.1186/17541611-9-2.
Gelman, A., Hwang, J., Vehtari, A., 2014. Understanding predictive information criteria
for Bayesian models. Comput. Stat. https://doi.org/10.1007/s11222-013-9416-2.
Gelman, A., Meng, X.L., 1998. Simulating normalizing constants: from importance sampling to bridge sampling to path sampling. Sci. Stat. https://doi.org/10.1214/ss/
1028905934.
Goodman, J., Weare, J., 2010. Ensemble samplers with aﬃne invariance. Commun. Appl.
Math. Comput. Sci. 5, 65–80. https://doi.org/10.2140/camcos.2010.5.65.
Green, P.J., 1995. Reversible jump Markov chain monte carlo computation and Bayesianmodel determination. Biometrika. https://doi.org/10.1093/biomet/82.4.711.

17

Journal of Hydrology 584 (2020) 124681

G. Brunetti, et al.

2008WR006836.
Schwarz, G., 1978. Estimating the dimension of a model. Ann. Stat. https://doi.org/10.
1214/aos/1176344136.
Sharma, B., Ahlert, R.C., 1977. Nitriﬁcation and nitrogen removal. Water Res. https://doi.
org/10.1016/0043-1354(77)90078-1.
Shaw, J.R., Bridges, M., Hobson, M.P., 2007. Eﬃcient Bayesian inference for multimodal
problems in cosmology. MNRAS 378, 1365–1370. https://doi.org/10.1111/j.13652966.2007.11871.x.
Šimůnek, J., van Genuchten, M.T., Šejna, M., 2016. Recent developments and applications of the HYDRUS computer software packages. Vadose Zo. J. 15, 25. https://doi.
org/10.2136/vzj2016.04.0033.
Skilling, J., 2006. Nested sampling for general Bayesian computation. Bayesian Anal.
https://doi.org/10.1214/06-BA127.
Thoma, M.J., Barrash, W., Cardiﬀ, M., Bradford, J., Mead, J., 2014. Estimating unsaturated hydraulic functions for coarse sediment from a ﬁeld-scale inﬁltration. Exp.
Vadose Zo. J. https://doi.org/10.2136/vzj2013.05.0096.
van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils1. Soil Sci. Soc. Am. J. 44, 892–898. https://doi.org/10.
2136/sssaj1980.03615995004400050002x.
van Middelkoop, J.C., Holshof, G., 2017. Nitrogen fertilizer replacement value of concentrated liquid fraction of separated pig slurry applied to grassland. Commun. Soil
Sci. Plant Anal. 48, 1132–1144. https://doi.org/10.1080/00103624.2017.1323101.
Vanderborght, J., Vereecken, H., 2007. Review of dispersivities for transport modeling insoils. Vadose Zo. J. https://doi.org/10.2136/vzj2006.0096.
Volpi, E., Schoups, G., Firmani, G., Vrugt, J.A., 2017. Sworn testimony of the model
evidence: Gaussian Mixture Importance (GAME) sampling. Water Resour. Res.
https://doi.org/10.1002/2016WR020167.
Vrugt, J.A., 2016. Markov chain Monte Carlo simulation using the DREAM software
package: theory, concepts, and MATLAB implementation. Environ. Model. Softw. 75,
273–316. https://doi.org/10.1016/j.envsoft.2015.08.013.
Wagner, B.J., Gorelick, S.M., 1987. Optimal groundwater quality management underparameter uncertainty. Water Resour. Res. 23, 1162–1174. https://doi.org/10.1029/
WR023i007p01162.
Watanabe, S., 2013. A widely applicable bayesian information criterion. J. Mach. Learn.
Res.
Wehrer, M., Jaesche, P., Totsche, K.U., 2012. Modeling the kinetics of microbial degradation of deicing chemicals in porous media under ﬂow conditions. Pollut.
Environ. https://doi.org/10.1016/j.envpol.2012.04.016.
Wei, G., Zhang, X., Ye, M., Yue, N., Kan, F., 2019. Bayesian performance evaluation of
evapotranspiration models based on eddy covariance systems in an arid region. Earth
Syst. Sci. Hydrol. https://doi.org/10.5194/hess-23-2877-2019.
Weinberg, M.D., 2012. Computing the bayes factor from a markov chain monte carlo
simulation of the posterior distribution. Bayesian Anal. https://doi.org/10.1214/12BA725.
Wöhling, T., Vrugt, J.A., 2011. Multiresponse multilayer vadose zone model calibration
using Markov chain Monte Carlo simulation and ﬁeld water retention data. Water
Resour. Res. 47. https://doi.org/10.1029/2010WR009265.
Zeng, X., Ye, M., Wu, J., Wang, D., Zhu, X., 2018. Improved nested sampling and surrogate-enabled comparison with other marginal likelihood estimators. Water Resour.
Res. https://doi.org/10.1002/2017WR020782.

Millington, R.J., Quirk, J.P., 1960. Transport in porous media. Trans 7th int Congr. Soil
Sci. https://doi.org/10.1016/B978-0-08-044490-1.50005-3.
Mukherjee, P., Parkinson, D., Liddle, A.R., 2006. A nested sampling algorithm for cosmological model selection. Astrophys J. 638, L51–L54. https://doi.org/10.1086/
501068.
Myung, I.J., 2000. The importance of complexity in model selection. J. Math. Psychol.
https://doi.org/10.1006/jmps.1999.1283.
Neuman, S.P., 2003. Maximum likelihood Bayesian averaging of uncertain model predictions. Stoch. Environ. Res. Risk Assess. https://doi.org/10.1007/s00477-0030151-7.
Pajor, A., 2017. Estimating the marginal likelihood using the arithmetic mean identity.
Bayesian Anal. https://doi.org/10.1214/16-ba1001.
Pande, S., Arkesteijn, L., Savenije, H., Bastidas, L.A., 2015. Hydrological model parameter
dimensionality is a weak measure of prediction uncertainty. Hydrol. Earth Syst. Sci.
Discuss. 12, 3945–4004. https://doi.org/10.5194/hessd-12-3945-2015.
Parton, W.J., 1996. The CENTURY model. In: Evaluation of Soil Organic Matter Models,
https://doi.org/10.1007/978-3-642-61094-3_23.
Paulson, N.H., Jennings, E., Stan, M., 2019. Bayesian strategies for uncertainty quantiﬁcation of the thermodynamic properties of materials. J. Eng. Sci. Int. https://doi.
org/10.1016/j.ijengsci.2019.05.011.
Pedretti, D., Fernàndez-Garcia, D., Bolster, D., Sanchez-Vila, X., 2013. On the formation of
breakthrough curves tailing during convergent ﬂow tracer tests in three-dimensional
heterogeneous aquifers. Water Resour. Res. https://doi.org/10.1002/wrcr.20330.
Peña-Haro, S., Pulido-Velazquez, M., Llopis-Albert, C., 2011. Stochastic hydro-economic
modeling for optimal management of agricultural groundwater nitrate pollution
under hydraulic conductivity uncertainty. Environ. Model. Softw. 26, 999–1008.
https://doi.org/10.1016/J.ENVSOFT.2011.02.010.
Peters, K., Hjorth, M., Jensen, L.S., Magid, J., 2011. Carbon, nitrogen, and phosphorus
distribution in particle size-fractionated separated pig and cattle slurry. J. Environ.
Qual. https://doi.org/10.2134/jeq2010.0217.
Ramos, T.B., Šimůnek, J., Gonçalves, M.C., Martins, J.C., Prazeres, A., Pereira, L.S., 2012.
Two-dimensional modeling of water and nitrogen fate from sweet sorghum irrigatedwith fresh and blended saline waters. Agric. Water Manage. https://doi.org/10.
1016/j.agwat.2012.05.007.
Razavi, S., Tolson, B.A., Burn, D.H., 2012. Review of surrogate modeling in water resources. Water Resour. Res. 48. https://doi.org/10.1029/2011WR011527.
Robinson, T.D., Eldred, M.S., Willcox, K.E., Haimes, R., 2008. Surrogate-based optimization using multiﬁdelity models with variable parameterization and corrected space
mapping. AIAA J. https://doi.org/10.2514/1.36043.
Salazar, O., Nájera, F., Tapia, W., Casanova, M., 2017. Evaluation of the DAISY model for
predicting nitrogen leaching in coarse-textured soils cropped with maize in the
Mediterranean zone of Chile. Agric. Water Manage. 182, 77–86. https://doi.org/10.
1016/J.AGWAT.2016.12.005.
Sato, T., Tanahashi, H., Loáiciga, H.A., 2003. Solute dispersion in a variably saturated
sand. Water Resour. Res. https://doi.org/10.1029/2002WR001649.
Schöniger, A., Wöhling, T., Samaniego, L., Nowak, W., 2014. Model selection on solid
ground: rigorous comparison of nine ways to evaluate Bayesian model evidence.
Water Resour. Res. 50, 9484–9513. https://doi.org/10.1002/2014WR016062.
Schoups, G., van de Giesen, N.C., Savenije, H.H.G., 2008. Model complexity control forhydrologic prediction. Water Resour. Res. 44. https://doi.org/10.1029/

18