The role of morphology in the spatial distribution of short-duration rainfall extremes in Italy

The dependence of rainfall on elevation has frequently been documented in the scientific literature and may be relevant in Italy, due to the high degree of geographical and morphological heterogeneity of the country. However, a detailed analysis of the spatial variability of short-duration annual maximum rainfall depths and their connection to the landforms does not exist. Using a new, comprehensive and position-corrected rainfall extreme dataset (I²-RED, the Improved Italian-Rainfall Extreme Dataset), we present a systematic study of the relationship between geomorphological forms and the average annual maxima (index rainfall) across the whole of Italy. We first investigated the dependence of sub-daily rainfall depths on elevation and other landscape indices through univariate and multivariate linear regressions. The results of the national-scale regression analysis did not confirm the assumption of elevation being the sole driver of the variability of the index rainfall. The inclusion of longitude, latitude, distance from the coastline, morphological obstructions and mean annual rainfall contributes to the explanation of a larger percentage of the variance, even though this was in different ways for different durations (1 to 24 h). After analyzing the spatial variability of the regression residuals, we repeated the analysis on geomorphological subdivisions of Italy. Comparing the results of the best multivariate regression models with univariate regressions applied to small areas, deriving from morphological subdivisions, we found that “local” rainfall–topography relationships outperformed the country-wide multiple regressions, offered a uniform error spatial distribution and allowed the effect of morphology on rainfall extremes to be better reproduced.

1 Introduction and background

The spatial patterns of rainfall depth statistics are known to be affected by the geomorphological setting (Smith, 1979; Basist et al., 1994; Prudhomme and Reed, 1998, 1999; Faulkner and Prudhomme, 1998). The impact of orography on daily, multi-daily and annual precipitation events can generally be attributed to the so-called “orographic enhancement of precipitation”, i.e., an increase in rainfall depth along the windward slope of a relief and a decrease on the lee side, due to the lifting and the consequent drying of the air mass (Smith, 1979; Daly et al., 1994; Frei and Schär, 1998; Napoli et al., 2019). In a complex landscape, this effect can also entail significant precipitation values on the lee side, due to landforms that cause a delay in the hydrometeorological formation of precipitation and falling raindrops (Smith, 1979).

The impact of the orography on extreme rainfall depths and the complicated atmosphere–orography interactions for large areas are still not sufficiently understood for sub-daily rainfall events. In a country like Italy, characterized by a high degree of morphological heterogeneity (Fig. 1), these relations assume an evident importance, considering the significant exposure to Mediterranean storms (Claps and Siccardi, 2000). The focus of this study is the entire Italian territory (≈ 300 000 km²), considered as a representative case, both in terms of variety of landforms and in terms of variability of the rainfall extremes, as will be seen in the following.

Figure 1 Elevation data with the boundaries of the 20 Italian administrative regions. Source: Shuttle Radar Topography Mission (Farr et al., 2007).

Most of the existing studies in Italy have focused on limited areas (Allamano et al., 2009; Caracciolo et al., 2012; Pelosi and Furcolo, 2015; Furcolo et al., 2016; Furcolo and Pelosi, 2018; Libertino et al., 2018; Formetta et al., 2022), and the only attempt to deal with sub-daily data covering the entire nation was made by Avanzi et al. (2015).

These studies suffered from a lack of a comprehensive and quality-assessed national database for sub-daily extremes. Several of them analyzed the Italian Alpine area. For instance, Frei and Schär (1998) focused on the entire European Alps region and showed that foothills enhance monthly and seasonal precipitation, while inner valleys produce an orographic shielding effect on rainfall. Nevertheless, they did not find a unique precipitation depth–elevation relationship that could be considered valid for the entire Alps and attributed the observed variability to the effects of slope and shielding. Allamano et al. (2009) investigated the dependence of sub-daily annual maximum rainfall depths on elevation over the Italian Alpine region. They found a significant decreasing trend for increasing elevations and a nonuniform slope coefficient over the longitude range. The slope of the rainfall depth–elevation regression was shown to decrease for event durations that increased from 1 to 24 h. Libertino et al. (2018) showed, in the western sector of the Italian Alps (Fig. 1), that shorter durations (1–3 h) are characterized by a negative slope coefficient as a function of elevation (statistically significant at a 5 % level), while longer durations (12–24 h) show a positive slope and also a significant correlation, and the trend of the extremes over 6 h loses significance with elevation. Over Trento province, Formetta et al. (2022) identified a reverse orographic effect for hourly and sub-hourly durations and an orographic enhancement for a duration of about 8 h (or longer).

Other regional works that attempted to identify orographic effects in the Mediterranean part of Italy are available for Campania and Sicily. Pelosi and Furcolo (2015) and Furcolo et al. (2016) analyzed the daily annual maximum rainfall depths over Campania (see Fig. 1 for the geographical location) and attempted to explain systematic variations as being the result of the presence of orographic barriers, identified through the application of an automatic geomorphological procedure (Cuomo et al., 2011). Their results showed a link between orographic elements and a local increase in rainfall depths and allowed orographic elements that produced enhanced variability of extreme rainfall to be identified. The same group later worked on sub-daily annual maximum rainfall depths (Furcolo and Pelosi, 2018) and proposed a power-law amplification factor of rainfall over three mountainous systems.

Caracciolo et al. (2012) found, in Sicily, that the longitude, latitude, distance from the sea and a concavity index are the variables that govern the spatial variability of rainfall depths. However, these authors found that no linear relationship between sub-daily annual maximum rainfall depths and elevation was significant at a 5 % level over the entire island of Sicily.

