Published August 25, 2016

Original Research

A Comprehensive Analysis of
the Variably Saturated HydraulicBehavior of a Green Roof in a
Mediterranean Climate
Giuseppe Brunetti,* Jirka Šimůnek, and Patrizia Piro

Core Ideas• A comprehensive numerical analysis
of hydrological processes in a green
roof was conducted.
• The substrate of the green roof demonstrated bimodal behavior.
• The validated model accurately
describes the hydraulic behavior of
the green roof.
• The soil moisture governs the performance of the green roof during
precipitation.

Low-impact developments (LIDs), such as green roofs, have proven to bevaluable alternatives for stormwater management and hydrological restoration. Mechanistic models are reliable and accurate tools for analysis of
the hydrologic behavior of LIDs, yet only a few studies provide a comprehensive numerical analysis of the hydrological processes involved and test
their model predictions against field-scale data. Moreover, more research isneeded to determine the unsaturated hydraulic properties of the substrates
used in LIDs. For these reasons, the aim of this study was to provide a comprehensive description of the hydrological behavior of an extensive green
roof installed at the University of Calabria. The soil hydraulic properties were
determined by using the simplified evaporation method. Both unimodal andbimodal soil hydraulic functions were used in the analysis. The estimatedparameters were then used in the HYDRUS-3D model to simulate a 2-molong period. Precipitation, irrigation, evaporation, and root water uptake
processes were included in the numerical analysis. The values of 0.74 and
0.8 of the Nash–Sutcliffe efficiency index for the model predictions using
unimodal and bimodal functions, respectively, confirmed the good agreement between the modeled and measured outflows. The bimodal model
was able to both accurately reproduce the hydrographs in both dry andwet periods and account for daily fluctuations of soil moisture. Finally, thevalidated model was used to carry out a hydrological analysis of the green
roof and its hydrological performance during the entire simulated period aswell as during single precipitation events.
Abbreviations: AIC, Akaike information criterion; LID, low-impact development; NSE,
Nash–Sutcliffe efficiency; RMSE, root mean square error; VGM, van Genuchten–Mualem.

During the last few decades, the area of impervious surfaces in urban areas has

G. Brunetti and P. Piro, Dep. of Civil
Engineering, Univ. of Calabria, Rende,
CS 87036, Italy; J. Šimůnek, Dep. of Environmental Sciences, Univ. of California,
Riverside, CA 92521. *Corresponding
author (giuseppe.brunetti@unical.it).

Vadose Zone J.
doi:10.2136/vzj2016.04.0032
Received 19 April 2016.
Accepted 14 June 2016.

© Soil Science Society of America

5585 Guilford Rd., Madison, WI 53711 USA.
All rights reserved.

exponentially increased as a consequence of demographic growth. This long-term process
has altered the natural hydrological cycle by reducing the infiltration and evaporation
capacity of urban catchments, while increasing surface runoff and reducing groundwater
recharge. Moreover, 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 global
warming (Kundzewicz et al., 2006; Min et al., 2011).
The combined effects of urbanization and climate change expose urban areas to anincreasing risk of flooding. In this context, urban drainage systems play a fundamentalrole in improving the resilience of cities. In recent years, an innovative approach to landdevelopment known as low-impact development (LID) has gained increasing popularity.
Low-impact development is a “green” approach to storm water management that seeks
to mimic the natural hydrology of a site using decentralized microscale control measures(Coffman, 2002). Low-impact development practices consist of bioretention cells, infiltration wells or trenches, stormwater wetlands, wet ponds, level spreaders, permeable
pavements, swales, green roofs, vegetated filter and buffer strips, sand filters, smaller culverts, and water harvesting systems. Low-impact developments are able to reduce runoff

Vadose Zone Journal | Advancing Critical Zone Sciencevolumes and pollutant loads and increase evapotranspiration.
Green roofs were able to significantly reduce peak rates of stormwater runoff (Getter et al., 2007) and retain rainfall volumes withretention 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 andNorth Carolina (Davis, 2008). Permeable pavements offered great
advantages in terms of runoff reduction (Carbone et al., 2014;
Collins et al., 2008), water retention, and water quality (Brattebo
and Booth, 2003). Considering that rooftops may represent asmuch as 40 to 50% of the total impervious surfaces in urban areas,green roofs are among the key choices for hydrologic restorationand stormwater management.
One of the key factors limiting the widespread use of LIDs is thelack of adequate modeling tools (Elliott and Trowsdale, 2007) that
could be used to design LIDs that function properly for particularclimatic conditions. Low-impact development modeling requires
an accurate description of the hydrological processes involved,
which are multiple and interacting. In recent years, researchers
have focused their attention on applying and developing empirical,
conceptual, and physically based models for LID analysis. In theirreview, Li and Babcock (2014) reported that there were >600 studies published worldwide involving green roofs, with a significant
portion of them related to modeling.
Zhang and Guo (2013) developed an analytical model to evaluate the long-term average hydrologic performance of green roofs.
Local precipitation characteristics were described using probabilistic methods, and the hydrological behavior of the system wasdescribed using mass balance equations. Kasmin et al. (2010)
developed a simple conceptual model of the hydrological behavior of green roofs during a storm event. The model input was thetime series of precipitation events and the output was runoff. The
water content in the green roof at any given time was between field
capacity and the residual water content. Evapotranspiration was
estimated using an empirical relationship accounting for the actual
water contents, the storm event’s characteristics, and the antecedent dry weather period. During a precipitation event, the porous
medium absorbed moisture until field capacity was reached. The
addition of further moisture produced subsurface flow.
She and Pang (2010) developed a physical model that combined an
infiltration module (based on the Green–Ampt equation) and a saturation module (SWMM). The model calculates the water content ina green roof in a stepwise manner from the initiation of precipitation
until saturation. In simulating the hydraulic response of green roofs
to precipitation, an infiltration module is used before field capacity isreached and when no drainage is produced, while a saturation module
is used after field capacity is reached and when drainage is produced.
However, because runoff and infiltration can occur simultaneously
during heavy precipitation, this stepwise approach may not be appropriate for a wide range of precipitation events.

VZJ | Advancing Critical Zone Science

Although analytical and conceptual models represent a viable
alternative to the numerical analysis of green roofs, their use suffers from several limitations. Conceptualization of the physicalprocesses involved often leads to simplification of the system anda reduction in numerical parameters. While in a physical model
each parameter has its own meaning, in conceptual models,
lumped parameters often incorporate different components of thedescribed process. These lumped parameters are case sensitive andneed to be calibrated against experimental data, implying a lack ofgenerality of the model itself. These drawbacks could represent abarrier to the use of modeling tools among practitioners who needreliable and generally applicable models.
Mechanistic models have proven to be a valid and reliable alternative to conceptual and analytical models for the analysis of green
roofs and LIDs in general. Carbone et al. (2015a) developed a onedimensional finite volume model for description of the infiltration
process during rainfall events in green roof substrates. The model
was based on the reduced advective form of the Richards equation, in which the soil water diffusivity was neglected. Metselaar(2012) used the SWAP software (van Dam et al., 2008) to simulate
the one-dimensional water balance of a substrate layer on a flat
roof with plants. Hilten et al. (2008) simulated peak flow andrunoff volume reduction of a 10-cm modular green roof (60 by
60 cm) using HYDRUS-1D (Šimůnek et al., 2008). In that study,
only the values of field capacity and wilting point were measured.
These parameters, in conjunction with the soil bulk density andparticle size distribution, were used to estimate the soil water retention curve using a pedotransfer function. Multiple 24-h storms
were used to generate precipitation data and simulate runoff todescribe the green roof ’s hydrologic response. Li and Babcock(2015) used HYDRUS-2D to model the hydrologic response ofa pilot 61- by 61-cm green roof system. Physical properties of thesubstrate were obtained using laboratory measurements on soilcores extracted from a green roof. The saturated hydraulic conductivity was measured using the falling-head method, while theresidual and saturated water contents were measured using thegravimetric method. The hanging water column method was used
to estimate the shape parameters of the unimodal van Genuchten
function (van Genuchten, 1980). The model was calibrated using
water content measurements obtained with time domain reflectometry sensors. The calibrated model was then used to simulate
the potential beneficial effects of irrigation management on the
reduction of runoff volumes.
Although physically based models have been widely and often
successfully used, very few studies have provided a comprehensive
analysis of the hydrological behavior of a green roof and validated
it against field-scale data. Moreover, studies that investigated theunsaturated hydraulic properties of green roof substrates have been
limited to the determination of some specific soil characteristics(e.g., field capacity, wilting point, or particle size distribution) andgenerally focused only on the soil water retention curves.
p. 2 of 17

