Generally, roofs are the best candidates for rainwater harvesting. In this context, the correct evaluation of the quantity and quality of runoff from roofs is essential to effectively design rainwater harvesting systems. This study aims to evaluate the performance of a kinematic wave based numerical model in simulating runoff on sloping roofs, by comparing the numerical results with the ones obtained from laboratory rainfall simulations on a real-scale Lusa ceramic tile roof. For all studied slopes, simulated discharge hydrographs had a good adjust to observed ones. Coefficient of determination and Nash–Sutcliffe efficiency values were close to 1.0. Particularly, peak discharges, times to peak, peak durations and runoff volumes were very well simulated.

INTRODUCTION

Rainwater harvesting is a water saving alternative strategy that presents many advantages and can provide solutions to address major water resources problems, such as fresh water scarcity, urban stream degradation and flooding (e.g. Leusbrock et al. 2015). In recent years, these problems have become global challenges, due to climatic change, population growth and increasing urbanisation.

Generally, roofs are the first to come into contact with rainwater; thus, they are the best candidates for rainwater harvesting. In this context, the correct evaluation of roof runoff quantity and quality is essential to effectively design rainwater harvesting systems (e.g. Farreny et al. 2011). Despite this, many studies usually focus on the qualitative aspects (e.g. Lye 2009; Egodawatta et al. 2012) in detriment of the quantitative aspects. Quantifying roof runoff depends not only on the climatic conditions (e.g. rainfall) but also on the roof type. In the Mediterranean region, ceramic tile roofs are a well-known trademark (Farreny et al. 2011). In Portugal, one of the most used is the Lusa ceramic tile (Sousa et al. 1998).

A reliable rainfall–runoff numerical model can provide an efficient and economic tool to predict beforehand possible impacts of climatic changes and land use changes (e.g. increasing urbanisation), allowing the development of alternative water management policies and water saving techniques (e.g. Borris et al. 2013; Jung et al. 2015). Despite the large amount of studies, rainfall–runoff modelling still involves several difficulties, given the high spatial and temporal rainfall variability (e.g. Notaro et al. 2013).

For decades, kinematic wave modelling has been used as a reliable tool in water resources (see the reviews in Singh 1996, 2001). Despite its wide acceptance, several authors continue investigating kinematic wave accuracy and applicability in a wide range of runoff studies (e.g. de Lima & Singh 2002; Deng et al. 2005; Isidoro & de Lima 2013; Abrantes et al. 2015; Mallari et al. 2015).

Laboratory and field studies using rainfall simulators have been widely used to investigate rainfall–runoff processes (e.g. Egodawatta & Goonetilleke 2008; Arguelles et al. 2013; de Lima et al. 2013; Montenegro et al. 2013; Ries et al. 2014). These studies enabled a detailed exploration and systematic replication of a large range of hydrologic conditions, such as rainfall spatial and temporal characteristics, providing for a fast way to obtain precise and consistent data that can be used to calibrate and validate numerical models (e.g. Deng et al. 2005; Abrantes et al. 2015; Mallari et al. 2015).

This study aims to evaluate the performance of a kinematic wave based numerical model in simulating runoff on sloping roofs, by comparing the numerical results with the ones obtained from laboratory rainfall simulations on a real-scale Lusa ceramic tile roof.

MATERIALS AND METHODS

Laboratory installation

Experiments were conducted in a real scale laboratory setup, schematised in Figure 1, comprising a rainfall simulator (used in, for example, de Lima et al. (2013), Montenegro et al. (2013), Carvalho et al. (2014) and Isidoro & de Lima (2015)), a Lusa ceramic tile roof and a continuous discharge measuring system.
Figure 1

Sketch of the experimental setup.

Figure 1

Sketch of the experimental setup.

The rainfall simulator (Figure 1) has a pressurised hydraulic system comprised of: (i) a steady downward-oriented full-cone nozzle (1/4-HH-14W FullJet from Spraying Systems Co., USA), with 3.58 mm orifice diameter, positioned 2.2 m above the geometric centre of the roof surface; (ii) a hydraulic system attached just in front of the nozzle to eliminate pressure fluctuations (more details in Isidoro & de Lima 2015); and (iii) a submerged pump (76.2 mm SQ from Grundfos Holding A/S, Denmark), installed in a constant head reservoir supplied with tap water. This system allowed a steady operating pressure of 1.4 bar at the nozzle, producing a total discharge of 7.5 L min−1, with a spray angle of 120 °.

