The significance of recent and short pluviometric time series for the assessment of flood hazard in the context of climate change: examples from some sample basins of the Adriatic Central Italy

Numerical hydrological models are increasingly a fundamental tool for the analysis of floods in a river basin. If used for predictive purposes, the choice of the “design storm” to be applied, once set other variables (as basin geometry, land use, etc.), becomes fundamental. All the statistical methods currently adopted to calculate the design storm, suggest the use of long rainfall series (at least 40–50 years). On the other hand, the increasingly high frequency of intense events (rainfalls and floods) in the last twenty years, also as a result of the ongoing climate change, testify to the need for a critical analysis of the statistical significance of these methods. The present work, by applying the Gumbel distribution (Generalized Extreme Value Type-I distribution) on two rainfall series (1951–2018 and 1998–2018) coming from the same rain gauges and the “Chicago Method” for the calculation of the design storm, highlights how the choice of the series may influence the formation of flood events. More in particular, the comparison of different hydrological models, generated using HEC-HMS software on three sample basins of the Adriatic side of central Italy, shows that the use of shorter and recent rainfall series results in a generally higher runoff, mostly in case of events with a return time equal or higher than 100 years.


Introduction
Climate is undoubtedly the principal driving force in hydrologic systems. Nevertheless, the significant increase of floods observed in many regions of the world is often attributed to land use change, while the climate forcing under the effects of increased atmospheric greenhouse gases is usually underestimated.
At the end of 2018 the Intergovernmental Panel on Climate Change [1] published a "Special Report" on the impacts of 1.5 °C global warming above pre-industrial levels and related global greenhouse gas emission pathways, contained in the Decision of the 21st Conference of Parties of the United Nations Framework Convention on Climate Change to adopt the Paris Agreement. According to this Report, climate change is expected to increase extreme precipitations over many countries of Europe and, consequently, the risk of flooding on urban areas or rivers.
Even the Directive 2007/60/EC of the European Parliament and of the Council of 23 October 2007 (Directive "Floods"-FD), adopted with the aim of assess and manage flood risk, aims at the reduction of the harmful impacts on human health, environment, cultural heritage and economic activities. In particular, the article 4, paragraph 2 of the Directive states that "Based on available or readily derivable information, such as records and studies on long term developments, in particular impacts of climate change on the occurrence of floods, a preliminary flood risk assessment shall be undertaken to provide an assessment of potential risks".
The preparation of flood hazard and flood risk maps play and a key role in the implementation of the Directive itself and all the Member States, before the deadline of 22 December 2013, had to make maps at the most appropriate scale for selected areas, according to different climatic probability scenarios. Nevertheless, even though many States completed the procedure in time, several studies carried out in recent years [2,3] evidenced that the proposed scenarios, in the light of the ongoing climate change, are in some cases unrealistic, especially with regard to very extreme events (associated to long return time). The cause is certainly to be found, in addition to the consequences linked to the land use changes on the slopes and to the anthropization of alluvial plains and riverbeds, to the increase of intense rainfall events during the last twenty years, sometimes related (in densely populated areas) to urban growth-driven, microclimatic change (urban heat islands).
The hydrological models currently used for the definition of flood hazard scenarios, involve the use of different formulas and algorithms for the calculation of the "net rain" and for the definition of the "rainfall-runoff transformation process". Nevertheless, if the choice of the method or algorithm can produce different but substantially comparable results, the definition of the best "design storm", (i.e. the rainfall event that is expected to occur in the considered site for a given probability of occurrence and a given duration of the storm) obtained through statistical-probabilistic approaches, results critical.
The present study, by applying the Chicago-type hyetograph (considered representative of the design storm for the study area) to 18 rain gauges belonging to three sample basins of the Adriatic side of central Italy, highlights how the calculation of the hyetograph itself carried out by analyzing historical rainfalls series of different length, gives rise to significantly different results when used as input in a hydrological model for the calculation of the peak of flow along a river basin. Comparing the peak flow rates generated using a Chicago hyetograph derived from a "long" time series (about 70 years, from 1951 to 2018) and from a more recent part of it (about 20 years, 1998-2018) emerges as the latter produces in most cases higher flow rates. Although the result may appear less reliable from a statistical point of view, this aspect is instead important in a context of climate change, such as the current one, in which extreme events occur with ever more frequency.