For these reasons, the aim of this study was an accurate andcomprehensive analysis of the hydrological behavior of green
roofs using the mechanistic model HYDRUS-3D to analyze an
extensive green roof installed at the University of Calabria. The
problem was addressed in the following way. First, the soil water
retention curve and the unsaturated hydraulic conductivity of thegreen roof substrate were measured using a simplified evaporation
method. The soil hydraulic parameters obtained were then used in
HYDRUS-3D numerical simulations of the green roof function
using precipitation, climate, and subsurface experimental data fora 2-mo-long period. The model was validated by comparing themodeled and measured subsurface flows using the Nash–Sutcliffe
efficiency index (Nash and Sutcliffe, 1970). Finally, the validated
model was used to evaluate the hydrologic behavior of the green
roof and its hydraulic response to single precipitation events.

66Materials and Methods
Green Roof and Site Description

The University of Calabria is located in the south of Italy, inthe vicinity of Cosenza (39°18¢ N, 16°15¢ E). The climate is
Mediterranean, with a mean annual temperature of 15.5°C andan average annual precipitation of 881.2 mm. The green roof ispart of the “Urban Hydraulic Park,” which also includes a permeable pavement, a bioretention system, and a sedimentation
tank connected to a treatment unit. An extensive green roof wasinstalled on the existing rooftop of the Department of MechanicalEngineering. The original impervious roof was divided into four
sectors. Two sectors are vegetated with native plants and differ
from each other by the drainage layer. Another sector is characterized by bare soil with only a few spontaneous plants. The last sector
is the original impervious roof. The maximum depth of the soil
substrate is 8 cm. This depth was selected to investigate both theenergetic (heat fluxes) and hydrologic (water fluxes) behavior of avery thin extensive green roof under the Mediterranean climate.
The soil substrate is composed of mineral soil with 74% gravel, 22%sand, and 4% silt and clay. The soil has a measured bulk densityof 0.86 g cm−3 and 8% organic matter, which was determined in
the laboratory using the Walkley–Black method. Three different

plant species were selected and planted. Cerastium tomentosum L.
and Dianthus gratianopolitanus Vill. are herbaceous plants suited
for well-drained soils; Carpobrotus edulis (L.) L. Bolus is a succulent plant characterized by a high drought tolerance, largely due
to the high leaf succulence and physiological adaptations such as
crassulacean acid metabolism (CAM) photosynthesis (Durhman
et al., 2006); CAM plants have greater water use efficiency thanC3 plants because their transpiration per unit of CO2 is reduced
due to the stomata opening at night for CO2 uptake (Sayed, 2001).
In this study, only one vegetated sector of the green roof was considered. Figure 1 displays a cross-section of the green roof; theconsidered sector has an area of 50 m2 and an average slope of 1%.
The green roof is divided into square elements of 50 by 50 cm (Fig.
2), with alternating vegetated and non-vegetated areas. The substrate has a maximum depth of 8 cm where plants are grown anda minimum depth of 4 cm where no vegetation is present (Fig. 1
and 2). This design was meant to minimize the weight on the green
roof support structure. A highly permeable geotextile is placed at
the bottom of the substrate to prevent soil from migrating into theunderlying layers. The drainage layer is composed of a polystyrene
foam and is characterized by a water storage capacity of 11 L m−2
and a drainage capacity of 0.46 L s−1 m−2 . Water accumulated
in the drainage layer can be transferred back up to the substrate
only by condensation on the geotextile. An anti-root layer and animpervious membrane complete the green roof.
A drip irrigation system was installed to provide water to the plants
during drought periods. The irrigation system is connected witha reuse system, which collects outflow from the green roof. Only
reused water was used for irrigating the green roof. The reusesystem is composed of a storage tank and a pump. When thestorage capacity of the tank (1.5 m3) is exceeded, water is directly
discharged into the drainage system. Drippers are located at thecenter of each square, and their distance from each other is approximately 50 cm. Drippers were also installed in non-vegetated areas
to utilize water from the storage tank by using the evapotranspiration capacity of the green roof. In this way, the volume of water
discharged into the drainage system is reduced, and the evaporative cooling effect of the green roof on the building is expected

Fig. 1. A typical cross-section of the green roof.

VZJ | Advancing Critical Zone Sciencep. 3 of 17

measured using a tipping bucket rain gauge with a resolution of
0.254 mm and an acquisition frequency of 1 min. Climatic data
were acquired with a frequency of 5 min. Data were processed andstored in a SQL database.

Fig. 2. A schematic of the green roof, showing both vegetated (gray)and non-vegetated (white) areas.

to increase. The irrigation system is activated at predefined times
by an electric valve, and the irrigation rate is measured by a water
counter with an acquisition frequency of 1 min. The total volumeof irrigation for the selected time period was 142 mm.
A weather station located directly at the site measured precipitation, velocity and direction of wind, air humidity, air temperature,atmospheric pressure, and global solar radiation. Precipitation was

A flux meter located at the base of the building, composed of a
polyvinyl chloride (PVC) pipe with a sharp-crested weir and apressure transducer, measured outflow from the green roof. The
pressure transducer (Ge Druck PTX1830) measured the waterlevel inside the PVC pipe and had a measurement range of 75 cm,
with an accuracy of 0.1% of the full scale. The pressure transducer
was calibrated in the laboratory using a hydrostatic water column,
linking the electric current intensity with the water level inside
the column. The exponential head-discharge equation for the flux
meters was obtained by fitting the experimental data. The subsurface flow data were acquired with a time resolution of 1 min andstored in a SQLITE database.
A 2-mo data set was selected for analysis (Fig. 3). This particular time period, which started on 1 Sept. 2015 and ended on
30 Oct. 2015, was selected because it involved highly variableclimatic conditions. Isolated precipitation occurred in September,which had a relatively high average temperature. These climaticconditions required irrigation of the green roof for 1 h during thenight. October was characterized by intense and frequent precipitation events. The total recorded precipitation for the whole period
was 431 mm, with an average air temperature of 20.2°C.

Fig. 3. Measured precipitation (black), irrigation (top), and subsurface (gray) fluxes for a selected time period in September and October 2015.

VZJ | Advancing Critical Zone Sciencep. 4 of 17

Hourly reference evapotranspiration was calculated using thePenman–Monteith equation (Allen et al., 1998). An average albedo
value of 0.2 was used in calculations of net shortwave radiation,
assuming that the albedo for vegetated areas was 0.23 (Lazzarin etal., 2005) and 0.17 for bare soil (Rosenberg et al., 1983).

Soil Hydraulic Properties
Evaporation Method

Modeling of water flow in unsaturated soils by means of the
Richards equation requires knowledge of the water retentionfunction, q(h), and the hydraulic conductivity function, K(h), foreach soil layer of the green roof, where q is the volumetric watercontent [L3 L−3], h is the pressure head [L], and K is the hydraulic conductivity [L T−1]. A broad range of methods exists for thedetermination of soil hydraulic properties in the field or in thelaboratory (Arya, 2002; Dane and Hopmans, 2002; Klute andDirksen, 1986). The numerical inversion of transient flow experiments represents one of the most accurate ways to determine soil
hydraulic properties (Šimůnek et al., 1998). Among these, thesimplified evaporation method (Schindler, 1980) is one of themost popular methods. Peters and Durner (2008) conducted acomprehensive error analysis of the simplified evaporation method
and concluded that it is a fast, accurate, and reliable method todetermine soil hydraulic properties in the measured pressure
head range and that the linearization hypothesis introduced by
Schindler (1980) causes only small errors. The evaporation method
was further modified by Schindler et al. (2010a, 2010b) to significantly extend the measurement range to higher pressures. For adetailed description of the modified evaporation method, pleasesee Schindler et al. (2010a, 2010b).
A drawback of the evaporation method is that it remains unsuitedfor direct determinations of conductivities near saturation(Wendroth et al., 1993). The determination of hydraulic conductivities remains reliable only in the dry range, in which hydraulicgradients are more pronounced. To improve the characterizationof the hydraulic conductivity function near saturation, alternative methods are required such as the multistep outflow method(Peters and Durner, 2008).
In this study, the simplified evaporation method with the extended
measurement range (down to −9000 cm) was used for the determination of the unsaturated hydraulic properties of the green roof
substrate. For a complete description of the system, please see UMS
GmbH (2015). The soil for the laboratory analysis was sampleddirectly from the green roof using a stainless-steel sampling ringwith a volume of 250 mL. The soil sample was saturated from thebottom before starting the evaporation test. The measurement unit
and tensiometers were degassed using a vacuum pump to reduce
the potential nucleation sites in the demineralized water. BecausePeters and Durner (2008) suggested a reading interval for structured soils of <0.1 d, the reading interval was set to 20 min to havehigh-resolution measurements. At the end of the experiment, theVZJ | Advancing Critical Zone Sciencesample was placed in an oven at 105°C for 24 h, and then the dryweight was measured.

