Journal of Hydrology 540 (2016) 1146–1161

Contents lists available at ScienceDirect

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

Research papers

A comprehensive numerical analysis of the hydraulicbehavior of a permeable pavement
Giuseppe Brunetti a,⇑, Jiří Šimůnek b, Patrizia Piro aabDepartment of Civil Engineering, University of Calabria, Rende, CS 87036, Italy
Department of Environmental Sciences, University of California, Riverside, CA 92521, USAa r t i c l e

i n f oArticle history:
Received 2 March 2016
Received in revised form 1 July 2016
Accepted 18 July 2016
Available online 20 July 2016
This manuscript was handled by PeterK. Kitanidis, Editor-in-Chief, with theassistance of Roseanna M. Neupauer,
Associate Editor
Keywords:
Permeable pavement
Urban drainage
Infiltration
Sensitivity analysis
Global optimization
Particle swarm

a b s t r a c t
The increasing frequency of flooding events in urban catchments related to an increase in impervious surfaces highlights the inadequacy of traditional urban drainage systems. Low Impact Development (LID)
techniques have proven to be a viable and effective alternative by reducing stormwater runoff andincreasing the infiltration and evapotranspiration capacity of urban areas. However, the lack of adequate
modeling tools represents a barrier in designing and constructing such systems. This paper investigatesthe suitability of a mechanistic model, HYDRUS-1D, to correctly describe the hydraulic behavior of permeable pavement installed at the University of Calabria. Two different scenarios of describing thehydraulic behavior of the permeable pavement system were analyzed: the first one uses a singleporosity model for all layers of the permeable pavement; the second one uses a dual-porosity modelfor the base and sub-base layers. Measured and modeled month-long hydrographs were compared using
the Nash-Sutcliffe efficiency (NSE) index. A Global Sensitivity Analysis (GSA) followed by a Monte Carlofiltering highlighted the influence of the wear layer on the hydraulic behavior of the pavement and identified the ranges of parameters generating behavioral solutions. Reduced ranges were then used in the calibration procedure conducted with the metaheuristic Particle swarm optimization (PSO) algorithm forthe estimation of hydraulic parameters. The best fit value for the first scenario was NSE = 0.43; for thesecond scenario, it was NSE = 0.81, indicating that the dual-porosity approach is more appropriate fordescribing the variably-saturated flow in the base and sub-base layers. Estimated parameters were validated using an independent, month-long set of measurements, resulting in NSE values of 0.43 and 0.86
for the first and second scenarios, respectively. The improvement in correspondence between measured
and modeled hydrographs confirmed the reliability of the combination of GSA and PSO in dealing withhighly dimensional optimization problems. Obtained results have demonstrated that PSO, due to its easiness of implementation and effectiveness, can represent a new and viable alternative to traditional optimization algorithms for the inverse estimation of unsaturated hydraulic properties. Finally, the results
confirmed the suitability and the accuracy of HYDRUS-1D in correctly describing the hydraulic behaviorof permeable pavements.
Ó 2016 Elsevier B.V. All rights reserved.

1. Introduction
Progressing urbanization, connected with the demographic
growth of the last decades, has led to an increase in impervious
surfaces in urban catchments at the expense of natural areas. This
long-term process has resulted in the alteration of the naturalhydrological cycle by reducing the infiltration and evaporation
capacity of urban catchments, increasing surface runoff, and reducing groundwater recharge. While some studies highlighted adecrease of recharge as a result of an increase of impervious⇑ Corresponding author.
E-mail address: giusep.bru@gmail.com (G. Brunetti).
http://dx.doi.org/10.1016/j.jhydrol.2016.07.030
0022-1694/Ó 2016 Elsevier B.V. All rights reserved.

surfaces, other studies identified an increase in recharge due tothe leakage of water from an urban infrastructure, such as sewer
and water supply systems. The effect of urbanization on groundwater recharge is discussed in detail by Price (2011).
Another important factor is that the frequency of extreme rainfall events, characterized by high intensity and short duration, is
expected to increase in the near future as a consequence of climatechange (Kundzewicz et al., 2006; Min et al., 2011). For example, arecent study of Wasko and Sharma (2015) identified a strong correlation between intense precipitation peaks and high temperatures. They concluded that the expected global warming could
lead to an increase of short-duration floods. The correlation
between atmospheric temperature and extreme rainfall intensities

G. Brunetti et al. / Journal of Hydrology 540 (2016) 1146–1161

was also confirmed in other studies (e.g., Westra et al., 2014). This
will be accompanied by a more frequent occurrence of flooding
events in urban areas (Carbone et al., 2015b).
The traditional approach to urban drainage systems focuses oncollecting stormwater in piped networks and transporting it offsite as quickly as possible. The increasing frequency of flooding
events proves that a new design paradigm for drainage systems
is needed. This approach must aim to restore the natural hydrological cycle of urban catchments by increasing their evapotranspiration and infiltration capacity. In recent years, Low Impact
Development (LID), an innovative approach to land development,
has gained increasing popularity. LID is a ‘green’ approach forstormwater management that seeks to mimic the natural hydrology of a site using decentralized micro-scale control measures(Coffman, 2002). LID practices consist of bioretention cells, infiltration wells/trenches, stormwater wetlands, wet ponds, level spreaders, permeable pavements, swales, green roofs, vegetated filter/
buffer strips, sand filters, smaller culverts, and water harvesting
systems. Several studies have evaluated the benefits of LIDs. For
example, Newcomer et al. (2014) used a numerical model todemonstrate the benefits of LIDs, in particular of an infiltration
trench, 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. Environmental benefits
included ecological preservation, mitigation of air and water pollution, enhancement of urban hydrology, a decrease of urban heatisland effects, a reduction of energy consumption, etc. Furthermore, green roofs were able to significantly reduce storm-waterrunoff and retain rainfall volume with retention efficiencies ranging from 40% to 80% (Bengtsson et al., 2004); bioretention cells
were shown to reduce average peak flows by at least 45% during
a series of rainfall events in Maryland and North Carolina (Davis,
2008). Even though the results of available studies are encouraging, more research is needed to precisely assess the impact of LIDs
on the hydrological cycle.
Most impervious surfaces in urban catchments consist of roofs,
roads, parking lots and road shoulders. The development of any
large impervious surface commonly leads to multiple impacts on
stream systems. These impacts include higher peak stream flows,which cause channel incision, bank erosion, and increased sediment transport (Trimble, 1997; Whipple et al., 1981). Another consequence of these impervious surfaces is the reduction of
infiltration, which lowers groundwater recharge (Rose and Peters,2001) and potentially also stream base flow (DeWalle et al.,
2000; Simmons and Reynolds, 1982). Permeable pavements represent one solution to the problem of increased stormwater runoff
and decreased stream water quality. They consist of a surface concrete layer, a filter layer made of sand and other materials, a stonybase, and sub-base layers. Permeable pavements offer great advantages in terms of runoff reduction (Collins et al., 2008), water
retention, and water quality (Brattebo and Booth, 2003).
In spite of many well-known benefits of permeable pavements
and other LID practices, the transition to sustainable urban drainage systems is very slow. One of the key limiting factors in thewidespread adoption of such systems is the lack of adequate analytical and modeling tools (Elliot and Trowsdale, 2007). The availability of an effective LID modeling software could encourage awider adoption of LID principles. Although several stormwatermodels can be applied to the LID analysis (Elliot and Trowsdale,2007), most of them do not incorporate accurate descriptions ofhydrological processes involved, which leads to inaccurate predictions. Moreover, existing tools do not incorporate automaticparameter optimization techniques and sensitivity analysis routines, which have proven to be fundamental when the model
includes multiple parameters. In recent years, researchers have

1147

focused their attention on applying and developing physicallybased models for the LID analysis (Carbone et al., 2015a), however
more research is still needed in this direction.
For example, the HYDRUS software suite (Šimůnek et al., 2008)
has been widely used in the literature for the description of thehydraulic behavior of green roofs (Hilten et al., 2008; Li andBabcock, 2015; Newcomer et al., 2014; Palla et al., 2009), withexcellent agreement between numerical simulations and experimental data. Newcomer et al. (2014) investigated the effects of
LIDs on recharge. In their study, the HYDRUS-2D software wasused to simulate flow from an infiltration trench and an irrigatedlawn installed at the San Francisco State University. While themodel was calibrated by comparing the simulated and measured
recharge, only few details were given about the calibration procedure. The calibrated model was then used to simulate the behaviorof LIDs for future precipitation scenarios. Hilten et al. (2008) used
HYDRUS-1D to study the effectiveness of green roofs in mitigatingstormwater. Simulations were run using HYDRUS-1D for a 24-h
design storm to determine peak flow, retention, and detention time
for runoff. Li and Babcock (2015) used HYDRUS-2D to model thehydrologic response of a pilot green roof system. The root-meansquare error deviation (RMSD) between the modeled water contents and field measurements ranged between 0.38 and 1.74%. Thissuggests that the use of mechanistic models, such as HYDRUS, represents one of the most valuable alternatives to empirical and conceptual models for the LID analysis.
Among all LID practices, permeable pavements are those thatlack modeling tools able to describe their hydraulic behavior most.
The heterogeneity of materials that compose a permeable pavement, together with the high infiltration rates (Brattebo and
Booth, 2003), which may lead to preferential flow and especiallyin the base and sub-base layers, pose complex problems in the
numerical modeling of these systems. Very few modeling tools
exist in the literature for permeable pavements. One of them isincluded in the Storm Water Management Model (SWMM)(Gironás et al., 2010). However, results obtained by SWMM have
proven to be inaccurate, especially in the description of the effectsof base and sub-base layers on the infiltration processes (Zhang
and Guo, 2015). HYDRUS has also been used for the descriptionof variably-saturated flow in permeable pavements. Illgen et al.
(2007) used HYDRUS-2D for the numerical analysis of a permeable
pavement and calibrated the model against experimental data collected at a laboratory test facility. The calibrated model was thenused to simulate different scenarios not investigated during thelaboratory campaign. The Illgen et al. (2007) study provided only
limited details about the calibration of soil hydraulic parameters
and their uncertainty and sensitivity. The occurrence of preferential flow in the permeable pavement was also not investigated.
Moreover, the model was used to simulate a laboratory test facility,
the behavior of which can differ from a field scale experimental
facility. On the other hand, Carbone et al. (2014) used HYDRUS1D to model a permeable pavement at the field scale. The
HYDRUS-1D model was calibrated against four different rainfall
events with optimal results. In this study, the permeable pavement
was modeled as a single homogeneous layer and the differences
between hydraulic properties of different layers were neglected.
Furthermore, the numerical simulations were event-based. In bothstudies, calibration of soil hydraulic properties was carried out
manually without taking advantage of more recent global optimization algorithms. This indicates that research in this direction
is limited, with only inconclusive results that need to be further
investigated.
The lack of studies that provide a comprehensive description of
the hydraulic behavior of a permeable pavement at the field scale
and that propose a general methodology for the estimation of its
hydraulic parameters suggests that research is particularly needed