Geological, geomorphological and climatic context
The study area (around 2600 km 2 ), which includes the Musone (651 km 2 ), Esino (1219 km 2 ) and Potenza (775 km 2 ) river basins, encompasses a wide territory of the Adriatic side of Central Italy, characterized by a typical mountain/high-hilly landscape in the central-western portion and few flat sectors located along the main alluvial plains and eastward, near the coast (Figure 1). -a coastal sector (C), between 0m and 537m a.s.l. characterized by a clayey-sandy-conglomeratic bedrock and mean annual precipitation between 600 mm and 850 mm; -an intermediate medium-low hilly sector (M), between 95m and 1475m a.s.l. with mainly arenaceous-conglomeratic and sandy-pelitic bedrock and mean annual precipitation between 850 mm and 1100 mm; -an inner high-hilly and mountain sector (H), between 271 m and 1692 m a.s.l with the presence of mainly calcareous lithotypes and mean annual rainfalls between 1100 mm, and 1700 mm. The trend of the average monthly temperatures and rainfalls, over the entire study area (calculated over the period 1950-2018) is shown in Figure 3.   Such climatic features characterize what is currently defined as "Adriatic-sublitorial" regime [4]. The number of rainy days in this sector ranges between 60 and 75 while the daily mean intensity between 10 and 12 mm; the absolute monthly maximum is usually in November (secondary maximum in spring). Summer is rather dry, particularly near the coast, where periods without rainfall exceed 40 days. The mean annual temperatures on the other hand vary between 12.5 °C and 15.5 °C; the annual temperature excursion is around 18 °C, while the daily one goes from 7 °C along the coast to 10 °C in the intramountain basins. The Aridity Index-AI values varies from 19.9 to 38, defining an area which presents an average, slight summer aridity in the littoral portion and the low hills. The Modified Fournier erosivity Index-MFI ranges from 8.2. to 9.6 indicating a low erosion potential.