Parameter Estimation
The numerical optimization procedure HYPROP-FIT (Pertassek
et al., 2015) was used to simultaneously fit retention and hydraulic conductivity functions to the experimental data obtained using
the evaporation method. HYPROP-FIT is a computer program
designed to fit unimodal and multimodal retention functions tomeasured water retention data and to compute the corresponding
relative hydraulic conductivity function. The fitting is accomplishedby a nonlinear optimization algorithm that minimizes the sum of
weighted squared residuals between model predictions and measurements. The software uses the Shuffled Complex Evolution
algorithm proposed by Duan et al. (1992), which is a global parameter estimation algorithm. The software includes a corrected fit of
the hydraulic functions by the “integral method” to avoid bias in thehydraulic properties near saturation (Peters and Durner, 2006), an
Hermitian spline interpolation of the raw measured data to obtainsmooth and continuous time series of measured data, and automatic
detection of the validity range of conductivity data near saturation,
where the hydraulic gradients become too small to yield reliable data.
Two different models were evaluated for the description of soil
hydraulic properties. The unimodal van Genuchten–Mualem(VGM) model (van Genuchten, 1980) was used first:
ìï
1
ïï
ï
n m
Q = ïí éê1 + (a h ) ùú
ïï ë
û
ïï1
ïî
Q=

if h £ 0[1]
if h > 0q-q rq s -q rìïïï K sQ L
K =ï
í
ï
ï
ï
ï
îK s

m = 1-{éêë1-(1-Q )ùúû }
1/ m

m 2

if h < 0[2]

if h > 01nwhere Q is the effective saturation (dimensionless), a is a parameter related to the inverse of the air-entry pressure head [L−1], q s
and q r are the saturated and residual water contents, respectively(dimensionless), n and m are pore-size distribution indices (dimensionless), Ks is the saturated hydraulic conductivity [L T−1], and L
is the tortuosity and pore-connectivity parameter (dimensionless).
Because the unimodal VGM model cannot always describe thefull complexity of measured data, the bimodal model of Durner(1994), which constructs the retention and hydraulic conductivity
functions by a linear superposition of two or more VGM functions,
was used next:
p. 5 of 17

ìï
üï
ïìï
ïï
ïï
ïï 2
1
ïý
ïï å wi ïí
ni ù mi ï
Q = í i=1 ïï éïïïï ê1 + (a i h ) ú ïïï
û þï
îï ë
ïï
ï1
ïî

if h < 0[3]2
ìïïï 2 w a é1-(a h )ni Q ù üïïïï æ 2
Lì
å
ö
ê
ú
i
i
i
i
1
i
=
ï
ï
ïï ç
ë
ûï
K ç w i Q i ÷÷÷ ïí
ý if h < 0 [4]2÷
K = ïí s çèç å
ï
ï
ø
w
a
i =1
ï
ï
å
i
i
ïï
1
i
=
ï
ï
ï
ï
î
þ
ïï
ïï K sif h > 0î

where w is a weighting factor and i refers to the ith pore system.
Although the Ks value is commonly fixed to the measured valueof the saturated hydraulic conductivity, some studies have showedthat this can introduce bias in the unsaturated hydraulic conductivity function when using the traditional VGM model. Schaapand Leij (2000) and Schaap et al. (2001) confirmed that fixing Ks
to a measured value of the saturated hydraulic conductivity led toa systematic overestimation of the hydraulic conductivity at mostpressure heads. Furthermore, Schaap et al. (2001) demonstrated
that the hydraulic conductivity estimated by fitting Ks provided amuch better description of the hydraulic conductivity at negative
pressure heads than fixing it at the measured saturated hydraulicconductivity. In addition, Schaap and Leij (2000) found that thefitted value of the tortuosity L was often negative, with an optimalvalue of −1. For these reasons, all the parameters were initially
included in the optimization.
The goodness-of-fit was evaluated in terms of the root mean squareerror (RMSE), while the Akaike information criterion (AIC) (Hu,
1987) was used to choose between different hydraulic conductivity functions with different numbers of optimized parameters.
The software also provides 95% confidence intervals to assess theuncertainty in parameter estimation.

Modeling Theory

Water Flow and Root Water Uptake
The HYDRUS-3D software (Šimůnek et al., 2008) was used todescribe the morphological complexity of the green roof, which
simultaneously includes multiple soil depths, both vegetated andnon-vegetated areas, and drip irrigation. The green roof consistsof four square elements, which are regularly repeated (Fig. 2). The
hydrologic response of the entire green roof can be well described
as a superposition of the behavior of these four elements.
HYDRUS-3D is a three-dimensional model for simulating the
movement of water, heat, and multiple solutes in variably saturated
porous media. HYDRUS-3D numerically solves the Richards
equation for multidimensional unsaturated flow:

VZJ | Advancing Critical Zone Science

[6]

S ( h ) = a( h ) S p

if h > 0¶q= Ñ éë K Ñ( h - z )ùû - S¶twhere S is a sink term [L3 L−3 T−1], defined as a volume of water
removed from a unit volume of soil per unit of time due to plantwater uptake. Feddes et al. (1978) defined S as

[5]

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, in
which water uptake is assumed to be zero close to soil saturation(h1) and for pressure heads higher (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 d−1 in the model
simulation), stomata start to close at lower pressure heads (h3) (inabsolute value) than at low potential transpiration rates (1 mm d−1)(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).
As explained above, green roof plants were selected to suit
Mediterranean climate conditions. Hanscom and Ting (1978) conducted a comprehensive experimental campaign on the behavior of
succulent plants under water stress. They observed that during time
periods with water and salt stress, plants closed their stomata and, asa consequence, little or no transpiration occurred even during daylight hours. Thus, the plants were capable of withstanding extended
periods of drought. In the same study, well-watered plants exhibited
normal C3 photosynthesis mechanisms, with the maximum CO2
uptake occurring during the day. This behavior was reported also
by Starry et al. (2014). Considering that the combined effects ofirrigation and precipitation limited the drought periods, it appears
reasonable to assume that a normal C3 mechanism occurred. For
these reasons, the parameters reported by Wesseling et al. (1991) for
pasture were slightly modified in this study. In particular, h1 and h2
were set to −1 and −10 cm, respectively, to increase actual transpiration for near-saturated conditions. Parameters used in the water
stress response function are reported in Table 1.
The local potential root water uptake Sp was calculated from thepotential transpiration rate Tp. Beer’s equation was first used
Table 1. Feddes’ parameters for the water stress response function usedin numerical simulations.
Feddes’ parameter

Pressure head

h1

−1

h2

−10

h3

−200

h4

−800

h5

−8000

cm

p. 6 of 17

to partition reference evapotranspiration, calculated using thePenman–Monteith equation (Allen et al., 1998), into potential
transpiration and potential soil evaporation fluxes (e. g., Ritchie,1972). 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 invegetated areas. For a detailed explanation of evapotranspiration
partitioning, please see Sutanto et al. (2012).
As described above, the vegetated and non-vegetated green roof elements alternate, while plants are located in the center of vegetatedareas. HYDRUS-3D allows the consideration of a spatially variable
root distribution. A cylinder with a radius of 20 cm and a depthof 8 cm, in the center of the vegetated area, was used to model theroot zone. The root density was assumed to be uniform inside thecylinder and zero in the remaining part of the numerical domain.
The total potential transpiration flux from a transport domain is,in HYDRUS, equal to potential transpiration Tp multiplied by
the surface area associated with vegetation. This total potential
transpiration flux is then distributed throughout the entire rootzone for the computation of the actual root water uptake.