1148

G. Brunetti et al. / Journal of Hydrology 540 (2016) 1146–1161

in the development and identification of accurate modeling tools
for the analysis of LID practices, especially for permeable pavements. The aim of this study is to investigate the suitability of
the HYDRUS mechanistic model to correctly describe unsaturated
flow in typical permeable pavement, installed at the experimentalsite of the University of Calabria. Multiple uniform and nonequilibrium flow models included in HYDRUS-1D, such as single- anddual-porosity models, are used to describe the hydraulic behaviorof the permeable pavement. The problem is addressed in the following way. First, a Global Sensitivity Analysis (GSA) is carriedout to prioritize hydraulic parameters and identify those that are
non-influential. Results of the GSA, combined with a Monte Carlofiltering approach, are used to investigate the parameter space
and identify behavioral regions. These regions are then used in
the calibration process conducted with the Particle Swarm Optimization (PSO) algorithm. The use of PSO for the determinationof unsaturated hydraulic properties represents a new important
application of this method. Finally, the calibrated model is validated on an independent set of measurements.
2. Materials and methods
2.1. 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 an averageannual precipitation of 881.2 mm. The permeable pavement is partof the ‘‘Urban Hydraulic Park,” which also includes an extensive
green roof, a bioretention system, and a sedimentation tank connected with a treatment unit. The permeable pavement has an areaof 154 m2, an average slope of 2%, and a total depth of the profile of
0.98 m. Fig. 1 shows a schematic of the permeable pavement, consisting of 5 layers.
The surface wear layer consists of porous concrete blocks characterized by high permeability. Base, sub-base and bedding layers
were constructed by following the suggestions of the Interlocking
Concrete Pavement Institute (ICPI), which recommends certain
ASTM stone gradations. The ASTM numbers and corresponding
gradations can be found in ASTM D 448, Standard Classificationfor Sizes of Aggregate for Road and Bridge Construction. The ASTM
No. 57, used for the base layer, is characterized by a porosity of
about 30–35%. The ASTM No. 2 is used in the sub-base layer for

its stability and a high volumetric porosity of about 40%. The ASTM
No. 8 is used for the bedding layer and the protection layer and has
a porosity of about 20% of volume. The bedding layer is composedof a mixture of sand, glass sand, and zeolite to improve the pollutant removal efficiency of the permeable pavement for typical contaminants of stormwater runoff. A high permeability geotextile
with a fiber area weight of 60 g/m2 is placed at the interface
between the bedding layer and the base layer to prevent sand from
migrating into the bottom layers. An impervious membrane isplaced at the bottom of the profile to prevent water from percolating into deeper horizons. The protection layer which is composedof coarse sand is placed between the sub-base layer and the impervious membrane. The baseflow is collected in a horizontal drain,
which consists of a perforated PVC pipe, and is conducted to amanhole for quantity and quality measurements.
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 by a tipping bucket rain gauge with a resolution of
0.254 mm and an acquisition frequency of 1 min. Climatic data
are acquired with a frequency of 5 min. Data are processed andstored in the SQL database.
Two flux meters, composed of a PVC pipe with a sharp-crested
weir and a pressure transducer, measure baseflow and runoff fromthe permeable pavement. The pressure transducer (Ge Druck
PTX1830) measures the water level inside the PVC pipe and hasa range of measurement of 75 cm with an accuracy of 0.1% of thefull scale. The pressure transducers were calibrated in the laboratory by using a hydrostatic water column, linking the electric current intensity with the water level inside the column. The
exponential head-discharge equations for the two PVC flux meters
were obtained by fitting the experimental data with a coefficient ofdetermination R2 = 0.999 for both devices. Runoff and baseflow
data were acquired with a time resolution of 1 min and stored in
the SQL database. No measurements of pressure heads or volumetric water contents inside the pavement were taken.
Two month-long data sets were selected for further analysis(Fig. 2). The first data set, which started on 2014-01-15 and ended
on 2014-02-15, was used for parameter optimization and sensitivity analysis. Total precipitation and total potential evapotranspiration for the first data set were 274 mm and 43 mm, respectively.
The second data set, which started on 2014-03-01 and ended on
2014-03-31, was used for model validation. Total precipitation

Fig. 1. A schematic of the permeable pavement.

1149

G. Brunetti et al. / Journal of Hydrology 540 (2016) 1146–1161

Fig. 2. Precipitation and subsurface flow during the optimization (top) and validation (bottom) time periods.

and total potential evapotranspiration for the second data set were
175 mm and 81 mm, respectively. The second data set was
selected so that it had significantly different meteorological data
than during the first period. The optimization set is characterizedby multiple rain events with few dry periods. The validation set
has fewer rain events, which are concentrated at the beginningand end of the time period and separated by a relatively long dry
period between. Surface runoff was not observed during these time
periods.
Potential evaporation was calculated using the PenmanMonteith equation (Allen et al., 1998). The permeable pavement
was installed in 2013 and has been constantly exposed to atmospheric conditions and traffic since then that has altered the surface roughness and color. For these reasons, an albedo of 0.25
was used as suggested by Levinson and Akbari (2002) for weathered gray cement.
2.2. Theory
Water flow simulations were conducted using the HYDRUS-1D
software (Šimůnek et al., 2008). HYDRUS-1D is a one-dimensional
finite element model for simulating the movement of water, heat,and multiple solutes in variably-saturated porous media. HYDRUS1D implements multiple uniform (single-porosity) and nonequilibrium (dual-porosity and dual-permeability) water flow models(Šimůnek and van Genuchten, 2008). In this study, two different
conceptual models were used to represent flow in the permeable
pavement (Table 1).
Scenario I assumed that water flow in all five soil layers of thepermeable pavement can be described using the classical single-

Table 1
Conceptual models representing water flow in the permeable pavement.
Soil layer

Scenario I

Scenario II

Wear
Bedding
Base
Sub-base
Protection

Single Porosity
Single Porosity
Single Porosity
Single Porosity
Single Porosity

Single Porosity
Single Porosity
Dual Porosity – MIM
Dual Porosity – MIM
Single Porosity

porosity approach (SPM). Unsaturated water flow is then describedusing the one-dimensional Richards equation:




@h
@
@h
¼
KðhÞ
þ1
@z @z
@z

ð1Þ