All of the previously mentioned analyses refer to an analytic relationship that connects annual maximum rainfall depths of various durations, i.e., the average depth–duration (ADD) curve of the simple-scaling approach, which is usually represented by a power law:

$$\begin{array}{}\text{(1)}& {\stackrel{\mathrm{‾}}{h}}_d=a\cdot d^n,\end{array}$$

where ${\stackrel{\mathrm{‾}}{h}}_d$ is the average of the annual maximum rainfall depths of duration d, a is a scale factor and n is a scaling exponent. Coefficient a represents the best unbiased linear estimation of the 1 h average rainfall depths, considering that ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{1}}\cong a$.

Avanzi et al. (2015) analyzed the spatial variability of the ADD curve parameters, a and n, at a national scale, as obtained from measurements of 1494 stations distributed throughout Italy. They referred to the so-called “reverse orographic effect”, i.e., the relationship found between parameter a and elevation, which shows a decreasing trend. On the other hand, the scaling exponent n appears to increase nonlinearly with elevation. More details are provided in the following section.

On the basis of the described background and the significant improvements offered by a new, up-to-date rainfall dataset, i.e., the Improved Italian-Rainfall Dataset, I²-RED (Mazzoglio et al., 2020a), the present study has considered more than 3700 stations with at least 10 years of data to relate the average rainfall depths in all the durations (index rainfall) to several morphological variables and investigate their dependency on elevation and on other geomorphological and climatological parameters.

Searching for models that allow the index rainfall to be estimated for various durations in any location in Italy is the first, important, necessary step toward addressing the building of depth–duration–frequency curves over the entire country. For this purpose, simple (Sect. 2) and multiple (Sects. 3 and 4) national-scale regression models were first investigated. We then introduced four geomorphological classifications to perform local-scale regression analysis in order to tackle the evident spatial clustering of the regression residuals (Sect. 5). The comparisons made between the results obtained from the wide-area and the local regressions allowed the role of the morphology in rainfall variability to be discussed, as shown in Sect. 6. Some conclusions are drawn in Sect. 7.

2 National-scale simple regression analysis

2.1 Methods

As the first step of the analysis, we investigated the influence of elevation on the spatial distribution of the average of annual maximum rainfall depths. We calculated the ADD curve parameters for all the stations of the I²-RED dataset (Mazzoglio et al., 2020a) to compare them with previous studies (mainly Avanzi et al., 2015). Parameters a and n (cf. Eq. 1) were initially obtained by means of a linear regression of the logarithm of the average of all the available extremes over the 1 to 24 h durations. We computed the median values of these parameters, for all over Italy, to compare them with those of Avanzi et al. (2015), who grouped the stations into elevation ranges of 50 up to 1000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ and into intervals of 100 m for higher elevations. We then plotted both series of medians (a and n) against the median elevation of each interval (to consider that the distribution of the rain gauges in each elevation interval was skewed) and fitted regression models.

We studied the differences between the measured and estimated rainfall statistics to assess the effectiveness of the regression models, considering the observed averages of the extremes over 1, 3, 6, 12 and 24 h. We obtained performance indices for each station using the estimation errors Δ_d:

$$\begin{array}{}\text{(2)}& {\mathrm{\Delta }}_d=h_{\text{avg}}\left(d\right)-\widehat{a}\cdot d^{\widehat{n}},\end{array}$$

where h_avg(d) is the sample average of the extreme rainfall depth for duration d, while $\widehat{a}$ and $\widehat{n}$ are the estimates of parameters a and n.

In this paper, we show and discuss only the results related to the shortest and the longest of the five durations (1 and 24 h), as they can be considered the most representative of the different classes of rainfall events (convective and stratiform, respectively). The corresponding dependent variables are called ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{1}}$ and ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{24}}$ in the following.

The error statistics that were computed are the bias, the mean absolute error (MAE), the root mean square error (RMSE) and the Nash–Sutcliffe model efficiency (NSE) coefficient (Nash and Sutcliffe, 1970; Wasserman, 2004). Among all the statistics, particular attention was dedicated to spatial bias, i.e., the bias evaluated as the difference between the spatial mean of the observations over a generic area, and the corresponding values predicted by the model.

2.2 Results

By applying the procedure described in Sect. 2.1, we obtained results that are in agreement with those of Avanzi et al. (2015):

1. Parameter a decreases linearly with elevation (R²=0.89), through the following equation:

   $$\begin{array}{}\text{(3)}& a=\mathrm{30.61}-\mathrm{0.0060}\cdot z,\end{array}$$

   which is comparable with the equation obtained in Avanzi et al. (2015):

   $$\begin{array}{}\text{(4)}& a=\mathrm{29.17}-\mathrm{0.0062}\cdot z.\end{array}$$

2. Parameter n increases nonlinearly with elevation (R²=0.86), through the following equation:

   $$\begin{array}{}\text{(5)}& n=\mathrm{0.54}-\exp [-\mathrm{0.000077}\cdot (z+\mathrm{1650}\left)\right].\end{array}$$

   For comparison purposes, Avanzi et al. (2015) obtained

   $$\begin{array}{}\text{(6)}& n=\mathrm{0.54}-\exp [-\mathrm{0.00086}\cdot (z+\mathrm{1452}\left)\right]\end{array}$$

   for the latter parameter, with only a slight difference in R² (0.89).

   The fitting of the four models is reported in Supplement no. 1.

As already mentioned, parameter a is roughly equivalent to ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{1}}$. Its overall inverse dependence on elevation is somewhat counterintuitive, even though other authors have confirmed this dependence (e.g., Allamano et al., 2009; Marra et al., 2021).