Fig. 4. Details of the drainage layer, where d is the thickness of theopen space between drainage holes, and the geotextile supporting thegreen roof substrate.

Numerical Domain and Boundary Conditions

layer has only a limited effect on green roof outflow. For thesereasons, only the soil substrate was modeled in this study.

The two main elements that form a green roof are the soil substrate and the drainage layer. While the role of the substrate iswell known because it governs the dynamics of infiltration andevapotranspiration, the importance of the drainage layer for thehydraulic behavior of a green roof is only partially describedin the literature, especially with respect to the modeling of its
function. The drainage layer is frequently modeled as an open
reservoir (e.g., Locatelli et al., 2014; Vesuviano et al., 2014). Oncethe drainage layer’s storage capacity is reached, the excess wateris drained through holes into outflow drains. This guarantees
a high permeability of the system and avoids the formation of
ponding on top of the substrate layer even for intense precipitation. An open space of 1 cm separates the soil substrate anddrainage holes (Fig. 1 and 4).

While precipitation and potential evaporation (different in vegetatedand bare areas; see Table 2) were uniformly distributed on the soil
surface, the drip irrigation was modeled at predefined surface points.
Drippers can be idealized as point sources with a specified irrigation
flux. However, if the irrigation flux is applied to a single boundary node and this flux exceeds the infiltration capacity of this node,
problems with numerical convergence can occur. To avoid such
numerical problems, the irrigation flux should be distributed across alarger surface area, which should ideally represent the wetting radius.
This area must be large enough to avoid surface ponding. In thisstudy, the irrigation flux was distributed across a circular area witha radius of 5 cm located in the center of each element. As a result, no
ponding was observed during numerical simulations.

Water accumulated in the drainage layer can return to the soil
substrate only by evaporation and subsequent condensation on
the geotextile at the bottom of the soil. In this small air space,potential evaporation is expected to be limited due to the microclimatic conditions to which water in the drainage layer is exposed.
The enclosed air space is expected to be characterized by relatively
high humidity, considering the combined effects of soil moistureand the vicinity of the water table of the drainage layer. Moreover,radiation and air turbulence can be considered negligible in thisenclosed air space. The only factors that can thus produce evaporation are the air temperature and air humidity. However, the above
considerations suggest that the effects of evaporation and microcondensation can be neglected, especially at the field scale. Thisimplies that variations in the water level in the drainage layer are
limited and, consequently, the storage capacity of the drainage

VZJ | Advancing Critical Zone ScienceThe surface of the green roof was thus exposed to precipitation,
evaporation, and irrigation. As a result, three different boundary
Table 2. Fluxes considered under different types of boundary conditions.
Boundary
condition

Flux†

Atmosphericprecipitation, potential evaporation (= ET0f ), and potential
transpiration [= ET0(1 − f )]Variable flux 1

precipitation and potential evaporation

Variable flux 2

precipitation, irrigation, potential evaporation (= ET0f ), andpotential transpiration [= ET0(1 − f )]Seepage face

seepageZero flux

no flux† ET0, reference evapotranspiration; f, distribution coefficient dependent on
the leaf area index (Ritchie, 1972).

p. 7 of 17

conditions were specified at the surface of the modeled domain
and two boundary conditions at its bottom (Fig. 5). Table 2 summarizes the various fluxes considered in the various types of usedboundary conditions.
The “atmospheric” boundary condition, which was assigned toareas under vegetation, can exist in three different states: (i) precipitation and/or potential evaporation fluxes, (ii) a zero pressure
head (full saturation) during ponding when both infiltrationand surface runoff occurs, and (iii) an equilibrium between thesoil surface pressure head and the atmospheric water vapor pressure head when atmospheric evaporative demand cannot be metby the substrate. The threshold pressure head, which was set
to −30,000 cm, divides the evaporation process from the soil
surface into two stages: (i) a constant-rate stage when actual
evaporation, equal to potential evaporation, is limited only by
the supply of energy to the surface, and (ii) the falling-rate stage
when water movement to the evaporating sites near the surface
is controlled by subsurface soil moisture and the soil hydraulicproperties and when actual evaporation, calculated as a resultof the numerical solution of the Richards equation, is smaller
than potential evaporation.
A special option of HYDRUS-3D was used to treat the “variableflux” boundary conditions as the “atmospheric” boundary conditions (i.e., with the limiting pressure heads described above). The
Variable Flux 1 boundary condition included precipitation andpotential evaporation and was assigned to bare soil areas. Because
no vegetation was present in these areas, the reference evapotranspiration was not partitioned as for the “atmospheric” boundary
condition but was fully assigned to potential evaporation. Thisapproach shares some similarities with the “dual” crop coefficient
introduced by Allen et al. (1998):
ETc = ET0 ( K cb + K e )[7]

where ETc is the actual crop evapotranspiration, ET0 is the reference evapotranspiration, Kcb is the basal crop coefficient, and Ke
is the empirical soil evaporation coefficient, which accounts for
multiple factors affecting soil evaporation, such as soil texture andavailable soil moisture. In the case of bare soil, Kcb becomes zero
because no vegetation is present and ETc is related only to the soil
evaporation coefficient (Torres and Calera, 2010). In HYDRUS,
soil evaporation is modeled using the two-stage model with thethreshold pressure head (described above), which directly accounts
for factors affecting soil evaporation and which thus does not
require the use of Ke.
The Variable Flux 2 boundary condition, which involved precipitation, irrigation, and evaporation, was applied to the circular areas
with a radius of 5 cm where drippers were located.
A seepage-face boundary condition was specified at the bottomof the soil substrate under vegetated areas because the geotextile
is exposed to atmospheric pressure. 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.
In non-vegetated elements, a zero flux boundary condition wasapplied, except in small circular areas that represented drainage
holes (Fig. 1 and 4). Three circular areas, each with a radius of
0.5 cm at the bottom of the non-vegetated elements, were modeled
as seepage faces. Considering the occurrence of high nonlinearities
and fluxes around these drainage holes, the finite element mesh
was refined here (to 0.5 cm) to guarantee a good accuracy of thenumerical solution. No-flux boundary conditions were used at theremaining boundaries.
The initial pressure head was assumed to be constant in the entire
domain and was set equal to −330 cm, which is usually assumed tobe the field capacity. The numerical model is expected to be sensitive to the initial condition only during the first few simulated days.

Fig. 5. Spatial distribution of consideredboundary conditions.

VZJ | Advancing Critical Zone Sciencep. 8 of 17

The three-dimensional simulated domain had a surface area of 1 m2 ,a maximum height of 8 cm, and a total volume of 0.06 m3. The
domain was discretized into three-dimensional prismatic elements
using the MESHGEN Plus tool of HYDRUS-3D. No mesh stretching was used, and the finite element (FE) mesh was isotropic. The
generated FE mesh had 10,709 nodes and 49,027 three-dimensional
elements. The quality of the FE mesh was assessed by checking themass balance error reported by HYDRUS-3D at the end of the simulation. Mass balance errors, which in this simulation were always<1%, are generally considered acceptable at these low levels.

Statistical Evaluation
The Nash–Sutcliffe efficiency (NSE) index (Nash and Sutcliffe,1970) was used to evaluate the agreement between measured andmodeled hydrographs:
é Tobs
mod 2 ùê å i=1( Q i - Q i ) úú
NSE = 1- êê Tobs
obs 2 úê å i=1( Q i - Q mean ) ú
ë
û

[8]

where T is the total number of observations, Qiobs is the ith meaobs
sured value, Qimod is the ith simulated value, and Q meanis themean value of the observed data. The NSE index ranges between−¥ and 1.0, is equal to 1 in the case of a perfect agreement, andgenerally, 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 indicator for evaluating the overall fit of ahydrograph (Servat and Dezetter, 1991).

66Results and Discussion
Soil Hydraulic Properties