where h is the volumetric water content [–], h is the soil water pressure head [L], K(h) is the unsaturated hydraulic conductivity [L T1],
t is time [T], and z is the soil depth [L]. The soil hydraulic propertiesare described by the van Genuchten – Mualem relation (vanGenuchten, 1980):(h¼

hs hrm þ hrð1þðajhjÞn Þ

if h 6 0hs

if h > 0h  hrhs  hr
8
h

 i2
1 m
<
K s SLe 1  1  Sme
K¼
:
Ks

ð2Þ

Se ¼

if h 6 0
if h > 0

ð3Þ

1150

m¼1

G. Brunetti et al. / Journal of Hydrology 540 (2016) 1146–1161

1nwhere hr [–] is the residual water content, hs [–] is the saturated
water content, Ks [L T1] is the saturated hydraulic conductivity, n
is a pore-size distribution index [–], a is a parameter related tothe inverse of the air-entry pressure [L1], L indicates the tortuosity
and is usually assumed to be 0.5 for many soils, and Se is the effective saturation [–]. In order to simplify the model (to lower thenumber of unknown parameters), the residual water content of alllayers was fixed. In particular, the residual water content for thewear and bedding layers was assumed to be 0.045 and 0.03, respectively, while the residual water content for both the base and subbase layers was assumed to be 0.0, considering that they were composed of crushed stones. Furthermore, considering that the bedding
layer and the protection layer had the same stone gradation, ASTM
No. 8, the same set of parameters was used for both. Despite of allthese considerations, this scenario still involves 16 parameters (hs,a, n, and Ks for 4 soil layers).
Scenario II assumes a single-porosity model for the wear layer,
the bedding layer, and the protection layer, and a dual-porosity
model for the base and sub-base layers. This configuration was
selected in order to consider the occurrence of preferential flowin the coarse layers of the pavement.
The base and sub-base layers are composed of crushed stones,
with particle size diameters ranging from 2.5 to 37 mm in the baselayer and from 20 to 75 mm in the sub-base layer. Crushed stones
were washed before installation in order to remove fine particles.
This narrow gradation provides a high volume of voids and increases
the water storage and infiltration capacities of these two layers.
From a physical point of view, the structure of the base and subbase materials closely resembles fractured aquifers (Barenblatt
et al., 1960). Fractured aquifers are represented by a blocky matrix
system intercepted by fractures. Open and well-connected fractures
represent high permeability pathways that are many orders of magnitude more permeable than the porous rock matrix. At the sametime, one of the characteristics of a fractured aquifer is that the fractures occupy a much smaller volume than the pores of the rockmatrix. Traditionally, fractured porous media are thus represented
by two separate flow domains: the high permeability (mobile)
domain, the network of connected fractures characterized by advective flow, and the low permeability (immobile) domain, dominatedby diffusion. The rock matrix also provides storage capacity becauseof its significantly larger volume than the fracture system. Typicalbreakthrough curves for a fractured aquifer are characterized byearly breakthrough and long tailing (Geiger et al., 2010). This isdue to the fact that the matrix has a delayed response to pressure
head changes that occur in the surrounding fractures. The resulting
pressure difference induces matrix-fracture interflow. This flow
takes place after initial fracture flow and before the matrix and fracture pressures equilibrate (Bai et al., 1994). Several studies havedemonstrated the long tailing from permeable pavements in discharge hydrographs (e.g., Brattebo and Booth, 2003; Fassman andBlackbourn, 2010) and attributed this effect to the storage and flow
through the base and sub-base layers.
The classical approach to model water flow in fractured porous
media is the so-called ‘‘dual-porosity” or ‘‘mobile-immobile water”(MIM) approach (Barenblatt et al., 1960; van Genuchten and
Wierenga, 1976; Warren and Root, 1963). This approach assumes
that flow occurs only in the mobile fracture domain, for whichan effective permeability must be known, while water in thematrix domain is immobile. Both domains are connected by a simple first-order transfer function, which accounts for the exchangeof fluid across the boundary of the two domains.
In the dual-porosity approach, the liquid phase is divided into
two domains:h ¼ hf þ hm

ð4Þ

where subscript f refers to the (mobile) fracture system, and subscript m refers to the immobile matrix domain. The dual-porosity
water flow formulation is based on a modified Richards equation
for flow in fractures and a mass balance equation for moisturedynamics in the matrix:




@hf
@
@h
KðhÞ
þ 1  Cw
¼
@z
@z
@z

ð5Þ

@hm
¼ Cw
@t

ð6Þ

where Cw is the mass transfer between two domains, which isassumed to be proportional to the difference in effective saturationsof the two regions (Šimůnek and van Genuchten, 2008; Simunek
et al., 2003):

im
Cw ¼ x  Sme  Se

ð7Þ

where x is a first-order coefficient [T1]. Compared to assuming apressure head based driving force for the mass transfer, the dualporosity model based on (7) requires significantly less parameters
since one does not need to know the retention function (and corresponding parameters) for the matrix region explicitly, but only its
residual and saturated water contents (Simunek et al., 2003). The
residual water content for the mobile domain of both the baseand sub-base layers is assumed to be 0.0 (Simunek et al., 2003).
The tortuosity factor, L, is again assumed to be 0.5 for all layers. Scenario II thus includes 20 parameters (additionally also x and hs of
the immobile domain for the base and subbase layers).
2.2.1. Numerical domain and boundary conditions
The numerical domain representing the stratigraphy of the permeable pavement was divided in 5 layers. The bedding layer andthe protection layer had the same properties since they were constructed using the same ASTM No. 8 stone gradation. A relatively
fine, finite element mesh with a constant element size of 0.5 cm
was used in order to minimize mass balance errors and avoid
non-convergent runs during sensitivity analysis and parameter
optimization. An atmospheric boundary condition was applied at
the pavement surface using (a) precipitation and potential evaporation fluxes, (b) a prescribed zero pressure head (saturation) during ponding, and (c) equilibrium between the pavement surface
water content and atmospheric water vapor when atmosphericevaporative demand could not be met by the wear layer. A seepage
face boundary condition was specified at the bottom of the protection layer. A seepage face boundary acts as a zero pressure headboundary when the bottom boundary node is saturated and as a
no-flux boundary when it is unsaturated. The initial conditions
were specified in terms of the soil water pressure head and were
set to linearly increase with depth, from 90 cm at the top of theflow domain (z = 0) to 0.5 cm at the bottom (z = 98). The surface
layers are assumed to be drier than the bottom layers since they
are directly exposed to evaporation.
2.2.2. Objective function
The Nash-Sutcliffe Efficiency (NSE) index (Nash and Sutcliffe,1970) is used for evaluating the agreement between hydrographs:2

2 3
PT  obs
mod
Q

Q
i
6 i¼1 i7
NSE ¼ 1  4
2 5PT  obs
obs
 Q meani¼1 Q ið8Þ

where Qobs
is the ith measured value, Qmod
is the ith simulatediivalue, and Qobs
mean is the mean value of observed data. The NSE coefficient ranges between 1 and 1.0, is equal to 1 in case of a perfect1151

G. Brunetti et al. / Journal of Hydrology 540 (2016) 1146–1161

agreement, and, generally, values between 0.0 and 1.0 are considered acceptable (Moriasi et al., 2007). The NSE has been usedbecause it is often reported to be the best measure for evaluating
the overall fit of a hydrograph (Sevat et al., 1991).

of the model output variance that would disappear when parameter
Xi is fixed. When the model is additive, i.e., when it does not include
interactions between input factors, then the first-order index is sufficient for decomposing the model’s variance. For additive models,the following relation is valid:

2.2.3. Global sensitivity analysis
Most existing environmental models include a high number ofparameters. This aspect creates a major problem in their application, as the parameter estimation becomes a high-dimensionaland mostly nonlinear problem. To solve this problem, several optimization algorithms were developed (Beven and Binley, 1992;
Duan et al., 1992; Poli et al., 2007; Vrugt et al., 2003). Moreover,environmental optimization studies are often affected by the equifinality problem (Beven, 2006) when multiple sets of parameters
can produce similar results. This problem is exacerbated when
the number of parameters is significant and only limited information about their interactions and their effects on the output isavailable. However, it is not always necessary to include all model
parameters in the optimization process because some of them
could be measured or estimated, and some may have negligible
effects on the output of the model for a particular application. A
sensitivity analysis (SA) can identify the most influential parameters and their interactions and how these parameters affect theoutput (Saltelli et al., 2005).
The principal steps of a SA are: Factors Prioritization (FP), Factors Fixing (FF), Variance Cutting (VC), and Factors Mapping (FM)(Saltelli and Tarantola, 2004). The aim of FP is to identify factorsthat one should measure in order to obtain the greatest reductionin the uncertainty of the output. Conversely, FF identifies factorsthat are non-influential. By applying these two settings, the modeler is able to reduce the dimension of the optimization problem
and have a complete appreciation of the parameters’ influences
and interactions.
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 of parameters one-at-a-time while
keeping the others constant (Cheviron and Coquet, 2009; Houska
et al., 2013; Rezaei et al., 2015). However, when the model includes
interactions between parameters, results of the OAT analysis areinaccurate because parameter interactions can be identified only
by changing multiple parameters simultaneously. For this reason,
when the property of a model is a priori unknown, a Global Sensitivity Analysis (GSA) is always preferred (Saltelli and Annoni,2010). Practitioners call this analysis a model-free setting.
One of the most widespread algorithms for the GSA is thevariance-based Sobol0 method (Sobol0 , 2001). Variance-based
methods aim to quantify the amount of variance that each parameter contributes to the unconditional variance of the model output.
For the Sobol0 method, these amounts are represented by Sobol’s
sensitivity indices (SI’s). These indices give quantitative information about the variance associated with a single parameter or
related to interactions of multiple parameters. For a more complete
explanation about the Sobol0 method, please refer to Sobol0 (2001).
Sobol’s sensitivity indices are expressed as follows:X
Si ¼ 1First Order Si ¼

Vi
V
V ij
V

ð10Þ

X
Sij þ   

ð11ÞSecond Order 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”. This index can be described as the fraction

ð12Þ

i

Even when the model includes interactions between parameters, the first-order index remains the measure to use for FP(Saltelli and Tarantola, 2004). On the other hand, the total effectindex, ST, gives information about a non-additive part of the model.
A significant difference between ST and Si indicates an importantrole of an interaction for the parameter considered. Essentially,
the total effect index, STi, gives a fraction of the total variance thatwould be left when all factors but Xi were fixed. STi = 0 is a condition necessary and sufficient for Xi to be non-influential. Therefore,
Xi can be fixed at any value within its range of uncertainty without
affecting the output unconditional variance. The total effect is themeasure to use for FF.
Considering that environmental models are generally highly
nonlinear, it is almost impossible to calculate the variances using
analytical integrals. Hence, Monte Carlo integrals are often applied,which are based on sampling the parameter space in q samples.
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. For an accurate description of the calculation of Sobol’s indices please refer to Saltelli et al.
(2010).
Basically, the calculation of Sobol’s indices requires q  (2p + 1)
model evaluations, where p is the number of input factors. However, Saltelli (2002) introduced a method that requires only q  (p+ 2) model evaluations. To sample the parameters’ space we used
Sobol’s quasi-random sampling technique (Sobol0 , 2001).
One of the most important aspects of the GSA is the choice ofthe number of samples, q. An increase in the number of samples
will increase the accuracy of Sobol’s indices. However, a high q
implies a higher number of model evaluations. The number of samples is case-sensitive; it depends on the structure of the model andon the type of simulations performed. A convergence analysis ofSobol’s indices is the recommended procedure for estimating q.
However, this approach is time consuming because it needs torepeat the GSA several times by increasing the number of samples
until the variability of indices between two consecutive analyses isbelow a threshold value for all parameters.
In a recent study, Nossent et al. (2011) gave a comprehensive
description of the influence of q on the accuracy of a GSA for anenvironmental model that included 26 parameters. Nossent et al.
(2011) reported that for most parameters, less than 5000 samples
were sufficient to reach a stable solution. An extensive review of
the GSA in hydrological models is reported in Song et al. (2015).
Here, we report the number of model runs for each GSA performed,
together with the type of GSA, the number of parameters of themodel, and the objective function used. For the GSA based on
Sobol’s method, the number of model runs rarely exceeds
100,000. Due to considerations discussed above, a value ofq = 5000 was chosen in our study. Table 2 summarizes the characteristics of the GSA for the two scenarios considered.

Table 2
Number of parameters and HYDRUS-1D runs for both scenarios.
Scenario

Number of parameters

Model runs

I
II

16
20

90,000110,000

1152

G. Brunetti et al. / Journal of Hydrology 540 (2016) 1146–1161

In order to assess the accuracy of estimations of the sensitivity
indices, the bootstrap confidence intervals (BCIs) (Efron andTibshirani, 1986) were estimated. The basic idea of the bootstrapping is that, in absence of any other information about the distribution, the sample contains all the 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
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 q samples used for the model evaluation
were sampled 1000 times with replacement, whereby Sobol’s
indices were calculated for each resampling. In this way, 95% confidence intervals are constructed by using the percentile method
and the moment method (Archer et al., 1997).
The sensitivity analysis was conducted using the programming
language Python and in particular, the Sensitivity Analysis Library(SALib) (Usher et al., 2015). An elaborated script overwrites the
input file containing the parameters for different materials at each
iteration. The script then executes HYDRUS-1D, which usually runs
less than one second. If the HYDRUS-1D run is not finished after
15 s, it is considered non-convergent; the script then terminates
the process and attributes a large negative value to the objectivefunction. The same negative value is attributed when the lengthof the modeled hydrograph is shorter than one month, which indicates that the run was unsuccessful. Values of the objective function are stored in a one-dimensional array for the subsequent
computation of sensitivity indices. Table 3 reports the initial rangeof all evaluated parameters in the two scenarios. The initial conditions were not included in the GSA because their effects on thehydrograph for a month-long simulation are assumed to be limited
to only the first few days.
2.2.4. Monte Carlo filtering
In the context of an optimization framework, results of the GSA
can be used to extract useful information about the problem structure. The GSA preliminarily identifies the subset of input factorsthat drive most of the variation in the model output; to establish
their optimal values, these sensitive parameters can be further
investigated by using a Monte Carlo filtering approach. Filtering

Table 3Ranges of parameters used in the GSA for both scenarios.
Parameter

Scenario I
Initial range

Scenario II

hs1 [–]
a1 [1/cm]
n1 [–]
Ks1 [cm/min]
hs2 [–]
a2 [1/cm]
n2 [–]
Ks2 [cm/min]
hs3 [–]
a3 [1/cm]
n3 [–]
Ks3 [cm/min]hs,im3 [–]x3 [1/min]
hs4 [–]
a4 [1/cm]
n4 [–]
Ks4 [cm/min]hs,im4 [–]x4 [1/min]

0.2–0.50.001–0.3
1.1–4.5
1.0–20.0
0.2–0.50.001–0.3
1.1–4.5
1.0–20.0
0.01–0.40
0.001–0.3
1.1–4.5
1.0–100.0––
0.01–0.4
0.001–0.3
1.1–4.5
1.0–100.0––

0.2–0.50.001–0.3
1.1–4.5
1.0–20.0
0.2–0.50.001–0.3
1.1–4.5
1.0–20.0
0.001–0.1
0.001–0.3
1.1–4.5
1.0–100.0
0.15–0.4
0.00001–0.009
0.001–0.1
0.001–0.3
1.1–4.5
1.0–100.0
0.15–0.4
0.00001–0.009

techniques are used to explore the parameter space pertaining to
the single or multiple optima. This is particularly relevant when
dealing with mechanistic models that almost always contain illdefined parameters and are thus referred to as overparameterized models (Draper and Smith, 1981).
The Monte Carlo filtering is often coupled with the regionalized
sensitivity analysis (RSA) (Hornberger and Spear, 1981). The RSA
generally requires two tasks: (a) a qualitative description of thesystem behavior and (b) a binary classification of the model output
that divides solutions into two behavioral and non-behavioral
groups. However, the main drawback of the RSA is that nohigher-order analysis is performed and thus interactions between
parameters are not investigated. In the GSA, a complete descriptionof main effects and interactions is given. The GSA has been combined effectively with the GLUE analysis (Beven and Binley,1992) in the context of the parameter optimization (Ratto et al.,2001). In Ratto et al. (2001), the sample generated for the GLUE
analysis is also used for the computation of variance-based sensitivity indices.
In this study, the GSA is coupled with a basic Monte Carlo filtering. The aim of this step is to identify behavioral regions in theparameter space and to reduce the uncertainty in the following
parameter estimation step by using the same sample and runs of
the GSA. For each parameter set used in the GSA, a value of theobjective function is calculated. Potential solutions are divided into
two groups: behavioral, solutions with NSE > 0.0, and non-behavioral, solutions with NSE 6 0.0. Two different types of analysis were
performed on the filtered sample: (a) Kernel density estimationand (b) correlation analysis.
2.2.5. Kernel Density Estimation (KDE)The KDE plots have been used to identify regions with a highdensity of behavioral solutions. The KDE is a non-parametric estimator of the probability density function (PDF) of a random variable (Silverman, 1981). A kernel is a special type of PDF with an
added property that it must be even. The KDE bi-variate plots havebeen used because they give a smooth qualitative representation of
PDFs in a bi-dimensional space. The uni-variate KDE has also been
computed for each parameter. The KDE plots have been calculatedusing a Gaussian kernel and the Scott procedure for the determination of a bandwidth (Scott, 1992).
2.2.6. Correlation analysis
The correlation analysis helps to identify particular interaction
structures between parameters. Detecting high values of correlation coefficients suggests a way to reduce the input factor space.
In particular, when the coefficient is positive, the couple of parameters acts in the model as a quotient/difference, and when it is negative, the parameters act as a product/sum.
2.2.7. Particle swarm optimization
Inverse modeling is a procedure to estimate unknown parameters of the model from experimental data. One of the major reasons
to apply inverse modeling is to estimate parameters that cannot bedirectly measured for various reasons. Numerous applications of
inverse modeling for the estimation of soil hydraulic properties
exist in the literature (Abbaspour et al., 2004; Hopmans et al.,
2002; Vrugt et al., 2008, 2004). The gradient methods(Marquardt, 1963) have been used most widely among hydrologists and soil scientists. However, these methods are sensitive tothe 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 com-1153

G. Brunetti et al. / Journal of Hydrology 540 (2016) 1146–1161
Table 4
Parameters used in the PSO optimization.
N

c1

c2

w

69

0.2673.395

0.444plex 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) is used. PSO
has been used in multiple studies involving inverse modeling withcomplex environmental models (Gill et al., 2006; Jiang et al., 2010;
Zambrano-Bigiarini and Rojas, 2013). However, so far it has notbeen used for the determination of unsaturated hydraulic properties. PSO is a relatively new algorithm for evolutionary computation methodology, but its performance has proven to becomparable to various other, more established methodologies(Kennedy and Spears, 1998; Shi et al., 1999). One of the main
advantages of PSO is the easiness of its implementation (Lianget al., 2006). PSO is characterized by an algorithm based on asocial-psychological metaphor involving individuals that interact
with each other in a social world. PSO was inspired by the behaviorof schools of fish or flocks of birds as they seek food or otherresources. In PSO, collections of ‘‘particles” explore the search
space, looking for a global or near-global optimum. Particles in
PSO keep track of their best positions thus far obtained in the
search space and the best positions obtained by their neighboring
particles. The global best position is what all particles tend to follow. A detailed description of the PSO algorithm is given in Shi andEberhart (1998).
The most important parameters in the PSO are: c1, c2, and w. c1
and c2 are constant parameters known as the cognitive and social
parameters, respectively, and w is the inertia-weight, which plays
a key role in the optimization process by providing balance
between exploration and exploitation. A large w facilitates a global
search while a small one facilitates a local search. The w parameteris very similar to the ‘‘temperature” parameter in the simulated
annealing algorithm. While several strategies have been used in
the literature for the inertia weight, in this study, a constant valueof w has been used (Shi and Eberhart, 1998).
In PSO, each particle is influenced by its r nearest neighbors.
The arrangement of neighbors that influence a particle is called
the topology of the swarm. Different types of neighborhoods arereported in the literature (Akat and Gazi, 2008). In this study, theall topology is used, in which the neighborhood encompasses theentire swarm. The PSO parameters used in this study for both scenarios are reported in Table 4 and are as suggested by Pedersen(2010).
A modified version of the PySwarm Python Library was used for
the PSO analysis. Similar to the GSA, a Python script has been written for the optimization process. The script overwrites the input
file of HYDRUS-1D containing the hydraulic parameters for the different layers, runs the executable module, and retrieves the valueof the objective function. A large negative value of NSE is attributed
to non-convergent runs, as defined above.
3. Results and discussion