Rainfall data analysis
Hydrological modelling, as known, requires the definition of design storms or precipitation hyetographs. Design storms are then used as inputs to hydrological models, while the resulting flows and flow rates in a river basin are calculated using rainfall-runoff and flow routing procedures.
The present study calculates the design storm following the "Chicago Method". This method establishes two different analytical equations for rainfall intensity over time, one valid before and another one valid after the peak rate, both derived from a Depth-Duration-Frequency (DDF) analytical expression, that preserves the same volumes of all rainfall intensities. Intensity-Duration-Frequency (DDF) curves, on the other hand, express in a synthetic way, for a given return time (T) (or its inverse, probability of exceedance) and a duration (t) of a rainfall event, the information on the maximum height of precipitation (h) and the maximum rainfall intensity (i). Generally, DDF curves can be described by the expression: in which a and n are parameters that have to be estimated through a probabilistic approach. The cumulative probability function F X (x) represents the probability of not exceeding the value of the rainfall height h by that random variable. In this study, the cumulative distribution function used was the classical Generalized Extreme Value (GEV) Type-I distribution (or Gumbel distribution) where X is a random variable, x is a possible value of X, ξ is the location parameter calculated using 0.5772 and α is the scale parameter calculated using /6, where is the mean and the variance of the data set [7,8]. Rainfall data used in this study come from 18 weather stations, managed by the Marche Region Multi-risk Functional Center and considered representatives of the area analyzed ( Figure 4). These weather station (whose characteristic parameters are given in Table 1) are distributed as follows: -10 in the coastal sector; -4 in the medium and low-hilly sector; -4 in the high-hilly and mountainous sector.  To calculate the design storms, the maximum rainfall values recorded each year for durations of 1, 3, 6, 12 and 24 hours were collected for each rain gauge and, where available, for the period 1950-2018.
To verify the existence of a positive or negative trend of rainfall data and their statistical significance, the non-parametric Mann-Kendall and Sen's slope estimator tests, very frequently used for the analysis of hydro-meteorological time series [9,10], were performed on both the "long" (1950-2018) and "short" (1998-2018) periods, to check also for possible changings.
The Mann-Kendall test S is generally calculated as: where is the number of data points, and are the data values in time series and ( > ), respectively and is the sign function as: The variance is calculated as: where n is the number of data, m is the number of tied groups (i.e. sets of sample data having the same value) and t i is the number of ties of extent i. In cases where the sample size n > 10, the standard normal test statistic Z is computed as follows: Positive values of indicate increasing trends while negative values show decreasing trends. Testing trends is usually done at a specific significance level. In this study, significance levels of 99% and 95% were used.
The Sen's slope estimator test (Q) describes the slope of trend in the sample of N pairs of data: where and are the data values at times j and ( > ), respectively. The Q sign evidences data trend reflection, while its value indicates the steepness of the trend. To determine whether the slope is statistically different than zero, one should calculate the confidence interval of Q at specific probability; same as for the Mann-Kendall test, significance levels of 99% and 95% have been chosen.
The results of the tests have shown in Table 2. From the analysis of results it is evident that only some weather stations register statistically significant trends; nevertheless, the latter always show a positive trend and the number of stations systematically increases (for all the durations analyzed) in the "short" period .
After this step, rainfall data were subsequently subjected to statistical analysis using the GEV Type I distribution to obtain the DDF curves. The spatialization of rainfall data over the whole study area functional to the use within the hydrologic model, has been finally implemented using the "Thiessen polygons" method [11]. Although other methods are also commonly used for the spatial distribution of rainfall data, the choice was here forced due to the software used for the hydrological modeling. As will be explained in more detail below, the association of specific design storm to each weather station within HEC-HMS, can only be achieved using the aforementioned tool; other methods only allow the insertion of a single maximum rainfall height, a procedure considered less suitable also for the successive calibration carried out with the observed river flow data.

Land use
Concerning the land use, vector data from the III rd level of the CORINE Land Cover 2012 inventory (CLC-2012) [12] have been collected (Figure 5a); the analysis of data, evidences that around 68% is characterized by agricultural areas, 27% by Forest and Seminatural areas and only 5% by Artificial surfaces. The characteristics of the soils, on the other hand, have been extrapolated generalizing and integrating data from the Soil Map of the Marche Region at 1: 250,000 scale. The different types of soil have been merged according to the World Reference Base for Soil Resources-WRB [13]. Four main groups of the WRB characterize the study area: Cambisols, Leptosols, Regosols, Calcisols.
Cambisols, very common in temperate and boreal regions, are brownish soils with weak horizon differentiation and characterized by the presence of a cambic horizon (Bw), below the organic-mineral one. These soils are predominant in the sedimentary basin placed within the calcareous ridges of the Apennine chain. Leptosols are very shallow soils with minimal development, formed typically on hard rock or highly calcareous materials; in the study area they are present along the Apennine ridges and are mainly developed over colluvial and slope quaternary deposits. Regosols, mostly developed along the alluvial plains, are usually not shallow enough to be Leptosols and consist of very weakly developed mineral soils in unconsolidated materials with only an ochric surface horizon. Calcisols, finally, are soils with a significant secondary accumulation of calcium carbonate and are common on the low-hilly areas of the region where a mainly pelitic and sandy-pelitic bedrock outcrops. All these groups are characterized by dominant calcareous horizons, in line with the geological features of the territory and are usually well-drained soils with fine to medium texture; given their high spatial variability, a further generalization was carried out in order to obtain more homogeneous areas, in agreement with the geomorphological characteristics of the territory (Figure 5b). According to the NRCS (Natural Resource Conservation Service) hydrologic method [14], for each type of soil have been assigned values as a percentage of the Hydrologic Soil Groups classified as follows on the base on their infiltration rate: -Group A: soils with low runoff potential and high infiltration rate even when thoroughly wetted. They consist mainly of deep, well drained sand or gravel and have a high rate of water transmission (greater than 7.6 mm/hr). -Group B: soils with moderate infiltration rate when thoroughly wetted and consist of moderately deep to deep, moderately well to well drained soils with moderately fine to moderately coarse textures. They have a moderate rate of water transmission (3.8-7.6 mm/hr). -Group C: soils with low infiltration rate when thoroughly wetted and consist chiefly of soils with a layer that inhibits downward movement of water and soils with moderately fine to fine texture. They have a low rate of water transmission (1.3-3.8 mm/hr). -Group D: soils with high runoff potential. They have very low infiltration rate when thoroughly wetted, they may contain a permanent water table and consist of clay soils with a high swelling potential, shallow soils over nearly impervious material and may contain a claypan or clay layer at or near the surface. These soils have a very low rate of water transmission (0-1.3 mm/hr).