Table 1 Comparison of national-scale error statistics related to the estimates performed with our data and those of Avanzi et al. (2015). The results were obtained with Eqs. (3) to (6).

The error statistics computed on the two sets of residuals, in this work and in that of Avanzi et al. (2015), are listed in Table 1. The results show that the increase in the number of stations and the recording length achieved in I²-RED have led to an improvement compared to the results of Avanzi et al. (2015). This result is not surprising, but more insights can be derived from the observation of the spatial distribution of the residuals, which were not discussed explicitly in the previous literature.

In this regard, we mapped differences Δ₁ and Δ₂₄ to investigate where the under- and overestimations show spatial coherence. The maps, reported in Fig. 2, clearly show that clusters of residuals with high residuals of the same sign emerge in various areas of the country: for instance, many coherent errors larger than 3 times the MAE are present in the Liguria region (see Fig. 1 for the geographic position) for both durations. Therefore, despite the high R² values, significant spatially correlated errors can undermine the practical validity of these general relationships.

Figure 2 Residuals of the estimations of the 1 h (a) and 24 h (b) durations performed using Eqs. (3) and (5). Source: I²-RED data (Mazzoglio et al., 2020a).

On the basis of these results, the need for a more detailed spatial analysis of these variables became evident. A set of new analyses, aimed at reducing the local bias and increasing the reliability of the results, was therefore introduced.

3 National-scale multiple regression analysis

3.1 Methods

In an attempt to improve the evaluation of the relationships between rainfall and topography, we undertook an analysis of the relationships between rainfall and several geomorphological (and climatological) parameters, which may complement the explanatory power of elevation. Unlike what was done in Avanzi et al. (2015), multivariate models were used in the literature to relate rainfall statistics and various morphological variables, both of which were evaluated at the same location. In these approaches, no aggregated or median spatial statistics of rainfall were considered. Prudhomme and Reed (1998, 1999), for instance, identified meaningful geographic and morphological attributes of each location as good explanatory variables of the daily rainfall maxima in Scotland. They showed that obstruction indices, derived from the orography, and the distance from the coastline are able to define how morphological barriers influence the characteristics of the extremes. These appear to work better than the EXPO variable used by Basist et al. (1994) and Konrad (1996).

Basist et al. (1994) defined EXPO as the distance between a rain gauge and an upwind barrier whose elevation is at least 500 m higher than the station. Konrad (1996) suggested an elevation of the barrier at least 150 m higher than the station. Prudhomme and Reed (1998) also tried to use this variable, setting the elevation difference at 200 m, but concluded that the definition of EXPO has several drawbacks, as it is based on arbitrary thresholds and is defined assuming a specific and subjective direction.

Introducing new variables with omnidirectional meaning, as the distance from the sea, the obstruction and the barrier, which are evaluated in the eight main directions, Prudhomme and Reed (1998) were able to explain a much larger percentage of variability in the annual maximum daily rainfall than that explained by the EXPO variable.

Caracciolo et al. (2012) applied this latter approach on the island of Sicily (south of Italy): they found that the longitude, elevation, a barrier obstruction index and the distance from the coastline are able to represent the spatial variability of parameter a for the whole island, while the longitude, elevation, a concavity index and the slope are able to satisfactorily describe the variability of exponent n. They also noticed that different descriptors became significant when analyzing smaller portions of the island.

Based on the above considerations, in this work we followed the approach suggested by Prudhomme and Reed (1998, 1999), considering two groups of variables computed for each station location:

• a.

  Geographic and climatic variables. These do not require computation and do not depend on the landscape forms, that is, longitude (LONG, expressed in the WGS84 UTM32N reference system, in m), latitude (LAT, expressed in the WGS84 UTM32N reference system, in m), elevation above sea level (z, in m), minimum distance from the coastline (C, in km) and mean annual rainfall (MAR, taken from Braca et al., 2021, in mm); the latter represents a very robust climatological variable, which is seldom used as ancillary information but easily available throughout the world thanks to the presence of various rainfall databases (Schneider et al., 2011; Fick and Hijmans, 2017; Muñoz Sabater, 2019).

• b.

  Morphological variables or descriptors. These are based on a digital elevation model (DEM), computed for each cell in a square grid. These variables are as follows:

  • slope (S, in degrees), defined as the angle of the inclination of the terrain to the horizontal, which is evaluated using the eight closest DEM cells;

  • obstruction (OBS, in degrees), defined as the maximum angle needed to overcome the highest orographic obstacles in the eight main cardinal directions (i.e., the maximum of the angles subtended by the line that connects the rain gauge with the highest orographic peak within a 15 km radius in the eight main directions; see Fig. 3);

  • barrier (BAR, in m), defined as the distance between the rain gauge and the highest orographic obstacle defined in OBS (Prudhomme and Reed, 1998, 1999; see Fig. 3);

  • maximum slope angle (MSA, in degrees), i.e., the angle with the greatest slope needed to overcome obstacles within a 15 km radius in the eight main directions (see Fig. 3);

  • maximum slope angle distance (MSAD, in m), defined as the equivalent of BAR but computed with respect to MSA (see Fig. 3);

  • openness (OP, in radians), defined as a mean angular measurement of the relationships between the surface of the relief and the horizontal distances, in the eight main directions (Yokoyama et al., 2002).

Figure 3 Representation of the MSA, MSAD, OBS and BAR morphological variables.