Table 5First-order (S1) and total (ST) effect indices (in decreasing order) with their bootstrap
confidence intervals (BCI) for parameters of Scenario I.
Parameter

S1

S1 (BCI)

ParameterST

ST (BCI)n1 [–]
a1 [1/cm]
Ks1 [cm/min]
hs3 [–]
a4 [1/cm]
a2 [1/cm]
n3 [–]
Ks4 [cm/min]n4 [–]
hs1 [–]
n2 [–]
Ks3 [cm/min]
hs4 [–]
a3 [1/cm]
hs2 [–]
Ks2 [cm/min]
Sum

0.298
0.102
0.051
0.023
0.020
0.017
0.014
0.009
0.009
0.009
0.007
0.004
0.001
0.001
0.004
0.005
0.563

0.054
0.040
0.040
0.024
0.023
0.022
0.029
0.025
0.035
0.028
0.023
0.022
0.022
0.024
0.019
0.016

n1 [–]
a1 [1/cm]
Ks1 [cm/min]
hs1 [–]
n4 [–]
Ks3 [cm/min]n3 [–]
a3 [1/cm]
a2 [1/cm]
n2 [–]
a4 [1/cm]
hs3 [–]
hs2 [–]
Ks2 [cm/min]
Ks4 [cm/min]
hs4 [–]

