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.