The values of all of these variables depend on the landscape forms and can vary according to the resolution of the used DEM. In our case, after thorough consideration, we adopted the Shuttle Radar Topography Mission (SRTM) DEM, which has a resolution of 30 m (Farr et al., 2007). However, the openness required to be evaluated on the SRTM DEM resampled at a resolution of 500 m due to computational limitations. This computation was conducted with the SAGA “Topographic openness” module, using a radial limit of 5 km. This value was obtained after testing different radial limits and selecting the one that presents the best correlation with the mean rainfall depth.

Multiple linear regression models were built, based on the following relationship:

$$\begin{array}{}\text{(7)}& Y=\mathbf{X}\cdot \mathit{\beta }+\mathit{\delta }=\sum _{i=\mathrm{1}}^N{\mathit{\beta }}_i{X}_i+\mathit{\delta },\end{array}$$

where the dependent variable Y is related to the matrix of the independent variables X or covariates. β in Eq. (7) is the vector that contains the regression model coefficients, and δ is the vector of the residuals.

In order to select the best model equation, the number i of covariates can be increased as necessary, according to the criteria of statistical significance of the estimated parameters. Caracciolo et al. (2012), for instance, used a stepwise regression approach. In this paper, we have preferred to use a generalized multiple regression approach, whereby increasing the number of covariates to i+1 does not necessarily preserve the descriptors that were the most significant at step i. This approach entails always considering all the possible combinations of two, three or four covariates until the “best” model is found. Other tests made using five or more variables did not lead to significantly higher $R_{\text{adj}}^{\mathrm{2}}$ values.

The best regression model was selected on the basis of an analysis of the regression residuals, favoring models with the highest adjusted coefficient of determination, $R_{\text{adj}}^{\mathrm{2}}$. Student's t test was used to quantify the significance of the independent variables. We checked for the possible presence of multicollinearity for each model in which all the covariates passed Student's t test, as this could lead to the formulation of an unstable model. Multicollinearity was measured using the variance inflation factor (VIF), which is determined by placing the jth independent variable as the dependent variable and calculating the coefficient of determination $R_{\text{adj}}^{\mathrm{2}}$ of the multiple regression performed on the remaining p−1 independent variables (Eq. 8).

$$\begin{array}{}\text{(8)}& {\text{VIF}}_j={\displaystyle \frac{\mathrm{1}}{\mathrm{1}-R_j^{\mathrm{2}}}}\end{array}$$

Values of VIF greater than 5 were associated with an unacceptable level of multicollinearity, and the corresponding model was discarded (Montgomery et al., 2012).

3.2 Results

The equations of the best regression models (built using two to four variables) are reported in Eqs. (9) to (11) for ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{1}}$ and in Eqs. (12) to (14) for ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{24}}$. The $R_{\text{adj}}^{\mathrm{2}}$ values in Eqs. (9) to (11) are 0.46, 0.52 and 0.54, respectively. The coefficients of determination are higher for ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{24}}$, i.e., 0.66, 0.67 and 0.68 (Eqs. 12 to 14, respectively).

$$\begin{array}{}\text{(9)}& {\displaystyle }{h}_{\mathrm{1}}=\mathrm{20.3163}-\mathrm{0.0080}\cdot z+\mathrm{0.0117}\cdot \text{MAR}\text{(10)}& {\displaystyle }\begin{array}{rl}{h}_{\mathrm{1}}=& -\mathrm{21.6293}-\mathrm{0.0061}\cdot z+\mathrm{0.0134}\cdot \text{MAR}\\ & +\mathrm{26.2682}\cdot \text{OP}\end{array}\text{(11)}& {\displaystyle }\begin{array}{rl}{h}_{\mathrm{1}}=& -\mathrm{10.6928}-\mathrm{0.0051}\cdot z-\mathrm{0.0273}\cdot C\\ & +\mathrm{0.0131}\cdot \text{MAR}+\mathrm{19.8449}\cdot \text{OP}\end{array}\text{(12)}& {\displaystyle }{h}_{\mathrm{24}}=\mathrm{16.1392}-\mathrm{0.0937}\cdot C+\mathrm{0.0712}\cdot \text{MAR}\text{(13)}& {\displaystyle }\begin{array}{rl}{h}_{\mathrm{24}}=& \mathrm{33.1529}-\mathrm{1.8574}\cdot {\mathrm{10}}^{-\mathrm{5}}\cdot \text{LONG}-\mathrm{0.1319}\cdot C\\ & +\mathrm{0.0701}\cdot \text{MAR}\end{array}\text{(14)}& {\displaystyle }\begin{array}{rl}{h}_{\mathrm{24}}& =\mathrm{127.2773}-\mathrm{2.9498}\cdot {\mathrm{10}}^{-\mathrm{5}}\cdot \text{LONG}\\ & -\mathrm{1.9130}\cdot {\mathrm{10}}^{-\mathrm{5}}\cdot \text{LAT}-\mathrm{0.0971}\cdot C+\mathrm{0.0735}\\ & \cdot \text{MAR}\end{array}\end{array}$$

Considering the three ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{1}}$ models (Eqs. 9 to 11), it is possible to notice the negative slope coefficient associated with elevation. This confirms what was discussed in the previous section. On the other hand, it is remarkable that the best models found for ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{24}}$ do not include the elevation: this outcome can be explained by considering the fact that MAR is always significant, regardless of the number of variables. Models in which MAR was excluded actually present z as a significant covariate but with less relevance than the regression models for ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{1}}$.