0.745
0.508
0.421
0.247
0.224
0.210
0.194
0.181
0.176
0.176
0.170
0.167
0.151
0.138
0.138
0.136
>1.0

0.042
0.032
0.032
0.025
0.146
0.127
0.035
0.024
0.024
0.028
0.033
0.031
0.030
0.023
0.038
0.022

variance, the pore-size distribution index n1 and the air-entry pressure parameter a1. The third most influential parameter, the saturated hydraulic conductivity Ks1, has the effect, which is only halfof the second most influential parameter, a1. Ten parameters havea first-order index lower than 1%, which indicates that their main
effect on the output variance is negligible. Table 5 also shows thatthe sum of all first-order indices is less than 1, which means thatthe model is non-additive. Only 56% of variance is attributable tothe first-order effects, which indicates that interactions between
parameters play a fundamental role.
The right part of Table 5 (ST) shows that almost 75% of variancein simulated outflow is caused by n1, either by the variation of theparameter itself (30%) or by interactions with other parameters.
Together with a1 (51%) and Ks1 (42%), it is the most influential
parameter for simulated flow. It can be noted that the saturated
hydraulic conductivity, Ks1, has a relatively low main effect but arelatively high total effect. That indicates that this parameter has
a limited direct effect on the variance of the objective function,
but it has an effect in interactions with other parameters.
The effect of the sub-base layer on the output is less significant,while the wear layer strongly conditions the output. That behavior
is in agreement with results reported in the literature. Illgen et al.
(2007), in his laboratory campaign, confirmed that the wear layer
has the major influence on the infiltration capacity of the permeable pavement, while the base and sub-base layers have a minor
impact and act as a storage tank. The total index is always greaterthan zero, which implies that all parameters influence the output
variance either directly or by their interactions, and thus noparameter can be fixed without affecting the uncertainty of theoutput.
Scatter plots for the plain Monte Carlo runs for the two mostsensitive parameters, a1 and n1, are displayed in Fig. 3. The scatterplots show that there is no clear pattern of factors driving bad solutions. Particular trends in the solutions were further identified by
the regression lines. They indicate that there is a trend for parameter a1, with admissible solutions in the left part of the plot. On theother hand, the distribution of values for parameter n1 is flat, andthus no conclusions can be made about the position of a denser
region of behavioral solutions in the high-dimensional space.

3.1. Sensitivity analysis – Scenario I
3.2. Monte Carlo filtering – Scenario I
As discussed above, the basic outcome of Sobol’s SA are thefirst-order (S1) and total (ST) sensitivity indices. Table 5 presents
these two indices and their relative bootstrap confidence intervals(BCI). In the left part of Table 5 (S1), it can be seen that only two
parameters exhibit a significant direct influence on the output’s

A Monte Carlo Filtering procedure was applied to the runs of the
GSA. The threshold value of NSE = 0.0 produced a filtered sample
composed of 1452 behavioral solutions. Fig. 4 shows the univariate
and bivariate KDE plots and the correlation plots for the wear layer.

1154

G. Brunetti et al. / Journal of Hydrology 540 (2016) 1146–1161

Fig. 3. Scatter plots for pair relations a1-NSE (left) and n1-NSE (right) for Scenario I. The red line is a regression line. (For interpretation of the references to color in this figurelegend, the reader is referred to the web version of this article.)

Fig. 4. Bivariate KDE plots (below diagonal), univariate KDE plots (diagonal), and correlation plots (above diagonal) for Scenario I.

The maximum Pearson correlation coefficient (in absolute values) was 0.42 between parameters a1 and n1. It is also evident from
Fig. 4 that a moderate negative correlation is present for parameters hs1-a1, and a positive correlation for parameters hs1- n1, while
for the other parameters, the correlation is negligible. The univariate KDEs for parameters hs1, n1, and Ks1 indicate a platykurtic distribution of behavioral solutions characterized by multimodality.
Parameter a1 exhibits a leptokurtic distribution, for which a denser
region of good solutions is clearly identifiable in the range of
0.001–0.1.

This behavior is more clear in the bivariate KDE plots. The
bivariate KDE for a1-n1 highlights the presence of a denser region
for values of n1 in the range 2.5–4.5, a behavior that was not evident from the univariate KDE. The comparison between univariate
and bivariate plots reveals that the latter gives a much more comprehensive description of the response surface. This aspect is exacerbated when the model is governed by interactions between
parameters, which is clearly highlighted by Sobol’s indices. In sucha case, the high-dimensional inspection of the parameter spaceprovides significant information.

G. Brunetti et al. / Journal of Hydrology 540 (2016) 1146–1161

The saturated hydraulic conductivity, Ks1, for which the univariate KDE indicates a multimodal behavior, exhibits a denser region
in the range of 10.0–20.0; this region is clearly identifiable in thebivariate plot of K1- n1.
3.3. Sensitivity analysis – Scenario II
Results of the GSA for Scenario II are reported in Table 6. Also
for Scenario II, parameters a1 and n1 exhibit the highest main
effects on the output’s variance (about 35%). For Scenario II, the differences are even more evident than for Scenario I. Parameters n1
and a1 have a first-order index of 30% and 5%, respectively, while
all remaining parameters remain well under 5%. Nine parameters
have a first-order index lower than 1%. The main effects represent
53% of the output variance, which clearly indicates both that themodel output is again (similarly as for Scenario I) partially driven
by interactions between parameters, and that the model is nonadditive.
The right part of Table 6 (ST) shows that the output variance islargely influenced by n1, either directly (30%) or by interactions
with other parameters (64%). Similar to Scenario I, parameters a1
Table 6First-order (S1) and total (ST) effect indices (in decreasing order) with their bootstrap
confidence intervals (BCI) for parameter of Scenario II.
Parameter

S1

S1 (BCI)

ParameterST

ST (BCI)n1 [–]
a1 [1/cm]
hs3 [–]
n3 [–]
Ks3 [cm/min]
a4 [1/cm]
hs2 [–]
hs4 [–]
a3 [1/cm]
Ks2 [cm/min]
Ks4 [cm/min]
hs1 [–]
hs,im3 [–]a2 [1/cm]
Ks1 [cm/min]
x4 [1/min]n2 [–]
x3 [1/min]hs,im4 [–]n4 [–]
Sum

0.302
0.054
0.030
0.026
0.018
0.018
0.017
0.014
0.013
0.012
0.011
0.007
0.006
0.005
0.001
0.001
0.001
0.003
0.004
0.006
0.534

0.026
0.029
0.045
0.024
0.022
0.020
0.018
0.022
0.025
0.026
0.031
0.023
0.017
0.016
0.022
0.027
0.026
0.021
0.020
0.017

n1 [–]
a1 [1/cm]
n3 [–]
hs3 [–]
a3 [1/cm]
hs1 [–]
a4 [1/cm]
Ks1 [cm/min]n4 [–]
a2 [1/cm]
Ks3 [cm/min]n2 [–]
hs4 [–]
Ks2 [cm/min]
Ks4 [cm/min]
hs2 [–]
hs,im3 [–]x3 [1/min]hs,im4 [–]x4 [1/min]

0.640
0.387
0.383
0.294
0.291
0.271
0.269
0.259
0.256
0.229
0.222
0.217
0.201
0.195
0.186
0.185
0.149
0.143
0.138
0.125
>1.0

0.023
0.027
0.020
0.027
0.022
0.019
0.019
0.018
0.013
0.017
0.017
0.022
0.017
0.023
0.021
0.018
0.016
0.016
0.013
0.020

1155

and n1 are the most influential parameters, and the model’s output
is mainly driven by the wear layer. Four of the first eight mostinfluential parameters are related to the wear layer. The main difference between Scenarios I and II is the influence of the base andsub-base layers on the model’s output. This is evident from Fig. 5,
in which the average ST for each layer is reported for both scenarios. For both scenarios, modeling results are most sensitive to thewear layer, which strongly influences the output’s variance. However in Scenario II, the influence of the wear layer is partiallyreduced and redistributed to other layers. It is evident that the
adoption of the dual-porosity model for the unsaturated hydraulicproperties significantly affects the influence of the base and subbase layers on the model’s output. The dynamics of sensitivity
indices between the two scenarios suggest that the physicaldescription of unsaturated flow in the sub-base layer is an important element in numerical simulations.
Similar to scenario I, all parameters influence the model’s output, either by the variation of the parameters themselves or by
their mutual interactions. The condition for FF is never achieved
for all parameters.
Scatter plots for the plain Monte Carlo runs for the two mostsensitive parameters, a1 and n1, are displayed in Fig. 6. It can beseen that there is again no clear pattern of factors driving bad solutions. The regression lines indicate that there is a slight trend,
which is higher for parameter a1, to have admissible solutions in
the left part of the plot. The optimum appears flat, however.
3.4. Monte Carlo filtering – Scenario II
A Monte Carlo Filtering procedure was again applied to the runsof the GSA. The filtered sample now consisted of 28,107 behavioral
solutions. The filtered sample of behavioral solutions for Scenario II
was considerably larger than for Scenario I. This indicates that theimplementation of the dual-porosity model leads to higher valuesof the objective function.
Fig. 7 shows the univariate and bivariate KDE plots as well as
the correlation plots for parameters of the wear layer. It is evident
that no clear correlation exists between various parameters(Fig. 7), except for a negative correlation trend between parameters a1 and n1, but only with a small magnitude. The maximum correlation coefficient, in absolute values, was 0.531 between
parameters a4 and n4.
The univariate KDE for parameters hs1-Ks1 indicates a platykurtic distribution of behavioral solutions without a clear identification of a denser region across the parameter space. On the other
hand, for parameters a1 and n1, the univariate KDEs indicate a more