Hydrological modeling
The hydrological model for the study area was realized through two distinct phases. In a former phase, data concerning the morphometry and the hydro-geomorphological characteristics of the river basins (land use, soil type etc) were preprocessed using the tools for ArcGIS "HEC-GeoHMS" version 10.2, developed by the U.S. Army Corps of Engineers [15]. This software allows to process spatial data to obtain dimensional, morphological and hydrological characteristics of the river basins. These data were subsequently used to perform hydrological simulations through HEC-HMS software v. 4.3, also developed by the U.S. Army Corps of Engineers [16].
The morphological characterization of the river basins through HEC-GeoHMS was carried out using a 1m high-resolution DTM (obtained processing LiDAR images provided by the Italian Ministry of Environment): this allows to disaggregate the three watersheds (Musone, Esino and Potenza) into a series of interconnected sub-basins (50 in total) and to divide the stream networks in reaches and junctions (Figure 6a). Physical characteristics such as the river length and slope, subbasin slope, subbasin centroid location and elevation, longest flow have been also extracted. The evaluation of the hydrological parameters (necessary for the subsequent calculation of infiltration and runoff coefficients in HEC-HMS) was then performed by applying the "Curve Number" method, developed by the Soil Conservation Service of the United States. The method provides the computation of the parameter Curve Number (CN), a dimensionless parameter that ranges between 0 and 100 and expresses the capacity of the soil to produce direct runoff: the higher the value the greater the runoff produced with the same rainfall. A grid file of the Curve Number was obtained by geoprocessing operation in ArcGIS, combining the previously created shapefiles of the land use and of the soil type ( Figure 6b). The variability in the CN depends on several factors as rainfall intensity and duration, soil moisture conditions, vegetation cover, temperature; these factors are collectively called Antecedent Runoff Conditions (ARCs). ARCs are divided into three classes: II for average conditions, I for dry conditions, and III for wet conditions [17]. Taking into account the purpose of the work, i.e. the runoff computation in critical conditions, the ARCs III were chosen and the initial grid file (corresponding to the ARCs II) was modified according to the following formula: The CN values thus obtained can be considered reliable, as they are in line with those obtained in previous studies through calibration procedures with real events (Materazzi, 2015).