Regardless of the number of the variables considered, and despite the marked increase in the corresponding value of $R_{\text{adj}}^{\mathrm{2}}$, compared to the simple regression, we found that the residuals of the multivariate regressions were still characterized by spatial clustering and high local errors, basically all in the same areas in Fig. 2. In other words, resorting to additional variables but keeping a uniform relationship between each variable and precipitation over all of Italy does not produce a decisive reduction in the bias for large areas of the country. Thus, we decided to deconstruct the modeling approach and to look for clues of distinct generating mechanisms in distinct areas of Italy.

4 Sub-national-scale multiple regression analysis

4.1 Methods

In this section, an additional paradigm is introduced into the models for the spatial variability of precipitation to reduce the spatial bias, namely the selection of limited areas to build “local” regression models, as an alternative to using data for the whole of Italy. Such an attempt was already made by Caracciolo et al. (2012), who borrowed the subdivision criterion from previous regional frequency analyses. In this work, we have focused on the role of geography and morphology in the spatial variability of annual maximum rainfall depths.

To better understand how to move from national-scale relationships to relationships valid for smaller areas, we started by considering the Alpine area separately from the Apennine region along the entire peninsula, and from the two main islands (Sicily and Sardinia; see Fig. 1 for the geographic positions), as a first approximation. We then built four different multivariate models: (1) the Alpine region (i.e., from Piedmont, including the western part of Liguria, eastward up to Friuli Venezia Giulia, delineated using the SOIUSA classification, as suggested by Accorsi, 2016); (2) the Apennine region, including peninsular Italy; (3) Sicily; and (4) Sardinia. We evaluated the best regression models for these four regions, as described in Sect. 3.1, using up to four covariates.

Table 2 Error statistics of the multiple regression models at a national scale and for the four macro-regions described in Sect. 4.1, for ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{1}}$ and ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{24}}$.

4.2 Results

The new set of models built for the four regions were tested by computing the error statistics over the entire country. The obtained results indicated that they provided higher $R_{\text{adj}}^{\mathrm{2}}$ than for the national case and better error statistics (see Table 2 for a comparison with the previous multivariate approach). The better results achieved in terms of RMSE, MAE and NSE at the national scale are due to the improvements obtained for the two main islands (Sicily and Sardinia). More insights are provided in Sect. 6.

It is interesting to compare the results obtained for the individual Alpine region with those of Allamano et al. (2009), who analyzed almost the same area. In that case, the ADD curve parameters appeared to be related to elevation and longitude. For the different durations Allamano et al. (2009) also estimated a regression model by linear regression between rainfall depth, elevation and longitude. The dependence of short-duration rainfall on elevation and longitude was found to be statistically significant for all the time intervals, except for the 1 h duration: in this case, the longitude was not statistically significant. In our application, the best relationships found for ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{1}}$ and ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{24}}$ are those of Eqs. (15) and (16) (characterized by $R_{\text{adj}}^{\mathrm{2}}=\mathrm{0.75}$ and 0.76, respectively):

$$\begin{array}{}\text{(15)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{‾}}{h}}_{\mathrm{1}}=& \mathrm{60.9365}-\mathrm{1.6664}\cdot {\mathrm{10}}^{-\mathrm{5}}\cdot \text{LAT}-\mathrm{0.0046}\cdot z\\ & +\mathrm{0.0148}\cdot \text{MAR}+\mathrm{25.1825}\cdot \text{OP}\end{array}\text{(16)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{‾}}{h}}_{\mathrm{24}}=& \mathrm{59.0632}-\mathrm{7.2955}\cdot {\mathrm{10}}^{-\mathrm{5}}\cdot \text{LONG}-\mathrm{0.2223}\cdot C\\ & +\mathrm{0.4306}\cdot \text{OBS}+\mathrm{0.0822}\cdot \text{MAR}.\end{array}\end{array}$$

As expected, the ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{1}}-z$ relationship has a negative slope, and the Eq. (15) does not include the longitude as covariate, in agreement with Allamano et al. (2009). The same negative relationship is found in a 24 h equation that includes z (which produces an $R_{\text{adj}}^{\mathrm{2}}=\mathrm{0.74}$, that is, lower than that of Eq. 16). These findings confirm the results of Allamano et al. (2009), who found a general decrease in rainfall depth for increases in elevation for all the durations (up to 24 h). Equation (16) also confirms that, although ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{1}}$ decreases systematically with elevation over the whole alpine region, the dependence of ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{24}}$ on z decreases as the longitude increases, i.e., moving westward.

The full set of equations used for the four regions is provided in Supplement no. 2, together with the $R_{\text{adj}}^{\mathrm{2}}$.

Although the improvements achieved with multivariate models over the simple regressions are evident, we found that they were not decisive in providing a homogeneous spatial distribution of the errors. We in fact observed that, even with the best model, we were not able to reduce the clustering effect shown in Fig. 2 for the peninsular region (see also Supplement no. 3). We believe that a model capable of describing the observed spatial variability of the index rainfall simultaneously at a national and a local level requires additional insights, which can be obtained using a finer spatial segmentation of Italy.

5 Local-scale simple regression analysis of morphological regions

5.1 Methods

On the basis of the considerations presented above pertaining to the spatial clustering of residuals, we examined the possibility of obtaining a meaningful segmentation of large areas in subdomains that could be used to obtain “local” relationships between annual maximum rainfall depths and terrain properties. The main reasoning behind the segmentation is that some macroscopic morphological differences can determine markedly different behaviors of the relationships between rainfall and elevation (or other local variables). One example concerns what happens in the windward and leeward sides of mountain ridges, which represent transversal obstacles to the humid masses coming from the sea. Accordingly, we considered some general geomorphological classifications of the landscape that delineate homogeneous areas based on the homogeneity of the macroscopic land properties, such as convexity and texture.