Fig. 5. The average total index, ST, for different layers for both scenarios.

1156

G. Brunetti et al. / Journal of Hydrology 540 (2016) 1146–1161

Fig. 6. Scatter plots for pair relations a1-NSE (left) and n1-NSE (right) for Scenario II. The red line is a regression line. (For interpretation of the references to color in this figurelegend, the reader is referred to the web version of this article.)

Fig. 7. Bivariate KDE plots (below diagonal), univariate KDE plots (diagonal), and correlation plots (above diagonal) for Scenario II.

leptokurtic distributions, especially for n1, for which a denser
region of solutions between 1.1 and 2.8 is identifiable.
The bivariate KDEs give a better description of the location of
behavioral regions in the bidimensional parameter space than theunivariate KDEs. The bivariate KDE for the two most sensitive
parameters, a1 and n1, indicate the presence of a denser region in
the range of n1 = (1.1, 2.8), and a1 = (0.01, 0.15). The bivariate plots,
hs1- a1 and hs1- n1, indicate the presence of a denser region in therange of hs1 = (0.25, 0.40), a region that was not clearly indicatedby the univariate plot for hs1. The saturated hydraulic conductivity,

Ks1, exhibits a multimodal behavior characterized by several
potential regions of interest. A potential behavioral region may be
identified in the range of Ks1 = (7.0, 15.0).
3.5. Particle swarm optimization
The results and conclusions from the coupled GSA-Monte Carlofiltering analysis were used to reduce the ranges of parameters for
the PSO. The reduction was applied only for parameters that exhibited well identifiable behavioral regions in multivariate plots. The

G. Brunetti et al. / Journal of Hydrology 540 (2016) 1146–1161

original ranges were kept for parameters that displayed high multimodality, in order to avoid the convergence of PSO to the local
optimum. Table 7 reports the new ranges for all parameters.
Fig. 8 compares measured and modeled hydrographs for thetwo scenarios. The PSO for Scenarios I and II resulted in NSE valuesof 0.43 and 0.81, respectively. Both NSE values of the objectivefunction are higher than zero and thus admissible (Moriasi et al.,2007). However, the implementation of the dual-porosity modelfor the base and sub-base layers in Scenario II provides a moreaccurate description of the hydraulic behavior of the permeable
pavement. In particular, the dual-porosity model is able to accurately reproduce the fast hydraulic response of the permeable

Table 7
Reduced ranges of optimized parameters for the optimization process.
Parameter

Scenario I
Reduced range

Scenario II

hs1 [–]
a1 [1/cm]
n1 [–]
Ks1 [cm/min]
hs2 [–]
a2 [1/cm]
n2 [–]
Ks2 [cm/min]
hs3 [–]
a3 [1/cm]
n3 [–]
Ks3 [cm/min]hs,im3 [–]x3 [1/min]
hs4 [–]
a4 [1/cm]
n4 [–]
Ks4 [cm/min]hs,im4 [–]x4 [1/min]

0.2–0.50.001–0.1
3.0–4.5
10.0–20.0
0.25–0.5
0.2–0.3
1.1–4.5
1.0–20.0
0.20–0.40
0.001–0.05
1.1–4.5
1.0–100.0––
0.01–0.2
0.15–0.3
2.0–4.0
1.0–100.0––

0.2–0.40.001–0.15
1.1–2.8
1.5–20.0
0.2–0.4
0.1–0.2
1.1–4.5
3.0–20.0
0.001–0.05
0.001–0.05
1.5–4.5
20.–100.0
0.2–0.40.00001–0.009
0.001–0.05
0.15–0.3
1.5–3.5
1.0–100.0
0.15–0.3
0.00001–0.009

1157

pavement and the long-tailing behavior of the measured hydrograph. The modeled hydrograph for Scenario I appears less accurate in reproducing the dynamics of the observed hydrograph,
especially the fast response of the pavements to precipitation.
Optimized parameters for the two scenarios are reported inTable 8. Significant differences emerge between the two scenarios
in terms of estimated values of the saturated water contents, hs1
and a1; differences between estimated values of the saturated
hydraulic conductivities, Ks1 and n1, are less pronounced. For layer
2, while estimated values of saturated water contents are very similar, huge differences arise between estimated pore-size distribution indices, n2, which for Scenario I is less than half of its value
for Scenario II. Also, Ks2 is considerably lower for Scenario II thanfor Scenario I. Estimated values of dual-porosity parameters confirm the assumptions made about the fractured nature of the baseand subbase layers. While the saturated water content for themobile domain is very low, the porous matrix possesses a highstorage capacity as indicated by the large value of the immobile
saturated water content. In particular, the overall porosity of thebase layer is about 40% and 30% for the subbase layer. The result
for the base layer is slightly higher than the prescriptions of ICPI,
which recommends a porosity of 30–35%. The estimated porosity
for the subbase layer is 30%, which is lower than the prescribed
porosity of about 40%. This difference can be related to the simplifications made in the mobile-immobile dual porosity model for thedescription of preferential flow and uncertainties related to theeffective graduation of the stone material used. However, the significant increase in the accuracy between the single-porosity
model and the dual-porosity model suggests that the hydraulicbehavior of the base and subbase layers is strongly affected by fast
preferential flows in interconnected fractures and the accumulation of water in the rock matrix. This behavior is in agreement withresults reported in the literature. For example, Illgen et al. (2007)
reported that the water contents in the base and sub-base layers

Fig. 8. Comparison between the modeled and measured hydrographs for Scenarios I (top) and II (bottom) for the optimization process.

1158

G. Brunetti et al. / Journal of Hydrology 540 (2016) 1146–1161

Table 8Optimized soil hydraulic parameters for both scenarios.
Layer

hr

hs

a

N

Ks

L

hr,im

hs,im

x

Scenario I
Wear
Bedding
Base
Sub-base
Protection

0.045
0.03
0
0
0.03

0.2
0.3
0.2
0.01
0.3

0.002
0.3
0.023
0.27
0.3

3.0
4.47
2.85
2.41
4.47

10
20
68.7
96.7
20

0.5
0.5
0.5
0.5
0.5

–
–
–
–
–

–
–
–
–
–

–
–
–
–
–

Scenario II
Wear
Bedding
Base
Sub-base
Protection

0.045
0.03
0
0
0.03

0.287
0.298
0.044
0.001
0.298

0.03
0.113
0.021
0.247
0.113

2.67
3.04
4.33
2.17
3.04

7.33
3.87
93.2
56.3
3.87

0.5
0.5
0.5
0.5
0.5

–
–
0
0
–

–
–
0.35
0.29
–

–
–
0.000170.0013–

only marginally increased during rainfall events, and that thelower layers act as a storage tank.
Both scenarios exhibit low values of porosity for the base andsub-base layers. For scenario II, the total porosity is divided
between the mobile and immobile domains. Flow is restricted only
to highly conductive and interconnected fractures, which represent
a relatively small part of the domain, while the immobile domain
provides the storage capacity. While Scenario II assumes overlapping and interacting continua, Scenario I assumes a single continuum approach for all layers. When the optimized value ofporosity is very low, such as for the sub-base layer, it is necessary
to interpret the optimized values differently than for typical
Richards’ type flow. In such case, especially for flow in crushedstones, the model tends to approximate a combination of film flow
and fingering that likely occur in this layer. This hydraulic behavior
is similar to the one reported, for example, by Hodnet and Bell(1990) for unsaturated flow in a medium composed largely of chalkcobbles. In their study, Tokunaga and Wan (1997) analyzed theinfluence of film flow on unsaturated flow in fractures. High velocities of film flow measured in their study suggested that film flow
contributed significantly to preferential flow in fractured rocks.
Our model, based on a macroscopic description of this fast unsaturated flow, shares some similarities with the active fracture model
proposed by Liu et al. (1998). This approach divides the pore space
into two parts, active and inactive. Flow and transport occurs only
within the active pore space, with the inactive part simply
bypassed. Liu et al. (1998) further assumed that van Genuchten(1980) relations are approximately valid for the active pore space.
In a separate study, Liu et al. (2003) reports values of porosity
between 0.01 and 0.03 for the pore space used with the active fracture model.

3.6. Confidence regions
Since parameter estimation involves a variety of possible errors,
including measurement errors, model errors, and numerical errors,
an uncertainty analysis of the optimized parameters constitutes animportant part of parameter estimation. In order to evaluate theuncertainty associated with the estimated parameters, a confidence region around the best solutions optimized with PSO were
calculated using HYDRUS-1D. HYDRUS-1D uses the linear approximation method to identify the confidence region around estimated parameters b, resulting in ellipsoid contours centered at b.
Although restrictive and only approximately valid for nonlinear
problems, an uncertainty analysis provides a means to compare
confidence intervals between parameters, thereby indicating
which parameters should be independently measured or estimated. Confidence intervals have been calculated using the Student’s t distribution with a confidence level of 95%. It is evident

from Table 9 that confidence intervals are narrower for Scenario
II, and that the most uncertain parameters are the saturated
hydraulic conductivities for different layers.