Results and discussion
The rainfalls series collected for each of the 18 rain gauges, as described in chapter 2b, have been processed for the construction of the DDF curves for different event durations, return times and length of the time series (Figure 7a). More specifically, for each rain gauge, the analysis was performed on: -two historical series (1951-2018; 1998-2018) -five durations (1, 3, 6, 12, 24 hours) -two return times: 100 yrs (corresponding to an Annual Exceedance Probability-AEP of 1%) and 200 yrs (corresponding to an Annual Exceedance Probability-AEP of 0.5%). The results of this analysis are globally shown in Table 3. For each rain gauge (and consequently for each return time and for both series analyzed) the Chicago-type hyetographs was computed; these hyetographs ( Figure 7b) were later used as input within each of the 12 HEC-HMS hydrological models (3 river basins × 2 return times × 2 rainfalls series). The regionalization of the rainfall data, as mentioned in paragraph 2b, was carried out using the Thiessen polygon method. Once the design storm was developed, the basin model was defined in HEC-HMS. Some data relating to the sub-basins (including the NC), previously calculated in HEC-GeoHMS, were automatically imported into the model. The time of corrective action was instead calculated using the Kirpich Method, developed for small drainage basins in Tennessee and Pennsylvania, with basin areas from 1 to 112 acres 0.0078 . . (9) where Tc is the time of concentration in minutes, L is the maximum hydraulic flow length in feet, and H is the difference in elevation in feet between the outlet of the watershed and the hydraulically most remote point in the watershed. Basing on the values obtained, a duration of 6 hours for the design storm of each rain gauge, in line with the time of concentration of all the sub-basins, was chosen. "Net rain" and "rainfall-runoff transformation process" were finally evaluated by applying the Soil Conservation Service methods, as the SCS Curve Number Loss and the SCS Unit Hydrograph Transform respectively.
Concerning the river reaches the choice of the routing procedure fell on the SCS Lag method, which is best suited for short stream segments with a predicable travel time that doesn't vary with flow depth. Once the data has been entered, the simulation has been started for each of the models provided. The results are shown in Table 4.

Continued on next page
For each basin analyzed, the table shows the list of sub-basins, the flow rate calculated for different return times (100 yrs or 200 yrs) and for the two historical series analyzed (1951-2018 or 1998-2018). The rightmost column, with different colors, indicates the percentage variations (positive or negative) of the flow rate calculated from the most recent rainfall series compared to the longer one. In particular, red color evidences flow rate increase, green color flow rate decrease, yellow color positive or negative deviations not exceeding 10%; this latter value takes into account the uncertainty and the number of rain gauges (not high) used for the statistical analysis.
The analysis of results, clearly shows that most of the Potenza and Musone sub-basins, if modeled using the rainfall series 1998-2018, increases more than 30% the flow rate for both 100 yrs and 200 yrs return times. The Musone river, on the other hand, has a more uncertain behavior. On the other hand, observing the same results in Figure 8a,b (where the areas of the Thiessen polygons are also reported) it is evident that most of the "uncertain" percentage variations (within a range of ±10%) are only associated with 3-4 rain gauges, while those negative essentially to only one. These apparent rainfalls anomalies, once the quality of data coming from these sensors has been verified, will require future climatological insights.

Conclusions
The present work, comparing the results of hydrological models derived from statistical analyses of different historical rainfalls series showed that: -Mann-Kendall and Sen's Slope estimator tests evidenced a positive, although weak and not evenly distributed, positive trend of rainfall starting from 1950; this trend is slightly more evident for the period 1998-2018; -the DDF curves of 18 pluviometers uniformly distributed in the study area, calculated using a recent rainfalls series , show heights of precipitation significantly higher than those derived from a longer rainfall series (1951-2018) for the return times considered (100 yrs and 200 yrs); -Chicago-type design storms, calculated from the recent rainfalls series, generate higher flood rates when entered as input data in a hydrological model; -the hydrological model developed for the study area shows how the flow rate generated by the use of the 1998-2018 pluviometric series is in many cases more than 30% higher than the other; -values apparently in contrast with those cited are to be associated with the response of few, single sensors whose reliability have to be verified or adequately justified from a climatological point of view; -the results obtained by the present study, although preliminary and limited to a small sector of central Italy, show that the ongoing ascertained climate change is associated with a change in the rainfall regime, turning towards an increase in rainfall intensity and frequency of extreme events; -the use of recent but short rainfall series, although less significant from a statistical point of view, can therefore be more reliable if included in a numerical modeling for the definition of future risk scenarios.