We considered four geomorphological classifications (GCs) and denominated them as GC1 to GC4, according to their diversity and success in the geomorphological literature (see the Data availability section for more information).

The first considered classification, called GC1, was proposed by Iwahashi and Pike (2007); they classified the Earth's surface into 16 topographic types, at a 1 km resolution, based on slope gradient, local convexity and surface texture. We vectorized the raster map, which is available on the European Soil Data Centre website, and then, to reduce the presence of small areas, which could have an extent of just a few square kilometers, all the areas covered by fewer than 10 pixels (10 km²) were merged with the adjacent class. Among the four different classifications that were used, this is the only one that has a worldwide coverage, as all the other classifications are available at a national scale. A detailed description of the methodologies used by the authors is available in the related references, thus allowing all the classifications to be reproduced over other nations.

The second classification – GC2 – is the Carta delle Unità Fisiografiche dei Paesaggi italiani (“Map of the physiographic units of Italian landscapes”) and is included in the Carta della Natura (“Map of nature”; Amadei et al., 2003). A vector map, which was obtained by means of a visual interpretation of satellite images, aided by the analysis of further land cover maps and morphological–lithological characteristics, was available at a 1:250 000 scale.

The third classification – GC3 – was proposed by Guzzetti and Reichenbach (1994). It was obtained, in vector format, by combining an unsupervised three-class cluster analysis of four properties of altitude (altitude itself, slope curvature, frequency of slope reversal and elevation–relief ratio) with a visual interpretation of morphometric maps and an inspection of geological and structural maps.

The fourth classification – GC4 – is the one that delineates areas with the greatest detail, as it is based on local morphometric properties of the landscape. It was proposed by Alvioli et al. (2020a), who considered a set of 439 watersheds, covering the whole of Italy, grouped into seven clusters on the basis of the various properties of the slope units within each basin, e.g., a distribution of slope unit sizes and aspects. In this work, adjacent watersheds of the same class were collapsed (GIS Dissolve), thus producing a total of 178 areas. Geomorphologically homogeneous terrain partitions were defined as “slope units” that were delimited by drainage and divide lines and delineated with a method that was first introduced by Alvioli et al. (2016) and which is widely used in the literature for geomorphological zonation purposes.

An additional geomorphological classification, which was proposed by Meybeck et al. (2001) and which has a worldwide coverage, was also considered. It is based on a combination of a relief roughness index and elevation and in principle could have been a good fifth candidate. However, it was not included in this analysis because, except for a very large geographical zones, the resulting delineated areas often contained very few rain gauges, which would have made it impossible to perform the desired statistical analyses.

Coherently with the aim of addressing connections between terrain properties and rainfall at a more local level, we built a set of linear regression models between elevation and index rainfall for all the classifications, considering an individual model for each outlined geomorphological zone. Only the internal rain gauges in each of these homogeneous areas with a minimum of five available stations that had to ensure at least 100 m of difference in elevation were considered for the regressions. There were four possible outcomes of the applications: (a) a positive and significant correlation (at the 5 % level), (b) a negative and significant correlation; (c) a nonsignificant correlation and (d) an insufficient number of stations or an insufficient difference in elevation.

5.2 Results

The results for ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{1}}$ are presented hereafter, while the details on ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{24}}$ are available in Supplement no. 4. The results obtained for each geomorphological zone are mapped in Fig. 4a–d. Blue areas denote geographical zones, where ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{1}}$ increases together with elevation, while the red palette applies to zones where rainfall decreases with elevation. The color intensity is proportional to the respective slopes. The light gray color denotes zones in which the linear regression is not statistically significant (at a 5 % level), while dark gray denotes insufficient data (case d). A comparison of the maps (Fig. 4a–d) clarifies that the more detailed the geomorphological zonation is, the less likely it is to satisfy the requirements necessary to build a significant regression. On the other hand, if one applies regression models to finer geomorphological classifications, it is possible to see that the regression sign is not uniform over the entire country. For example, with regard to 1 h data, it is possible to clearly recognize the presence of zones with a positive rainfall depth versus elevation trend for pre-hill/plain morphology in both GC3 (Fig. 4c) and GC4 (Fig. 4d).

Figure 4 Slope coefficients of the regression between the mean 1 h rainfall depth and elevation for GC1 (a), GC2 (b), GC3 (c) and GC4 (d) and the R² of the regression between the mean 1 h rainfall depth and elevation for GC1 (e), GC2 (f), GC3 (g) and GC4 (h). Geomorphological data source: Iwahashi and Pike (2007), Amadei et al. (2003), Guzzetti and Reichenbach (1994) and Alvioli et al. (2020a).

The spatial distribution of the light gray zones is an important piece of information: no trend can be assumed over these areas because the p value is greater than 0.05. Consequently, ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{1}}$ can be considered constant over these areas. Finally, the occurrence of the dark gray zones is directly connected to the kind of classification: the smaller the areas delineated by the classification are, the more likely it is that the requirement of having at least five rain gauges with at least 100 m difference in elevation is not satisfied. In this regard, we can observe that the above requirements are not met for the elevation difference, i.e., in plain areas, and it is necessary to assume a constant ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{1}}$ in the area as being the most reasonable value.

The maps in Fig. 4a–d show that the availability of more detail in the spatial analysis of the relationship between rainfall depth and elevation has a remarkable effect on both the sign of the regression and the slope of the regression line in several areas. In addition, even the quality of the relationship can improve, as can be seen from a comparison of Fig. 4e–h: far more areas with high R² can be seen in Fig. 4h than in Fig. 4e. This allows us to conclude that lower values of R² are obtained in wider areas.