Soil hydraulic properties measured using the evaporation method
are displayed in Fig. 6. The soil water retention curve is welldescribed across the entire water content range (Fig. 6). The
retention data point close to log(h) = 4 (h in cm) was obtained byusing the air-entry pressure head of the ceramic. At first inspection, the behavior of the retention curve appears not to be perfectly
sigmoidal, which may indicate the presence of a secondary poresystem (Durner, 1994). Measured points of the hydraulic conductivity function are concentrated in the dry range between 10 and
30% of the volumetric water content. This is common when theevaporation method is used to measure soil hydraulic properties ofcoarse-textured soils such as the substrate of the green roof.
The measured data were imported into the HYPROP-FIT software to fit the analytical hydraulic property functions. The
unimodal VGM model (van Genuchten, 1980) was fitted first.
The RMSE values for retention and conductivity functions were0.02 cm3 cm−3 and 0.13 (in log K, cm d−1), respectively. An AIC
of −874 was obtained when L was included in the optimization. The unimodal function introduced a high bias, especiallyin the hydraulic conductivity function. The bimodal Durner(1994) model (Eq. [3–4]) was fitted next. The RMSE values forthe retention and conductivity functions were 0.005 cm3 cm−3)
and 0.07 (in log K, cm d−1), respectively. An AIC of −1298 was
obtained when the value of L was fixed to 0.5, as this is the value
usually assumed in the literature for many soils. Figure 6 displays
a comparison between the measured data and their fit using theunimodal and bimodal retention functions. The estimated soil
hydraulic parameters with their confidence intervals are reported
in Table 3.
It is evident that the bimodal function provides a more accurate description of the retention curve. A significant difference
between unimodal and bimodal functions emerges in the dry
range between log(h) = 4 and 6, where the unimodal function estimates lower water contents. The unimodal function overestimatesFig. 6. Measured and modeled values of the retention curve q[log10(h)] (left) and the hydraulic conductivity functions K[log10(h)] (center) and K(q)(right). The measured values are scatter points; the solid and dashed lines are the fitted bimodal and unimodal functions, respectively.

VZJ | Advancing Critical Zone Sciencep. 9 of 17

Table 3. Estimated soil hydraulic parameters and their confidence intervals (CIs) for the unimodal and bimodal hydraulic functions.
Unimodal

BimodalParameter

Estimation

CI

Estimation

CI

Residual water content (q r)

0

0.05

0.070

0.007

Saturated water content (q s)

0.551

0.01

0.562

0.003

Air-entry pressure head index for the first pore system (a1), cm−1

0.13

0.03

0.843

0.07

Pore-size distribution index for the first pore system (n1)

1.25

0.06

1.24

0.04

Saturated hydraulic conductivity (Ks), cm d−1

4700

3500

12600

3700

Air-entry pressure index for the secondary pore system (a 2), cm−1

–

–

0.01

0.001

Pore-size distribution index for the secondary pore system (n2)

–

–

1.97

0.08

Weight coefficient (w2)

–

–

0.422

0.01

Tortuosity and pore connectivity parameter (L)

0.53

0.02

0.5

–

the water contents in the range between log(h) = 0 and 2 (h incm). This is confirmed by the RMSE values, which are higher for
the unimodal function, especially for the hydraulic conductivity
function, for which the RMSE value is twice as large as for thebimodal function. The AIC value is higher (in absolute values) forthe bimodal function, indicating that the model is better suited todescribe the soil hydraulic properties of the evaluated soil than theunimodal model. Figure 6 displays measured data and the fitted
multimodal retention and conductivity curves.

a generally well-posed problem. On the other hand, the maximum
correlation coefficient for the unimodal model was 0.97 for q r
and n, which indicates ill-posedness of the optimization problem.
Narrow confidence intervals for the parameters q r, q s, n1, a 2 , n2 ,and w2 indicate high confidence in their estimation. The fitted
saturated hydraulic conductivity Ks and the parameter a1 for the
first pore system exhibited the largest uncertainties. As explained
above, the evaporation method is not accurate for the determination of the hydraulic conductivity near saturation, and this fact
is reflected in the estimation of Ks. To improve the accuracy inthe estimation of the hydraulic conductivity near saturation, othermethods should be used.

The change in the slope of the retention curve near saturation isreflected in the hydraulic conductivity function. However, thisdecrease in the hydraulic conductivity is not highly pronounced.
Significant differences between the two hydraulic conductivity
models occur only near saturation and in the dry range. The unimodal model estimates significantly lower hydraulic conductivity
values close to saturation than the bimodal model. While the fitted
saturated water contents q s are almost the same for both models,significant differences exist in the residual water contents q r.
Fitted values of q r are 0.0 and 0.062 for the unimodal and bimodal
models, respectively. The value estimated by the unimodal modelseems unrealistically low, and it could indicate inaccuracy in thedescription of the soil hydraulic properties in the dry range.
Both models indicate a soil characterized by very high permeability, which corresponds well with the textural composition of thegreen roof substrate. This characteristic is well suited for green roof
substrates, which must guarantee fast drainage and avoid water
ponding on the surface even during intense precipitation. The
volumetric water contents corresponding to field capacity and thewilting point were 21 and 10%, respectively, for the bimodal modeland 21 and 8%, respectively, for the unimodal model.
The maximum correlation between optimized parameters forthe bimodal model was 0.88 for q r and n1, a result that is quite
common when only a few points are measured in the dry range.
However, only five correlation coefficients were >0.8, indicating

VZJ | Advancing Critical Zone ScienceThe above discussed analysis suggests that the bimodal model
could provide a slightly better description of the soil hydraulicproperties than the unimodal model.

Model Validation
Parameters obtained with the evaporation method were used in
HYDRUS-3D to describe the soil hydraulic properties of the green
roof substrate. Figure 7 shows a comparison between measured andsimulated outflows from the green roof when using both unimodaland bimodal functions of soil hydraulic properties.
The NSE indices for measured and simulated green roof outflows
were 0.74 and 0.8 when the unimodal and bimodal functions wereused, respectively. These values indicate that both models wereable to accurately describe the hydraulic behavior of the green roof,
while a higher accuracy was achieved by the bimodal model. It isevident from Fig. 7 that the unimodal model failed in reproducing small outflows after irrigation, which were observed at the
beginning and at the end of September. The inserts in Fig. 7 show
simulated vs. measured green roof outflows. The same plot also
shows a bisector line, which indicates a perfect agreement between
simulated and measured outflows, and a linear regression line. The
good performance of the models are confirmed by the coefficientof determination R2 = 0.85 and 0.82 of the linear regressions forp. 10 of 17

Fig. 7. A comparison between measured and simulated outflows vs. time for selected dates in September and October 2015 and against each other (inthe insert). The red and black lines in the insert are bisector and linear regression lines, respectively.

the bimodal and unimodal models, respectively. The comparison
between bisector and regression lines indicates that both models
slightly overestimated the outflow.
A further analysis of the fit between measured and simulated green
roof outflows was performed by analysis of the residuals, which
is displayed in Fig. 8 using a lag plot. The lag plot is constructedby comparing neighboring residuals with respect to time (i and i− 1), where i is time with a measured value. A lag plot is useful for
examining the dependency of the error terms. Any non-random
pattern in a lag plot suggests that the variance is not random. No
particular pattern emerges from analysis of the lag plot, suggestingthat the errors are random for both the unimodal and bimodal
functions (Fig. 8).
The bimodal model performed well during both dry and wet periods. In September, which was characterized by sparse precipitation,
significant evapotranspiration, and daily irrigation, the model

VZJ | Advancing Critical Zone Sciencewas able to reproduce both the fast response of the green roof toprecipitation and the small response to drip irrigation. This also
indicates good accuracy in the estimation of daily fluctuations insoil water contents due to the combined effects of both evaporation
and root water uptake.
The overestimation of the peak flux values (Fig. 7) can be related
to the combined effects of uncertainties in measured precipitationand estimated unsaturated hydraulic conductivity. During large
precipitation events, a significant part of the domain near the geotextile (a seepage face) becomes saturated or nearly saturated, andunder such conditions the hydraulic conductivity plays a fundamental role in the infiltration process. Considering the uncertainty
and a possible bias introduced by the evaporation method forhydraulic conductivity values near saturation, it is reasonable toassume that the errors in the predicted peak fluxes are related tothis uncertainty. A more accurate description of the hydraulic conductivity close to saturation should help in improving the accuracy
p. 11 of 17

Fig. 8. A lag plot of the residuals between measured and simulated outflows. Results are for the HYDRUS model with the unimodal (left) and bimodal(right) functions of soil hydraulic properties.

of the model. Moreover, some uncertainty also arises from theevaluation of the soil depth in non-vegetated areas. Although anaverage depth of 4 cm was assumed, a slightly higher value may
also be realistic considering the structure of the green roof (Fig. 1
and 4), and this would slightly modify the hydraulic response of
the simulated green roof.
It can be concluded that HYDRUS, with both unimodal andbimodal functions of soil hydraulic properties, can accurately
describe the hydraulic behavior of the considered green roof. The
NSE indices are high, and residuals are randomly distributed forboth models, with a slightly better performance with the bimodal
functions. Because the bimodal model has proven to be moreaccurate in reproducing the real behavior of the green roof undersmall inflows, such as irrigation (Fig. 7), it was selected for further
analysis.