The ceramic tile roof (Figure 1) consisted on an assembly of seven rows by seven lines of UM Lusa ceramic tiles (from Umbelino Monteiro SA, Portugal), covering a total area of 4.2 m2 (1.5 × 2.8 m2), supported by a metal structure with adjustable longitudinal slope. According to the manufacturer, UM Lusa ceramic tiles' dimensions are 0.455 × 0.254 m2, weight 3.15 kg, have low permeability, high resistance to frost and are produced according to Portuguese technical standards. Manning's roughness coefficient for ceramic tiles was considered 58.8 m1/3 s−1 (Chow 1959).

Experiments were conducted only considering the five most central tile rows of the roof, positioned at the slopes of 3.5, 23.0 and 40.5%. The two most peripheral rows worked as buffer zones in order to compensate for water ejected outside the other rows due to splash.

Rainfall characteristics

Rainfall intensity spatial distribution shown in Figure 2 was obtained capturing precipitated water, during 1 min, in plastic recipients placed over the roof in a uniformly arranged grid. Mean rainfall intensity was obtained by collecting runoff, during 1 min, from the roof covered with pre-wetted impermeable plastic sheet, to avoid tile absorption. Uniformity coefficients of rainfall intensity spatial distribution, calculated according Christiansen (1942), were 85.3%, 81.9% and 78.4%, respectively, for the slopes of 3.5%, 23.0% and 40.5%. Rainfall measurements were repeated three times.
Figure 2

Rainfall intensity spatial distribution at the tile roof level and mean rainfall intensity of the five most central rows, for the three defined slopes (S). Isohyets (obtained using Kriging) and mean rainfall intensity are in mm h−1. The arrow (→) represents the flow direction and the symbol ◙ represents the position of the nozzle. Values are averages from three repetitions.

Figure 2

Rainfall intensity spatial distribution at the tile roof level and mean rainfall intensity of the five most central rows, for the three defined slopes (S). Isohyets (obtained using Kriging) and mean rainfall intensity are in mm h−1. The arrow (→) represents the flow direction and the symbol ◙ represents the position of the nozzle. Values are averages from three repetitions.

Raindrop mean diameter and mean velocity, calculated from measurements with a laser precipitation monitor (Thies Laser Distrometer from Adolf Thies GmbH & Co., Germany), were, respectively, 0.52 mm and 1.41 m s−1 (more details in Carvalho et al. (2014)).

Tiles absorption intensity

Water absorption tests were conducted to characterise UM Lusa ceramic tiles absorption intensity, for the three defined slopes. The tiles, pre-dried in an oven during 24 h at 105 °C, were placed over high precision weighting devices and were subjected to simulated rainfall during 30 min. Monitoring the combined weight of the tiles and absorbed water, it was possible to characterise the tiles absorption intensity temporal evolution, as shown in Figure 3. These tests were repeated three times for each defined slope.
Figure 3

UM Lusa ceramic tiles water absorption intensity, for the three defined slopes (S). Values, averages from three repetitions, were adjusted to a power function where the absorption intensity (a) is in mm h−1 and the time (t) is in min. Coefficient of determination (r2) is also shown.

Figure 3

UM Lusa ceramic tiles water absorption intensity, for the three defined slopes (S). Values, averages from three repetitions, were adjusted to a power function where the absorption intensity (a) is in mm h−1 and the time (t) is in min. Coefficient of determination (r2) is also shown.

Experimental procedure

Rainfall simulations were conducted to produce 15 discharge hydrographs for the five most central rows of the roof positioned at the defined slopes of 3.5, 32.0 and 40.5%. Each simulation comprised one constant rainfall event, with the duration of 1 min, approximately equal to the time of concentration (time to peak) of the roof positioned at the lowest slope of 3.5% (the highest time of concentration of the three slopes). All simulations were repeated three times to assure their reliability.

At the beginning of the experiments, the roof was wetted for 5 min. After that, rainfall events were simulated consecutively with a 5–10 min interval between them. Thus, the tiles' initial moisture and absorption condition were approximately the same for all simulations.

Discharge hydrographs at the roof downstream end were obtained using the continuous discharge measuring system (Figure 1). This system consisted of plastic recipients placed over weighting devices with a precision of 0.01 g (one for each row), to monitor runoff during rainfall events and runoff recessions, in a total of 3 min, with a temporal resolution of 1 s. Runoff was conveyed to the plastic recipients using funnels, minimising runoff losses and the potential turbulence caused by falling water, thus obtaining smoother hydrographs.

KINEMATIC WAVE MODELLING OF ROOF RUNOFF