The same analysis was conducted on ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{24}}$, where all of the above outcomes were confirmed, except for the sign of the precipitation vs. elevation relationship (Supplement no. 4).

Table 3 National-scale error statistics for the 1 h interval. Statistics for GC1, GC2 and GC3 were only evaluated over areas where the regression was statistically significant at a 5 % level, while GC4 was tested on both statistically significant areas and over the entire nation, using the mean rainfall, where there was a p value > 0.05 or where the requirement of at least five rain gauges with at least 100 m difference in elevation was not satisfied.

Table 4 National-scale error statistics for the 24 h interval. Statistics for GC1, GC2 and GC3 were only evaluated over areas where the regression was statistically significant at a 5 % level, while GC4 was tested on both statistically significant areas and over the entire nation, using the mean rainfall, where there was a p value > 0.05 or where the requirement of at least five rain gauges with at least 100 m difference in elevation was not satisfied.

5.3 Error analysis

To test the reliability of the regression models built over the GCs, the linear equations found in each geomorphological zone were applied to all the rain gauge positions, to obtain errors that could be examined at the country scale. The global indices computed for the GC areas in which the regressions were statistically significant are reported in Table 3 for ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{1}}$ and in Table 4 for ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{24}}$. These results clearly show a lower performance of the GC1 than the national-scale regression model. On the other hand, the error statistics in Tables 3 and 4 show that GC4 produces the smallest errors, and this geomorphological subdivision therefore presents the best performances. It is possible to understand this result by considering that GC4 uses watershed units, while the other classifications are based on the automatic processing of digital terrain data.

A constant value of the index rainfall computed as the spatial average of ${\stackrel{\mathrm{‾}}{h}}_d$ was adopted over any areas where the morphological regression model was not statistically significant. In this way, it was possible to compute the statistical indexes at the whole country scale. The error statistics obtained with this application are reported in the last rows of Tables 3 and 4.

6 Discussion

The different regression models used in this work to investigate the role of morphology in the spatial distribution of sub-daily annual maximum rainfall depths produced results deserving some comments. First of all, it must be mentioned that a nationwide multiple regression model that includes morpho-climatic attributes represents a significant step forward with respect to the simple regression model, as the error statistics show. In this approach, working at a national scale and given the elongated shape of the Italian peninsula, geographic location was expected to play a major role in the spatial distribution of extremes, even though this evidence was not mentioned in similar national-scale analyses (see, e.g., Faulkner and Prudhomme, 1998, for the UK, and Avanzi et al., 2015, for Italy). The role of geography progressively weakened while seeking further improvements, in terms of MAE and RMSE, through the application of distinct multiple regressions to four macro-regions, i.e., the Alps, peninsular Italy and the main islands (Supplement no. 5).

Our findings show that while the 24 h index rainfall exhibits a clear overall dependence on the geographic location at a full national scale (Eq. 14), the same does not apply to 1 h extremes (Eq. 11). In an area with a lesser span in latitude (the Italian Alps), instead, the 1 h extremes curiously show some dependence on latitude (Eq. 15). Formetta et al. (2022) followed a similar reasoning, recognizing the role of geography and elevation as they partitioned by longitude and elevation even a small area (the province of Trento) before applying their statistical analyses.

While the multivariate regression can be a good tool to express geographic dependence, and on 24 h extremes the national scale helps in drawing some general findings, the residual errors in large clustered areas are still very significant. Therefore, geographic attributes seem not to drive uniformly the variability of rainfall extremes all over Italy, as the high residuals of the multiple regression over these areas do not apparently follow any latitudinal/longitudinal gradient. These findings can only derive from a national-scale analysis.

The better suitability of the application of multiple regressions on four regions is confirmed by the increase of the adjusted coefficient of determination ($R_{\text{adj}}^{\mathrm{2}}$), as reported in Sect. 3.2 and in Supplement no. 2. Moreover, while the national-scale multiple regression model provides high negative residuals over Sardinia and high positive residuals over Sicily, the four-region multiple regression model significantly improves this result (see Supplement no. 3–5 for more details). However, similar improvements were not achieved in the peninsular and alpine areas of the country.

Table 5 Error statistics for the 1 and 24 h intervals at a national scale. Average spatial values are used for the gray areas in Fig. 4d. The bias of the national simple regression is different from zero, being evaluated as ${\text{bias}}_d=\frac{\mathrm{1}}{n}\cdot \sum h_{\text{avg}}\left(d\right)-\widehat{a}\cdot d^{\widehat{n}}$.

The subsequent investigations undertaken in Sect. 5 descend from the above considerations; i.e., the building of regressions in morphological regions that are a fraction of the whole area is an attempt to overcome the highlighted lack of regularity in the dependence between rainfall and geography. Among all the considered geomorphological classifications, the selection of rain gauges for the model application is more effective in the case of GC4 (Alvioli et al., 2020a), which also embeds hydrographic information. The GC4 model behaves reasonably well for both the 1 and 24 h durations, compared to the multiple regression models, as far as the national scale is considered. Table 5 summarizes all the previously mentioned statistics.