3.7. Model validation
In order to evaluate the reliability of the estimated parameters,
the model has been validated on another independent set of experimental data. Fig. 9 shows a comparison between measured andmodeled hydrographs for the two scenarios during the validation
period.
The value of the objective functions are NSE = 0.43 for Scenario I
and NSE = 0.86 for Scenario II. For Scenario I, the value of the objective function remains the same, which confirms the reliability of
the calibrated model. Although the simulated hydrograph provides
an overall sufficiently accurate description of the hydraulic behavior of the pavement, it is less accurate during rainfall events, which
may be a time period of main interest. For Scenario II, the value ofthe objective function actually increased and reached the value
NSE = 0.86, which is very high and reflects the accuracy of the modeled hydrograph. Also the description of the hydraulic behavior of
the pavement during rainfall events is optimal. This capability of
the calibrated model is important when dealing with the analysis
of combined traditional drainage systems and LID techniques. A
correct description of the hydrograph during precipitation gives
information about the lag time and the intensity of peak flow,

Table 9Confidence intervals (CI) for optimized parameters for both scenarios.
Parameter

hs1 [–]
a1 [1/cm]
n1 [–]
Ks1 [cm/min]
hs2 [–]
a2 [1/cm]
n2 [–]
Ks2 [cm/min]
hs3 [–]
a3 [1/cm]
n3 [–]
Ks3 [cm/min]hs,im3 [–]x3 [1/min]
hs4 [–]
a4 [1/cm]
n4 [–]
Ks4 [cm/min]hs,im4 [–]x4 [1/min]

Scenario I

Scenario II

Value

CI

Value

CI

0.2
0.002
3
10
0.3
0.3
4.47
20
0.2
0.023
2.85
68.73
–
–
0.01
0.27
2.41
96.7
–
–

0.057
0.0006
0.783
4.4
0.072
0.062
1.08
6.5
0.048
0.003
0.537
19.7
–
–
0.002
0.018
0.121
9.2
–
–

0.287
0.029
2.67
7.33
0.29
0.11
3.04
3.87
0.044
0.021
4.33
93.2
0.35
0.00017
0.001
0.247
2.17
56.3
0.288
0.00130.007
0.0008
0.058
0.272
0.009
0.002
0.052
0.084
0.001
0.0005
0.139
3.172
0
0.000003
0.00003
0.004
0.039
1.051
0
0.00002

G. Brunetti et al. / Journal of Hydrology 540 (2016) 1146–1161

1159

Fig. 9. Comparison between the modeled and measured hydrograph for the two scenarios for the validation period.

which are fundamental for both a comprehensive hydraulic analysis of drainage systems, and for the evaluation of benefits of LIDs
implementation. The initial part of the hydrograph appears to beunderestimated, which may be related to the influence of theunknown initial conditions. The model was not able to reproduce
outflow induced by the precipitation event on March 15. Thismay be related to an overestimation of potential evaporation calculated using a literature value of albedo, which could result inan overestimation of the storage capacity of the pavement at the
beginning of the precipitation event, which had a total volume of
6 mm. As a result, the model predicted that the pavement retained
all the precipitation volume. A better characterization of evaporation could help in increasing the accuracy of the model, which isalready high.
Fig. 10 directly compares the measured outflows with those calculated by the two modeling scenarios. The red1 bisector line represents conditions when modeled and measured outflows are
perfectly matched. Linear regression lines are reported for both scenarios. Since the Scenario I tends to overestimate the outflow fluxes,
the difference between the bisector and the linear regression line
(gray) for scenario I is substantial. On the other hand, Scenario II
tends to only slightly underestimate the outflow fluxes, and thus
the slopes of the bisector and the linear regression line (black) for
Scenario II are similar. The simulated hydrographs for both scenarios
tend to introduce some bias in the estimation of peak flows. This
aspect is related to the choice of the NSE as the objective function
for the optimization. The NSE is focused on the general behavior of
the hydrograph rather than on particular components such as peakflows. A multi-objective optimization that would include an objective function targeted to peak flow estimates could represent a more
appropriate approach if estimates of peak flows were the main goalof calibration. However, even of great interest, the multi-objective

1
For interpretation of color in Fig. 10, the reader is referred to the web version ofthis article.

Fig. 10. Comparison between the modeled and measured outflows for the two
scenarios for the validation period.

optimization is out of the scope of this paper. Overall, the validation
process demonstrated the reliability of the calibrated models forboth scenarios.

4. Conclusions
In this paper, we investigated the suitability of the mechanistic
model, HYDRUS-1D, to correctly describe the hydraulic behavior ofa permeable pavement installed at the University of Calabria. We
considered two different scenarios in describing the system. In Scenario I, we assumed that flow on all layers can be described using a

1160

G. Brunetti et al. / Journal of Hydrology 540 (2016) 1146–1161

single-porosity model, while in Scenario II, we assumed that adual-porosity mobile-immobile model is needed to describe flowin the base and subbase layers. The widely used Nash-Sutcliffe efficiency index was used to assess the models. A Global Sensitivity
Analysis, coupled with a Monte Carlo filtering procedure, was carried out before the model calibration. Sensitivity analysis results
suggested that the model is non-additive and mainly driven byparameter interactions in both scenarios. The first-order effects
only accounted for 56% of output variance for Scenario I and 53%
for Scenario II. Sensitivity analysis also revealed that the wear layer
mainly influenced the hydraulic behavior of the pavement. A subsequent Monte Carlo filtering procedure was applied to the runs
performed during the sensitivity analysis in order to identify the
behavioral regions and to reduce parameter uncertainty. Both univariate and bivariate Kernel Density Estimation plots were used toinspect the response surfaces and identify the behavioral regions.
This analysis revealed the high multimodality of the response surfaces, which suggested the use of a global optimization algorithm
for parameter estimation. Correlation coefficients of the filteredsample were also computed, indicating a general low correlation
between parameters. Based on the results of the Monte Carlo filtering, a heuristic global optimization method based on the Particle
Swarm algorithm was used for parameter estimation. The calibrated model for Scenario I exhibited an optimum NSE = 0.43,
while for Scenario II, it reached NSE = 0.81. The optimized parameters were then validated against an independent set of experimental data, resulting in NSE = 0.43 for Scenario I and NSE = 0.86 for
Scenario II. The results of optimization and validation clearly indicated that the implementation of the dual-porosity model for thebase and subbase layers produced more accurate results than the
single-porosity model and described much better the hydraulicbehavior of pervious pavement. Results also confirmed the validityof the assumption that the hydraulic behavior of the base and subbase layers was similar to the behavior of a fractured rock, which ischaracterized by the highly permeable interconnected fractures
and the highly storative rock matrix. The main advantage in using
a simple, dual-porosity, mobile-immobile model with a saturationbased mass transfer is that this model requires only two additionalparameters compared to the single-porosity model. Further significant improvements could be obtained by characterizing thehydraulic properties of the wear layer in the laboratory, as suggested by the sensitivity analysis.
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” Action II: ‘‘Interventions to support industrial research”.
References
Abbaspour, K.C., Johnson, C.A., van Genuchten, M.T., 2004. Estimating uncertain
flow and transport parameters using a sequential uncertainty fitting procedure.
Vadose Zone J. 3, 1340–1352. http://dx.doi.org/10.2113/3.4.1340.
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.
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.

Bai, M., Ma, Q., Roegiers, J.-C., 1994. A nonlinear dual-porosity model. Appl. Math.
Model. 18, 602–610. http://dx.doi.org/10.1016/0307-904X(94)90318-2.
Barenblatt, G., Zheltov, I., Kochina, I., 1960. Basic concepts in the theory of seepageof homogeneous liquids in fissured rocks [strata]. J. Appl. Math. Mech. 24,
1286–1303. http://dx.doi.org/10.1016/0021-8928(60)90107-6.
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.
Beven, K., 2006. A manifesto for the equifinality thesis. J. Hydrol. 320, 18–36. http://dx.doi.org/10.1016/j.jhydrol.2005.07.007.
Beven, K., Binley, A., 1992. The future of distributed models: model calibration anduncertainty prediction. Hydrol. Process. 6, 279–298. http://dx.doi.org/10.1002/
hyp.3360060305.
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.
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.
Carbone, M., Brunetti, G., Piro, P., 2015a. Modelling the hydraulic behaviour ofgrowing media with the explicit finite volume solution. Water 7, 568–591.
http://dx.doi.org/10.3390/w7020568.
Carbone, M., Turco, M., Brunetti, G., Piro, P., 2015b. A cumulative rainfall function
for subhourly design storm in Mediterranean urban areas. Adv. Meteorol. 2015,
1–10. http://dx.doi.org/10.1155/2015/528564.
Cheviron, B., Coquet, Y., 2009. Sensitivity analysis of transient-MIM HYDRUS-1D:case study related to pesticide fate in soils. Vadose Zone 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).
Davis, A.P., 2008. Field performance of bioretention: hydrology impacts. J. Hydrol.
Eng. 13, 90–95. http://dx.doi.org/10.1061/(ASCE)1084-0699(2008) 13:2(90).
DeWalle, D.R., Swistock, B.R., Johnson, T.E., McGuire, K.J., 2000. Potential effects ofclimate change and urbanization on mean annual streamflow in the UnitedStates. Water Resour. Res. 36, 2655–2664. http://dx.doi.org/10.1029/
2000WR900134.
Draper, N., Smith, H., 1981. Applied regression analysis. Technometrics, second ed.
http://dx.doi.org/10.1198/tech.2005.s303.
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. Modell. 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.
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.
Gill, M.K., Kaheil, Y.H., Khalil, A., McKee, M., Bastidas, L., 2006. Multiobjective
particle swarm optimization for parameter estimation in hydrology. Water
Resour. Res. 42. http://dx.doi.org/10.1029/2005WR004528, n/a–n/a.
Gironás, J., Roesner, L.a., Rossman, L.a., Davis, J., 2010. A new applications manual
for the Storm Water Management Model (SWMM). Environ. Model. Softw. 25,
813–814. http://dx.doi.org/10.1016/j.envsoft.2009.11.009.
Hilten, R.N., Lawrence, T.M., Tollner, E.W., 2008. Modeling stormwater runoff fromgreen roofs with HYDRUS-1D. J. Hydrol. 358, 288–293. http://dx.doi.org/
10.1016/j.jhydrol.2008.06.010.
Hodnet, M.G., Bell, J.P., 1990. Processes of water movement through a chalk coombe
deposit in Southeast England. Hydrol. Process. 4, 361–372. http://dx.doi.org/
10.1002/hyp.3360040406.
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.
Hornberger, G.M., Spear, R.C., 1981. An approach to the preliminary analysis ofenvironmental systems. J. Environ. Manage. 12, 7–12. http://dx.doi.org/
10.1016/S0272-4944(89)80026-X.
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.
Illgen, M., Harting, K., Schmitt, T.G., Welker, A., 2007. Runoff and infiltration
characteristics of pavement structures – review of an extensive monitoring
program. Water Sci. Technol. 56, 133–140. http://dx.doi.org/10.2166/
wst.2007.750.