Experimental results were compared to the ones simulated with a one-dimensional numerical model based on the kinematic wave approximation of the Saint-Venant shallow water equations. These are approximations of the shallow water mass and momentum conservation laws (Singh 1996) and can be expressed as: 
formula
1
where h(x,t) is the runoff depth [L], Q(x,t) is the runoff discharge [L3 T−1], q(x,t) is the rainfall excess intensity [L T−1], W is the runoff width [L], t is the simulation time [T] and x is distance along the flow direction [L].
The kinematic wave model assumes that the hydrostatic pressure distribution is valid across the runoff depth, the surface tension forces are negligible and the momentum coefficient variation along the flow direction is negligible. By also assuming that the bed slope (S) equals the friction slope (Sf) and using existing open channel flow friction equations, we can express the runoff discharge at any point and time as function of the runoff depth, as follows: 
formula
2
where α is the kinematic wave resistance parameter [L1/3 T−1] and β is the Bakmeteff dimensionless exponent [-]. Since runoff was assumed as turbulent, Manning's formula can be used to estimate α and β, as and , where kM is the Manning's roughness coefficient [L1/3 T−1] and S is the bed slope [L L−1].
Substituting Equation (2) in Equation (1), the kinematic wave equation can be written as: 
formula
3
In this formulation, rainfall excess intensity (q) was estimated as the rainfall intensity (p) subtracted by the tiles absorption intensity (a), as follows: 
formula
4
 
formula
5
where p is the rainfall intensity [L T−1] and a is the tiles absorption intensity [L T−1].
Rainfall was inserted in the model as a spatial varied pattern of discrete rainfall intensity blocks kept constant during the rainfall simulation time. Tiles absorption was considered constant in space and was inserted in the model as a power curve in function of the rainfall simulation time (t) adjusted by an initial absorption time (ta) that represents the tiles absorption intensity at the beginning of the simulation, as follows: 
formula
6
where ta is the initial absorption time [T], A is a constant of the power curve [L T−2] and B is an dimensionless exponent of the power curve [-]. The power function was adjusted from the absorption tests (Figure 3). Since the tiles' initial moisture and absorption condition were approximately the same for all simulations, simulated hydrographs were obtained considering only three initial absorption times (one for each slope), representing the tiles' absorption intensity at the beginning of the simulation. Therefore, according to the rainfall characteristics (e.g. spatial pattern) and the absorption characteristics of the tiles (e.g. initial moisture condition), rainfall excess intensity varied in space and time.
Equation (3) can be solved using the second-order single-step Lax–Wendroff numerical scheme, which can be expressed in finite difference form as (Singh 1996): 
formula
7
where i identifies the time (t), j identifies the space (x), Δt is the time discretisation [T] and Δx is the space discretisation [L].
For the downstream boundary , the following first-order scheme is employed (Singh 1996): 
formula
8
The solution of the kinematic wave equations (Equations (7) and (8)) requires the initial condition for and the upstream boundary condition for , where L is the total length of the surface along the flow direction [L].
To guarantee the numerical stability of the kinematic wave equations, Courant–Friedrichs–Lewy condition for linear numerical stability was satisfied (Singh 1996): 
formula
9
The performance of the numerical model was evaluated using the coefficient of determination (r2), the Nash–Sutcliffe efficiency (NSE) and the Nash–Sutcliffe efficiency with logarithmic values (lnNSE), comparing observed to simulated discharge at the downstream boundary, as follows (e.g. Krause et al. 2005): 
formula
10
 
formula
11
 
formula
12
where r2 is the Determination coefficient [-], NSE is the Nash-Sutcliffe efficiency [-], lnNSE is the Nash-Sutcliffe efficiency with logarithmic values [-], is the observed discharge at the downstream boundary and instant i [L3 T−1], is the simulated discharge at the downstream boundary and instant i [L3 T−1], is the average observed discharge at the downstream boundary [L3 T−1], is the average simulated discharge at the downstream boundary [L3 T−1] and n is the total runoff discharge data. The range of r2 lies between 0 and 1. A value of 0 means no correlation at all whereas a value of 1 means that the dispersion of the simulation is equal to that of the observation. The range of NSE and lnNSE lies between −∞ and 1. An efficiency of 1 means a perfect fit and an efficiency lower than 0 indicates that the mean value of the observed time series would have been a better predictor than the model.

RESULTS AND DISCUSSION

Experimental results

Discharge hydrographs observed from the experiments, for the three slopes, are shown in Figure 4, where it is possible to observe the effect of slope on runoff.
Figure 4

Discharge hydrographs observed from the experiments, for the roof positioned at the three defined slopes (S). Average values and standard deviation bars from the five most central rows of the tile roof and three repetitions.

Figure 4