an increasing trend, caused by variations in actual evapotranspiration, as confirmed by the model. Figure 10 shows simulated actualroot water uptake and evaporation from vegetated areas.
The first part of September (Fig. 10) was characterized by relatively high evapotranspiration, which lowered water contents
in the soil substrate. As a consequence, only negligible outflow
was produced by irrigation. The evapotranspiration rate dropped
during the rain events between 9 and 11 Sept. 2015 due to the
combined effects of high air humidity and low solar radiation.

Hydrological Analysis of the Green Roof
Cumulative inflow and outflow fluxes of the green roof arereported in Fig. 9. The green roof, coupled with the reuse system,
was able to reduce the runoff volume by 25% during the considered
period. Considering that the volume of water stored in the green
roof substrate was only 16 mm (3% of the total inflow), evapotranspiration was the main factor in reducing the runoff volume.
The steep gradients in cumulative outflow (Fig. 9) indicate thatthe green roof had a fast response to precipitation, with a negligible delay in the hydrograph. This aspect is directly related to
the limited thickness of the green roof, which reduces the possible
delay effect. On the other hand, cumulative outflow appeared flat
when irrigation was applied. This behavior is particularly evident
in September, when only negligible outflow fluxes were observed.
At the beginning of October, cumulative outflow started to exhibit

VZJ | Advancing Critical Zone ScienceFig. 9. A comparison between cumulative inflow and outflow from thegreen roof, simulated by HYDRUS-3D.

p. 12 of 17

Fig. 10. Simulated actual root water uptake (top)and evaporation (bottom) from vegetated areasof the green roof.

At the end of September, due to the irrigation and lower evapotranspiration rates, the soil water contents were higher and actual
transpiration lower. This behavior is shown in detail in Fig. 11.
Between 25 and 30 Sept. 2015, the water content in vegetatedareas ranged between 0.45 and 0.50, which corresponded to apressure head range between −10 and −3 cm. In this pressure
head range, plant transpiration is limited because of anaerobic
conditions induced by high soil water contents (see the Feddes
parameters in Table 3). Under these conditions, significant outflow was measured after irrigation (Fig. 3). This indicates thatthe model reproduced the dynamics of evapotranspiration in thegreen roof with good accuracy and that the parameters chosen
for the water stress function are reasonable.
The comparison between pressure heads at the bottom of vegetatedand non-vegetated areas of the green roof helps in identifying thedifferent hydraulic responses of these two sections (Fig. 12). Before
irrigation started, the bottom pressure head in the non-vegetatedsection was higher (in absolute value) than in the vegetated section.
This indicates that the non-vegetated section dried out more andfaster than the vegetated section due to its very low thickness. As
irrigation started, the pressure head quickly increased in the nonvegetated section, reaching the seepage condition and producing
outflow. On the other hand, in the vegetated section, the pressure
head only approached the seepage condition without generating
outflow. This indicates that only the non-vegetated section of thegreen roof was responsible for outflow after irrigation and thatthe irrigation system could be optimized to avoid the formationof outflow after irrigation.

VZJ | Advancing Critical Zone Science

Hydrological Performance during
Precipitation Events
The hydrological response to single precipitation events is an
important characteristic in the analysis of green roofs and LIDs
in general. The measures such as peak flow reduction Pred (%), thepeak flow delay tdel (min), and the volume reduction Vred (%) provide information about the hydrological benefits of the LID system
to the urban drainage system. Four distinct rainfall events of a significant magnitude, one in September and three in October, wereidentified for further analysis of the hydrological performance ofthe green roof during single precipitation events. Figure 13 shows
precipitation and measured and modeled outflow fluxes for eachprecipitation event, as well as the total volume of the entire precipitation event Vprec (mm) and measures Pred and Vred.
It is evident from Fig. 13 that the model exhibited good accuracy in reproducing the green roof hydraulic response also forsingle precipitation events. Modeled Vred and Pred are in agreement with information reported above. The model tended tooverestimate the peak flow, especially for precipitation events
characterized by an early peak (7 and 10 Oct. 2015), while theaccuracy was very high for precipitation events with a delayed
peak (9 Sept. and 21 Oct. 2015). Modeled volume reductions
were sufficiently accurate. Only the analysis of the precipitationeven on 21 October exhibited a significant deviation between
the measured and modeled volume reductions.
The hydrological response of the green roof varied considerably for
different precipitation events. Delay time tdel was on the order of
p. 13 of 17

Fig. 11. Precipitation (dark area) and modeled (gray area) and measured (red line) outflow for four selected rainfall events in the analysis of the hydrological performance of the green roof during single precipitation events; Predmeas and Predmod are the measured and modeled peak flow reductions,
respectively, and Vredmeas and Vredmod are the measured and modeled volume reductions, respectively.

15 min in all events, indicating that the delay effect was limited dueto the limited thickness of the substrate. The largest Pred and Vred
were observed for the event of 7 October. During this event, theprecipitation peak occurred at the beginning of the event, which
was preceded by more than 2 d of a dry period since irrigation
was stopped on 5 October. Under these conditions, peak flow waspartially attenuated by the relatively low initial water content of
the soil substrate (Fig. 11). The effective volume reduction of the
precipitation event (Veff ) was 13 mm.
Although the precipitation event of 10 October was also characterized by an early peak, the measured peak reduction of 45% was
lower than that measured on 7 October. This can be attributed tothe higher water content at the beginning of the event. While on 7
October the initial water content was about 0.32, it was about 0.39
on 10 October. The higher water content on 10 October resulted

VZJ | Advancing Critical Zone Science

in a lower attenuation capacity of the substrate. For the events of
9 and 21 September, Pred was only 7%. For both events, peak flow
was preceded by low-intensity precipitation, which increased thesoil water content and, as a consequence, lowered the capacity of
the green roof to attenuate the peak of precipitation. The effective
volume reduction (Veff ) for the event of 21 October was 20.4 mm,
which was the highest among all of the evaluated precipitation
events. Before this event started, the water content of the substratereached the lowest value of 0.2, significantly increasing the soil
storage capacity.
It can be concluded that the green roof response to precipitation,
during the analyzed period, was influenced primarily by the antecedent substrate moisture and secondarily by the precipitation
pattern. Positively skewed and leptokurtic precipitation distributions tended to exhibit the largest peak flow reduction, while
p. 14 of 17

Fig. 12. Simulated water contents (q) in the vegetated section of the green roof at a depth of 4 cm.

negatively skewed and platykurtic precipitation distributions
were likely to have lower peak reductions. The volume retainedby the green roof was mainly determined by the initial substratemoisture and was independent of the precipitation pattern. It
must be emphasized that these conclusions are not general but are
restricted to the specific investigated green roof. The validated
model could be used in a sensitivity analysis, in which, together
with other factors (e.g., substrate moisture, slope, substrate depth,etc.), the effect of the precipitation pattern on the green roof
response can be investigated using synthetically designed storms(Carbone et al., 2015b).

66Conclusions

The goal of this study was to provide a comprehensive numerical analysis of the hydrological behavior of a green roof. The
widely used mechanistic model HYDRUS-3D was used to model
the hydrological responses to precipitation, evapotranspiration,and irrigation fluxes of green roof installed at the University ofCalabria. The green roof was characterized by a high degree ofcomplexity; it included different soil depths, vegetated and nonvegetated areas, and non-uniformity in the boundary conditions.
Moreover, considering the lack of studies with regard to the unsaturated hydraulic properties of soil substrates in green roofs, thesimplified evaporation method was used to determine both the soil
water retention curve and the unsaturated hydraulic conductivity function of the soil used in the green roof. The soil exhibiteda bimodal pore structure, characterized by a weakly pronounced
secondary pore system and a high hydraulic conductivity. The
bimodal functions proposed by Durner (1994) best represented
the hydraulic properties of the substrate. The unimodal VGM
relationships were also tested, but their description of unsaturated
hydraulic properties was less accurate, especially for the unsaturated hydraulic conductivity.
Soil hydraulic parameters obtained from the evaporation method
were then used in HYDRUS-3D to model the hydrologicalbehavior of the green roof during a 2-mo time period. Both theunimodal and bimodal functions of soil hydraulic properties