Analyzing the error statistics computed globally at the national scale, it seems that the four-region multiple regression approach is the most precise. However, this is not necessarily true at a local scale. In order to clarify the drawbacks that large-scale regression models can produce, for the 1 h case, we compared the residuals obtained from the four-region multiple regression model (Fig. 5a) in the areas identified by GC4 with the residuals of the GC4 regression model by selecting: (1) the GC4 areas that were statistically significant (Fig. 5b) and (2) the entire nation (Fig. 5c). The mean rainfall depths were considered over nonstatistically significant areas (i.e., the gray areas visible in Fig. 5b). The GC4 regression models were statistically significant for ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{1}}$ for 45 % of the Italian area, for a total of 31 different areas, while the GC4 model for ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{24}}$ were statistically significant in 49 % of the area, for a total of 47 different zones (figures dealing with the 24 h case are available in the Supplement). From a comparison of the maps in Fig. 5, it is possible to note that the multiple regression model has a spatially nonuniform bias while the average bias obtained from the individual models in the zones selected by GC4 is zero all over Italy. Maps of all the other statistics are reported in Supplement no. 5 (Fig. S4 in the Supplement). This outcome is also evident for the 24 h case (see, for example, the maps of the bias – Fig. S5a–e in the Supplement).

Figure 5 The 1 h case. Local bias for the four-region multiple regression model (a), the GC4 simple linear regression model over statistically significant areas (b) and the GC4 simple linear regression model over all the areas (c). Geomorphological data source: Alvioli et al. (2020a).

An additional comparison was undertaken to investigate the local bias. In this case, we computed the bias for each subdomain of GC4. We compared the bias values obtained using the following four conditions: (1) the national-scale simple regression model (Eqs. 3 and 5), (2) the national-scale multiple linear regression model (Eqs. 11 and 14), (3) the four-region multiple regression model (Eqs. S3, S6, S9, S12, S15, S18, S21 and S24 in Supplement no. 2) and (4) the GC4 simple regression model (Sect. 5).

The results are illustrated in the maps of Fig. 6, which shows the best regression model for each area in different colors: Fig. 6a and b are related to ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{1}}$, while Fig. 6c and d refer to ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{24}}$; Fig. 6a and c only highlight the situations where significant regressions were found. The results in Fig. 6b and d include the bias calculated in nonstatistically significant areas with respect to the spatial average of the rainfall depths. The good results obtained in the areas where the spatial mean values are adopted can be seen by comparing the borders of the GC4 areas with the clusters of the residuals of the multiple regression model (see Supplement no. 3). A dedicated multiple regression model was built for the island of Sardinia: nevertheless, the bias all over the GC4 areas is smaller when the local spatial average is used. A good correspondence between the residual clusters and the GC4 borders is evident.

From Fig. 6, it is possible to conclude that the morphological subdivisions allow a set of simple linear regression models to be built that can perform better almost everywhere than the other wide-area models in terms of local bias.

Figure 6 Absolute bias assessment for all the regression models used for the 1 h case (a, b) and 24 h case (c, d). The color refers to the model that provides the lowest absolute value of the bias. The GC4 model bias in cases (a) and (c) was only evaluated for statistically significant areas, while it was evaluated over every area in (b) and (d). Geomorphological data source: Alvioli et al. (2020a).

7 Conclusions

In this paper, we have analyzed the role of orography and morphology in short-duration annual maximum rainfall depths, taking advantage of a new and comprehensive database for Italy, I²-RED (Mazzoglio et al., 2020a). The approach finds its relevance in the first use of the most complete and updated data collection of short-duration annual maxima available for the whole Italian territory. As regards to the previous knowledge on the topic, our analyses allowed us to better understand, confirm and extend previous results from the literature.

The results described in this paper show that a national-scale simple regression model of the precipitation vs. elevation presents some weaknesses (high residual values, high local- and national-scale bias, and low NSE coefficient, etc.) and therefore needs to be improved.

The use of multiple regression models introduces some benefits, such as a reduction of MAE and RMSE at the national scale; nevertheless they were not successful in reducing the local bias.

Considering the necessity of working on smaller domains, we analyzed several geomorphological classifications which are able to preserve the intrinsic value of the statistically significant landscape variables that emerge in regression models. Four different geomorphological classifications available in literature were used to provide criteria for the identification of homogeneous regions. We applied simple linear regression models over these homogeneous domains and compared the performances at both a national and a local level. Among all the considered classifications, the selection of rain gauges for the model application was found to be more effective in the case of GC4 (Alvioli et al., 2020a), which embeds hydrographic information.

The best approach was selected by evaluating the error statistics for the bias at both a national and a local scale and at a national scale for MAE, RMSE and NSE. The obtained results have shown that using simple linear regression applied to the GC4 model produces better results than all the others, in the areas in which the GC4 model is statistically significant, in terms of bias. As far as national statistics are concerned, considering the mean rainfall depths in the gray areas in Fig. 5b does not significantly affect the performance of GC4, in terms of MAE, RMSE and NSE, in particular for the 1 h duration. In short, we propose using the GC4 model where possible and adopting the (spatial) mean value of the rainfall depths in case of a nonstatistically significant relationship.

This work has led to the following conclusions. The relationship between precipitation and elevation is not meaningful in all the areas in Italy, as already pointed out by Caracciolo et al. (2012) for the island of Sicily. In this work, this concept has systematically been extended to the whole country, and significant relationships have only been obtained for 45 % of the area for ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{1}}$ and 49 % for ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{24}}$. As far as the model that we suggest using is concerned, that is GC4, we are aware that improvements are possible, considering that no significant regressions were found over 55 % (${\stackrel{\mathrm{‾}}{h}}_{\mathrm{1}}$) and 51 % (${\stackrel{\mathrm{‾}}{h}}_{\mathrm{24}}$) of the territory. However, it should be pointed out that the rainfall station density is not sufficient for the application of the method proposed here over 9 % of the territory.

Details regarding the model based on GC4 and numerical values of the regression parameters are provided in the Data availability section.