Discharge hydrographs observed from the experiments, for the roof positioned at the three defined slopes (S). Average values and standard deviation bars from the five most central rows of the tile roof and three repetitions.

Considering the standard deviation, both runoff initiations, peak discharges and runoff volumes were very similar for all roof slopes, being, respectively, 13.3 ± 0.8 s, 7.2 ± 0.5 mL s−1 and 424.8 ± 20.2 mL. For all roof slopes, runoff coefficient of the roof was approximately 0.87. Furthermore, small differences on these parameters can be attributed to differences of rainfall intensity on each slope. Nevertheless, the gentle slope of 3.5% clearly produced a delay on the rising and recession limbs, when compared to the other slopes. Time to peak was 57.4 s, 23.0 s and 20.2 s, respectively, for the slopes of 3.5%, 33.0% and 40.5%.

Numerical results

Discharge hydrographs simulated with the kinematic wave numerical model and observed from the experiments are shown in Figure 5, using the power function for the tile absorption intensity adjusted from the absorption tests. For all studied slopes, simulated discharge hydrographs had a good adjust to observed ones. Particularly, peak discharges, times to peak, peak durations and runoff volumes were very well simulated, as shown in Figure 6. However, time to runoff always started earlier in numerical simulations.
Figure 5

Comparison between observed and simulated (kinematic wave model) discharge hydrographs for the five most central rows of the tile roof positioned at the three defined slopes (S). Coefficient of determination (r2), NSE and lnNSE are also shown. Observed values are averages from three repetitions.

Figure 5

Comparison between observed and simulated (kinematic wave model) discharge hydrographs for the five most central rows of the tile roof positioned at the three defined slopes (S). Coefficient of determination (r2), NSE and lnNSE are also shown. Observed values are averages from three repetitions.

Figure 6

Comparison between observed and simulated time to peak, peak duration, peak discharge and runoff volume for the roof positioned at the three defined slopes (S). For each slope, five values are presented corresponding to the five most central rows of the tile roof. Values are the averages from three repetitions.

Figure 6

Comparison between observed and simulated time to peak, peak duration, peak discharge and runoff volume for the roof positioned at the three defined slopes (S). For each slope, five values are presented corresponding to the five most central rows of the tile roof. Values are the averages from three repetitions.

Coefficient of determination (r2), NSE and lnNSE, comparing the observed to simulated discharge (Equations (10)–(12)) are also shown in Figure 5. The numerical model had a good performance on simulating runoff for all three slopes, since in r2 and NSE values were close to 1.0. However, lnNSE values were very low. While NSE strongly overestimates peak discharge values whereas lower values are neglected, with lnNSE the influence of the low discharge values is increased in comparison to the peak ones. Therefore, it is clear that the numerical model better estimated peak discharge values than the lower one. This happened due to the already noticed differences in time to runoff.

On average, the best model performance was obtained for the middle slope of 23.0%, with an r2 and an NSE of 0.98 and an lnNSE of 0.46. The worst adjustment was obtained for the highest slope of 40.5%, with an r2 and an NSE of 0.96 and an lnNSE of 0.42. However, differences in r2, NSE and lnNSE values for all three slopes were negligible.

CONCLUSIONS

A kinematic wave based numerical model was used to runoff on a Lusa ceramic tile roof, with good adjustment for all studied slopes. Coefficient of determination and NSE values were close to 1.0. Simulated peak discharges, times to peak and peak durations were close to the observations from the rainfall experiments. However, lower discharge values at the beginning of runoff were not simulated with a good performance.

Future work should include other types of ceramic tiles (e.g. Marseille, Canudo), other types of roofs, initial moisture conditions and absorption conditions and rainfall intensities (e.g. low rainfall intensity below absorption intensity).

ACKNOWLEDGEMENTS

The first author acknowledges CNPq, Brazil for the financial support through the Post Doctoral Grant 206872/2014-3. The second and third authors acknowledge FCT, Portugal for the financial support through the Doctoral Grant SFRH/BD/103300/2014 and the Project PTDC/ECM/HID/4259/2014. The authors also acknowledge the contribution of Antonio Mousinho Fernandes and Gabriel Silveira for their help during the experimental work.

The laboratory experiments described in this study were conducted at the Laboratory of Hydraulics, Water Resources and Environment of the Department of Civil Engineering of the Faculty of Science and Technology of the University of Coimbra, Portugal.

REFERENCES

