- Email: [email protected]
Monte Carlo analysis of probability of inundation of Rome
Environmental Modelling & Software 22 (2007) 1409e1416 www.elsevier.com/locate/envsoft
Monte Carlo analysis of probability of inundation of Rome L. Natale a, F. Savi b,* b
a Department of Hydraulic and Environmental Engineering, University of Pavia, Italy Department of Hydraulics, Transportation and Highways, University of Rome ‘‘La Sapienza’’, Via Eudossiana 20, I-00184 Rome, Italy
Received 14 October 2005; received in revised form 28 November 2006; accepted 11 December 2006 Available online 26 March 2007
Abstract Rome’s monumental centre has often been inundated by Tiber River. In the last decades of the 19th century, river walls were erected to protect Rome from floods so that the last significant flood, which occurred in 1937, caused only marginal damages. Although the probability of inundation of the city seems to be now substantially reduced, the evaluation of the residual risk is still worthwhile. With this aim, rainfall, rainfallerunoff, river flood propagation and street flooding processes are simulated in detail to produce the inundation scenarios analysed by the Monte Carlo method. The study shows that severe floods, having a return period greater than 180 years, overtop both the left and right river banks and inundate the northern outskirts of Rome, while extreme events, with 1000 years return period, submerge large parts of the monumental centre of Rome. Ó 2007 Elsevier Ltd. All rights reserved. Keywords: Flooding of urban areas; Semi distributed rainfall-runoff model; Flood routing model; Regionalisation; Monte Carlo method
1. Introduction The monumental centre of Rome was frequently inundated by the Tiber River until river walls, called muraglioni, were erected between the end of the 19th and the beginning of the 20th century (Fig. 1). They protect Rome from floods comparable to the catastrophic 1870 event that submerged the centre of the city. Later on, in the mid of 20th century, several dams were built along the river that store a significant portion of the flood wave coming from the upper part of the watershed. Corbara reservoir, the most important of these new structures, provides 165 hm3 of active storage. In the last century only three e 1900, 1915, 1937 e significant flood occurred, producing marginal damages. During the 1937 flood, water flowed over both the left and right banks near ponte Milvio, the most ancient roman bridge in Rome. It is located upstream of Rome’s muraglioni which the 1937
* Corresponding author. Tel.: þ39 0644585491; fax: þ39 0644585065. E-mail addresses: [email protected] (L. Natale), [email protected] (F. Savi). 1364-8152/$ - see front matter Ó 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.envsoft.2006.12.004
flood did not overtop. Since a relatively small area was flooded, the common belief has taken hold that Rome is quite safe from future floods. But, in the last decades, the town has significantly expanded along the river course and new districts have developed in the poorly protected floodplain, upstream and downstream of the muraglioni. Calenda et al. (1997) estimated the inundation probability of the modern city using the standard procedure based on evaluation of the probability distribution function of the maximum annual discharge at the Ripetta reference cross section (Fig. 1). They computed water levels along the river by integrating the 1D steady-state flow equation with different discharges to evaluate the probability of overflow. But the 75-year long Ripetta record suffers some limitations with respect to its usefulness for flood-frequency analysis. During the recording period, only one inundation occurred that of 1937. Ripetta measurements are influenced by outflows that occur upstream, where no gauging stations are available. The sample is not statistically homogeneous due to flood routing in the Corbara reservoir that came into operation only in 1965.
1410
L. Natale, F. Savi / Environmental Modelling & Software 22 (2007) 1409e1416
Fig. 1. Map of monumental centre of Rome.
Since Ripetta observations are not truly representative of the flood conditions in Rome, we extended the analysis by generating synthetic samples of flood discharges using a suite of mathematical models and Monte Carlo procedures. 2. The mathematical models The following suite of mathematical models generated the set of synthetic scenarios of inundation of Rome (Natale and Savi, 2004). The KORNA stochastic rainfall model produces time series of daily/hourly rainfall distributed over the Tiber basin. The TEVERE BASIN conceptual hydrologic model transforms the synthetic rainfall into hourly runoff from 43 Tiber sub-basins. The TEVERE RIVER model simulates flood wave propagation in the lower reach of Tiber River from Corbara Lake to Tyrrhenian Sea. The URBE model simulates the propagation of the inundation in Rome.
2.1. Rainfall stochastic model Since the shape of the flood wave that arrives in Rome is due to the superimposition of flood waves coming from the upper basins and from the downstream tributaries, the model used in the Monte Carlo procedure should reproduce, at least approximately, the spatial distribution of the rainfall. For this reason, Kottegoda et al. (2003) developed the KORNA model that simulates statistically realistic daily rainfall runs through a two-station model at a key station and satellite stations. The method maintains the first three moments of the station rainfall time series and the spatial correlation structure of the rainfall series. To maintain periodicity and persistence of daily rainfall, a two-state first-order Markov model is adopted with a two-harmonic representation of the seasonal variation in the Markov parameters (Kottegoda et al., 2004). The model of daily rainfall was calibrated by using historical series recorded at 608 gauging stations. The time series varied from 2 to 74 years, with an average length of about 20 years.
L. Natale, F. Savi / Environmental Modelling & Software 22 (2007) 1409e1416
Daily rainfall was desegregated into hourly values through dimensionless accumulated hourly amounts generated by a beta distribution. It is postulated that the occurrence process of hourly rainfall has a geometric distribution, conditioned on the total daily rainfall (Kottegoda et al., 2003). The model used to desegregate daily into hourly rainfalls was calibrated on the basis of 38 storms recorded at six hourly-recording rainfall stations. With this model, 100 series of daily rainfalls, averaged over the nine rainfall regions of the Tiber watershed were generated (Fig. 2). Each synthetic series is 1000 years long. 2.2. Rainfallerunoff model This model belongs to the class of semi distributed rainfallrunoff models able to represent spatial variability of runoff generation and simulate both quick and slow flow components (Croke et al., 2006). The rainfallerunoff mathematical model is described in Fig. 3, where the meaning of the 13 parameters is explained. Rainfall intensity r is divided into surface runoff q and infiltration rate f. The infiltration rate is computed by
q
IUH sur
surface runoff
r rain
1411
q
f
+
+ f
SOIL infiltration rate
p
+
αp
αp
S
channelized flow
IUH sub subsurface runoff
(1-α)p
(1-α)p
deep percolation
r
= rainfall intensity
f
= infiltration rate
q
= surface runoff q= max ( 0,r-f )
S
= soil water content
p
= subsurface runoff
f
p
fi
pf
ff
pi Sif
Sff
S
Sip
Sfp
S
Fig. 3. Scheme of the rainfallerunoff model.
applying the mass balance equation to a control volume of soil of unit area. The water balance equation of the conceptual reservoir is: dS ¼f p dt
ð1Þ
where S is the water level in the reservoir and p is the subsurface flow. The surface runoff that results is represented as: q ¼ maxð0; r f Þ:
Fig. 2. Rainfall regions of Tiber watershed.
The fraction a p of the subsurface runoff p contributes to the surface runoff, while the remaining fraction (1 a) p is lost in deep percolation: the variables f and p depend on soil water content S. Separate convolution of the two runoff components q and p to surface and subsurface Gamma IUH gives the surface runoff Q. The model does not take snowfall and snowmelt into account. The Tiber watershed, shown in Fig. 4, is divided into 40 small sub-basins located downstream of Corbara dam and three larger sub-basins, upstream of the dam. Hourly and daily rainfall hyetographs of 12 storms and corresponding hourly stage hydrographs recorded at the closure section of 11 sub-basins were used to calibrate the model. Since the number of hourly recorded series was low, daily amounts of rainfall were desegregated at hourly scale with reference to the most similar recorded hourly hyetograph chosen by multi-criteria analysis (Natale and Savi, 2004). Average hourly rainfall on each sub-basin was obtained by Kriging spatial interpolation (Bastin et al., 1984). On the whole, 66 floods were used to calibrate the rainfalle runoff model on 11 sub-basins (Gallagher and Doherty, 2007). The values of the model parameters were calibrated (Conti et al., 2002) by means of the genetic algorithm proposed by Wang (1991, 1997). The Curve Number CNe for each gauged sub-basin, estimated from maps of soil type and soil use (SCS, 1985), was compared to the average of the observed values of CNo
L. Natale, F. Savi / Environmental Modelling & Software 22 (2007) 1409e1416
1412
where X is the generic parameter of subsurface flow model shown in Fig. 2 (a, fi, ff, Sif, Sff, pi, pf, Sip, Sfp). Coefficients bj and gj of regressions Eq. (2) were estimated using the data of gauged sub-basins. Since the value of single storm CN depends on antecedent soil moisture content and evapotranspiration (Young and Carleton, 2006), the ratio CN/CNm was related to the antecedent precipitation index API (Heggen, 2001) and to the evapotranspiration ratio ETp0 according to the linear relationship:
N
Corbara reservoir
CN=CNm ¼ a þ bAPI þ cETp0
ð3Þ
Coefficients a, b and c were calibrated using the data of gauged sub-basins. The mean seasonal variation of daily temperature, averaged over 30 years, is considered in Thornthwaite (1946) formula: ETp0 ¼ cT
a
ð4Þ
where T is the mean temperature averaged over the 15 days antecedent the storm, assumed uniform over the Tiber watershed. The coefficients a and c are: a ¼ 0:016 J; a
Ty an
eni
rrh
Rome sea
0 5 10
20
30
40
50
60
75
100Km
Fig. 4. Sub-basins of Tiber watershed considered in the rainfallerunoff model.
directly deduced from the recorded data; estimated values match the observed ones for the upper sub-basins (upstream and just downstream of Corbara), yielding the ratio rCN ¼ CNe/CNo ¼ 1.05 on average, and are slightly greater than observed ones for downstream sub-basins, with rCN ¼ 1.23. At the end of the calibration procedure, we estimated the values of 13 model parameters for the 11 gauged sub-basins and the mean values of the Curve Number CNm for both gauged and ungaged sub-basins. For gauged sub-basins CNm ¼ CNo, for ungaged sub-basins CNm ¼ CNe rCN. The four parameters of the Gamma IUH (both for surface and subsurface flow) were correlated to area and average elevation of the gauged sub-basins by linear regressions in order to estimate the parameters of the rainfallerunoff model of ungaged sub-basins. The remaining nine parameters, all of which refer to subsurface flow model were linearly correlated to the CN value of the single storm: Xj ¼ bj CN þ gj j ¼ 1; .; 9
ð2Þ
c ¼ 1:6ð10=JÞ P 1:514 where J ¼ 12 and Ti is the monthly mean temi¼1 ðTi =5Þ perature in the ith month. In applying Monte Carlo procedure, continuous synthetic rainfall runs are generated as described in Section 2.1. Since in the Monte Carlo procedure we consider only the floods generated from the main storms of each synthetic run, the values of the model parameters, referring to subsurface flow are calculated only for these storms. For each storm and for each subbasin: the values of API and ETp0 are evaluated to compute CN by Eq. (3); then the values of the nine parameters are evaluated by Eq. (2); the flood hydrograph is produced by the calibrated model. 2.3. Flood routing model The flood routing model in the looped network of eight reaches, schematised in Fig. 5, includes two channels crossing Tiber island (Isola Tiberina), crossed by four bridges and two sills, which strongly influence flow, and two channels, Fiumara Grande and Canale Fiumicino, through which the Tiber flows out to the Tyrrhenian sea. In each reach, 1D equations of free-surface boundary conditions and gradually varied flow were considered (Cunge et al., 1980): vY vQ B þ þq ¼0 vt vs
ð5aÞ
vQ v Q2 vY þ þ gA þ gASf ¼ 0 vt vs A vs
ð5bÞ
L. Natale, F. Savi / Environmental Modelling & Software 22 (2007) 1409e1416
1413
street/channel, flow is mostly one-dimensional and only in small areas, mainly at crossroads and squares, is it likely that two-dimensional effects are important. This suggests that a significant simplification may be introduced by modelling the complex hydraulic system as 1D open channel network, provided that some hydraulic conditions are specified at the junctions (Braschi et al., 1991). Further simplifications are introduced: the roads were modelled as channels connecting nodes located in crossroads and squares. The storage capacity in each node depends on the area of the open spaces, as streets, gardens, courtyard, etc., in the area of influence of the node. Therefore, this storage capacity is computed by multiplying the area of influence of each node by the density of the buildings. The contribution of local inertia in the flow equation of channels is disregarded and the continuity equation of channels is included in the continuity equations of the nodes. According to these assumptions, the mass balance equation for each of the M nodes of the network reads:
Fig. 5. Hydrographic network considered in the flood routing model.
where s is the spatial coordinate along the reach, Y is the water surface elevation, Q is the discharge, B is the surface width, A is the wetted cross section, g is the gravitational acceleration, Sf is the slope friction and q is the unit discharge flowing over the banks. For each node, where M reaches are converging/diverging, mass and total head balance were imposed: M X
Qk ¼ 0
N vh X ¼ Qj þ Qe with i ¼ 1; .; M vt j¼1
ð7Þ
where S is the influence area of each node that can be flooded, Qj are the discharges in each of the Ni channels connected to the node i and Qe is the discharge lost in basements, cellars and underground storages. The discharge in the channel is computed using the analytical Bresse’s solution of the 1D steady-state flow equation (Henderson, 1966). The system of M nonlinear algebraic equations is solved by iteration with the over-relaxation method.
ð6aÞ 3. Monte Carlo procedure and definition of the scenarios
k¼1
H1 þ DH1 ¼ H2 þ DH2 ¼ / ¼ HM þ DHM
Si
ð6bÞ
where Qk and Hk are discharges and total heads in the upstream/downstream sections of the reaches converging/diverging at the node, and DHk are minor losses at the node. Eqs. (5a) and (5b) were discretised according to Preissmann’s implicit finite difference scheme (Cunge et al., 1980). The resulting system of equations, including boundary conditions, was solved by using the condensation method proposed by Uan (1984). The model simulates flow through hydraulic singularities along the river as well as bridges (65), dams (4), sills and drop structures (12), by imposing internal boundary conditions, i.e. rating curves. In the same way, inflows from tributaries are simulated. The computational domain was discretised into more than 670 cross section. The model was calibrated by reproducing six recent floods and 11 historical floods. 2.4. Urban inundation model The street network can by conceptually assimilated to a system of channels supplied by the main river channel. In each
First of all, 100 continuous series of 1000 years of daily rainfalls were synthetically generated; daily rainfalls of major storms were desegregated into hourly rainfall and flood hydrographs from 43 Tiber tributaries were obtained, running the rainfallerunoff model. Combining and routing the flood waves from tributaries in the propagation model, the hydrograph of the maximum flood, in each year of each series, in any cross section was computed. Time series M ¼ 2, 3, 4, 6, 10, 20, 31, 50, 100, 200, 500, 1000 years long were randomly extracted from each 1000 years long series. In any cross section, a return period of M years has been assigned to the value of the highest peak discharge in M years: in this way, 100 samples of peak discharge were obtained for each M years return period. The average of the values in the sample having M year return period, is the estimate of the expected M years return period value (discharge, water level, etc.). To evaluate the influence of Corbara reservoir, Monte Carlo procedure was performed twice, both (a) considering and (b) ignoring routing in the reservoir. In case (a), it is assumed that the reservoir is half full at the beginning of the flood, as estimated by Calenda and Mancini
L. Natale, F. Savi / Environmental Modelling & Software 22 (2007) 1409e1416
1414
(1999). Constant outflow of 300 m3/s is maintained until water level in the reservoir reaches normal pool elevation: for higher water elevation, outflow equals inflow. In case (b) reservoir outflow always equals inflow.
10 F censored GEV
8
Monte Carlo
6
The sample of 75 annual maximum peak flows, observed at Ripetta gauging station from 1922 to 2002, is labelled OS in the following discussion. The sample of 62 naturalised annual maxima, computed by Calenda et al. (1997) by removing Corbara storage routing from measured hydrographs is labelled NS in the following discussion. For the purposes of our analysis 1870, 1878, 1900 and 1915 flood peak discharges were estimated from historical measures. Nonparametric Wilcoxon and ManneWhitney tests (A€ıvazian, 1970; McCuen, 2003) indicated that the homogeneity hypothesis should be rejected for OS sample and two sub-samples e OS1 e ante-Corbara (1922e1964) and e OS2 e post-Corbara (1965e2002) should be used instead of the complete sample. Consequently, probability distribution of present time maxima was estimated from OS2 sample e for this step EV1 and GEV (Hosking et al., 1985) functions were considered. Pearson and KolmogoroveSmirnov goodness of fit tests (Kottegoda and Rosso, 1997; McCuen, 2003) rejected the GEV probability law, but the EV1 was accepted. The theoretical PDF of peak discharges were compared to the corresponding Monte Carlo frequency distribution, shown in Fig. 6. The sample of naturalised peak discharge NS at Ripetta was extended to include historical records of 1870, 1878, 1900 and 1915 floods. Since low floods, occurred in very dry years, tend to bias the estimates drawn from extended NS sample,
Y Gumbel
4. Frequency analisys of maximum annual peak discharge 4
2
0
-2
-4 0
1000
3000
2000
4000
Q (m3/s) Fig. 7. Probability distribution and Monte Carlo frequency distribution of maximum annual naturalised flood.
especially if large recurrence periods are considered, we discarded values lower than the median of the systematic sample. Following the procedure suggested by Lian (1990) and Hirsch and Stedinger (1987), the plotting positions of NS integrated sample were estimated by separating systematic records from historical floods. The parameters of the GEV probability distribution, fitted to the extended and censored sample were estimated by Partial Probability Weighted Moments Method e PPWM e (Wang, 1990). In Fig. 7, the empirical sample of the maximum annual peak of naturalised floods at Ripetta is plotted on probabilistic paper and compared to censored GEV and to Monte Carlo frequency distribution.
8
3
F GEV EV1 Monte Carlo
2
6
T=10 years 1
Y EV1
Y Gumbel
4 0
T=200 years
2 -1
0
-2 Ferrovia Roma Nord bridge Ripetta
-3
-2 0
1000
2000
3000
4000
Q (m3/s) Fig. 6. Probability distribution and Monte Carlo frequency distribution of maximum annual regulated flood.
0
1000
2000
3000
4000
5000
6000
7000
Q (m3/s) Fig. 8. Frequency distribution of maximum peak discharges in Rome for T ¼ 20 and T ¼ 200 years (EV1 probabilistic paper).
L. Natale, F. Savi / Environmental Modelling & Software 22 (2007) 1409e1416
1415
For 100 to 1000 years return periods, the Monte Carlo values are about 15% lower than the GEV ones and about 5% lower than Peak Over Threshold estimates. Statistical analysis overestimates peak discharges for high values of return period because this method does not consider that overtopping of the banks a short distance upstream tends to reduce the highest peak discharges at the Ripetta river section. Fig. 8 shows the empirical distribution of two samples of peak discharges, calculated at the Ferrovia Roma Nord and Ripetta cross sections. For T ¼ 10 years, the frequency distributions of peak discharges are quite similar; for T ¼ 200 years, it is observed that peak discharges at Ripetta section are limited as water flows over upstream banks. The volumes of water flowing over right banks are significantly greater than the ones passing over left banks.
In the past, overflow from sewers directly connected to the river and the associated water-table rise produced only minor inundations in Rome, typically resulting in flooded cellars, basements and low lands. The floodwater flowed through the centre of the city with negligible velocities. But catastrophic floods, with many casualties, could occur if banks were to be overtopped upstream of ponte Milvio and waters were to run along via Flaminia and other historical roads with potentially dangerous velocities. Results from the Monte Carlo procedure indicate that these calamities could occur under present conditions if the flood discharge were to exceed 3200 m3/s at Ferrovia Roma Nord, located about 10 km upstream of ponte Milvio (Fig. 1).
Fig. 9. Flooded area for T ¼ 200 years.
Fig. 10. Flooded area for T ¼ 1000 years.
5. Present time inundation scenarios
1416
L. Natale, F. Savi / Environmental Modelling & Software 22 (2007) 1409e1416
For the time being, only inundations due to overtopping of levees and river walls are possible because Rome is now protected by muraglioni. Monte Carlo analysis assumes that inundation is to be expected for a given return period when overflow occurs in more than 50 over the 100 simulation cases considered in the study. Corbara reservoir significantly affects the shape of the flood hydrograph, so that the return period of inundation increases from 65 without Corbara to 180 years with Corbara. Since storage capacity of depressed riverside areas near ponte Milvio is about 1.0 hm3, on the left bank, and 0.7 hm3, on the right bank, two inundation scenarios can be identified. When inundation volume is lower than the storage capacity of depressed areas, then only the northern outskirts of the city are inundated. On other hand, for large outflow, inundation propagates downstream and the monumental centre of Rome is submerged. According to the Monte Carlo computations, the first case occurs for something less than a 200 years recurrence period, while the second case occurs, on average, every 1000 years. Applying the URBE model, both inundation scenarios were simulated: flooded areas are contoured in Figs. 9 and 10, respectively. The analysis of inundation scenarios emphasises the importance of several interrelated factors, natural and human, which affect the onset and the development of the flood: - the monumental centre of Rome is now protected by muraglioni and inundation must be considered as extremely rare; - the Corbara reservoir retards the arrival of the flood wave from the upper Tiber and avoids superposition of flood waves originating from the upper and lower part of the watershed so that the resulting flood wave is smoothed and peak discharges in Rome are substantially reduced; - water volumes flowing on the streets are likely to be small and only very rare floods can now reach the monumental centre of Rome. These conclusions are consistent with the present situation of the Tiber valley upstream of Rome. In the future, the hydraulic safety of Rome depends on the continued storage capacity of these areas. If the storage capacity in these areas were to decrease, the hydraulic risk for the city would increase significantly. References A€ıvazian, S., 1970. E´tude statistique des de´pendences. MIR, Moscow. Bastin, G., Lorent, B., Duque´, C., Gevers, M., 1984. Optimal estimation of the average areal rainfall and optimal selection of rain gauge locations. Water Resources Research 20 (4), 463e470.
Braschi, G., Gallati, M., Natale, L., 1991. Modelling floods in urban areas. In: ASFPM 15th Annual Conference, Denver. Calenda, G., Mancini, P.C., 1999. The role of the Corbara reservoir on the Tiber River in the flood protection of the town of Rome. In: Waterpower 99. ASCE, Las Vegas. Calenda, G., Di Malta, L., Mancini, C.P., Ubertini, L., 1997. Distribuzione di probabilita` dei colmi di piena del Tevere a Roma. L’Acqua 5, 13e22. Conti, M., Ferrante, M., Savi, F., 2002. A genetic algorithm to calibrate rainfallerunoff model of Tiber Catchment. In: IASTED International Conference on Applied Simulation and Modelling, Creta. Croke, B.F.W., Andrews, F., Jakerman, A.J., Cuddy, S.M., Luddy, A., 2006. IHACRES classic plus: a redesign of the IHACRES rainfall-runoff model. Environmental Modelling & Software 21 (3), 426e427. Cunge, J., Holly, F.M., Verwey, A., 1980. Practical Aspects of Computational River Hydraulics. Pittman Publ., Boston. Gallagher, M., Doherty, J., 2007. Parameter estimation and uncertainty analysis for a watershed model. Environmental Modelling & Software 22, 1000e1020. Heggen, R.J., 2001. Normalized antecedent precipitation index. Journal of Hydrologic Engineering, ASCE 6 (5), 377e381. Henderson, F.M., 1966. Open Channel Flow. MacMillan Pub., New York. Hirsch, R.M., Stedinger, J.M., 1987. Plotting positions for historical floods and their precision. Water Resources Research 23 (4), 715e727. Hosking, J.R.M., Wallis, J.R., Wood, E.F., 1985. Estimation of the generalised extreme value distribution by the method of probability weighted moments. Technometrics 27 (3), 251e261. Kottegoda, N.T., Natale, L., Raiteri, E., 2003. A parsimonious approach to stochastic modelling and disaggregation of daily rainfall. Journal of Hydrology 274, 47e61. Kottegoda, N.T., Natale, L., Raiteri, E., 2004. Some considerations of periodicity and persistence in daily rainfalls. Journal of Hydrology 296, 23e37. Kottegoda, N.T., Rosso, R., 1997. Probability, Statistics, and Reliability for Civil and Environmental Engineers. McGraw-Hill Co. Inc., New York. Lian, G.S., 1990. New plotting position rule for flood records considering historical data and palaeologic information. In: International Workshop on New Uncertainty Concepts in Hydrology and Water Resources, Warsaw. McCuen, R.H., 2003. Modeling Hydrologic Change: Statistical Methods. Lewis Publishers, CRC Press, New York. Natale, L., Savi, F., 2004. Analisi Monte Carlo degli scenari di inondazione di Roma. Atti dei Convegni Lincei, Roma. S.C.S., 1985. National Engineering Handbook, Section 4 e Hydrology. Water Resources Publication, Littleton. Thornthwaite, C.W., 1946. An approach toward a rational classification of climate. Geographical Review 38, 55e94. Uan, M., 1984, Principe de re´solution de mode´lisation filaires des e´coulements a surface libre: application a un ressue maille, EDF HE/43-80-40 Report, Grenoble. Wang, Q.J., 1990. Estimation of the GEV distribution from censored samples by method of partial probability weighted moments. Journal of Hydrology 120, 103e114. Wang, Q.J., 1991. The genetic algorithm and its application to calibrating conceptual rainfallerunoff models. Water Resources Research 27 (9), 2467e 2471. Wang, Q.J., 1997. Using genetic algorithms to optimise model parameters. Environmental Modelling & Software 12 (1), 27e34. Young, D.F., Carleton, J.N., 2006. Implementation of a probabilistic curve number method in the PRZM runoff model. Environmental Modelling & Software 21, 1172e1179.
Monte Carlo analysis of probability of inundation of Rome L. Natale a, F. Savi b,* b
a Department of Hydraulic and Environmental Engineering, University of Pavia, Italy Department of Hydraulics, Transportation and Highways, University of Rome ‘‘La Sapienza’’, Via Eudossiana 20, I-00184 Rome, Italy
Received 14 October 2005; received in revised form 28 November 2006; accepted 11 December 2006 Available online 26 March 2007
Abstract Rome’s monumental centre has often been inundated by Tiber River. In the last decades of the 19th century, river walls were erected to protect Rome from floods so that the last significant flood, which occurred in 1937, caused only marginal damages. Although the probability of inundation of the city seems to be now substantially reduced, the evaluation of the residual risk is still worthwhile. With this aim, rainfall, rainfallerunoff, river flood propagation and street flooding processes are simulated in detail to produce the inundation scenarios analysed by the Monte Carlo method. The study shows that severe floods, having a return period greater than 180 years, overtop both the left and right river banks and inundate the northern outskirts of Rome, while extreme events, with 1000 years return period, submerge large parts of the monumental centre of Rome. Ó 2007 Elsevier Ltd. All rights reserved. Keywords: Flooding of urban areas; Semi distributed rainfall-runoff model; Flood routing model; Regionalisation; Monte Carlo method
1. Introduction The monumental centre of Rome was frequently inundated by the Tiber River until river walls, called muraglioni, were erected between the end of the 19th and the beginning of the 20th century (Fig. 1). They protect Rome from floods comparable to the catastrophic 1870 event that submerged the centre of the city. Later on, in the mid of 20th century, several dams were built along the river that store a significant portion of the flood wave coming from the upper part of the watershed. Corbara reservoir, the most important of these new structures, provides 165 hm3 of active storage. In the last century only three e 1900, 1915, 1937 e significant flood occurred, producing marginal damages. During the 1937 flood, water flowed over both the left and right banks near ponte Milvio, the most ancient roman bridge in Rome. It is located upstream of Rome’s muraglioni which the 1937
* Corresponding author. Tel.: þ39 0644585491; fax: þ39 0644585065. E-mail addresses: [email protected] (L. Natale), [email protected] (F. Savi). 1364-8152/$ - see front matter Ó 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.envsoft.2006.12.004
flood did not overtop. Since a relatively small area was flooded, the common belief has taken hold that Rome is quite safe from future floods. But, in the last decades, the town has significantly expanded along the river course and new districts have developed in the poorly protected floodplain, upstream and downstream of the muraglioni. Calenda et al. (1997) estimated the inundation probability of the modern city using the standard procedure based on evaluation of the probability distribution function of the maximum annual discharge at the Ripetta reference cross section (Fig. 1). They computed water levels along the river by integrating the 1D steady-state flow equation with different discharges to evaluate the probability of overflow. But the 75-year long Ripetta record suffers some limitations with respect to its usefulness for flood-frequency analysis. During the recording period, only one inundation occurred that of 1937. Ripetta measurements are influenced by outflows that occur upstream, where no gauging stations are available. The sample is not statistically homogeneous due to flood routing in the Corbara reservoir that came into operation only in 1965.
1410
L. Natale, F. Savi / Environmental Modelling & Software 22 (2007) 1409e1416
Fig. 1. Map of monumental centre of Rome.
Since Ripetta observations are not truly representative of the flood conditions in Rome, we extended the analysis by generating synthetic samples of flood discharges using a suite of mathematical models and Monte Carlo procedures. 2. The mathematical models The following suite of mathematical models generated the set of synthetic scenarios of inundation of Rome (Natale and Savi, 2004). The KORNA stochastic rainfall model produces time series of daily/hourly rainfall distributed over the Tiber basin. The TEVERE BASIN conceptual hydrologic model transforms the synthetic rainfall into hourly runoff from 43 Tiber sub-basins. The TEVERE RIVER model simulates flood wave propagation in the lower reach of Tiber River from Corbara Lake to Tyrrhenian Sea. The URBE model simulates the propagation of the inundation in Rome.
2.1. Rainfall stochastic model Since the shape of the flood wave that arrives in Rome is due to the superimposition of flood waves coming from the upper basins and from the downstream tributaries, the model used in the Monte Carlo procedure should reproduce, at least approximately, the spatial distribution of the rainfall. For this reason, Kottegoda et al. (2003) developed the KORNA model that simulates statistically realistic daily rainfall runs through a two-station model at a key station and satellite stations. The method maintains the first three moments of the station rainfall time series and the spatial correlation structure of the rainfall series. To maintain periodicity and persistence of daily rainfall, a two-state first-order Markov model is adopted with a two-harmonic representation of the seasonal variation in the Markov parameters (Kottegoda et al., 2004). The model of daily rainfall was calibrated by using historical series recorded at 608 gauging stations. The time series varied from 2 to 74 years, with an average length of about 20 years.
L. Natale, F. Savi / Environmental Modelling & Software 22 (2007) 1409e1416
Daily rainfall was desegregated into hourly values through dimensionless accumulated hourly amounts generated by a beta distribution. It is postulated that the occurrence process of hourly rainfall has a geometric distribution, conditioned on the total daily rainfall (Kottegoda et al., 2003). The model used to desegregate daily into hourly rainfalls was calibrated on the basis of 38 storms recorded at six hourly-recording rainfall stations. With this model, 100 series of daily rainfalls, averaged over the nine rainfall regions of the Tiber watershed were generated (Fig. 2). Each synthetic series is 1000 years long. 2.2. Rainfallerunoff model This model belongs to the class of semi distributed rainfallrunoff models able to represent spatial variability of runoff generation and simulate both quick and slow flow components (Croke et al., 2006). The rainfallerunoff mathematical model is described in Fig. 3, where the meaning of the 13 parameters is explained. Rainfall intensity r is divided into surface runoff q and infiltration rate f. The infiltration rate is computed by
q
IUH sur
surface runoff
r rain
1411
q
f
+
+ f
SOIL infiltration rate
p
+
αp
αp
S
channelized flow
IUH sub subsurface runoff
(1-α)p
(1-α)p
deep percolation
r
= rainfall intensity
f
= infiltration rate
q
= surface runoff q= max ( 0,r-f )
S
= soil water content
p
= subsurface runoff
f
p
fi
pf
ff
pi Sif
Sff
S
Sip
Sfp
S
Fig. 3. Scheme of the rainfallerunoff model.
applying the mass balance equation to a control volume of soil of unit area. The water balance equation of the conceptual reservoir is: dS ¼f p dt
ð1Þ
where S is the water level in the reservoir and p is the subsurface flow. The surface runoff that results is represented as: q ¼ maxð0; r f Þ:
Fig. 2. Rainfall regions of Tiber watershed.
The fraction a p of the subsurface runoff p contributes to the surface runoff, while the remaining fraction (1 a) p is lost in deep percolation: the variables f and p depend on soil water content S. Separate convolution of the two runoff components q and p to surface and subsurface Gamma IUH gives the surface runoff Q. The model does not take snowfall and snowmelt into account. The Tiber watershed, shown in Fig. 4, is divided into 40 small sub-basins located downstream of Corbara dam and three larger sub-basins, upstream of the dam. Hourly and daily rainfall hyetographs of 12 storms and corresponding hourly stage hydrographs recorded at the closure section of 11 sub-basins were used to calibrate the model. Since the number of hourly recorded series was low, daily amounts of rainfall were desegregated at hourly scale with reference to the most similar recorded hourly hyetograph chosen by multi-criteria analysis (Natale and Savi, 2004). Average hourly rainfall on each sub-basin was obtained by Kriging spatial interpolation (Bastin et al., 1984). On the whole, 66 floods were used to calibrate the rainfalle runoff model on 11 sub-basins (Gallagher and Doherty, 2007). The values of the model parameters were calibrated (Conti et al., 2002) by means of the genetic algorithm proposed by Wang (1991, 1997). The Curve Number CNe for each gauged sub-basin, estimated from maps of soil type and soil use (SCS, 1985), was compared to the average of the observed values of CNo
L. Natale, F. Savi / Environmental Modelling & Software 22 (2007) 1409e1416
1412
where X is the generic parameter of subsurface flow model shown in Fig. 2 (a, fi, ff, Sif, Sff, pi, pf, Sip, Sfp). Coefficients bj and gj of regressions Eq. (2) were estimated using the data of gauged sub-basins. Since the value of single storm CN depends on antecedent soil moisture content and evapotranspiration (Young and Carleton, 2006), the ratio CN/CNm was related to the antecedent precipitation index API (Heggen, 2001) and to the evapotranspiration ratio ETp0 according to the linear relationship:
N
Corbara reservoir
CN=CNm ¼ a þ bAPI þ cETp0
ð3Þ
Coefficients a, b and c were calibrated using the data of gauged sub-basins. The mean seasonal variation of daily temperature, averaged over 30 years, is considered in Thornthwaite (1946) formula: ETp0 ¼ cT
a
ð4Þ
where T is the mean temperature averaged over the 15 days antecedent the storm, assumed uniform over the Tiber watershed. The coefficients a and c are: a ¼ 0:016 J; a
Ty an
eni
rrh
Rome sea
0 5 10
20
30
40
50
60
75
100Km
Fig. 4. Sub-basins of Tiber watershed considered in the rainfallerunoff model.
directly deduced from the recorded data; estimated values match the observed ones for the upper sub-basins (upstream and just downstream of Corbara), yielding the ratio rCN ¼ CNe/CNo ¼ 1.05 on average, and are slightly greater than observed ones for downstream sub-basins, with rCN ¼ 1.23. At the end of the calibration procedure, we estimated the values of 13 model parameters for the 11 gauged sub-basins and the mean values of the Curve Number CNm for both gauged and ungaged sub-basins. For gauged sub-basins CNm ¼ CNo, for ungaged sub-basins CNm ¼ CNe rCN. The four parameters of the Gamma IUH (both for surface and subsurface flow) were correlated to area and average elevation of the gauged sub-basins by linear regressions in order to estimate the parameters of the rainfallerunoff model of ungaged sub-basins. The remaining nine parameters, all of which refer to subsurface flow model were linearly correlated to the CN value of the single storm: Xj ¼ bj CN þ gj j ¼ 1; .; 9
ð2Þ
c ¼ 1:6ð10=JÞ P 1:514 where J ¼ 12 and Ti is the monthly mean temi¼1 ðTi =5Þ perature in the ith month. In applying Monte Carlo procedure, continuous synthetic rainfall runs are generated as described in Section 2.1. Since in the Monte Carlo procedure we consider only the floods generated from the main storms of each synthetic run, the values of the model parameters, referring to subsurface flow are calculated only for these storms. For each storm and for each subbasin: the values of API and ETp0 are evaluated to compute CN by Eq. (3); then the values of the nine parameters are evaluated by Eq. (2); the flood hydrograph is produced by the calibrated model. 2.3. Flood routing model The flood routing model in the looped network of eight reaches, schematised in Fig. 5, includes two channels crossing Tiber island (Isola Tiberina), crossed by four bridges and two sills, which strongly influence flow, and two channels, Fiumara Grande and Canale Fiumicino, through which the Tiber flows out to the Tyrrhenian sea. In each reach, 1D equations of free-surface boundary conditions and gradually varied flow were considered (Cunge et al., 1980): vY vQ B þ þq ¼0 vt vs
ð5aÞ
vQ v Q2 vY þ þ gA þ gASf ¼ 0 vt vs A vs
ð5bÞ
L. Natale, F. Savi / Environmental Modelling & Software 22 (2007) 1409e1416
1413
street/channel, flow is mostly one-dimensional and only in small areas, mainly at crossroads and squares, is it likely that two-dimensional effects are important. This suggests that a significant simplification may be introduced by modelling the complex hydraulic system as 1D open channel network, provided that some hydraulic conditions are specified at the junctions (Braschi et al., 1991). Further simplifications are introduced: the roads were modelled as channels connecting nodes located in crossroads and squares. The storage capacity in each node depends on the area of the open spaces, as streets, gardens, courtyard, etc., in the area of influence of the node. Therefore, this storage capacity is computed by multiplying the area of influence of each node by the density of the buildings. The contribution of local inertia in the flow equation of channels is disregarded and the continuity equation of channels is included in the continuity equations of the nodes. According to these assumptions, the mass balance equation for each of the M nodes of the network reads:
Fig. 5. Hydrographic network considered in the flood routing model.
where s is the spatial coordinate along the reach, Y is the water surface elevation, Q is the discharge, B is the surface width, A is the wetted cross section, g is the gravitational acceleration, Sf is the slope friction and q is the unit discharge flowing over the banks. For each node, where M reaches are converging/diverging, mass and total head balance were imposed: M X
Qk ¼ 0
N vh X ¼ Qj þ Qe with i ¼ 1; .; M vt j¼1
ð7Þ
where S is the influence area of each node that can be flooded, Qj are the discharges in each of the Ni channels connected to the node i and Qe is the discharge lost in basements, cellars and underground storages. The discharge in the channel is computed using the analytical Bresse’s solution of the 1D steady-state flow equation (Henderson, 1966). The system of M nonlinear algebraic equations is solved by iteration with the over-relaxation method.
ð6aÞ 3. Monte Carlo procedure and definition of the scenarios
k¼1
H1 þ DH1 ¼ H2 þ DH2 ¼ / ¼ HM þ DHM
Si
ð6bÞ
where Qk and Hk are discharges and total heads in the upstream/downstream sections of the reaches converging/diverging at the node, and DHk are minor losses at the node. Eqs. (5a) and (5b) were discretised according to Preissmann’s implicit finite difference scheme (Cunge et al., 1980). The resulting system of equations, including boundary conditions, was solved by using the condensation method proposed by Uan (1984). The model simulates flow through hydraulic singularities along the river as well as bridges (65), dams (4), sills and drop structures (12), by imposing internal boundary conditions, i.e. rating curves. In the same way, inflows from tributaries are simulated. The computational domain was discretised into more than 670 cross section. The model was calibrated by reproducing six recent floods and 11 historical floods. 2.4. Urban inundation model The street network can by conceptually assimilated to a system of channels supplied by the main river channel. In each
First of all, 100 continuous series of 1000 years of daily rainfalls were synthetically generated; daily rainfalls of major storms were desegregated into hourly rainfall and flood hydrographs from 43 Tiber tributaries were obtained, running the rainfallerunoff model. Combining and routing the flood waves from tributaries in the propagation model, the hydrograph of the maximum flood, in each year of each series, in any cross section was computed. Time series M ¼ 2, 3, 4, 6, 10, 20, 31, 50, 100, 200, 500, 1000 years long were randomly extracted from each 1000 years long series. In any cross section, a return period of M years has been assigned to the value of the highest peak discharge in M years: in this way, 100 samples of peak discharge were obtained for each M years return period. The average of the values in the sample having M year return period, is the estimate of the expected M years return period value (discharge, water level, etc.). To evaluate the influence of Corbara reservoir, Monte Carlo procedure was performed twice, both (a) considering and (b) ignoring routing in the reservoir. In case (a), it is assumed that the reservoir is half full at the beginning of the flood, as estimated by Calenda and Mancini
L. Natale, F. Savi / Environmental Modelling & Software 22 (2007) 1409e1416
1414
(1999). Constant outflow of 300 m3/s is maintained until water level in the reservoir reaches normal pool elevation: for higher water elevation, outflow equals inflow. In case (b) reservoir outflow always equals inflow.
10 F censored GEV
8
Monte Carlo
6
The sample of 75 annual maximum peak flows, observed at Ripetta gauging station from 1922 to 2002, is labelled OS in the following discussion. The sample of 62 naturalised annual maxima, computed by Calenda et al. (1997) by removing Corbara storage routing from measured hydrographs is labelled NS in the following discussion. For the purposes of our analysis 1870, 1878, 1900 and 1915 flood peak discharges were estimated from historical measures. Nonparametric Wilcoxon and ManneWhitney tests (A€ıvazian, 1970; McCuen, 2003) indicated that the homogeneity hypothesis should be rejected for OS sample and two sub-samples e OS1 e ante-Corbara (1922e1964) and e OS2 e post-Corbara (1965e2002) should be used instead of the complete sample. Consequently, probability distribution of present time maxima was estimated from OS2 sample e for this step EV1 and GEV (Hosking et al., 1985) functions were considered. Pearson and KolmogoroveSmirnov goodness of fit tests (Kottegoda and Rosso, 1997; McCuen, 2003) rejected the GEV probability law, but the EV1 was accepted. The theoretical PDF of peak discharges were compared to the corresponding Monte Carlo frequency distribution, shown in Fig. 6. The sample of naturalised peak discharge NS at Ripetta was extended to include historical records of 1870, 1878, 1900 and 1915 floods. Since low floods, occurred in very dry years, tend to bias the estimates drawn from extended NS sample,
Y Gumbel
4. Frequency analisys of maximum annual peak discharge 4
2
0
-2
-4 0
1000
3000
2000
4000
Q (m3/s) Fig. 7. Probability distribution and Monte Carlo frequency distribution of maximum annual naturalised flood.
especially if large recurrence periods are considered, we discarded values lower than the median of the systematic sample. Following the procedure suggested by Lian (1990) and Hirsch and Stedinger (1987), the plotting positions of NS integrated sample were estimated by separating systematic records from historical floods. The parameters of the GEV probability distribution, fitted to the extended and censored sample were estimated by Partial Probability Weighted Moments Method e PPWM e (Wang, 1990). In Fig. 7, the empirical sample of the maximum annual peak of naturalised floods at Ripetta is plotted on probabilistic paper and compared to censored GEV and to Monte Carlo frequency distribution.
8
3
F GEV EV1 Monte Carlo
2
6
T=10 years 1
Y EV1
Y Gumbel
4 0
T=200 years
2 -1
0
-2 Ferrovia Roma Nord bridge Ripetta
-3
-2 0
1000
2000
3000
4000
Q (m3/s) Fig. 6. Probability distribution and Monte Carlo frequency distribution of maximum annual regulated flood.
0
1000
2000
3000
4000
5000
6000
7000
Q (m3/s) Fig. 8. Frequency distribution of maximum peak discharges in Rome for T ¼ 20 and T ¼ 200 years (EV1 probabilistic paper).
L. Natale, F. Savi / Environmental Modelling & Software 22 (2007) 1409e1416
1415
For 100 to 1000 years return periods, the Monte Carlo values are about 15% lower than the GEV ones and about 5% lower than Peak Over Threshold estimates. Statistical analysis overestimates peak discharges for high values of return period because this method does not consider that overtopping of the banks a short distance upstream tends to reduce the highest peak discharges at the Ripetta river section. Fig. 8 shows the empirical distribution of two samples of peak discharges, calculated at the Ferrovia Roma Nord and Ripetta cross sections. For T ¼ 10 years, the frequency distributions of peak discharges are quite similar; for T ¼ 200 years, it is observed that peak discharges at Ripetta section are limited as water flows over upstream banks. The volumes of water flowing over right banks are significantly greater than the ones passing over left banks.
In the past, overflow from sewers directly connected to the river and the associated water-table rise produced only minor inundations in Rome, typically resulting in flooded cellars, basements and low lands. The floodwater flowed through the centre of the city with negligible velocities. But catastrophic floods, with many casualties, could occur if banks were to be overtopped upstream of ponte Milvio and waters were to run along via Flaminia and other historical roads with potentially dangerous velocities. Results from the Monte Carlo procedure indicate that these calamities could occur under present conditions if the flood discharge were to exceed 3200 m3/s at Ferrovia Roma Nord, located about 10 km upstream of ponte Milvio (Fig. 1).
Fig. 9. Flooded area for T ¼ 200 years.
Fig. 10. Flooded area for T ¼ 1000 years.
5. Present time inundation scenarios
1416
L. Natale, F. Savi / Environmental Modelling & Software 22 (2007) 1409e1416
For the time being, only inundations due to overtopping of levees and river walls are possible because Rome is now protected by muraglioni. Monte Carlo analysis assumes that inundation is to be expected for a given return period when overflow occurs in more than 50 over the 100 simulation cases considered in the study. Corbara reservoir significantly affects the shape of the flood hydrograph, so that the return period of inundation increases from 65 without Corbara to 180 years with Corbara. Since storage capacity of depressed riverside areas near ponte Milvio is about 1.0 hm3, on the left bank, and 0.7 hm3, on the right bank, two inundation scenarios can be identified. When inundation volume is lower than the storage capacity of depressed areas, then only the northern outskirts of the city are inundated. On other hand, for large outflow, inundation propagates downstream and the monumental centre of Rome is submerged. According to the Monte Carlo computations, the first case occurs for something less than a 200 years recurrence period, while the second case occurs, on average, every 1000 years. Applying the URBE model, both inundation scenarios were simulated: flooded areas are contoured in Figs. 9 and 10, respectively. The analysis of inundation scenarios emphasises the importance of several interrelated factors, natural and human, which affect the onset and the development of the flood: - the monumental centre of Rome is now protected by muraglioni and inundation must be considered as extremely rare; - the Corbara reservoir retards the arrival of the flood wave from the upper Tiber and avoids superposition of flood waves originating from the upper and lower part of the watershed so that the resulting flood wave is smoothed and peak discharges in Rome are substantially reduced; - water volumes flowing on the streets are likely to be small and only very rare floods can now reach the monumental centre of Rome. These conclusions are consistent with the present situation of the Tiber valley upstream of Rome. In the future, the hydraulic safety of Rome depends on the continued storage capacity of these areas. If the storage capacity in these areas were to decrease, the hydraulic risk for the city would increase significantly. References A€ıvazian, S., 1970. E´tude statistique des de´pendences. MIR, Moscow. Bastin, G., Lorent, B., Duque´, C., Gevers, M., 1984. Optimal estimation of the average areal rainfall and optimal selection of rain gauge locations. Water Resources Research 20 (4), 463e470.
Braschi, G., Gallati, M., Natale, L., 1991. Modelling floods in urban areas. In: ASFPM 15th Annual Conference, Denver. Calenda, G., Mancini, P.C., 1999. The role of the Corbara reservoir on the Tiber River in the flood protection of the town of Rome. In: Waterpower 99. ASCE, Las Vegas. Calenda, G., Di Malta, L., Mancini, C.P., Ubertini, L., 1997. Distribuzione di probabilita` dei colmi di piena del Tevere a Roma. L’Acqua 5, 13e22. Conti, M., Ferrante, M., Savi, F., 2002. A genetic algorithm to calibrate rainfallerunoff model of Tiber Catchment. In: IASTED International Conference on Applied Simulation and Modelling, Creta. Croke, B.F.W., Andrews, F., Jakerman, A.J., Cuddy, S.M., Luddy, A., 2006. IHACRES classic plus: a redesign of the IHACRES rainfall-runoff model. Environmental Modelling & Software 21 (3), 426e427. Cunge, J., Holly, F.M., Verwey, A., 1980. Practical Aspects of Computational River Hydraulics. Pittman Publ., Boston. Gallagher, M., Doherty, J., 2007. Parameter estimation and uncertainty analysis for a watershed model. Environmental Modelling & Software 22, 1000e1020. Heggen, R.J., 2001. Normalized antecedent precipitation index. Journal of Hydrologic Engineering, ASCE 6 (5), 377e381. Henderson, F.M., 1966. Open Channel Flow. MacMillan Pub., New York. Hirsch, R.M., Stedinger, J.M., 1987. Plotting positions for historical floods and their precision. Water Resources Research 23 (4), 715e727. Hosking, J.R.M., Wallis, J.R., Wood, E.F., 1985. Estimation of the generalised extreme value distribution by the method of probability weighted moments. Technometrics 27 (3), 251e261. Kottegoda, N.T., Natale, L., Raiteri, E., 2003. A parsimonious approach to stochastic modelling and disaggregation of daily rainfall. Journal of Hydrology 274, 47e61. Kottegoda, N.T., Natale, L., Raiteri, E., 2004. Some considerations of periodicity and persistence in daily rainfalls. Journal of Hydrology 296, 23e37. Kottegoda, N.T., Rosso, R., 1997. Probability, Statistics, and Reliability for Civil and Environmental Engineers. McGraw-Hill Co. Inc., New York. Lian, G.S., 1990. New plotting position rule for flood records considering historical data and palaeologic information. In: International Workshop on New Uncertainty Concepts in Hydrology and Water Resources, Warsaw. McCuen, R.H., 2003. Modeling Hydrologic Change: Statistical Methods. Lewis Publishers, CRC Press, New York. Natale, L., Savi, F., 2004. Analisi Monte Carlo degli scenari di inondazione di Roma. Atti dei Convegni Lincei, Roma. S.C.S., 1985. National Engineering Handbook, Section 4 e Hydrology. Water Resources Publication, Littleton. Thornthwaite, C.W., 1946. An approach toward a rational classification of climate. Geographical Review 38, 55e94. Uan, M., 1984, Principe de re´solution de mode´lisation filaires des e´coulements a surface libre: application a un ressue maille, EDF HE/43-80-40 Report, Grenoble. Wang, Q.J., 1990. Estimation of the GEV distribution from censored samples by method of partial probability weighted moments. Journal of Hydrology 120, 103e114. Wang, Q.J., 1991. The genetic algorithm and its application to calibrating conceptual rainfallerunoff models. Water Resources Research 27 (9), 2467e 2471. Wang, Q.J., 1997. Using genetic algorithms to optimise model parameters. Environmental Modelling & Software 12 (1), 27e34. Young, D.F., Carleton, J.N., 2006. Implementation of a probabilistic curve number method in the PRZM runoff model. Environmental Modelling & Software 21, 1172e1179.