Journal of Hydrology 548 (2017) 263–277

Contents lists available at ScienceDirect

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

Research papers

On the use of surrogate-based modeling for the numerical analysis of
Low Impact Development techniques
Giuseppe Brunetti a,⇑, Jirka Šimůnek b, Michele Turco a, Patrizia Piro aabDepartment of Civil Engineering, University of Calabria, Rende, CS 87036, Italy
Department of Environmental Sciences, University of California, Riverside, Riverside, CA 92521, USAa r t i c l e

i n f oArticle history:
Received 19 October 2016
Received in revised form 28 February 2017
Accepted 6 March 2017
Available online 8 March 2017
This manuscript was handled by Corrado
Corradini, Editor-in-Chief, with theassistance of Philip Brunner, Associate
Editor
Keywords:
Urban hydrology
Surrogate-based modeling
Sensitivity analysis
Infiltration
LIDs

a b s t r a c t
Mechanistic models have proven to be accurate tools for the numerical analysis of the hydraulic behaviorof Low Impact Development (LIDs) techniques. However, their widespread adoption has been limited bytheir computational cost. In this view, surrogate modeling is focused on developing and using a computationally inexpensive surrogate of the original model. While having been previously applied to various
water-related and environmental modeling problems, no studies have used surrogate models for the
analysis of LIDs. The aim of this research thus was to investigate the benefit of surrogate-based modelingin the numerical analysis of LIDs. The kriging technique was used to approximate the deterministic
response of the widely used mechanistic model HYDRUS-2D, which was employed to simulate thevariably-saturated hydraulic behavior of a contained stormwater filter. The Nash-Sutcliffe efficiency(NSE) index was used to compare the simulated and measured outflows and as the variable of interestfor the construction of the response surface. The validated kriging model was first used to carry out a
Global Sensitivity Analysis of the unknown soil hydraulic parameters of the filter layer, revealing that
only the shape parameter a and the saturated hydraulic conductivity Ks significantly affected the model
response. Next, the Particle Swarm Optimization algorithm was used to estimate their values. The NSE
value of 0.85 indicated a good accuracy of estimated parameters. Finally, the calibrated model was validated against an independent set of measured outflows with a NSE value of 0.8, which again corroborated the reliability of the surrogate-based optimized parameters.
Ó 2017 Elsevier B.V. All rights reserved.

1. Introduction
During the last few decades, stormwater management has
become a major component of the prevention of floods in urban
areas and for the preservation of water resources. An increase ofimpervious surfaces, connected with demographic growth, hasaltered the natural hydrological cycle by reducing the infiltrationand evaporation capacity of urban catchments while also increasing surface runoff. In their report, the Organization for Economic
Co-operation and Development (OECD) (2013) identified an
expected increase in flash and urban floods in large parts of Europe
as one of the major issues for the future.
In this context, urban drainage systems play a fundamental rolein improving the resilience of cities. In recent years, an innovative
approach to land development known as a Low Impact Development (LID) has gained increasing popularity. A LID is a ‘green’
approach to storm water management that seeks to mimic the nat⇑ Corresponding author.
E-mail address: giusep.bru@gmail.com (G. Brunetti).
http://dx.doi.org/10.1016/j.jhydrol.2017.03.013
0022-1694/Ó 2017 Elsevier B.V. All rights reserved.

ural hydrology of a site using decentralized micro-scale control
measures (Coffman, 2002). LID practices consist of bioretention
cells, infiltration wells/trenches, storm water wetlands, wet ponds,
level spreaders, permeable pavements, swales, green roofs, vegetated filter/buffer strips, sand and gravel filters, smaller culverts,
and water harvesting systems. Several studies have evaluated thebenefits of LIDs. For example, Newcomer et al. (2014) used anumerical model to demonstrate the benefits of LIDs, and an infiltration trench in particular, on recharge and local groundwater
resources for future climate scenarios. In another paper, Berardi
et al. (2014) demonstrated how green roofs may contribute to
the development of more sustainable buildings and cities. Green
Roofs (GR) were able to significantly reduce peak rates of storm
water runoff (Getter et al., 2007) and retain rainfall volumes withretention efficiencies ranging from 40% to 80% (Bengtsson et al.,2004). Permeable pavements offered great advantages in terms ofrunoff reduction (Carbone et al., 2014; Collins et al., 2008), water
retention, and water quality (Brattebo and Booth, 2003). Even
though the results of available studies are encouraging, more

264

G. Brunetti et al. / Journal of Hydrology 548 (2017) 263–277

research is needed to precisely assess the impact of LIDs on thehydrological cycle.
As pointed out by several authors (e.g., Elliot and Trowsdale,
2007; Wong et al., 2006), there is a strong demand for predictivemodels that can be applied across a range of locations and conditions to predict the general performance of a range of stormwater
treatment measures. In recent years, researchers have focused
their attention on applying and developing empirical, conceptual,and physically-based models for LIDs analysis. In their review article, Li and Babcock (2014) reported that there were >600 papers
published worldwide involving green roofs, with a significant portion of them related to modeling. Several studies demonstrated
that physically-based models can provide a rigorous descriptionof various relevant processes such as variably-saturated water
flow, evaporation and root water uptake, solute transport, heattransport, and carbon sequestration. Brunetti et al. (2016a,
2016b) used a mechanistic model, HYDRUS-3D (Šimůnek et al.,2016; Šimůnek et al., 2008), to analyze an extensive green roof in
a Mediterranean climate. The model, previously validated against
field scale measurements, was used to investigate the hydraulicresponse of a green roof to single precipitation events and its
hydrological behavior during a two-month period. Metselaar(2012) used the SWAP model (van Dam et al., 2008) to simulate
the one-dimensional water balance of a substrate layer on a flat
roof with plants. Li and Babcock (2015) used HYDRUS-2D to model
the hydrologic response of a pilot green roof system. The model
was calibrated using water content measurements obtained with
TDR (Time Domain Reflectometer) sensors. The calibrated model
was then used to simulate the potentially beneficial effects of irrigation management on the reduction of runoff volumes. The
VFSMOD model (Munoz-Carpena and Parsons, 2004) was extensively used for the analysis of the hydraulic behavior and solutetransport of vegetated filter strips (Abu-Zreig et al., 2001;
Dosskey et al., 2002).
However, physically-based modeling often involves highly nonlinear, partial, differential equations that are solved using various
numerical approximation methods, requiring a high computationalcost. Moreover, a comprehensive simulation framework includes
model calibration, sensitivity analysis, and uncertainty quantification aimed at enhancing confidence in the model and its ability todescribe real world systems. These tasks require running the simulation model hundreds or thousands of times and thus the computational cost exponentially increases.
Surrogate modeling focuses on developing and using a computationally inexpensive surrogate of the original model. The main aim
is to approximate the response of an original simulation model,
which is typically computationally intensive, for various quantities
of interest (Razavi et al., 2012). Surrogate models have been widely
applied in various water-related and environmental modelingproblems. Khu and Werner (2003) used artificial neural networks(ANN) in conjunction with genetic algorithms (GA) to reduce thecomputational budget required in the uncertainty quantification
framework of the rainfall-runoff model SWMM. The GA was firstused to identify the areas of higher importance in the parameter
space and ANNs were then used to approximate the response surface in these areas (Khu and Werner, 2003). Borgonovo et al. (2012)
tested a surrogate model for the estimation of the sensitivity
indices of an environmental model. Zhang et al. (2009) evaluated
ANN and Support Vector Machine (SVM) for approximating the Soiland Water Assessment Tool (SWAT) model in two watersheds.
Keating et al. (2010) used a surrogate model to carry out a comparison between the null-space Monte Carlo sampling (NSMC) and theDiffeRential Evolution Adaptive Metropolis (DREAM) algorithm forparameter estimation and uncertainty quantification. In anotherstudy, Laloy et al. (2013) used Polynomial Chaos Expansion (PCE)to emulate the output of a large-scale flow model. The surrogate

model was used in a Bayesian analysis framework to derive theposterior distribution of different parameters. In their study,
Younes et al. (2013) used a surrogate model to estimate three soil
hydraulic parameters from a drainage experiment. In particular,
PCE was used to run a Monte Carlo Markov Chain (MCMC) analysis.
However, although the widespread diffusion of surrogate modelingtools could drastically reduce computational budgets, their use forphysically-based modeling of LIDs is still unexploited.
The primary objective of this paper is to investigate the suitability of surrogate modeling for the numerical analysis of LIDs techniques by analyzing data from a real case study. The mechanistic
model HYDRUS-2D is first used to simulate the hydraulic behaviorof a Stormwater Filter (SF) at the University of Calabria, Italy. The
surrogate model, based on kriging, is then used to carry out a Global Sensitivity Analysis (GSA) and a Global Optimization of soil
hydraulic parameters. The use of a surrogate model for the sensitivity analysis of model outputs to soil hydraulic properties represents a new application of this technique that can provide asignificant contribution in this field.
The problem is addressed in the following way. First, the evaporation method is used to measure the soil hydraulic properties ofthe vegetated substrate above the gravel filter, for which thehydraulic properties were unknown. The measured soil hydraulicproperties of the vegetated substrate and the selected ranges ofparameters of the filter layer are then used in HYDRUS-2D to set
up the model. A Latin Hypercube Sampling (LHS) plan is used tobuild a first trial of the surrogate model. Before continuing withthe other tasks, the surrogate model is validated and improvedby using specific infill criteria. Once validated, the surrogate model
is first used for the GSA based on Sobol’s method to compute thesensitivity measures, and then for the inverse parameter estimation carried out using the Particle Swarm Optimization (PSO) algorithm. Finally, estimated parameters are used in the original
mechanistic model for the validation purpose.

2. Materials and methods
2.1. Stormwater filter and site description
The University of Calabria is located in the south of Italy, in thevicinity of Cosenza (39°180 N 16°150 E). The climate is Mediterranean with a mean annual temperature of 15.5 °C and averageannual precipitation of 881.2 mm. The stormwater filter (SF) has
a surface area of 125 m2, an average slope of 2%, and a total profile
depth of 0.75 m. Fig. 1 shows a schematic of the SF.
The filter layer is covered by a vegetated soil substrate with ameasured bulk density of 1.59 g/cm3. A high permeability geotextile with a fiber area weight of 60 g/m2 is placed at the interface
between the soil substrate and the filter layer to prevent fine particles from migrating into the underlying layer. The filter layer iscomposed of a gravelly material characterized by a high permeability. An impervious membrane is placed at the bottom of the profile
to prevent water from percolating into deeper horizons.
The SF is used to treat stormwater runoff from the adjoining
impervious parking lot, which is characterized by an area of
220 m2. Stormwater runoff from the parking lot is first conveyed
into a manhole and then to an instrumented channel where theflow rate is measured by a flux meter composed of a rectangular,
sharp crested weir coupled with a pressure transducer. The pressure transducer (Ge Druck PTX1830) measures the water level
inside the channel and has a range of measurements of 75 cm withan accuracy of 0.1% of the full scale. The pressure transducer wascalibrated in the laboratory using a hydrostatic water column, linking the electric current intensity with the water level inside theG. Brunetti et al. / Journal of Hydrology 548 (2017) 263–277