G. Brunetti et al. / Journal of Hydrology 540 (2016) 1146–1161
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.
Jiang, Y., Liu, C., Huang, C., Wu, X., 2010. Improved particle swarm algorithm forhydrological parameter optimization. Appl. Math. Comput. 217, 3207–3215.
http://dx.doi.org/10.1016/j.amc.2010.08.053.
Kennedy, J., Eberhart, R., 1995. Particle swarm optimization. Eng. Technol., 1942–
1948
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.
Kundzewicz, Z., Radziejewski, M., Pínskwar, I., 2006. Precipitation extremes in thechanging climate of Europe. Clim. Res. 31, 51–58. http://dx.doi.org/10.3354/
cr031051.
Levinson, R., Akbari, H., 2002. Effects of composition and exposure on the solarreflectance of Portland cement concrete. Cem. Concr. Res. 32, 1679–1698.
http://dx.doi.org/10.1016/S0008-8846(02)00835-9.
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.
Liu, H.H., Doughty, C., Bodvarsson, G.S., 1998. An active fracture model forunsaturated flow and transport in fractured rocks. Water Resour. Res. 34,
2633–2646. http://dx.doi.org/10.1029/98WR02040.
Liu, H.-H., Haukwa, C.B., Ahlers, C.F., Bodvarsson, G.S., Flint, A.L., Guertal, W.B., 2003.
Modeling flow and transport in unsaturated fractured rock: an evaluation of thecontinuum approach. J. Contam. Hydrol. 62–63, 173–188. http://dx.doi.org/
10.1016/S0169-7722(02)00170-5.
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.
Min, S.-K., Zhang, X., Zwiers, F.W., Hegerl, G.C., 2011. Human contribution to moreintense precipitation extremes. Nature 470, 378–381. http://dx.doi.org/
10.1038/nature09763.
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.
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.
Palla, A., Gnecco, I., Lanza, L.G., 2009. Unsaturated 2D modelling of subsurface waterflow in the coarse-grained porous matrix of a green roof. J. Hydrol. http://dx.doi.
org/10.1016/j.jhydrol.2009.10.008.
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.
Poli, R., Kennedy, J., Blackwell, T., 2007. Particle swarm optimization. Swarm Intell.
1, 33–57. http://dx.doi.org/10.1007/s11721-007-0002-0.
Price, K., 2011. Effects of watershed topography, soils, land use, and climate on
baseflow hydrology in humid regions: a review. Prog. Phys. Geogr. 35, 465–492.
http://dx.doi.org/10.1177/0309133311402714.
Ratto, M., Tarantola, S., Saltelli, A., 2001. Sensitivity analysis in model calibration:
GSA-GLUE approach. Comput. Phys. Commun. 136, 212–224. http://dx.doi.org/
10.1016/S0010-4655(01)00159-X.
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.
Rose, S., Peters, N.E., 2001. Effects of urbanization on streamflow in the Atlanta area
(Georgia, USA): a comparative hydrological approach. Hydrol. Process. 15,
1441–1457. http://dx.doi.org/10.1002/hyp.218.
Saltelli, A., 2002. Making best use of model evaluations to compute sensitivity
indices. Comput. Phys. Commun. 145, 280–297. http://dx.doi.org/10.1016/
S0010-4655(02)00280-1.
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.
Saltelli, A., Annoni, P., Azzini, I., Campolongo, F., Ratto, M., Tarantola, S., 2010.
Variance based sensitivity analysis of model output. Design and estimator for1161

the total sensitivity index. Comput. Phys. Commun. 181, 259–270. http://dx.doi.
org/10.1016/j.cpc.2009.09.018.
Saltelli, A., Tarantola, S., 2004. Sensitivity Analysis in Practice: A Guide to Assessing
Scientific Models 232.
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.
Scott, D.W., 1992. Frontmatter. Wiley Series in Probability and Mathematical
Statistics Applied Probability and Statistics Section. http://dx.doi.org/10.1002/
9780470316849.fmatter.
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. 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. Proc.
1999
Congr.
Evol.
Comput.,1945–1950
http://dx.doi.org/10.1109/
CEC.1999.785511.
Silverman, B., 1981. Using kernel density estimates to investigate multimodality. J.
R. Stat. Soc. 43, 97–99.
Simmons, D.L., Reynolds, R.J., 1982. Effects of urbanization on base flow of selectedsouth-shore streams, Long Island, New York. J. Am. Water Resour. Assoc. 18,
797–805. http://dx.doi.org/10.1111/j.1752-1688.1982.tb00075.x.
Simunek, J., Jarvis, N.J., van Genuchten, M.T., Gardenas, A., 2003. Review andcomparison of models for describing non-equilibrium and preferential flow andtransport in the vadose zone. J. Hydrol. 272, 14–35. http://dx.doi.org/10.1016/
S0022-1694(02)00252-4.
Šimůnek, J., van Genuchten, M.T., 2008. Modeling nonequilibrium flow andtransport processes using HYDRUS. Vadose Zone J. 7, 782. http://dx.doi.org/
10.2136/vzj2007.0074.
Šimůnek, J., van Genuchten, M.T., Šejna, M., 2008. Development and applications of
the HYDRUS and STANMOD software packages and related codes. Vadose Zone
J. 7, 587. http://dx.doi.org/10.2136/vzj2007.0077.
Sobol0 , 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.
Song, X., Zhang, J., Zhan, C., Xuan, Y., Ye, M., Xu, C., 2015. Global sensitivity analysis
in hydrological modeling: review of concepts, methods, theoretical framework,and applications. J. Hydrol. 523, 739–757. http://dx.doi.org/10.1016/j.
jhydrol.2015.02.013.
Tokunaga, T.K., Wan, J., 1997. Water film flow along fracture surfaces of porous rock.
Water Resour. Res. 33, 1287–1295. http://dx.doi.org/10.1029/97WR00473.
Trimble, S.W., 1997. Contribution of stream channel erosion to sediment yield froman urbanizing watershed. Science 278, 1442–1444. http://dx.doi.org/
10.1126/science.278.5342.1442 (80).
Usher, W., Xantares, Hadka, D., Bernardoct, Fernando, Herman, J., Mutel, C., 2015.
SALib: New Documentation, Doc Strings and Installation Requirements. http://dx.doi.org/10.5281/zenodo.31316.
Van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulicconductivity of unsaturated soils 1. Soil Sci. Soc. Am. J. 44, 892–898. http://dx.
doi.org/10.2136/sssaj1980.03615995004400050002x.
Van Genuchten, M.T., Wierenga, P.J., 1976. Mass transfer studies in sorbing porous
media I. Analytical solutions. Soil Sci. Soc. Am. J. 40, 473–480.
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 Zone J. 7, 843. http://dx.doi.org/10.2136/vzj2007.0078.
Warren, J.E., Root, P.J., 1963. The behavior of naturally fractured reservoirs. Soc. Pet.
Eng. J. 3, 245–255. http://dx.doi.org/10.2118/426-PA.
Wasko, C., Sharma, A., 2015. Steeper temporal distribution of rain intensity athigher temperatures within Australian storms. Nat. Geosci. 8, 527–529. http://dx.doi.org/10.1038/ngeo2456.
Westra, S., Fowler, H.J., Evans, J.P., Alexander, L.V., Berg, P., Johnson, F., Kendon, E.J.,
Lenderink, G., Roberts, N.M., 2014. Future changes to the intensity andfrequency of short-duration extreme rainfall. Rev. Geophys. 52, 522–555.
http://dx.doi.org/10.1002/2014RG000464.
Whipple, W., DiLouie, J.M., Pytlar, T., 1981. Erosional potential of streams inurbanizing areas. J. Am. Water Resour. Assoc. 17, 36–45. http://dx.doi.org/
10.1111/j.1752-1688.1981.tb02586.x.
Zambrano-Bigiarini, M., Rojas, R., 2013. A model-independent particle swarm
optimisation software for model calibration. Environ. Model. Softw. 43, 5–25.
http://dx.doi.org/10.1016/j.envsoft.2013.01.004.
Zhang, S., Guo, Y., 2015. SWMM simulation of the storm water volume control
performance of permeable pavement systems. J. Hydrol. Eng. 20, 06014010.
http://dx.doi.org/10.1061/(ASCE)HE.1943-5584.0001092.