Fig. 13. Pressure heads at the bottom of the vegetated (black) and non-vegetated (red) sectionsof the green roof simulated by HYDRUS-3D.
The yellow rectangular area in the top figure
is expanded in the bottom figure. The dashedline represents a seepage condition.

VZJ | Advancing Critical Zone Sciencep. 15 of 17

were used in the analysis. The numerical simulation considered
precipitation, evaporation, root water uptake, and irrigation. The
Feddes model was used to represent the water stress response
function of plants installed on the green roof. Feddes’ parameters
were adjusted by taking into account the behavior of succulent
plants such as Carpobrotus edulis under different soil watercontent regimes. The NSE index was used to compare the simulated and measured outflows. The NSE indices were 0.74 and
0.8 for the model predictions with the unimodal and bimodal
functions, respectively, indicating a good accuracy. Both models
slightly overestimated some outflow fluxes. The randomness of
the residuals was confirmed by the absence of evident patterns in
the lag plot. The main difference between the two models was in
the description of the hydraulic behavior of the green roof after
irrigation. While the unimodal model failed to reproduce small
measured outflows, the bimodal model gave a more accurate
description. This further confirmed the accuracy of this model.
The validated model with bimodal functions was used to analyze
the hydrological performance of the green roof during the entire
simulation period and to investigate its hydrological response to
single precipitation events. About 25% of the total inflow volume
was returned to the atmosphere by evaporation and transpiration. The analysis of simulated pressure heads revealed that thenon-vegetated section of the green roof was responsible for thesmall outflow fluxes following irrigation. This suggests that theirrigation system could be further optimized. The analysis of thehydrological response of the green roof to single precipitation
events highlighted the importance of the initial soil moisturefor volume and peak flow reductions. Furthermore, the analysis
revealed that positively skewed and leptokurtic precipitation distributions are likely to exhibit the largest peak flow reductions.
The volume reduction is mainly influenced by the initial soil
moisture and is insensitive to the precipitation pattern.
The performance of the model can be improved by a better
description of the soil hydraulic properties, especially nearsaturation, considering that during intense precipitation events,
significant portions of the domain become saturated. Moreover,uncertainties in rainfall measurements and in geometric characteristics of the green roof can introduce a further bias tothe simulated results. A sensitivity analysis followed by anuncertainty analysis can help in identifying the most sensitive
parameters and address the source of uncertainty (e.g., Brunetti
et al., 2016). Nevertheless, it can be concluded that the use ofmeasured soil hydraulic properties with a mechanistic model canbe a valuable tool for the analysis of green roofs and other LIDs
and can boost the widespread adoption of such systems as a viable
alternative to traditional urban drainage systems.

Acknowledgments

The study was co-funded by the Italian National Operative Project (PON)–Research andCompetitiveness for the convergence regions 2007/2013-I Axis “Support to structuralchanges,” 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.”

VZJ | Advancing Critical Zone Science

References

Allen, R.G., L.S. Pereira, D. Raes, and M. Smith. 1998. Crop evapotranspiration. Irrig. Drain. Pap. 56. FAO, Rome.
Arya, L.M. 2002. Simultaneous determination of water transmission
and retention properties: Wind and hot-air methods. In: J.H. Daneand G.C. Topp, editors, Methods of soil analysis. Part 4. Physical methods. SSSA Book Ser. 5. SSSA, Madison, WI. p. 916–926.
doi:10.2136/sssabookser5.4.c38
Bengtsson, L., L. Grahn, and J. Olsson. 2004. Hydrological function of a thin
extensive green roof in southern Sweden. Nord. Hydrol. 36:259–268.
Blanusa, T., M.M. Vaz Monteiro, F. Fantozzi, E. Vysini, Y. Li, and R.W.F. Cameron. 2013. Alternatives to Sedum on green roofs: Can broad leaf perennial plants offer better “cooling service”? Build. Environ. 59:99–106.
doi:10.1016/j.buildenv.2012.08.011
Brattebo, B.O., and D.B. Booth. 2003. Long-term stormwater quantity andquality performance of permeable pavement systems. Water Res.
37:4369–4376. doi:10.1016/S0043-1354(03)00410-X
Brunetti, G., J. Šimůnek, and P. Piro. 2016. A comprehensive numericalanalysis of the hydraulic behavior of permeable pavement. J. Hydrol.
540:1146–1161. doi:10.1016/j.jhydrol.2016.07.030
Carbone, M., G. Brunetti, and P. Piro. 2014. Hydrological performanceof a permeable pavement in Mediterranean climate. In: 14th
SGEM GeoConference on Water Resources: Forest, Marine andOcean Ecosystems, Albena, Bulgaria. 17–26 June 2014. Book
3. Vol. 1. Int. Multidiscip. Sci. GeoConf., Sofia, Bulgaria. p. 381–
388. doi:10.5593/SGEM2014/B31/S12.050.
Carbone, M., G. Brunetti, and P. Piro. 2015a. Modelling the hydraulic behaviour of growing media with the explicit finite volume solution. Water 7:568–591. doi:10.3390/w7020568
Carbone, M., M. Turco, G. Brunetti, and P. Piro. 2015b. A cumulative rainfall function for subhourly design storm in Mediterranean urban areas.
Adv. Meteorol. 2015:1–10. doi:10.1155/2015/528564
Coffman, L.S. 2002. Low-impact development: An alternative stormwater
management technology. In: R.L. France, editor, Handbook of water
sensitive planning and design. Lewis Publ., Boca Raton, FL. p. 97–123.
Collins, K.A., W.F. Hunt, and J.M. Hathaway. 2008. Hydrologic comparison of four types of permeable pavement and standard asphalt in eastern North Carolina. J. Hydrol. Eng. 13:1146–1157.
doi:10.1061/(ASCE)1084-0699(2008)13:12(1146)
Dane, J.H., and J.W. Hopmans. 2002. Water retention and storage: Pressure plate extractor. In: J.H. Dane and G.C. Topp, editors, Methods ofsoil analysis. Part 4. Physical methods. SSSA Book Ser. 5. SSSA, Madison,
WI. p. 688–690. doi:10.2136/sssabookser5.4.c25
Davis, A.P. 2008. Field performance of bioretention: Hydrology impacts. J.
Hydrol. Eng. 13:90–95. doi:10.1061/(ASCE)1084-0699(2008)13:2(90)Duan, Q., S. Sorooshian, and V. Gupta. 1992. Effective and efficient global optimization for conceptual rainfall–runoff models. Water Resour.
Res. 28:1015–1031. doi:10.1029/91WR02985
Durhman, A.K., D.B. Rowe, and C.L. Rugh. 2006. Effect of watering regimen on chlorophyll fluorescence and growth of selected green roof
plant taxa. HortScience 41:1623–1628.
Durner, W. 1994. Hydraulic conductivity estimation for soils withheterogeneous pore structure. Water Resour. Res. 30:211–223.
doi:10.1029/93WR02676
Elliott, A., and S. Trowsdale. 2007. A review of models for low impact
urban stormwater drainage. Environ. Modell. Softw. 22:394–405.
doi:10.1016/j.envsoft.2005.12.005
Feddes, R.A., P.J. Kowalik, and H. Zaradny. 1978. Simulation of field water
use and crop yield. PUDOC, Wageningen, the Netherlands.
Getter, K.L., D.B. Rowe, and J.A. Andresen. 2007. Quantifying the effect ofslope on extensive green roof stormwater retention. Ecol. Eng. 31:225–
231. doi:10.1016/j.ecoleng.2007.06.004
Hanscom, Z., and I.P. Ting. 1978. Responses of succulents to plant water
stress. Plant Physiol. 61:327–330. doi:10.1104/pp.61.3.327
Hilten, R.N., T.M. Lawrence, and E.W. Tollner. 2008. Modeling stormwater runoff from green roofs with HYDRUS-1D. J. Hydrol. 358:288–293.
doi:10.1016/j.jhydrol.2008.06.010
Hu, S. 1987. Book review: Akaike information criterion statistics. Math.
Comput. Simul. 29:452. doi:10.1016/0378-4754(87)90094-2
Kasmin, H., V.R. Stovin, and E.A. Hathway. 2010. Towards a generic rainfall–runoff model for green roofs. Water Sci. Technol. 62:898–905.
doi:10.2166/wst.2010.352