265

Fig. 1. A schematic of the experimental site (top) and a typical cross-section (bottom) of the stormwater filter.

column. An exponential head-discharge equation for the flux meter
was obtained by fitting the experimental data.
Measured runoff is next conveyed into a 14 m long, horizontal
perforated pipe where it is distributed on the top of the filter layer(Fig. 1). As shown in Fig. 1, the soil substrate is not used to treat
stormwater runoff, which is directly routed into the filter, but only
to increase the retention and evapotranspiration capacity of thesystem itself. The baseflow is collected in a horizontal drain, which
consists of a perforated PVC pipe, and is conducted to a manhole
for quantity and quality measurements. A second flux meter, composed of a PVC pipe with a sharp-crested weir and a pressure transducer, measures the flow rate. Runoff and baseflow data wereacquired with a time resolution of one minute and stored in a
SQL database. No measurements of pressure heads or volumetricwater contents inside of the filter were taken.
A weather station located directly at the site measures precipitation, wind velocity and direction, air humidity, air temperature,atmospheric pressure, and global solar radiation. Rain data aremeasured using a tipping bucket rain gauge with a resolution of
0.254 mm and an acquisition frequency of one minute. Climatic
data are acquired with a frequency of five minutes. Data are processed and stored in the SQL database.
Two month-long data sets were selected for the analysis (Fig. 2).
The first data set, which started on 2014-01-15 and ended on
2014-02-15, was used for obtaining the surrogate model. The second data set, which started on 2014-03-01 and ended on 2014-0331, was used for model validation. The precipitation totals for the
first and second data sets were 274 and 174 mm, respectively.
The second data set was selected because it had significantly different meteorological dynamics than during the first period. The optimization set is characterized by multiple rain events with few dry

periods. The validation set has fewer rain events, which are concentrated at the beginning and end of the time period and separated by a relatively long dry period.
Hourly reference evapotranspiration was calculated using thePenman-Monteith equation (Allen et al., 1998). Considering thatvegetation mainly consisted of herbaceous plants, an average valueof albedo of 0.23 was assumed in calculations of net short-wave
radiation (Breshears et al., 1997).
2.2. Evaporation method and parameter estimation
2.3.1. Evaporation method
Modeling of water flow in unsaturated soils by means of the
Richards equation requires knowledge of the water retention function, h(h), and the hydraulic conductivity function, K(h), for eachsoil layer of the SF, where h is the volumetric water content[L3L3], h is the pressure head [L], and K is the hydraulic conductivity [LT1]. In order to reduce the dimensionality of the optimization
problem, the soil hydraulic properties of the soil substrate were
measured in the laboratory using a simplified evaporation methodwith an extended measurement range (down to 9000 cm), as
proposed by Schindler et al. (2010a, 2010b). For a detailed description of the modified evaporation method, please refer to Schindler
et al. (2010a, 2010b).
Peters and Durner (2008) conducted a comprehensive error
analysis of the simplified evaporation method and concluded that
it is a fast, accurate, and reliable method to determine soil hydraulic properties in the measured pressure head range, and that thelinearization hypothesis introduced by Schindler (1980) causes
only small errors. The above cited method has already been usedin the LIDs analysis for the determination of the unsaturated soil

266

G. Brunetti et al. / Journal of Hydrology 548 (2017) 263–277

Fig. 2. Precipitation (black line) and subsurface flow (grey line) for the optimization (top) and validation (bottom) periods, respectively.

hydraulic properties of a green roof substrate (Brunetti et al.,2016a). In that study, the measured soil hydraulic properties wereused in HYDRUS-3D to simulate the hydraulic behavior of a green
roof and validated by providing optimal correspondence between
simulated and measured outflows. The simplified evaporation
method was similarly used in this study for the determination of
the unsaturated hydraulic properties of the soil substrate. For acomplete description of the system, please refer to UMS GmbH(2015).
The soil for the laboratory analysis was directly sampled fromthe SF using a stainless-steel sampling ring with a volume of
250 ml. The soil sample was saturated from the bottom before
starting the evaporation test. The measurement unit and tensiometers were degassed using a vacuum pump, in order to reduce thepotential nucleation sites in the demineralized water. Since
Peters and Durner (2008) suggested a reading interval for structured soils of less than 0.1 day, the reading interval was set to
20 min in order to have high resolution measurements. At theend of the experiment, the sample was placed in an oven at
105 °C for 24 h, and then the dry weight was measured.

2.4. Parameter estimation
The numerical optimization procedure, HYPROP-FIT (Pertassek
et al., 2015), was used to simultaneously fit retention and hydraulicconductivity functions to experimental data obtained using theevaporation method. Fitting was accomplished using a non-linear
optimization algorithm that minimizes the sum of weighted
squared residuals between model predictions and measurements.
The software uses the Shuffled Complex Evolution (SCE) algorithm
proposed by Duan et al. (1992), which is a global parameter estimation algorithm. The goodness-of-fit was evaluated in terms of
the Root Mean Square Error (RMSE), while the Akaike information
criterion (AIC) (Hu, 1987) was used to choose between different
hydraulic conductivity functions. The software also provides 95%

confidence intervals to assess the uncertainty in parameter estimation. In order to calculate the parameter uncertainties a linear
approximation of the covariance matrix for each estimated parameter is calculated. The confidence interval for the i-th parameter isthen computed by combining the covariance matrix and the upper
a/2 quantile of the Students t-distribution, where a is set to 0.05 forthe computation of the 95% confidence intervals.
2.5. Modeling theory
2.5.1. Water flow and root water uptake
The HYDRUS-2D software (Šimůnek et al., 2008) was used tomodel the hydraulic behavior of the SF. HYDRUS-2D is a twodimensional model for simulating the movement of water, heat,and multiple solutes in variably-saturated porous media.
HYDRUS-2D numerically solves the Richards equation for multidimensional unsaturated flow:@h¼ O½k  Oðh  zÞ  S@tð1Þ

where t is time (T), z is the vertical coordinate (L), and S is a sinkterm (L3L3T1), defined as a volume of water removed from a unit
volume of soil per unit of time due to plant water uptake. The unimodal van Genuchten–Mualem (VGM) model (van Genuchten,1980) was used to describe the soil hydraulic properties of thetwo layers:

