Figure 1 . a)Location and orography of the study area and location of the rain gauges used in this study; b) Decorrelation time of the highest 25% ordinary events organized by season. The red dots indicate the median values; bars indicates percentiles: 25-75th, 5-95th, 1-99th. The number of storms occurred across the stations in each season is reported.
2.2 Definition of the ordinary events
Ordinary events are all the independent realizations of a process of interest, in our case precipitation intensities at multiple durations. The here presented analysis is based on the storm-based identification of ordinary events proposed by Marra et al. (2020), in which “storms” are defined as independent meteorological objects, and “ordinary events” of each duration are extracted from the storms. For each station, storms are defined as wet periods separated by dry hiatuses of predefined length. We define as wet all the 5 min time intervals reporting at least 0.1 mm of precipitation, and separate storms using 24 hr dry hiatuses. A minimum duration of 30 min for a single storm is set to avoid individual tips to be considered as storms. Ordinary events are then defined as the maximum intensities observed over the duration of interest in each storm (details in Marra et al., 2020). Durations between 15 min to 24 hr are explored: 15 min, 30 min, 45 min, 1 h, 2 h, 3 h, 6 h, 12 h, 24h.
2.3 Tail of the ordinary events distribution
Previous studies show that subdaily precipitation intensities require three- (or more) parameters distributions (Papalexiou et al., 2018). However, their right tails can be well approximated using a two-parameter distribution which, in many cases, is found to be a Weibull distribution (e.g., Zorzetto et al. 2016; Marra et al., 2020). This means that a portion of their distribution including the extremes, which is here termed “tail”, can be approximated as\(\text{\ F}\left(x;\lambda,\kappa\right)=1-e^{-\left(\frac{x}{\lambda}\right)^{\kappa}}\), where λ is a scale parameter and κ is a shape parameter which determines the tail heaviness. Larger shape parameters are associated to lighter tails, and vice versa (see Figure S1). In particular, the tail is sub-exponential for κ > 1, exponential for κ=1, and heavier than exponential for κ < 1.
The choice of the left-censoring threshold follows the test described in Marra et al. (2020): the distribution parameters are estimated for different thresholds by censoring the values below the left-censoring threshold as well as the observed annual maxima. The maxima are then compared to the sampling confidence interval from the estimated distribution to assess whether they could be likely samples. Following the method suggested in Marra et al. (2019), we select the 75th percentile of the ordinary events for the left-censoring. This is in line with previous findings in areas dominated by convective processes (Marra et al., 2019; Marra et al., 2020). It should be recalled that the selection method implies a low sensitivity of the results to this threshold.
2.4 Extreme value model
The cumulative distribution \(\zeta(x)\) of extreme return levels \(x\)emerging from the underlying distribution of ordinary events with tail\(F\left(x;\lambda,\kappa\right)\) can be written as\(\zeta\left(x\right)=\ F\left(x;\lambda,\kappa\right)^{n}\), where \(n\) is the average number of ordinary events per year (Marra et al., 2019; Serinaldi et al., 2020). When one considers the j-th year of data, this formalism allows us to quantify return levels from individual years by inverting\(\zeta_{j}\left(x\right)=F\left(x;\lambda_{j},\kappa_{j}\right)^{n_{j}}\), where \(\lambda_{j}\) and \(\kappa_{j}\) are the parameters describing the ordinary events tail at the j-th year and \(n_{j}\) is the number of ordinary events in the year.
The parameters describing the ordinary events distribution tail are computed at each station, duration and year by left-censoring the lowest 75% of the ordinary events and using the least-squares method in Weibull‑transformed coordinates (Marani and Ignaccolo, 2015). After left-censoring, an average of ~14 ordinary events per year (including annual maxima) are used for parameter estimation. Yearly return levels are obtained by inverting the equation for\(\zeta_{j}\left(x\right)\). In this way, we obtain, for each station, yearly series of scale parameter, shape parameter, number of ordinary events, and return levels. Annual maxima (AM) series are also extracted.
2.5 Temporal trends analysis
We investigate the presence of monotonic trends in these quantities using the Regional Mann-Kendall test at the 0.05 significance level (Mann, 1945; Kendall, 1975; Helsel & Frans, 2006), and we quantify the average rate of change using the nonparametric Sen’s slope estimator (Sen, 1968). Serial correlation in the series was tested and found negligible. In case trends within the region are heterogeneous, the slope and significance estimated by the Regional Mann-Kendall test could be misleading (Gilbert, 1987). We verify the homogeneity of the trends at the different sites in the area by applying the Van Belle and Hughes test (1984). We find that homogeneity is verified for all the investigated variables. As spatial correlation among nearby stations could decrease the power of regional test, we include the correction proposed by Hirsch and Slack (1984).
If the null hypothesis of the Mann-Kendall test is true (i.e., no trend) about half of the pair comparisons between ordered data points is concordant and half discordant. . Considering that 2 yr return levels correspond to the theoretical median of the AM, we consider the estimated trend on the 2 yr return levels as our model quantification of the trend in the AM.
2.6 Validation of the statistical model
The ability of our statistical model to reproduce observed trends in AM is verified by accounting for stochastic uncertainty in a Monte Carlo framework. For each station i , year j and durationd , nijd Weibull-distributed ordinary events are generated according to the distribution parametersλijd and κijd , and the AM are extracted. The procedure is iterated 1000 times (which was found to provide coherent estimates of the 90% confidence interval), to obtain 1000 synthetic regional sets of AM series for each duration. The Regional Mann-Kendall test is then performed on these sets to obtain 1000 slopes estimates for each duration, which provide a quantification of the stochastic uncertainty in the trends of the modelled AM. It is worth noting that this confidence interval is obtained by neglecting spatial correlation in the local exceedance probability of the events, and it is thus to be considered as a lower limit to the true confidence interval. In fact, such a correlation would cause a loss of information in the regional pooling of the trend test, inflating the stochastic uncertainty in the outcome.
2.7 Differential impact of ordinary events change on annual maxima changes
The relative impact of trends in the ordinary events characteristics and frequency on the emerging trend in the AM is evaluated. For each station and duration, the trends on modelled AM are computed using different combinations in which inter-annual variability in the parameters is either considered or ignored. In the latter case, the median parameter is used. We thus obtain the following cases: one case with 3 time-varying parameters (real case), 3 combinations of 2 varying and 1 constant parameter, 3 combinations of 1 varying and 2 constant parameters, and one case of 3 constant parameters (no-change). Then the Regional Mann-Kendall test is applied to the resulting series.
2.8 Changes in the proportion of convective-like and other types of storms
Changes in the seasonal proportions between convective-like and other event types in different seasons are explored to investigate the seasonal and physical mechanisms underlying the observed trends. Events exceeding the left-censoring threshold at any of the durations are organized by seasons. The temporal decorrelation of the rain intensity timeseries is used as a proxy for broadly distinguishing between convective-like and other types of storms. The decorrelation time (Figure 1 b) is taken equal to the scale parameter of the exponential fitting of the temporal autocorrelation. This is thus the time lag at which the temporal autocorrelation drops to e-1. For each station and season, the yearly number of storms belonging to the two groups is calculated, and the significance and slope of the regional trend is estimated using the Regional Mann-Kendall test (p=0.05) and the Sen’s slope estimator. This shows if temporal changes in the proportion of different event types in the seasons emerged. A 2 hr threshold is found to optimally describe (that is, optimize the statistical significance) the temporal changes in our data and is therefore used as a proxy for distinguishing between convective‑like (decorrelation time ≤ 2 hr) and other event types (> 2 hr). Qualitatively analogous outcomes are obtained with thresholds between 1 and 3 hr.
3 Results and discussion
3.1 Regional trends on multi-duration extremes
Slopes for the regional trends for the nine investigated durations are reported in Figure 2 a. Hereinafter, slopes are normalized over the median value of each variable and expressed as percent change per year. As expected (Libertino et al., 2019), observed AM show positive trends at all durations. Statistically significant trends are observed for durations up to 6 hours and stronger increases for hourly and sub-hourly durations. The slopes estimated using the model (“modelled AM” in Figure 2 ) lie within the 90% confidence interval due to stochastic uncertainty (grey area), with the exception of the longest durations (12 and 24 h). Since at longer durations, the confidence interval is likely underestimated due to a larger correlation in the severity of the storms, this indicates that they are likely samples from our model. This means that the model well reproduces the trends in the observed AM.
The annual number of storms, uniquely defined for all durations (Marra et al., 2020), shows an increase (0.4% yr-1) (Figure 2 b). Trends in the scale parameter of the intensity distributions are always positive, indicating a general increase in the intensity of the largest 25% of the ordinary events, with larger and significant increases (up to 1.0% yr-1) for multi-hour durations (Figure 2 b). The shape parameter shows negative trends for sub-hourly durations and positive trends for longer durations (Figure 2 b), indicating that the proportion between heavy and mild events changed in different ways for short and long durations: increased tail heaviness is reported for sub-hourly durations and decreased tail heaviness for multi-hour durations (see Figure S1 for a visual interpretation of the effect of the shape parameter on tail-heaviness). At short durations the changes in the two parameters have a synergistic impact on extremes. Although the trend in individual parameters is not significant, observed and modelled AM experience stronger and significant changes. In contrast, at longer durations the changes in the parameters have opposing impact on extremes, and AM exhibit weaker increases, despite the increase of both scale parameter (significant) and yearly number of storms. In particular, where tail-heaviness has its strongest decrease (increase in the shape parameter), trends in extremes are at a minimum and are not significant.
These findings indicate that in the examined period (1991-2020) and area, AM exhibit significant changes, in particular for short-duration intensities, in agreement with previous studies (Libertino et al., 2019). Overall, our statistical model reproduces these trends accurately, and allows us to investigate the underlying statistical mechanisms. Changes in AM seem to be mostly influenced by changes in the tail-heaviness of the ordinary events, although trends in the shape parameter itself are not statistically significant.
3.2 Differential impact of ordinary events change on annual maxima changes
We investigate the impact of the trends in the individual model parameters on the trends in AM (Figure 2 c). The ‘real’ case in which all parameters change with time reproduces the trends in the modelled AM (line 1 in Figure 2 c). The other lines are a combination of varying and constant (median) parameters. Notably, the increase (+0.4% yr-1) in the number of yearly storms only has a marginal impact on the overall trends in extremes (same-color pairs of lines). Synergistic and opposing impacts of the other parameters are mostly evident by comparing the constant scale-parameter case (line 2) with the constant tail-heaviness case (line 3). When no changes in tail-heaviness are considered, AM show increasing trends whose magnitude can even increase with duration, instead of decrease (lines 3, 7). This analysis shows that little changes in the tail-heaviness (shape parameter) turn into large changes in extreme intensities, suggesting this is an important parameter explaining the observed AM trends in the region. Crucially, without considering changes in tail heaviness it is not possible to explain the large observed increase in short-duration AM, as well as the different response of short and long duration extremes. This has profound implications for change-permitting extreme value models in which tail heaviness is often assumed to remain constant.
3.3 Regional trends of extreme return levels
Our statistical model allows to directly quantify changes on specific rare return levels. In general, slopes are always significantly positive for sub-hourly durations and decrease with increasing duration until they lose significance for durations above 2-3 hr (Figure 2 d). For higher return levels, this behavior is enhanced: higher positive slopes are estimated for sub-hourly durations and lower not significant slopes for multi-hour durations. There is a duration interval between 1 and 2 hr where the trends don’t depend on return period, closely following the change in regime in which the trend in the shape parameter crosses zero, that is no change in tail heaviness.
The here adopted statistical framework gives the opportunity to quantify and evaluate the statistical significance of trends in rare return levels of interest for hydrological design and risk management. It could be argued that estimating rare return levels on a yearly basis should lead to unberable uncertainties. We showed here that the statistical significance of trends in yearly-modelled return levels as high as the 50 yr events is comparable to the statistical significance of trends in AM, suggesting a similar signal to noise ratio. Trends on extreme return levels estimated on yearly basis from our model are thus characterized by stochastic uncertainties comparable to the ones of AM.