p. 16 of 17

Klute, A., and C. Dirksen. 1986. Hydraulic conductivity and diffusivity. Laboratory methods. In: Methods of soil analysis. Part 1. Physical and mineralogical methods. 2nd ed. SSSA Book Ser. 5. SSSA and ASA, Madison,
WI. p. 687–734. doi:10.2136/sssabookser5.1.2ed.c28
Kundzewicz, Z., M. Radziejewski, and I. Pínskwar. 2006. Precipitation
extremes in the changing climate of Europe. Clim. Res. 31:51–58.
doi:10.3354/cr031051
Lazzarin, R.M., F. Castellotti, and F. Busato. 2005. Experimental measurements and numerical modelling of a green roof. Energy Build.
37:1260–1267. doi:10.1016/j.enbuild.2005.02.001
Li, Y., and R.W. Babcock. 2014. Green roof hydrologic performance and modeling: A review. Water Sci. Technol. 69:727–738.
doi:10.2166/wst.2013.770
Li, Y., and R.W. Babcock. 2015. Modeling hydrologic performance ofa green roof system with HYDRUS-2D. J. Environ. Eng. 141:04015036.
doi:10.1061/(ASCE)EE.1943-7870.0000976
Locatelli, L., O. Mark, P.S. Mikkelsen, K. Arnbjerg-Nielsen, M.B. Jensen,
and P.J. Binning. 2014. Modelling of green roof hydrological performance for urban drainage applications. J. Hydrol. 519:3237–3248.
doi:10.1016/j.jhydrol.2014.10.030

Schindler, U., 1980. Ein Schnellverfahren zur Messung der Wasserleitfähigkeit im teilgesättigten Boden an Stechzylinderproben. Arch. AckerPflanzenbau Bodenkd. 24:1–7.
Schindler, U., W. Durner, G. von Unold, L. Mueller, and R. Wieland. 2010a.
The evaporation method: Extending the measurement range of soil
hydraulic properties using the air-entry pressure of the ceramic cup. J.
Plant Nutr. Soil Sci. 173:563–572. doi:10.1002/jpln.200900201
Schindler, U., W. Durner, G. von Unold, and L. Müller. 2010b. Evaporation
method for measuring unsaturated hydraulic properties of soils: Extending the measurement range. Soil Sci. Soc. Am. J. 74:1071–1083.
doi:10.2136/sssaj2008.0358
Servat, E., and A. Dezetter. 1991. Selection of calibration objective functions in the context of rainfall–runoff modelling in a Sudanese savannah area. Hydrol. Sci. J. 36:307–330. doi:10.1080/02626669109492517
She, N., and J. Pang. 2010. Physically based green roof model. J. Hydrol.
Eng. 15:458–464. doi:10.1061/(ASCE)HE.1943-5584.0000138
Šimůnek, J., M.Th. van Genuchten, and M. Šejna. 2008. Development andapplications of the HYDRUS and STANMOD software packages andrelated codes. Vadose Zone J. 7:587–600. doi:10.2136/vzj2007.0077

Metselaar, K. 2012. Water retention and evapotranspiration of green
roofs and possible natural vegetation types. Resour. Conserv. Recycl.
64:49–55. doi:10.1016/j.resconrec.2011.12.009

Šimůnek, J., M.Th. van Genuchten, and O. Wendroth. 1998. Parameter estimation analysis of the evaporation method for determining soil hydraulic properties. Soil Sci. Soc. Am. J. 62:894­–905.
doi:10.2136/sssaj1998.03615995006200040007x

Min, S.-K., X. Zhang, F.W. Zwiers, and G.C. Hegerl. 2011. Human contribution to more-intense precipitation extremes. Nature 470:378–381.
doi:10.1038/nature09763

Starry, O., J.D. Lea-Cox, J. Kim, and M.W. van Iersel. 2014. Photosynthesis
and water use by two Sedum species in green roof substrate. Environ.
Exp. Bot. 107:105–112. doi:10.1016/j.envexpbot.2014.05.014

Moriasi, D.N., J.G. Arnold, M.W. Van Liew, R.L. Binger, R.D. Harmel, andT.L. Veith. 2007. Model evaluation guidelines for systematic quantification of accuracy in watershed simulations. Trans. ASABE 50:885–900.
doi:10.13031/2013.23153

Sutanto, S.J., J. Wenninger, A.M.J. Coenders-Gerrits, and S. Uhlenbrook.
2012. Partitioning of evaporation into transpiration, soil evaporation and interception: A comparison between isotope measurements and a HYDRUS-1D model. Hydrol. Earth Syst. Sci. 16:2605–2616.
doi:10.5194/hess-16-2605-2012

Nash, J.E., and J.V. Sutcliffe. 1970. River flow forecasting through conceptual models: I. A discussion of principles. J. Hydrol. 10:282–290.
doi:10.1016/0022-1694(70)90255-6
Pertassek, T., A. Peters, and W. Durner. 2015. HYPROP-FIT software user’s
manual, Version 3.0. UMS GmbH, Munich, Germany.
Peters, A., and W. Durner. 2006. Improved estimation of soil water retention characteristics from hydrostatic column experiments. Water Resour. Res. 42:W11401. doi:10.1029/2006WR004952

Taylor, S.A., and G.L. Ashcroft. 1972. Physical edaphology: The physics ofirrigated and nonirrigated soils. W.H. Freeman, San Francisco.
Torres, E.A., and A. Calera. 2010. Bare soil evaporation under high evaporation demand: A proposed modification to the FAO-56 model. Hydrol. Sci. J. 55:303–315. doi:10.1080/02626661003683249
UMS GmbH. 2015. Manual HYPROP, Version 2015-01. UMS GmbH, Munich, Germany.

Peters, A., and W. Durner. 2008. Simplified evaporation methodfor determining soil hydraulic properties. J. Hydrol. 356:147–162.
doi:10.1016/j.jhydrol.2008.04.016

van Dam, J.C., P. Groenendijk, R.F.A. Hendriks, and J.G. Kroes. 2008. Advances of modeling water flow in variably saturated soils with SWAP.
Vadose Zone J. 7:640–­653. doi:10.2136/vzj2007.0060

Ritchie, J.T. 1972. Model for predicting evaporation from a row
crop with incomplete cover. Water Resour. Res. 8:1204–1213.
doi:10.1029/WR008i005p01204

van Genuchten, M.Th. 1980. A closed-form equation for predicting thehydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–
898. doi:10.2136/sssaj1980.03615995004400050002x

Rosenberg, N.J., L.B. Blaine, and B.V. Shashi. 1983. Microclimate: The biological environment. John Wiley & Sons, New York.

Vesuviano, G., F. Sonnenwald, and V. Stovin. 2014. A two-stage storage
routing model for green roof runoff detention. Water Sci. Technol.
69:1191–1197. doi:10.2166/wst.2013.808

Sayed, O.H. 2001. Crassulacean acid metabolism 1975–2000, a check list.
Photosynthetica 39:339–352. doi:10.1023/A:1020292623960
Schaap, M.G., and F.J. Leij. 2000. Improved prediction of unsaturated hydraulic conductivity with the Mualem–van Genuchten model. Soil Sci.
Soc. Am. J. 64:843–851. doi:10.2136/sssaj2000.643843x
Schaap, M.G., F.J. Leij, and M.Th. van Genuchten. 2001. Rosetta:
A computer program for estimating soil hydraulic parameterswith hierarchical pedotransfer functions. J. Hydrol. 251:163–176.
doi:10.1016/S0022-1694(01)00466-8

VZJ | Advancing Critical Zone Science

Wendroth, O., W. Ehlers, H. Kage, J.W. Hopmans, J. Halbertsma, andJ.H.M. Wösten. 1993. Reevaluation of the evaporation method for determining hydraulic functions in unsaturated soils. Soil Sci. Soc. Am. J.
57:1436–1443. doi:10.2136/sssaj1993.03615995005700060007x
Wesseling, J., J. Elbers, P. Kabat, and B. Van den Broek. 1991. SWATRE: Instructions for input. Winand Staring Ctr., Wageningen, the Netherlands.
Zhang, S.H., and Y.P. Guo. 2013. Analytical probabilistic model for evaluating the hydrologic performance of green roofs. J. Hydrol. Eng.
18:19–28. doi:10.1061/(ASCE)HE.1943-5584.0000593

p. 17 of 17