(

H¼

H¼

if h 6 0

1

if h > 0ð2Þh  hrhs  hr(

K¼

1m
ð1þðajhjÞn Þ1

m 2

K s HL ½ð1  ð1  Hm ÞÞ 

if h < 0Ks

if h > 0ð3Þ

G. Brunetti et al. / Journal of Hydrology 548 (2017) 263–277

m¼1

1
n

where H is the effective saturation (), a is a shape parameter
related to the inverse of the air-entry pressure head (L1), hs andhr are the saturated and residual water contents, respectively (),n and m are pore-size distribution indices (), Ks is the saturated
hydraulic conductivity (LT1), and L is the tortuosity and poreconnectivity parameter ().
While the soil hydraulic properties of the soil substrate were
determined using the simplified evaporation method, those of
the filter were optimized in the surrogate analysis framework.
However, not all parameters were included in the optimization
process. The residual water content hr was fixed to 0, consideringthat the filter is composed of coarse gravel, and the tortuosity L
was set to 0.5, which is a common value in the literature. The initial
range of the investigated parameters is reported in Table 1.
Feddes et al. (1978) defined S as:SðhÞ ¼ aðhÞ  Sp

ð4Þ

where a(h) is a dimensionless water stress response function that
depends on the soil pressure head h and has a range of values
between 0 and 1, and Sp is the potential root water uptake rate.
Feddes et al. (1978) proposed a water stress response function inwhich water uptake is assumed to be zero close to soil saturation(h1) and for pressure heads larger (in absolute values) than the wilting point (h5). Water uptake is assumed to be optimal between two
specific pressure heads (h2, h3 or h4), which depend on a particularplant. At high potential transpiration rates (5 mm/day in the model
simulation) stomata start closing at lower pressure heads (h3) (inabsolute value) than at low potential transpiration rates (1 mm/d)(h4). Parameters of the stress response function for a majority ofagricultural crops can be found in various databases (e.g., Taylorand Ashcroft, 1972; Wesseling et al., 1991). Considering that thevegetation cover was mainly constituted of herbaceous plants,
parameters reported for grass in Wesseling et al. (1991) were usedin this study.
The local potential root water uptake Sp was calculated from thepotential transpiration rate Tp. Beer’s equation was first used topartition reference evapotranspiration, calculated using thePenman-Monteith equation (Allen et al., 1998), into potential transpiration and potential soil evaporation fluxes (e.g., Ritchie, 1972).
The partitioning of evapotranspiration into potential transpirationand potential evaporation allows the computation of different
actual fluxes in the soil-vegetation system. The Leaf Area Index(LAI) is needed to partition evaporation and transpiration fluxes.
In this study, a LAI value of 2.29 as reported by Blanusa et al.
(2013) for a sedum mix was used, considering its similarity with
the installed vegetation. For a detailed explanation of evapotranspiration partitioning, please refer to Sutanto et al. (2012).
HYDRUS-2D allows for the consideration of a spatially variable
root distribution. In this study, a homogeneous root zone within adepth of 15 cm was defined. The root density was assumed to beuniform inside the root zone and zero in the remaining part of
the numerical domain. The total potential transpiration flux froma transport domain is equal in HYDRUS to potential transpiration
Tp, multiplied by the surface length associated with vegetation.
Table 1
Ranges of investigated parameters for the surrogate-based analysis.
ParameterRange

hs []
a [1/cm]
n1 []
Ks1 [cm/min]

0.1–0.3
0.001–0.3
3.0–7.0
30.0–100.0267This total potential transpiration flux is then distributed over theentire root zone for the computation of the actual root water
uptake.

2.6. Numerical domain and boundary conditions
The two-dimensional domain had a length of 8.0 m and a depthof 0.75 m. The geotextile was not included in the model considering its negligible thickness, its limited hydraulic effect due to its
high permeability, and that its sole function was to separate thesoil substrate from the filter layer. The domain was discretized into
two-dimensional triangular elements using the MESHGEN tool of
HYDRUS-2D. The mesh was refined in the right part of the domain,
where the effect of the surface runoff from the parking lot was simulated. This refinement was necessary in order to numericallyaccommodate the significant pressure head gradients generated
by infiltration of runoff, and thus to reduce the mass balanceerrors. The generated FE mesh had 736 nodes and 1350 twodimensional elements. The quality of the FE mesh was assessedby checking the mass balance error reported by HYDRUS-2D atthe end of the simulation. Mass balance errors, which in this simulation were always below 1%, are generally considered acceptable
at these low levels.
The surface of the SF was exposed to precipitation, evapotranspiration, and surface runoff from the impervious parking lot. As
a result, in HYDRUS, two different boundary conditions were specified at the top of the modeled domain, as well as at its bottom(Fig. 3).
The ‘‘Atmospheric” boundary condition, which was assigned on
the surface of the soil substrate (green line in Fig. 3), can exist in
three different states: (a) precipitation and/or potential evaporation fluxes, (b) a zero pressure head (full saturation) during ponding when both infiltration and surface runoff occurs, and (c) an
equilibrium between the soil surface pressure head and the atmospheric water vapor pressure head when atmospheric evaporative
demand cannot be met by the substrate. The threshold pressure
head, which was set to 10,000 cm, divides the evaporation process from the soil surface into two stages: (1) a constant rate stage
when actual evaporation, equal to potential evaporation, is limited
only by the supply of energy to the surface, and (2) the falling rate
stage, when water movement to the evaporating sites near the surface is controlled by subsurface soil moisture and the soil hydraulicproperties. In such conditions, actual evaporation, calculated as aresult of the numerical solution of the Richards equation, is smaller
than potential evaporation.
The ‘‘Variable Flux” boundary condition, which included bothprecipitation and measured surface runoff, was used in the area
under the perforated pipe (red line in Fig. 3). Evaporation was
excluded since most of the surface was covered by the perforated
pipe, which reduced the exposure of the surface to wind and solar
radiation.
A seepage face boundary condition (brown line in Fig. 3) was
specified at the bottom left corner of the numerical domain to simulate the effect of the horizontal drain. A seepage face boundary
acts as a zero pressure head boundary when the boundary nodeis saturated and as a no-flux boundary when it is unsaturated. A
zero flux boundary condition (black line in Fig. 3) was applied to
all remaining boundaries of the domain to simulate the effect of
the impervious membrane placed at the bottom and on the sidesof the SF.
The initial conditions were specified in terms of the soil waterpressure head and were set to linearly increase with depth, from90 cm at the top of the flow domain (z = 0) to 0.5 cm at the bottom (z = 75). The surface layers are assumed to be drier than thebottom layers since they are directly exposed to the atmosphere.

268

G. Brunetti et al. / Journal of Hydrology 548 (2017) 263–277

Fig. 3. The spatial distribution of applied boundary conditions.

The numerical model is expected to be sensitive to the initial conditions only during the first few simulated days.

2.7. Surrogate based model
2.7.1. Kriging
There are two broad families of surrogate modeling techniques:
Response Surface Surrogates (RSS) and Lower-Fidelity Surrogates(LFS). While RSS are data-driven techniques for approximating
the response surface of high-fidelity (original) models based on alimited number of original model evaluations, LFS are essentiallycheaper-to-run, alternative simulation models with different levelsof accuracy (Razavi et al., 2012). As pointed out by Razavi et al.
(2012) in their review paper, LFS outperforms RSS when thedimensionality of the problem is high and the response surface
landscape is characterized by multimodality. In such circumstances, RSS would need a higher number of original model runs
to correctly approximate the response surface. O’Hagan (2006)
highlighted how the same number of design sites can lead to different parameters space coverages depending on the dimensionality of the problem. However, when the problem is low-dimensional
and the response surface is characterized by a low or moderate
multimodality, RSS are preferred since a limited number of highfidelity model runs is required to build a reliable surrogate model.
The present study involves a SF model with four parameters to
be investigated. Considering the low dimensionality of the problem, the kriging response surface approximation technique wasused. Unlike other RSS, kriging models have their origins in mining
and geostatistical applications involving spatially and temporally
correlated data. The kriging technique has also been referred toin the literature as a Gaussian Process (GP) prediction(Rasmussen and Williams, 2006; Sacks et al., 1989). A GaussianProcess is formally defined as being a probability distribution overa (possibly infinite) number of variables, such that the distribution
over any finite subset of them is a multi-variate Gaussian. As theGaussian distribution is fully specified by its mean and covariancematrix, the GP is specified by a mean and a covariance function(Mackay, 1998). The mean is usually assumed to be zero, and, insuch circumstances, the covariance function completely describes
the GP behavior. One of the most attractive features of GP is thatit treats the deterministic response of a computer model as therealization of a stochastic process, in particular a Gaussian random
process, thereby providing a statistical basis for fitting. This capability provides a first approximation of uncertainty associated witheach value predicted by the surrogate. Another advantage of kriging against other RSS techniques such as Artificial Neural Networks(ANN) or Support Vector Machines (SVM) is that kriging is an exact
emulator. An exact emulator precisely predicts all design sites used
to build the surrogate, while inexact emulators can introduce bias
in such sites. As described by Razavi et al. (2012), an exact emulation is recommended for approximating the deterministic responseof computer simulation models. Inexact emulators have smoothing
capabilities that can help when the response surface is noisy (e.g.,physical experiments), however this feature can lead to poor

approximation of the response surface when it is characterizedby multiple local minima.
The kriging model is a combination of a polynomial model and alocalized deviation model, which is based on a spatial correlationof samples (Eq. (4)):

yðxÞ ¼ f ðxÞ þ ZðxÞð5Þ

where y(x) is an unknown function of interest, f(x) is an approximation function, and Z(x) is the realization of a stochastic process withzero mean, the variance r2, and nonzero covariance. While f(x)globally approximates the response surface through design sites, Z(x) creates localized deviations. The covariance matrix of Z(x) is
given by Eq. (6):Cov ½Zðxi Þ; Zðx j Þ ¼ r2 Wð½Rðxi ; x j ÞÞð6Þ

where W is the p x p symmetric correlation matrix and R(x , x ) is thecorrelation function between two of the p sampled data points. R(xi,xj) can assume different forms and is specified by the user. In thisstudy, the Gaussian correlation function has been used (Eq. (7)):
iRðx ; x Þ ¼ exp 
i

j

N
X

s

j

!

i
k 2
j jxk  xk j

ð7Þ

k¼1

where N is the number of parameters, sj are the unknown correlation parameters used to fit the model, and xik and xjk are the kth components of the sample points xi and xj. Correlation parameters sj are
estimated using the maximum likelihood methodology. The ‘‘best”
kriging model is found by solving a j-dimensional, unconstrained,
nonlinear optimization problem. In this study, the PSO global optimization algorithm has been used to identify kriging parameters.
2.8. Design of experiments
The first step in the generation of a surrogate model is to sample
the response surface at some specific design sites. This procedure isusually referred to in the literature as the Design of Experiments(DoEs). As pointed out by Razavi et al. (2012), a sufficiently large
and well-distributed set of initial design sites is crucial for a successful application of a metamodeling framework. There are several DoEs methods available in the literature. Factorial design(Gutmann, 2001), Latin Hypercube Sampling (LHS) (McKay et al.,1979), and Symmetric Latin Hypercube Sampling (SLHS) (Yeet al., 2000) are the most commonly used. In this study, the LHS
has been used. The size of the DoEs sample is strongly dependent
on the complexity of the original response surface and computational budget available. The kriging model requires at least N+1
design sites to fit, while additional sites will improve the accuracyof the surrogate. Several relations were proposed in the literature
to choose the size p of the initial sample (e.g., Gutmann, 2001;
Regis and Shoemaker, 2004). In this study, the relation proposed
by Jones et al. (1998) has been used (Eq. (8)):

p ¼ 10N

ð8Þ

Considering that the number of investigated parameters was 4,
40 sampling points were generated using the LHS.

G. Brunetti et al. / Journal of Hydrology 548 (2017) 263–277

2.9. Approximation uncertainty framework
The surrogate model fitted on a DoEs sample is only the first
approximation of the original response function. Its accuracy canbe improved using further original model runs (infill points) inaddition to the initial sampling plan. The distribution of infill points
strongly depends on the complexity of the original response surface and on the type of analysis conducted. When the purposesof the analysis of the surrogate model includes uncertainty andsensitivity analyses, it is necessary to have an accurate global
approximation of the response surface. The adaptive-recursive
framework relies on the assumption that the optimal solution of
the surrogate represents well the original model. However, this is
not always true, especially when the response surface is characterized by high multimodality. When the approximation of the
response surface by the surrogate model is limited, this framework
is not recommended for the uncertainty and/or sensitivity analysis.
Jones (2001) concluded that the adaptive-framework is helpful atbest for local optimization.
In this study, the approximation uncertainty framework wasapplied to address the shortcoming of the adaptive-recursive
framework, which considers uncertainties associated with theapproximation. RSS, such as kriging, explicitly provide a measureof uncertainty since they treat the deterministic response of a computer simulation as the realization of a stochastic process. The
model can then be evaluated at points with the highest uncertainty, which can then be included among the design sites.
Although globally convergent, such an approach requires animpractically large number of original function evaluations. An
effective uncertainty based framework should balance exploitation(i.e., fine tuning of a good solution) and exploration (i.e., reducing
the overall uncertainty of the surrogate). The expectedimprovement approach (Schonlau, 1997) was used in this study.
An expected-improvement is a measure that statistically quantifies
the obtained improvement when a given point is evaluated by theoriginal model and added to design sites. For a complete description of the expected-improvement approach please refer to
Schonlau (1997).
In this study, the expected-improvement approach was used to
add 15 infill points to the initial design sites. In order to have good
accuracy, the surrogate was refitted after each new original model
evaluation (Razavi et al., 2012).
2.10. Surrogate validation
Validation of the RSS model is important for evaluating the reliability of the surrogate. Although exact emulators such as kriging
exactly interpolate the response surface at design sites, their accuracy in unexplored regions of the parameter space must be evaluated. The validation can be conducted by evaluating theagreement between values of the variable of interest predicted byboth the surrogate and original models on an independent set ofsample points. Cross validation strategies such as k-fold andleave-one-out cross validation have also been used in the literature(Wang and Shan, 2007). In the present study, an independent set ofsample points was generated using the LHS and used to validate themodel. Ten points were used to carry out the validation process. The
Pearson coefficient R2 was used to assess the agreement between
predicted and modeled values. As suggested by Forrester et al.
(2008), a value of the correlation coefficient higher than 0.8 indicates a surrogate model with good predictive capabilities.
2.11. Global Sensitivity Analysis (GSA)
A sensitivity analysis (SA) can identify the most influential
parameters, their interactions, and how these parameters affect269

the output (Saltelli et al., 2005). Most SAs performed in the literature of environmental sciences are the so-called ‘one-at-a-time’(OAT) sensitivity analyses, performed by changing the value ofparameters one-at-a-time while keeping the other parametersconstant (Cheviron and Coquet, 2009; Houska et al., 2013; Rezaei
et al., 2015). However, when the model includes interactions
between multiple parameters, results of the OAT analysis are inaccurate because parameter interactions can be globally identified
only by simultaneously changing multiple parameters. For this reason, when the property of a model is a priori unknown, a Global
Sensitivity Analysis (GSA) is always preferred (Saltelli andAnnoni, 2010). Practitioners call this analysis a model-free setting,which means that a particular application does not depend on particular assumptions regarding the behavior of the model, such as
linearity, monotonicity, etc (Saltelli and Annoni, 2010).
Variance-based methods aim to quantify the amount of variance that each parameter contributes to the unconditional variance of the model output. In Sobol’s method, these measures are
represented by Sobol’s sensitivity indices (SIs). These indices givequantitative information about the variance associated with a single parameter or with interactions of multiple parameters. For a
more complete explanation about Sobol’s method, please refer toSobol’ (2001). Sobol’s sensitivity indices are expressed as follows:FirstOrder Si ¼

Vi
V
V ij
V

ð10Þ

X
Sij þ :::

ð11ÞSecondOrder Sij ¼
Total ST ¼ Si þ

ð9Þ

j–i

where Vi is the variance associated with the ith parameter and V isthe total variance. The first-order index, Si, is denoted in the literature as the ‘‘main effect.” When the model is additive, i.e., when it
does not include interactions between input factors, the firstorder index is sufficient for decomposing the model’s variance.
For additive models, the following relation is valid:X
Si ¼ 1

ð12Þ

i

On the other hand, the total effect index, ST, gives information
about a non-additive part of the model. STi = 0 is a condition necessary and sufficient for Xi to be non-influential. For an accurate
description of the calculation of Sobol’s indices please refer to
Saltelli et al. (2010).
When the model is nonlinear, as most environmental models
are, Sobol’s indices are calculated using Monte Carlo integrals.
Obviously, the accuracy in the estimation of integrals becomes
more accurate as the number of samples increases, which also
increases the computational cost of the SA. However, this limitation is avoided when using a surrogate model since the computational cost associated with the evaluation of a large number ofsamples is very low (O’Hagan, 2006; Oakley and O’Hagan, 2004).
For this reason, 1000 samples for a total of 30,000 surrogate modelruns were used in this study. To sample the parameters’ space weused Sobol’s quasi-random sampling technique (Sobol’, 2001).
In order to assess the accuracy of the estimations of the sensitivity indices, the bootstrap confidence intervals (BCIs) (Efron andTibshirani, 1986) were estimated. The basic idea of the bootstrapping is that the sample contains all available information about theunderlying distribution. In our particular case, we were interested
in computing the uncertainty of estimated sensitivity indices.
However, since their distribution is unknown it is not possible tocompute the confidence intervals analytically. The rationale ofthe bootstrap method is to replace the unknown distribution with

270

G. Brunetti et al. / Journal of Hydrology 548 (2017) 263–277

its empirical distribution and to compute the sensitivity indicesusing a Monte Carlo simulation approach where samples are generated by resampling the original sample used for the sensitivity
analysis. In our case, the samples used for the GSA were sampled
1000 times with replacement, whereby Sobol’s indices were calculated for each resampling. In this way, 95% confidence intervals areconstructed using the percentile method and the moment method(Archer et al., 1997).
2.12. Particle Swarm Optimization
Numerous applications of inverse modeling for the estimationof soil hydraulic properties exist in the literature (Abbaspour
et al., 2004; Hopmans et al., 2002; Vrugt et al., 2008; Vrugt et al.,2004). The gradient methods (Marquardt, 1963) have been mostwidely used among hydrologists and soil scientists. However, thesemethods are sensitive to the initial values of optimized parameters
and the algorithm often remains trapped in local minima, especially when the response surface exhibits a multimodal behavior.
These considerations inspired researchers to develop and use global optimization techniques such as the annealing-simplex method(Pan and Wu, 1998), genetic algorithms (Ines and Droogers, 2002),
shuffled complex methods (Vrugt et al., 2003), and ant-colony optimization (Abbaspour et al., 2001), among many others.
In this paper, a global search method based on Particle Swarm
Optimization (PSO) (Kennedy and Eberhart, 1995) was used. Thismethod simulates the behavior of a flock of birds collectively foraging for food (i.e., searching for the optimum of the objective function). In their recent study, Brunetti et al. (2016b) used PSO to
estimate the soil hydraulic properties of a permeable pavementwith satisfactory results. PSO is a relatively new algorithm for evolutionary computation methodology, but its performance has proven to be comparable to various other more establishedmethodologies (Kennedy and Spears, 1998; Shi et al., 1999). Oneof the main advantages of PSO is the easiness of its implementation(Liang et al., 2006). A detailed description of the PSO algorithm isgiven in Shi and Eberhart (1998).
In PSO, each particle represents a possible solution of the problem. First, particles are placed in the search space, each characterized by a particular value of the objective function. Each particle
then changes its position after exchanging information about its
own current and best positions with other members of the swarm.
The next iteration starts after all particles have changed their positions. The most important parameters in the PSO are c1, c2, and w.
c1 and c2 are constant parameters known as the cognitive andsocial parameters, respectively, which drive the search behaviorof the algorithm. Depending on the values of c1 and c2 the PSO
can be more or less ‘‘responsive”. However, their values should
be selected carefully because large values of these parameters
can lead to instabilities in the algorithm. w is the inertia-weight,
which plays a key role in the optimization process by providingbalance between exploration and exploitation. In PSO, each particle
is influenced by its nearest neighbors. The arrangement of neighbors that influence a particle is called the topology of the swarm.
Different types of neighborhoods are reported in the literature(Akat and Gazi, 2008). In this study, the all topology is used, inwhich the neighborhood encompasses the entire swarm. The PSO
parameters used in this study for both scenarios are reported inTable 2 and are as suggested by Pedersen (2010).

2.13. Objective function
The Nash-Sutcliffe efficiency (NSE) index (Nash and Sutcliffe,1970) was used to evaluate the agreement between measured
and modeled hydrographs and as the variable of interest in the surrogate analysis:

2 XT

3
2
ðQ obs Q mod
Þ
i
i
i¼1
4
5
NSE ¼ 1  XT
2
obs
obs
ðQ

Q
Þ
i
mean
i¼1

ð13Þ

where Qobs
is the ith measured value, Qmod
is the ith simulatediivalue, and Qobs
mean is the mean value of observed data. The NSE index
ranges between 1 and 1.0, is equal to 1 in case of perfect agreement, and generally, values between 0.0 and 1.0 are considered
acceptable (Moriasi et al., 2007). The NSE index was used because
it is often reported to be a valid measure for evaluating the overall
fit of a hydrograph (Sevat et al., 1991).
It is important to emphasize that subsurface outflow from a LID
system is among the most important outputs in the analysis ofurban drainage systems. In our case, the stormwater filter was
impervious at the bottom and hydraulically connected with thesewer system. An accurate numerical reconstruction of the subsurface hydrograph is thus fundamental in order to quantify its effecton the drainage system.
3. Results and discussion
3.1. Evaporation method
Soil hydraulic properties measured using the evaporation
method are displayed in Fig. 4. The retention data point close tolog (|h|) = 4 (h in cm) was obtained by using the air-entry pressure
head of the ceramic cup of the tensiometers. Measured retention
points were not available in the very dry range, between 2.7 and
3.8, since cavitation occurred in the tensiometers. The behaviorof the retention curve appears to be sigmoidal and characterized
by a clearly identifiable air-entry pressure head at log (|h|) = 2 (hin cm). Below this pressure head, the soil quickly desaturated, indicating a narrow pore-size distribution. Measured points of thehydraulic conductivity function were sparser and concentrated in
the dry range between volumetric water contents of 0.10 and
0.20. The measured soil porosity and bulk density were 0.44 and1.59 g/cm3, respectively.
The unimodal van Genuchten-Mualem model (van Genuchten,1980) was fitted to measured points using the HYPROP-FIT software. The RMSE values for retention and conductivity functions
were 0.04 (cm3cm3) and 0.6 (in log K, cm/day), respectively. The
VGM function (full lines in Fig. 4) described the retention data well,especially in the dry and medium-wet regions (volumetric watercontents of 0.5–0.3), while it introduced some bias near saturation
where it poorly described the sharp increase in water retention at
the air-entry pressure head. However, the low RMSE value wasconsidered acceptable for the purposes of the present study. The
RMSE value for the hydraulic conductivity was higher, indicatinga slightly worse performance of the VGM function in describing
the hydraulic conductivity of the substrate. The estimated soil
hydraulic parameters (reported with their confidence intervals in
Table 3) were used in HYDRUS-2D to model the hydraulic behaviorof the soil substrate.
3.2. Kriging approximation of the response surface

Table 2
Parameters used in the PSO optimization.
Swarm size

c1

c2

w

63

0.73

2.02

0.36

The DoEs sample, generated with the LHS technique, was used
to build the first approximation of the response surface for theinvestigated soil hydraulic parameters. First, the HYDRUS-2D

G. Brunetti et al. / Journal of Hydrology 548 (2017) 263–277

271

Fig. 4. Measured values and modeled functions of soil water retention, h (log10(|h|)) (left) and the unsaturated hydraulic conductivity, K(log10(|h|)) (center) and K(h) (right).
Symbols represent the measured values, and full lines the fitted VGM functions.

Table 3
Estimated soil hydraulic parameters and their confidence intervals for the bimodal
hydraulic function.
Parameter

2.5%

Estimated value

97.5%

hr ()
hs ()a (1/cm)n ()Ks (cm/day)
L ()

0.02
0.43
0.021
1.84
200
0.63

0.03
0.44
0.025
1.97
260
0.44

0.04
0.45
0.029
2.10
320
0.25

model was executed 40 times (Eq. (8)), and the NSE index wascomputed for each run and stored in a 1D array. A single run ofthe original HYDRUS-2D model required almost 1 min of CPU time
on a laptop equipped with a CPU IntelÒ Core i7-4700 MQ 2.40 GHz
processor and 8 GB of RAM. Next, the LHS sample and the NSE array
were used in the PSO optimization framework to estimate the kriging parameters. To check its accuracy, the obtained kriging model
was validated on another independent sample generated with
the LHS. As shown in Fig. 5, the validation sample covered valuesof NSE ranging from 0.2 to 0.8, providing information about the
response surface for both less and more accurate portions of theparameters space.
At the first inspection, the kriging model based on the initialsample exhibited a moderate accuracy. The determination coefficient R2 for the initial kriging model was 0.91, which already indicated an overall accuracy of the surrogate model (Forrester et al.,2008). This confirmed the good coverage of the DoEs sample. However, as shown in Fig. 5, the surrogate, while being highly accurate
for low values of the objective function, introduced a significant
bias for high values of NSE, which are those of interest in an optimization framework. The regression line for the initial surrogate (adashed grey line in Fig. 5) almost overlapped the bisector (a blackline in Fig. 5), which indicates a perfect agreement between thesurrogate and the original model in the range of 0.2–0.4, while it
underestimated values of the response surface in the region
around 0.8. The underestimation of the response surface values
in the region where optimal parameter values are likely locatedcould influence the next surrogate-based optimization of soil
hydraulic properties.
To increase the accuracy of the kriging model, the approximation uncertainty framework was used next. As described in themethodology section, 15 infilled points were added to the initialdesign sites using the expected improvement approach. As shownin Fig. 5 and as expected, the infilled kriging model outperformed

Fig. 5. Comparison between the HYDRUS-2D and kriging-predicted values of the
NSE for the validation sample. The initial (grey diamonds) and infilled (red circles)
kriging models are compared. A bisector (a black line) and regression lines for the
initial (a dashed grey line) and infilled (a dashed red line) kriging models arereported.

the initial kriging model. The determination coefficient R2
increased to 0.98, indicating that the infilled kriging improved
the description of the response surface. This behavior was confirmed by the regression line (a dashed red line in Fig. 5), whichalmost overlapped a bisector line. Moreover, the accuracy of the
kriging model improved for high values of NSE, which are thoseof interest in an optimization process, while remaining similarlyhigh for low values of NSE. This global accuracy of the surrogate
is fundamental for the GSA, which explores the response surface
landscape. However, it should be noted that even when the surrogate accuracy improves after applying the approximation uncertainty framework, a certain degree of bias remains. A possible
solution to further investigate this bias would be to increase thevalidation sample size. However, there are no clear indications in
the literature regarding the size of the validation sample in similar
problems (i.e. optimization and sensitivity analysis of soil hydraulic properties), making the choice rather subjective. Moreover,272

G. Brunetti et al. / Journal of Hydrology 548 (2017) 263–277

increasing the validation sample size would significantly increase
the computational cost of the analysis. In our particular case, inview of all the possible sources of uncertainty that affect the analysis, the validation of the surrogate-based optimized parametersagainst an independent set of data, performed towards to end of
the study, was specifically intended to provide a final check on
the accuracy of the estimated parameters. Considering the relatively high accuracy of the surrogate, no additional points werethus added to the design sites considering the relatively high accuracy of the surrogate, and the final surrogate was used for the GSA
and optimization.

3.3. Global Sensitivity Analysis
The validated kriging model was next used in the GSA. Sobol’s
sensitivity indices, with their confidence intervals for each parameter, are reported in Table 4.
The sensitivity analysis revealed that the model was additive.
This was confirmed by the sum of the first-order indices, which
was almost 1, and by negligible differences between the firstorder and total-effect indices for each parameter. An additivemodel Y = f(X1, X2,. . ., XN) can be decomposed into a sum of N functions, where N is the number of parameters. This means that theeffects of interactions between model parameters on model results
were negligible.
As shown in Table 4, the most influential parameter was theshape parameter a. Its first-order (S1) and total effect (ST) indiceswere more than an order of magnitude higher than corresponding
indices for the second most influential parameter Ks. The pore-size
distribution index n and the saturated water content hs exhibitedthe lowest sensitivity, indicating their marginal role on the output’s variance. Moreover, since their total effects were almost zero,these parameters can be fixed to any feasible value in the parameter space without affecting the value of the objective function,
reducing the dimensionality of the inverse problem to only two
parameters, a and Ks. Such results are very useful in an optimization framework, since they can simplify the parameter estimation
procedure. Some of the total effect indices were only slightly larger
than the first order indices. This is mainly due to approximation in
the numerical integration of the total unconditional variance(Sobol’, 2001).
It must be emphasized that additivity is quite unusual for environmental models, which are generally characterized by high nonlinearity and interactions between parameters (Brunetti et al.,2016b; Nossent et al., 2011). In our particular case, the dominant
effect of the a parameter on the hydrograph makes it difficult to
identify interactions between parameters. This may also indicatethat water flow in the filter layer deviates from the traditional Darcian behavior and involves other physical processes, such as preferential and film flows. In such circumstances, some parameters
can exhibit a negligible effect on the model’s response. This may
indicate that the proposed modeling approach could be affectedby model discrepancy, a concept introduced by Kennedy andO’Hagan (2001) to assess sources of uncertainty due to underlying
missing physics, numerical approximations, and other inaccuracies
Table 4
The first-order (S1) and total (ST) effect indices (in a decreasing order) with theirbootstrap confidence intervals (BCI) for the soil hydraulic parameters.
Parameter

S1

S1 (BCI)ST

ST (BCI)a [1/cm]

0.93
0.04
0.01
0.002
1.0

0.24
0.08
0.02
0.01

0.94
0.05
0.008
0.001
1.0

0.05
0.007
0.001
0.001

Ks [cm/min]n []hs []
Sumof the model. When the model discrepancy is not directly accounted
for, the model parameters can be treated as simple tuning parameters, which often act as a simplified surrogate for some more complex process that is not modeled in the simulator (Brynjarsdottir
and O’Hagan, 2014). In our case, it is plausible that the shapeparameter a acted as a tuning parameter, condensing and representing more complex phenomena, namely preferential and film
flows. It must be emphasized that the surrogate-based sensitivity
analysis can be extremely useful in detecting different sources of
uncertainties. Future researches should explicitly account formodel discrepancy in the analysis (Brynjarsdottir and O’Hagan,2014) to examine its impact on the model’s behavior, or performthe same type of analysis using more complex models which
account for preferential flows in the filter layer.
As mentioned in the methodology section, the GSA required
30,000 evaluations of the surrogate model. The computational costof the kriging-based sensitivity analysis was limited to 1–2 s on alaptop equipped with a CPU IntelÒ Core i7-4700 MQ 2.40 GHz processor and 8 GB of RAM. On the other hand, since a single HYDRUS2D model run required approximately 1 min, the same type of GSA
performed using the original HYDRUS-2D model would haverequired approximately 21 days of continuous computation. This
clearly represents one of the main advantages of surrogate-basedmodeling: performing the same type of analysis with negligible
computation time and a similarly good level of accuracy.
3.4. Kriging-Based optimization
Using the results of the GSA, only shape parameter a and thesaturated hydraulic conductivity Ks were optimized. The saturated
water content hs and the pore-size distribution index n wereassumed to be 0.15 and 3.2, respectively, considering that the filter
layer consisted of coarse gravel, usually characterized by reduced
porosity and narrow pore-size distribution. The soil hydraulicparameters of the filter layer, including the optimized parameters,
are summarized in Table 5. The filter layer exhibited both a highvalue of the saturated hydraulic conductivity Ks (90 cm/min) anda very low value of the shape parameter a (0.001 1/cm).
The estimated parameters indicated that the hydraulic behaviorof the filter layer was characterized by high flow rates and negligible retention capacity, which are both typical for coarse textured
media. The optimized parameter values are similar to thosereported in Brunetti et al. (2016b) for the base layer of a permeable
pavement. In that study, the base layer consisted of crushed stones
and was modeled using either the classical VGM function or thedual-porosity approach to account for preferential flow. Specifically, for the unimodal VGM function, the authors reported a valueof 0.023 1/cm for the shape parameter a and a saturated hydraulicconductivity Ks of 68.7 cm/min. Moreover, the plausible occurrenceof film flow in the filter layer, which can support very high flow
rates, especially at near-zero matrix potential (Tokunaga, 2009),
needs to be contemplated. Under such circumstances, the hydraulic behavior of the material tends to deviate from the typical
Richard’s type flow, and the optimized parameters attempt toapproximate a combination of fingering and film flow that likely
occur in this layer.
Fig. 6 shows a comparison between the measured and modeledhydrographs for the optimization period. The PSO resulted in the
NSE value of 0.85, which confirmed the accuracy of the measured
and estimated parameters. As reported by Moriasi et al. (2007), values of NSE between 0.75 and 1.0 indicate a very good agreement
between hydrographs, and an adequate model calibration. The
model was able to correctly reproduce the fast hydraulic responseof the hydrograph during precipitations and to reasonably estimate
peak flows. The insert of Fig. 6 shows the simulated against measured SF outflows. The same plot also shows a bisector line, which

273

G. Brunetti et al. / Journal of Hydrology 548 (2017) 263–277
Table 5
Unimodal VGM parameters for the filter layer. The shape parameter a and the saturated hydraulic conductivity Ks were estimated using the PSO algorithm.
Soil hydraulic parameters
Layer
Filterhr ()0hs ()
0.15

a (1/cm)
0.001

n ()
3.2

Ks (cm/min)
90

L ()
0.5

Fig. 6. A comparison between measured and simulated outflows versus time and against each other (in the insert) for the calibration period. The full and dashed lines in theinsert are a bisector and linear regression line, respectively.

indicates a perfect agreement between simulated and measured
outflows, and a linear regression line. The good performance ofthe model is confirmed by the determination coefficient R2 = 0.85.
The comparison between the bisector and regression lines indicates
that, in general, the model slightly overestimated the outflow.
A careful inspection of simulated fluxes revealed that the modeltends to underestimate low flows (Fig. 6). This behavior could be
related to the coarse nature of the filter layer, which closely resembles ‘‘fractured aquifers”. Typical breakthrough curves for fractured
aquifers are characterized by early breakthrough and long tailing(Geiger et al., 2010). This is due to a delayed response in the matrix
to pressure head changes that occur in the surrounding fractures.
Brattebo and Booth (2003), Brunetti et al. (2016b), and Fassman
et al. (2010) observed a similar behavior in permeable pavements,
the base and sub-base layers of which are composed of crushedstones. In such circumstances, more complex models are needed
to simultaneously describe fast preferential flow and the matrixfracture interactions in the filter layer (Brunetti et al., 2016b).
Since the computational time associated with a single surrogate
model evaluation was negligible, the approximate response surface
was investigated. Fig. 7 shows the a-Ks response surface obtainedusing a regular grid of 40,000 points. The darkest areas represent
regions of the parameter space with high values of NSE. At the firstinspection, the response surface was characterized by a moderate
multimodality. As shown in Fig. 7, optimal solutions were concentrated in the left part of the plot for values of a between 0.001 and0.05. For a values > 0.05, the objective function dropped quickly to

Fig. 7. The a-Ks response surface obtained using a regular grid of 40,000 points. The
red lines indicate the cross sections reported in Fig. 8.

274

G. Brunetti et al. / Journal of Hydrology 548 (2017) 263–277

0.3–0.4. The NSE values slowly increased for values of a > 0.2. Thiscan be expected since the a parameter for gravels usually has relatively high values. However, in our case, high values of a were notable to accurately reproduce the hydraulic behavior of the SF. For
large a values, the filter layer quickly desaturates and increases
its retention capacity. This effect delays the hydraulic response of
the system to precipitation events. On the contrary, the use of
low a values delays the desaturation of the filter layer, thus reducing its retention capacity (Cey et al., 2006). This makes the filter
layer more responsive to precipitation events, which was the
behavior observed in the measured hydrograph. On the other hand,
the response surface exhibited limited variability in the Ks direction, since high values of NSE were guaranteed for a broad rangeof saturated hydraulic conductivity values.
These results confirmed the findings of the GSA, which clearly
indicated that the variance of the objective function was mainly
driven by the shape parameter a, with only a limited influence of
Ks. This behavior is shown in detail in Fig. 8.
Fig. 8 shows horizontal (Ks = 90 cm/min) and vertical(a = 0.001 1/cm) cross-sections (red lines in Fig. 7) through the
response surface. Yellow rectangular areas are expanded in theright part of the plot. With respect to Ks, the optimum was notclearly identifiable at the first inspection. The values of the objective function ranged between 0.8 and 0.85 in the entire range of thesaturated hydraulic conductivity, indicating its limited effect on
the response surface. The area around the PSO-optimized valueof Ks (90 cm/min) is expanded in the bottom-right corner ofFig. 8. From this plot, it is evident how the PSO algorithm well identified the optimal value Ks, even for a flat profile.
Conversely, Fig. 8 shows a completely different behavior for thea-NSE profile, for which an optimal region was identified in the leftpart of the plot for low values of the a parameter. While the gradient of the curve seemed to approach a maximum, it was not possible to clearly identify the optimum, which may have been outsideof the range imposed on the a parameter. As previously discussed,
this behavior could be related to a deviation from the Darcian flowin the filter layer, which requires a more complex modelingapproach. The further analyzes of the response surface indicatedthat values of a over 0.01 1/cm corresponded to a marked decrease
of NSE. This is evident from the expanded area in the top-right corner of Fig. 8. This finding is in agreement with the results of the
GSA, which identified a as the most influential parameter.
In order to verify whether additional soil hydraulic parameters,such as the saturated water content hs and the pore-size distribution index n, influenced the optimum, a surrogate-based optimization, which considered four soil hydraulic parameters was carriedout. The results of the optimization are listed in Table 6.
As shown in Table 6, the newly estimated parameters are very
similar to those estimated when hs and n were fixed. The saturated
water content was 0.17, which is slightly higher than the previously fixed value. Conversely, the pore-size distribution index n
was slightly lower, but the difference was again very small. The
two most sensitive parameters a and Ks had some small changes,
however overall results are in a good agreement with thosereported previously. Again, the filter layer exhibited a relatively
high permeability and a low value of the a parameter, which is
no longer situated at the border of the evaluated parameter space.
Since a is related to the inverse of the air-entry pressure head, thissuggests that its value is finite and could potentially be measured.
These numerical differences could be related to the effects of interactions between parameters, which, even if small in magnitude as
demonstrated by the GSA, can affect the parameters estimation.
Including these parameters in the optimization process can potentially reduce the parameters uncertainty and help in better defining the problem.
3.5. Model validation
In order to evaluate the reliability of the estimated parameters,
the model was validated using another independent set of experimental data. Fig. 9 shows a comparison between measured andmodeled hydrographs during the validation period.
The value of the objective function was NSE = 0.8, which again
confirmed the adequacy of the estimated parameters. The descrip-

Fig. 8. Horizontal [Ks = 90 cm/min] (top), and vertical [a = 0.001 1/cm] (bottom) response surface cross-sections. The yellow rectangular areas for both plots are expanded on
the right.

275

G. Brunetti et al. / Journal of Hydrology 548 (2017) 263–277

Table 6
Unimodal VGM parameters for the filter layer. The saturated water content hs, the shape parameters a and n, and the saturated hydraulic conductivity Ks were estimated using the
PSO algorithm.
Soil hydraulic parameters
Layer
Filterhr ()
0hs (-)
0.17

a(1/cm)
0.003

n ()
3.1

Ks (cm/min)
89

l ()
0.5

Fig. 9. A comparison between measured and simulated outflows versus time and against each other (in the insert) during the validation period. The full and dashed lines in
the insert are a bisector and linear regression line, respectively.

tion of the hydraulic behavior of the stormwater filter during rainfall events was satisfactory. This capability of the calibrated model
to correctly describe the hydraulic behavior of the system is important when dealing with the analysis of combined traditional drainage systems and LID techniques. A correct description of thehydrograph during precipitation events gives information about
the lag time and the intensity of peak flow, which are fundamental
for both a comprehensive hydraulic analysis of drainage systems,
and for the evaluation of benefits of LIDs implementations. The
model was not able to reproduce outflow induced by the precipitation event on March 15. This may be related to an overestimationof actual evapotranspiration calculated using values from the literature for albedo, LAI, and a hypothetical roots distribution, which
could result in an overestimation of the storage capacity of the
SF at the beginning of the precipitation event, which had a total
volume of 6 mm. As a result, the model predicted that the SF
retained all the precipitation volume. A better characterization of
evapotranspiration could help in increasing the accuracy of themodel, which was already high. As for the optimization data set,the model poorly described the long tailing behavior of the hydrograph after precipitations. This behavior is particularly evident forthe precipitation event on March 28.

4. Conclusions and summary
The aim of this study was to demonstrate the benefit ofsurrogate-based modeling in the numerical analysis of Low Impact
Development techniques. In particular, the unsaturated hydraulicproperties of a contained stormwater filter installed at the University of Calabria were evaluated in the study. The kriging technique
was used to approximate the deterministic response of the widelyused mechanistic model HYDRUS-2D, which was used to simulate
the variably-saturated hydraulic behavior of the filter. In order toreduce the dimensionality of the inverse problem, the simplified
evaporation method was used to determine the unsaturated soil
hydraulic properties of the soil substrate placed on the top of thefilter layer. The Nash-Sutcliffe efficiency index was used both tocompare the simulated and measured outflows, and as the variableof interest for the construction of the response surface.
The PSO heuristic algorithm was used to estimate the kriging
parameters based on an initial set of design sites obtained using
Latin Hypercube Sampling. The approximation uncertainty framework improved the accuracy of the surrogate model by using the
expected-improvement approach to select additional points tothe initial design sites. The kriging model was validated against
an unexplored set of points with satisfactory results. The obtained
surrogate was then used to perform a Global Sensitivity Analysis of
the hydraulic parameters of the filter layer based on Sobol’s
method, with a negligible computational cost. The sensitivity analysis revealed that the model is additive, and that two soil hydraulicparameters, the shape parameter a and the saturated hydraulicconductivity Ks, mainly affect the hydraulic response of the filter
layer. These two parameters were estimated using the PSO algorithm with a NSE value of 0.85, which indicated an good accuracyof the model. Moreover, the analysis of the response surface confirmed the results of the GSA, identifying a as the most influential
parameter. The reliability of the surrogate-based analysis was evaluated by validating the optimized parameters on an independent
dataset of measured outflows. A NSE value of 0.8 confirmed thereliability of the HYDRUS-2D model calibrated using the kriging
technique. In both the optimization and validation, the model
poorly described the long-tailed behavior of the hydrograph after
precipitations. This could be related to the model inadequacy in276

G. Brunetti et al. / Journal of Hydrology 548 (2017) 263–277

describing both the film and preferential flows, and the matrixfracture interactions that occurs in the filter layer. This suggeststhat more complex models are needed to better describe thehydraulic processes involved.
One of the most widespread criticism against the use of mechanistic models is their computational cost, which limits their adoption. This study has demonstrated how a surrogate-based analysis
can provide an effective solution in overcoming this limitation. In
this paper, the kriging technique was used for highly expensive
computational analyses, such as the GSA and the PSO, with goodresults. The sufficiently accurate reproduction of measured hydrographs, which is extremely important in the analysis of LIDs, confirmed the benefit and reliability of surrogate-based analysis.
This novel study represents the first contribution towards theuse of surrogates for LIDs analysis, which does not need to be limited only to the investigation of soil hydraulic properties. For
example, potential applications can be also targeted to the optimization of the morphological characteristics of the LID itself(depth, slope, plants, etc) for a particular climate, to the optimization of the adsorption properties of filter layers, or for specificdesign aims. Such studies could help in providing a better understanding of LID techniques while promoting the widespread adoption of such systems.
Acknowledgements
The study was co-funded by the Italian National Operative Project (PON): Research and Competitiveness for the convergenceregions 2007/2013-I Axis ‘‘Support to structural changes”, operative objective 4.1.1.1: ‘‘Scientific-technological generators of transformation processes of the productive system and creation of new
sectors”, and Action II: ‘‘Interventions to support industrial
research”. The surrogate-based analysis was conducted using the
Python language (Python Software Foundation, 2013).
References
Abbaspour, K.C., Schulin, R., van Genuchten, M.T., 2001. Estimating unsaturated soil
hydraulic parameters using ant colony optimization. Adv. Water Resour. 24,
827–841. http://dx.doi.org/10.1016/S0309-1708(01)00018-5.
Abbaspour, K.C., Johnson, C.A., van Genuchten, M.T., 2004. Estimating Uncertain
Flow and Transport Parameters Using a Sequential Uncertainty Fitting
Procedure. Vadose Zo. J. 3, 1340–1352. http://dx.doi.org/10.2113/3.4.1340.
Abu-Zreig, M., Rudra, R.P., Whiteley, H.R., 2001. Validation of a vegetated filter strip
model (VFSMOD). Hydrol. Process. 15, 729–742. http://dx.doi.org/10.1002/
hyp.101.
Akat, S.B., Gazi, V., 2008. Particle swarm optimization with dynamic neighborhood
topology: Three neighborhood strategies and preliminary results. In: 2008 IEEE
Swarm Intelligence Symposium. pp. 1–8. http://dx.doi.org/10.1109/SIS.2008.
4668298.
Allen, R.G., Pereira, L.S., Raes, D., Smith, M., 1998. FAO Irrigation and Drainage Paper
No. 56: Crop Evapotranspiration. FAO, Rome.
Archer, G.E.B., Saltelli, A., Sobol, I.M., 1997. Sensitivity measures, anova-like
Techniques and the use of bootstrap. J. Stat. Comput. Simul. 58, 99–120.
http://dx.doi.org/10.1080/00949659708811825.
Bengtsson, L., Grahn, L., Olsson, J., 2004. Hydrological function of a thin extensive
green roof in southern Sweden. Nord. Hydrol. 36, 259–268.
Berardi, U., GhaffarianHoseini, A., GhaffarianHoseini, A., 2014. State-of-the-art
analysis of the environmental benefits of green roofs. Appl. Energy 115, 411–
428. http://dx.doi.org/10.1016/j.apenergy.2013.10.047.
Blanusa, T., Vaz Monteiro, M.M., Fantozzi, F., Vysini, E., Li, Y., Cameron, R.W.F., 2013.
Alternatives to sedum on green roofs: can broad leaf perennial plants offer
better ‘‘cooling service”? Build. Environ. 59, 99–106. http://dx.doi.org/10.1016/
j.buildenv.2012.08.011.
Borgonovo, E., Castaings, W., Tarantola, S., 2012. Model emulation and momentindependent sensitivity analysis: an application to environmental modelling.
Environ.
Model.
Softw.
34,
105–115.
http://dx.doi.org/10.1016/j.
envsoft.2011.06.006.
Brattebo, B.O., Booth, D.B., 2003. Long-term stormwater quantity and quality
performance of permeable pavement systems. Water Res. 37, 4369–4376.
http://dx.doi.org/10.1016/S0043-1354(03)00410-X.
Breshears, D.D., Rich, P.M., Barnes, F.J., Campbell, K., 1997. Overstory-imposed
heterogeneity in solar radiation and soil moisture in a semiarid woodland. Ecol.
Appl. 7, 1201–1215. http://dx.doi.org/10.1890/1051-0761(1997) 007[1201:
OIHISR]2.0.CO;2.

Brunetti, G., Simunek, J., Piro, P., 2016a. A comprehensive analysis of the variablysaturated hydraulic behavior of a green roof in a mediterranean climate. Vadose
Zo. J. 15. http://dx.doi.org/10.2136/vzj2016.04.0032. In press.
Brunetti, G., Simunek, J., Piro, P., 2016b. A comprehensive numerical analysis of thehydraulic behavior of permeable pavement. J. Hydrol. 540, 1146–1161. doi:
10.2136/vzj2016.04.0032.
Brynjarsdottir, J., O’Hagan, A., 2014. Learning about physical parameters: theimportance of model discrepancy. Inverse Probl. 30, 114007. http://dx.doi.org/
10.1088/0266-5611/30/11/114007.
Carbone, M., Brunetti, G., Piro, P., 2014. Hydrological Performance of a Permeable
Pavement in Mediterranean Climate. In: 14th SGEM GeoConference on Water
Resources. Forest, Marine and Ocean Ecosystems. pp. 381–388. http://dx.doi.
org/10.5593/SGEM2014/B31/S12.050.
Cey, E., Rudolph, D., Therrien, R., 2006. Simulation of groundwater recharge
dynamics in partially saturated fractured soils incorporating spatially variablefracture apertures. Water Resour. Res. 42. http://dx.doi.org/10.1029/
2005WR004589.
Cheviron, B., Coquet, Y., 2009. Sensitivity analysis of transient-MIM HYDRUS-1D:case study related to pesticide fate in soils. Vadose Zo. J. 8, 1064. http://dx.doi.
org/10.2136/vzj2009.0023.
Coffman, L.S., 2002. Low-impact development: an alternative stormwater
management technology. In: France, R.L. (Ed.), Handbook of Water Sensitive
Planning and Design. Lewis Publishers Inc., pp. 97–123.
Collins, K.A., Hunt, W.F., Hathaway, J.M., 2008. Hydrologic comparison of four typesof permeable pavement and standard asphalt in eastern north Carolina. J.
Hydrol. Eng. 13, 1146–1157. http://dx.doi.org/10.1061/(ASCE)1084-0699(2008)
13:12(1146).
Dosskey, M.G., Helmers, M.J., Eisenhauer, D.E., Franti, T.G., Hoagland, K.D., 2002.
Assessment of concentrated flow through riparian buffers. J. Soil Water Conserv.
57, 336–343. VL - 57.
Duan, Q., Sorooshian, S., Gupta, V., 1992. Effective and efficient global optimization
for conceptual rainfall-runoff models. Water Resour. Res. 28, 1015–1031. http://dx.doi.org/10.1029/91WR02985.
Efron, B., Tibshirani, R., 1986. Bootstrap methods for standard errors, confidence
intervals, and other measures of statistical accuracy. Stat. Sci. 1, 54–75.
Elliot, A., Trowsdale, S., 2007. A review of models for low impact urban stormwaterdrainage. Environ. Model. Softw. 22, 394–405. http://dx.doi.org/10.1016/j.
envsoft.2005.12.005.
Fassman, E.A., Blackbourn, S., 2010. Urban runoff mitigation by a permeable
pavement system over impermeable soils. J. Hydrol. Eng. 15, 475–485. http://dx.doi.org/10.1061/(ASCE)HE.1943-5584.0000238.
Feddes, R.A., Kowalik, P.J., Zaradny, H., 1978. Simulation of Field Water Use and Crop
Yield. PUDOC, Wageningen.
Forrester, A.I.J., Sóbester, A., Keane, A.J., 2008. Engineering Design via Surrogate
Modelling: A practical Guide. J. Wiley. http://dx.doi.org/10.1002/
9780470770801.
Geiger, S., Cortis, A., Birkholzer, J.T., 2010. Upscaling solute transport in naturally
fractured porous media with the continuous time random walk method. Water
Resour. Res. 46. http://dx.doi.org/10.1029/2010WR009133. n/a–n/a.
Getter, K.L., Rowe, D.B., Andresen, J.A., 2007. Quantifying the effect of slope on
extensive green roof stormwater retention. Ecol. Eng. 31, 225–231. http://dx.
doi.org/10.1016/j.ecoleng.2007.06.004.
UMS GmbH, 2015. UMS (2015): Manual HYPROP, Version 2015-01.
Gutmann, H.M., 2001. A radial basis function method for global optimization. J.
Glob. Optim. 19, 201–227. http://dx.doi.org/10.1023/A:1011255519438.
Hopmans, J.W., Šimůnek, J., Romano, N., Durner, W., 2002. Inverse modeling oftransient water flow. In: Dane, J.H., Topp, G.C. (Eds.), Methods of Soil Analysis,
Part 4, Physical Methods. SSSA, Madison, WI, pp. 963–1008.
Houska, T., Multsch, S., Kraft, P., Frede, H.-G., Breuer, L., 2013. Monte Carlo based
calibration and uncertainty analysis of a coupled plant growth and hydrologicalmodel. Biogeosci. Discuss. 10, 19509–19540. http://dx.doi.org/10.5194/bgd-1019509-2013.
Hu, S., 1987. Akaike information criterion statistics. Math. Comput. Simul. 29, 452.
http://dx.doi.org/10.1016/0378-4754(87)90094-2.
Ines, A.V.M., Droogers, P., 2002. Inverse modelling in estimating soil hydraulicfunctions: a genetic algorithm approach. Hydrol. Earth Syst. Sci. Discuss. 6, 49–66.
Jones, D.R., 2001. A taxonomy of global optimization methods based on response
surfaces. J. Glob. Optim. 21, 345–383. http://dx.doi.org/10.1023/A:1012771025575.
Jones, D.R., Schonlau, M., William, J., 1998. Efficient global optimization of
expensive black-box functions. J. Glob. Optim. 13, 455–492. http://dx.doi.org/
10.1023/a:1008306431147.
Keating, E.H., Doherty, J., Vrugt, J.A., Kang, Q., 2010. Optimization and uncertainty
assessment of strongly nonlinear groundwater models with high parameter
dimensionality.
Water
Resour.
Res.
46.
http://dx.doi.org/10.1029/
200009WR8584. n/a-n/a.
Kennedy, J., Eberhart, R., 1995. Particle swarm optimization. Eng. Technol., 1942–
1948
Kennedy, J., O’Hagan, A., 2001. Bayesian calibration of computer models. J. R. Stat.
Soc. Ser. B (Statistical Methodol), 425–464. http://dx.doi.org/10.1111/14679868.00294.
Kennedy, J., Spears, W.M., 1998. Matching algorithms to problems: an experimental
test of the particle swarm and some genetic algorithms on the multimodal
problem generator. In: 1998 IEEE International Conference on Evolutionary
Computation Proceedings. IEEE World Congress on Computational Intelligence(Cat. No.98TH8360). IEEE, pp. 78–83. http://dx.doi.org/10.1109/ICEC.1998.
699326.

G. Brunetti et al. / Journal of Hydrology 548 (2017) 263–277
Khu, S.-T., Werner, M.G.F., 2003. Reduction of Monte-Carlo simulation runs for
uncertainty estimation in hydrological modelling. Hydrol. Earth Syst. Sci. 7,
680–692. http://dx.doi.org/10.5194/hess-7-680-2003.
Laloy, E., Rogiers, B., Vrugt, J.A., Mallants, D., Jacques, D., 2013. Efficient posterior
exploration of a high-dimensional groundwater model from two-stage Markov
chain Monte Carlo simulation and polynomial chaos expansion. Water Resour.
Res. 49, 2664–2682. http://dx.doi.org/10.1002/wrcr.20226.
Li, Y., Babcock, R.W., 2014. Green roof hydrologic performance and modeling: areview. Water Sci. Technol. 69, 727–738. http://dx.doi.org/10.2166/
wst.2013.770.
Li, Y., Babcock, R.W., 2015. Modeling hydrologic performance of a green roof systemwith HYDRUS-2D. J. Environ. Eng. 141, 04015036. http://dx.doi.org/10.1061/(ASCE)EE.1943-7870.0000976.
Liang, J.J., Qin, A.K., Suganthan, P.N., Baskar, S., 2006. Comprehensive learning
particle swarm optimizer for global optimization of multimodal functions. IEEE
Trans.
Evol.
Comput.
10,281–295.
http://dx.doi.org/10.1109/
TEVC.2005.857610.
Mackay, D.J.C., 1998. Introduction to Gaussian processes. Neural Networks Mach.
Learn. 168, 133–165. http://dx.doi.org/10.1007/s10067-003-0742-1.
Marquardt, D.W., 1963. An algorithm for least-squares estimation of nonlinear
parameters. J. Soc. Ind. Appl. Math. 11, 431–441. http://dx.doi.org/10.1137/
0111030.
McKay, M.D., Beckman, R.J., Conover, W.J., 1979. Comparison of three methods forselecting values of input variables in the analysis of output from a computer
code.
Technometrics
21,
239–245.
http://dx.doi.org/10.1080/
00401706.1979.10489755.
Metselaar, K., 2012. Water retention and evapotranspiration of green roofs andpossible natural vegetation types. Resour. Conserv. Recycl. 64, 49–55. http://dx.
doi.org/10.1016/j.resconrec.2011.12.009.
Moriasi, D.N., Arnold, J.G., Van Liew, M.W., Binger, R.L., Harmel, R.D., Veith, T.L.,
2007. Model evaluation guidelines for systematic quantification of accuracy inwatershed simulations. Trans. ASABE 50, 885–900. http://dx.doi.org/10.13031/
2013.23153.
Munoz-Carpena, R., Parsons, J.E., 2004. A design procedure for vegetative filter strips
using VFSMOD-W. Trans. Asae 47, 1933–1941.
Nash, J.E., Sutcliffe, J.V., 1970. River flow forecasting through conceptual models:Part I – a discussion of principles. J. Hydrol. 10, 282–290. http://dx.doi.org/
10.1016/0022-1694(70)90255-6.
Newcomer, M.E., Gurdak, J.J., Sklar, L.S., Nanus, L., 2014. Urban recharge beneath low
impact development and effects of climate variability and change. Water
Resour. Res. 50, 1716–1734. http://dx.doi.org/10.1002/2013WR014282.
Nossent, J., Elsen, P., Bauwens, W., 2011. Sobol’ sensitivity analysis of a complex
environmental model. Environ. Model. Softw. 26, 1515–1525. http://dx.doi.org/
10.1016/j.envsoft.2011.08.010.
O’Hagan, A., 2006. Bayesian analysis of computer code outputs: a tutorial. Reliab.
Eng. Syst. Saf. 91, 1290–1300. http://dx.doi.org/10.1016/j.ress.2005.11.025.
Oakley, J.E., O’Hagan, A., 2004. Probabilistic sensitivity analysis of complex models:a Bayesian approach. J. R. Stat. Soc. Ser. B Stat. Methodol. 66, 751–769. http://dx.
doi.org/10.1111/j.1467-9868.2004.05304.x.
OECD, 2013. Water and climate change adaptation: policies to navigate uncharted
waters. OECD Studies on Water. OECD Publishing, Paris. 10.1787/
9789264200449-en.
Pan, L., Wu, L., 1998. A hybrid global optimization method for inverse estimation of
hydraulic parameters: annealing-simplex method. Water Resour. Res. 34,
2261–2269. http://dx.doi.org/10.1029/98WR01672.
Pedersen, M.E.H., 2010. Good parameters for particle swarm optimization,
Technical Report HL1001, Hvass Laboratories.
Pertassek, T., Peters, A., Durner, W., 2015. HYPROP-FIT Software User’s Manual, V.
3.0.
Peters, A., Durner, W., 2008. Simplified evaporation method for determining soil
hydraulic properties. J. Hydrol. 356, 147–162. http://dx.doi.org/10.1016/j.
jhydrol.2008.04.016.
Python Software Foundation, 2013. Python Language Reference, version 2.7. Python
Softw. Found.
Rasmussen, C.E., Williams, C.K.I., 2006. Gaussian Processes for Machine Learning.
MIT Press, Cambridge, MA.
Razavi, S., Tolson, B.A., Burn, D.H., 2012. Review of surrogate modeling in water
resources. Water Resour. Res. 48. http://dx.doi.org/10.1029/2011WR011527. n/a–n/a.
Regis, R.G., Shoemaker, C.A., 2004. Local function approximation in evolutionary
algorithms for the optimization of costly functions. IEEE Trans. Evol. Comput. 8,
490–505. http://dx.doi.org/10.1109/TEVC.2004.835247.
Rezaei, M., Seuntjens, P., Joris, I., Boënne, W., Van Hoey, S., Campling, P., Cornelis, W.
M., 2015. Sensitivity of water stress in a two-layered sandy grassland soil tovariations in groundwater depth and soil hydraulic parameters. Hydrol. Earth
Syst. Sci. Discuss. 12, 6881–6920. http://dx.doi.org/10.5194/hessd-12-68812015.
Ritchie, J.T., 1972. Model for predicting evaporation from a row crop withincomplete cover. Water Resour. Res. 8, 1204–1213. http://dx.doi.org/
10.1029/WR008i005p01204.
Sacks, J., Welch, W.J., Mitchell, T.J., Wynn, H.P., 1989. Design and analysis ofcomputer experiments. Stat. Sci. 4, 409–435. doi:10.1214/ss/1177012413.
Saltelli, A., Annoni, P., 2010. How to avoid a perfunctory sensitivity analysis.
Environ. Model. Softw. 25, 1508–1517. http://dx.doi.org/10.1016/j.
envsoft.2010.04.012.

277

Saltelli, A., Tarantola, S., Saisana, M., Nardo, M., 2005. What is sensitivity analysis?
In: II Convegno Della Rete Dei Nuclei Di Valutazione E Verifica, Napoli 26, 27
Gennaio 2005, Centro Congressi Universitá Federico II, Via Partenope 36.
Saltelli, A., Annoni, P., Azzini, I., Campolongo, F., Ratto, M., Tarantola, S., 2010.
Variance based sensitivity analysis of model output. Design and estimator forthe total sensitivity index. Comput. Phys. Commun. 181, 259–270. http://dx.doi.
org/10.1016/j.cpc.2009.09.018.
Schindler, U., 1980. Ein Schnellverfahren zur Messung der Wasserleitfähigkeit imteilgesättigten Boden an Stechzylinderproben. Arch. für Acker- und Pflanzenbau
und Bodenkd. 24, 1–7.
Schindler, U., Durner, W., von Unold, G., Mueller, L., Wieland, R., 2010a. The
evaporation method: extending the measurement range of soil hydraulicproperties using the air-entry pressure of the ceramic cup. J. Plant Nutr. Soil Sci.
173, 563–572. http://dx.doi.org/10.1002/jpln.200900201.
Schindler, U., Durner, W., von Unold, G., Muller, L., 2010b. Evaporation method formeasuring unsaturated hydraulic properties of soils: extending themeasurement range. Soil Sci. Soc. Am. J. 74, 1071–1083. http://dx.doi.org/
10.2136/sssaj2008.0358.
Schonlau, M., 1997. Computer experiments and global optimization (PhD thesis).
University of Waterloo, Waterloo, Canada.
Sevat, E., Dezetter, A., Servat, E., 1991. Selection of calibration objective functions in
the context of rainfall-runoff modelling in a sudanese savannah area. Hydrol.
Sci. J. Des. Sci Hydrol. 36, 307–330. http://dx.doi.org/10.1080/
02626669109492517.
Shi, Y., Eberhart, R., 1998. A modified particle swarm optimizer. In: 1998 IEEE Int.
Conf. Evol. Comput. Proceedings. IEEE World Congr. Comput. Intell. (Cat.
No.98TH8360), pp.69–73. http://dx.doi.org/10.1109/ICEC.1998.699146.
Shi, Y., Eberhart, R.C., 1999. Empirical study of particle swarm optimization. In:
Proc. 1999 Congr. Evol. Comput. pp. 1945–1950. http://dx.doi.org/10.1109/CEC.
1999.785511.
Šimůnek, J., van Genuchten, M.T., Šejna, M., 2008. Development and applications of
the HYDRUS and STANMOD software packages and related codes. Vadose Zo. J.
7, 587. http://dx.doi.org/10.2136/vzj2007.0077.
Šimůnek, J., van Genuchten, M.T., Šejna, M., 2016. Recent DEVELOPMENTS andapplications of the HYDRUS computer software packages. Vadose Zo. J. 15, 25.
http://dx.doi.org/10.2136/vzj2016.04.0033.
Sobol’, I., 2001. Global sensitivity indices for nonlinear mathematical models andtheir Monte Carlo estimates. Math. Comput. Simul. 55, 271–280. http://dx.doi.
org/10.1016/S0378-4754(00)00270-6.
Sutanto, S.J., Wenninger, J., Coenders-Gerrits, A.M.J., Uhlenbrook, S., 2012.
Partitioning of evaporation into transpiration, soil evaporation andinterception: a comparison between isotope measurements and a HYDRUS1D model. Hydrol. Earth Syst. Sci. 16, 2605–2616. http://dx.doi.org/10.5194/
hess-16-2605-2012.
Taylor, S.A., Ashcroft, G.L., 1972. Physical Edaphology: The Physics of Irrigated andNonirrigated Soils. Freeman, San Francisco, W.H.
Tokunaga, T.K., 2009. Hydraulic properties of adsorbed water films in unsaturated
porous media. Water Resour. Res. 45, 1287–1295. http://dx.doi.org/10.1029/
2009WR007734.
van Dam, J.C., Groenendijk, P., Hendriks, R.F.A., Kroes, J.G., 2008. Advances ofModeling Water Flow in Variably Saturated Soils with SWAP. Vadose Zo. J. 7,640. http://dx.doi.org/10.2136/vzj2007.0060.
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. http://dx.
doi.org/10.2136/sssaj1980.03615995004400050002x.
Vrugt, J.A., Gupta, H.V., Bouten, W., Sorooshian, S., 2003. A Shuffled Complex
Evolution Metropolis algorithm for optimization and uncertainty assessment ofhydrologic model parameters. Water Resour. Res. 39. http://dx.doi.org/10.1029/
2002WR001642. n/a–n/a.
Vrugt, J.A., Schoups, G., Hopmans, J.W., Young, C., Wallender, W.W., Harter, T.,
Bouten, W., 2004. Inverse modeling of large-scale spatially distributed vadosezone properties using global optimization. Water Resour. Res. 40. http://dx.doi.
org/10.1029/2003WR002706. n/a–n/a.
Vrugt, J.A., Stauffer, P.H., Wöhling, T., Robinson, B.A., Vesselinov, V.V., 2008. Inverse
modeling of subsurface flow and transport properties: a review with newdevelopments. Vadose Zo. J. 7, 843. http://dx.doi.org/10.2136/vzj2007.0078.
Wang, G.G., Shan, S., 2007. Review of metamodeling techniques in support ofengineering design optimization. J. Mech. Des. 129, 370. http://dx.doi.org/
10.1115/1.2429697.
Wesseling, J., Elbers, J., Kabat, P., Broek, B. Van den, 1991. SWATRE: instructions for
input. Intern. Note, Winand Star. Cent.
Wong, T.H.F., Fletcher, T.D., Duncan, H.P., Jenkins, G.A., 2006. Modelling urban
stormwater treatment-A unified approach. Ecol. Eng. 27, 58–70. http://dx.doi.
org/10.1016/j.ecoleng.2005.10.014.
Ye, K.Q., Li, W., Sudjianto, A., 2000. Algorithmic construction of optimal symmetric
Latin hypercube designs. J. Stat. Plan. Inference 90, 145–159. http://dx.doi.org/
10.1016/S0378-3758(00)00105-1.
Younes, A., Mara, T.A., Fajraoui, N., Lehmann, F., Belfort, B., Beydoun, H., 2013. Use of
global sensitivity analysis to help assess unsaturated soil hydraulic parameters.
Vadose Zo. J. 12. http://dx.doi.org/10.2136/vzj2011.0150.
Zhang, X., Srinivasan, R., Van Liew, M., 2009. Approximating SWAT Model using
artificial neural network and support vector machine. JAWRA J. Am. Water
Resour. Assoc. 45, 460–474. http://dx.doi.org/10.1111/j.1752-1688.2009.00302.
x.