REFERENCES
Arguelles
A. C. C.
Jung
M.
Pak
G.
Aksoy
H.
Kavvas
M. L.
Yoon
J.
2013
Evaluation of overland flow model for a hillslope using laboratory flume data
.
Water Science & Technology
68
(
5
),
1188
1194
.
Borris
M.
Viklander
M.
Gustafsson
A.-M.
Marsalek
J.
2013
Simulating future trends in urban stormwater quality for changing climate, urban land use and environmental controls
.
Water Science & Technology
68
(
9
),
2082
2089
.
Carvalho
S. C. P.
de Lima
J. L. M. P.
de Lima
M. I. P.
2014
Using meshes to change the characteristics of simulated rainfall produced by spray nozzles
.
International Soil and Water Conservation Research
2
(
2
),
67
78
.
Chow
V. T.
1959
Open-Channel Hydraulics
.
McGraw-Hill Book Inc.
,
New York, NY
,
USA
.
Christiansen
J. E.
1942
Irrigation by Sprinkling, California Agricultural Experiment Station Bulletin 670
.
University of California
,
Berkeley, CA
,
USA
.
de Lima
J. L. M. P.
Singh
V. P.
2002
The influence of the pattern of moving rainstorms on overland flow
.
Advances in Water Resources
25
(
7
),
817
828
.
de Lima
J. L. M. P.
Carvalho
S. C. P.
de Lima
M. I. P.
2013
Rainfall simulator experiments on the importance of when rainfall burst occurs during storm events on runoff and soil loss
.
Zeitschrift für Geomorphologie
57
(
1
),
91
109
.
Deng
Z.-Q.
de Lima
J. L. M. P.
Singh
V. P.
2005
Transport rate-based model for overland flow and solute transport: parameter estimation and process simulation
.
Journal of Hydrology
315
(
1–4
),
220
235
.
Egodawatta
P.
Miguntanna
N. S.
Goonetilleke
A.
2012
Impact of roof surface runoff on urban water quality
.
Water Science & Technology
66
(
7
),
1527
1533
.
Farreny
R.
Morales-Pinzón
T.
Guisasola
A.
Tayà
C.
Rieradevall
J.
Gabarrell
X.
2011
Roof selection for rainwater harvesting: quantity and quality assessments in Spain
.
Water Research
45
(
10
),
3245
3254
.
Isidoro
J. M. G. P.
de Lima
J. L. M. P.
2013
An analytical closed form solution for 1D kinematic overland flow under moving rainstorms
.
Journal of Hydrological Engineering
18
(
9
),
1148
1156
.
Jung
M.
Kim
H.
Mallari
K. J. B.
Pak
G.
Yoon
J.
2015
Analysis of effects of climate change on runoff in an urban drainage system: a case study from Seoul, Korea
.
Water Science & Technology
71
(
5
),
653
660
.
Krause
P.
Boyle
D. P.
Bäse
F.
2005
Comparison of different efficiency criteria for hydrological model assessment
.
Advances in Geosciences
5
,
89
97
.
Leusbrock
I.
Nanninga
T. A.
Lieberg
K.
Agudelo-Vera
C. M.
Keesman
K. J.
Zeeman
G.
Rijnaarts
H. H. M.
2015
The urban harvest approach as framework and planning tool for improved water and resource cycles
.
Water Science & Technology
72
(
6
),
998
1006
.
Lye
D. J.
2009
Rooftop runoff as a source of contamination: a review
.
Science of The Total Environment
407
(
21
),
5429
5434
.
Mallari
K. J. B.
Kim
H.
Pak
G.
Aksoy
H.
Yoon
J.
2015
A comparison of two infiltration models applied to simulation of overland flow over a two-dimensional flume
.
Water Science & Technology
71
(
9
),
1325
1332
.
Montenegro
A. A. A.
Abrantes
J. R. C. B.
de Lima
J. L. M. P.
Singh
V. P.
Santos
T. E. M.
2013
Impact of mulching on soil and water dynamics under intermittent simulated rainfall
.
Catena
109
,
139
149
.
Notaro
V.
Fontanazza
C. M.
Freni
G.
Puleo
V.
2013
Impact of rainfall data resolution in time and space on the urban flooding evaluation
.
Water Science & Technology
68
(
9
),
1984
1993
.
Singh
V. P.
1996
Kinematic Wave Modelling in Water Resources: Surface-Water Hydrology
.
John Wiley and Sons Ltd
,
Chichester
,
UK
.
Sousa
A. V. S.
Silva
M.
de Moura
G. I.
Abrantes
V.
da Silva
J. R. M.
Freitas
V.
de Sousa
H.
1998
Manual de aplicação de telhas cerâmicas (Ceramic Tiles Application Manual)
.
APICER
,
Coimbra
,
Portugal
.