Stationary and Non-Stationary Frameworks for Extreme Rainfall Time Series in Southern Italy

Abstract: This study tests stationary and non-stationary approaches for modelling data series of hydro-meteorological variables. Specifically, the authors considered annual maximum rainfall accumulations observed in the Calabria region (southern Italy), and attention was focused on time series characterized by heavy rainfall events which occurred from 1 January 2000 in the study area. This choice is justified by the need to check if the recent rainfall events in the new century can be considered as very different or not from the events occurred in the past. In detail, the whole data set of each considered time series (characterized by a sample size N > 40 data) was analyzed, in order to compare recent and past rainfall accumulations, which occurred in a specific site. All the proposed models were based on the Two-Component Extreme Value (TCEV) probability distribution, which is frequently applied for annual maximum time series in Calabria. The authors discussed the possible sources of uncertainty related to each framework and remarked on the crucial role played by ergodicity. In fact, if the process is assumed to be non-stationary, then ergodicity cannot hold, and thus possible trends should be derived from external sources, different from the time series of interest: in this work, Regional Climate Models’ (RCMs) outputs were considered in order to assess possible trends of TCEV parameters. From the obtained results, it does not seem essential to adopt non-stationary models, as significant trends do not appear from the observed data, due to a relevant number of heavy events which also occurred in the central part of the last century.


Introduction
Modelling of hydro-meteorological variables is an important topic in the hydrology field, and it plays a crucial role in estimation of design values, referred to assigned return periods [1].In this context, the statistical approaches, which are usually adopted for this modelling, can be regrouped into two main classes: stationary and non-stationary models [2].Both classes constitute a wide ensemble of choices for a user, as they allow for taking into account uncertainty sources and/or possible trends, which may characterize the specific time series under investigation [2].
In all the models belonging to the stationary class, the values of a random variable X are assumed as independent and identically distributed (iid) with a probability distribution F X (x|Θ), where Θ is the parameters vector, invariant in time.Concerning the calculation of F X (x|Θ), the main source of uncertainty is constituted by extrapolation of the probability law beyond the range of observed values for X [3,4].
As regards the non-stationary class, the values of a random variable are independent but not identically distributed (i/nid), and the probability law F X (x|Θ(t)) is characterized by a parameters vector which is not constant, but it is a function of some covariates, usually only the time t [5][6][7].
Consequently, in this case, a further source of uncertainty is introduced, related to the assumed function linking parameters and the covariate t [8].In fact, evaluation of design values for X also implies that the function Θ(t) remains valid for the entire future design period.Moreover, another important aspect is the ergodicity [2,9], that is the possibility to estimate the properties of a stochastic process from only one realization (the observed time series); if the process is non-stationary, then ergodicity cannot hold.Thus, possible trends should be derived from external sources, different from the time series of interest [2].
In any case, data integration (i.e., using information from other sources) is necessary [10].In fact, the sample size N is usually insufficient in order to obtain a robust estimation of quantities related to high values of F X (.) [11].This insufficiency is even more penalizing for non-stationary approaches, where it is necessary to extrapolate not only towards values of F X (.) greater than the observed frequencies, but also in time, along the future design period, on the basis of the assumed parametric trends [12].
Possible external sources of information can be the data in neighbor sites (by using regionalization techniques), or future scenarios derived from climate models [8].However, this latter type of information should be used with great care, as it does not derive from observed series, but from outputs of models [2].
In this general context, the estimation of hydrological extremes, in future temporal horizons, presents high levels of uncertainty [3,4,13].
In the observed time series, the presence of some outliers can constitute an indicator of a non-stationary model, but these outliers sometimes occurred not only in the more recent years, thus making weaker the hypothesis of a presence of trends for some sites [14].
As an alternative to the hypothesis of a non-stationary model, the presence of outliers could be modelled by considering the imperfect knowledge of parameters (i.e., parameters are assumed as random variables with assigned probability distributions, [15]) and preserving the hypothesis of a stationary process.In this way, only the observed data of the specific site are used and, in a parsimonious framework, no further source of uncertainty is introduced if a parametric trend is not considered.
From all these briefly discussed aspects, the importance to carry out an investigation/strategy emerges, like that illustrated in this work.The procedure should firstly assess the applicability of models belonging to both classes for the case study of interest.Then, according to the principle of parsimony and the satisfaction of hypothesis tests, the second step is the choice of the "operational" model (which could not necessarily be the more complex one) for successive practical applications (as example, estimations of quantiles) [16].
In this work, the authors investigated both hypotheses (stationary and non-stationary) for some rainfall time series related to gauges located in the Calabria region (southern Italy).The authors selected the Calabria region as a study area, because from a previous work regarding this zone [17] in which only continuous daily data series of the period 1916-2006 were analyzed, no significant trends in number and magnitude of extreme events emerged, unlike other neighboring parts (Central Alps and Northern Italy [18,19]).This specific behavior, found for Calabria, was also detected for other parts in the Mediterranean area (for example, the Iberian Peninsula [20]), while a general decrease of annual and monthly precipitation can be asserted for the whole Mediterranean area [21].In this work, the authors extended the investigation to annual maximum time series with sub-daily durations (1-24 h), and they carried out an analysis which is different from the application of non-parametric tests like Mann-Kendall, adopted in [17].The adopted methodologies, for stationary and non-stationary frameworks, were based on the Two-Component Extreme Value (TCEV, [22]) probability distribution, which is the most applied for annual maximum data series in Calabria.

Study Area
The analyzed time series (detailed in Section 2.2) are related to the rain gauge network located in the Calabria region (southern Italy), characterized by an area of about 15,000 km 2 and a perimeter of about 800 km.It is surrounded by the Ionian (eastern coast) and Tyrrhenian (western coast) seas.Calabria presents a typical summer subtropical climate [17], also known as the Mediterranean climate.In detail, the Ionian part is mainly affected by warm currents from Africa, which induce high temperatures and heavy rainfall events with short duration; western air currents usually affect the Tyrrhenian part, where milder temperatures and many orographic rainfall events mainly occur; the inland zone is characterized by colder winters with snow and fresher summers with some rainfall events [23].
The mean annual precipitation (MAP) over the whole region is about 1000 mm and more than 80% of this amount is recorded in the season from October to April [24].MAP varies from west to east side and between mountains and lowlands.Along the west coast, it averages about 900 mm and decreases to 700 mm along the northern part of the Ionian coast, and it is usually greater than 1100 mm above elevations of 500-800 m (but also below 500 m a.s.l. for some stations located in the southern part of the region), while it is less than 900 mm for valleys and for coastal areas [24].
In general, the western part of Calabria presents greater values of MAP, while the most extreme events with short durations occur in the eastern part, as it is exposed to more intense (but less frequent) cyclones from North Africa and the Balkans [24].

Case Studies
The authors focused attention on annual maximum rainfall series.In particular, daily and sub-daily (1-24 h) durations were investigated, and time series with recent heavy rainfall events, which occurred from 1 January 2000, were analyzed.In order to make reliable comparisons with heavy rainfall of the last century, the authors only considered time series with a sample size N > 40 data (with N clearly related to the whole data set for specific stations, and not only to the data from 2000 onwards).The annual maximum data series were downloaded from the website of the Multi Risks Center of Calabria Region (www.cfd.calabria.it).
The list of the examined rain gauge time series is detailed in Table 1 (see also Figure 1 for the geographical location), in which the mean annual precipitation (MAP), the elevation, the sample size N, the analyzed durations, and the intensity and occurrence date of heavy rainfall events before and after 1 January 2000 are indicated for each investigated time series.Fourteen time series from 11 rain gauges were analyzed (for 3 stations, the authors investigated time series related to two different durations); 10 rain gauges out of 11 are located in the southern part of Calabria, from which 5 are in the Ionian zone, 4 in the central part, and 1 in the Tyrrhenian zone.
It can be observed that all the heavy rainfall accumulations, selected from 1 January 2000 for the 11 investigated rain gauges, belong to nine spatial storm events that occurred in Calabria Concerning rainfall intensity, 130.2 mm in 1 h registered in Vibo Valentia (ID code 2800) on 3 July 2006 represents the highest value among the recent events.Focusing on percentages of MAP, 430.2 mm in 24 h occurred in Ardore Superiore (ID code 2210), constitutes the maximum observed rate (about 46% of MAP = 928.7 mm), and it is also the highest value registered in 24 h from 1 January 2000.Moreover, within the analyzed time series (specifically, in 9 out of 14) there are rain data, observed in the last century, with a magnitude that is comparable with the recent heavy rainfall values, and they belong to 14 spatial storm events that occurred in the region (10 November 1932, 5 December 1933, 21 October 1935, 21 November 1935, 2 December 1938, 9 September 1939, 22 January 1946, 16-21 October 1951, 3 November 1953, 12 September 1958, 23 November 1959, 22 September 1965, 25 November 1993, 3 October 1996).Eleven events out of 14 occurred in SON season, while 2 events occurred in DJF season.Furthermore, for three time series, the maximum value observed in the last century is even higher than that recorded from 2000 to date: 325.1 mm in 1 day in Vibo Valentia (2 December 1938), 512 mm in 24 h in Chiaravalle Centrale (16 October 1951), and 637.9 mm in 24 h in Santa Cristina d'Aspromonte (16 October 1951).
From analysis of Table 1, it is clear that the almost all the heavy rainfall events (both recent and past) occurred during the SON season, in which there is mainly the development of large-scale storms coming from the lee of the Alps or the western Mediterranean for the Tyrrhenian zone, and from North Africa or the Balkans for the Ionian zone [24].Λ are reported in Table 2 [34].
. In this case, the authors assumed ( ) and the first level of regionalization for parameter estimation was used (as only time series with N > 40 data were analyzed).

Theoretical Background of the Adopted Frameworks
A Cumulative Density Function (CDF) of a random variable X can be usually indicated as F X (x) or F X (x|Θ) (or F X (x|Θ(t)) when parametric trends are introduced), in an indistinct way.
However, a user should pay attention when Θ is not assumed as deterministic but as a random vector, with an assigned joint density distribution g(Θ).In this case, for a stationary approach, the following relationship is obtained by applying the total probability law: and, in a similar way, for the non-stationary case: Concerning Equation (2), the authors only considered a deterministic parametric trend Θ(t) in this work, in order to obtain a truly non-stationary model.In fact, "fluctuations of the parameters should generate stationary compound or mixed distributions" [2,[25][26][27].Consequently, for a deterministic Θ(t) the joint probability g(.) corresponds to the Dirac's function δ(.), and Equation (2) can be simplified as: Obviously, a similar simplification can be obtained for Equation (1), when Θ is considered as deterministic and then g(.) is also in this case equal to δ(.), that means F X (x) = F X (x|Θ).
Overall, starting from the general Equations ( 1) and (3), it is possible to obtain several models, depending on the assumed expression for the probability distribution F X (x|Θ), on the deterministic or stochastic nature of parameters, and on the eventual presence of the parameters' trend (which implies other assumptions, as explained in Section 2.3.3.).

The TCEV Distribution
In this work, the Two-Component Extreme Value (TCEV, [22]) distribution was adopted, which is very suitable for a regionalization approach and considers two independent processes, both having a Poisson rate of occurrence and an exponential distribution of magnitudes above a threshold.From these assumptions, it is possible to demonstrate that the CDF of the annual maximum X, deriving from these two processes, is [22]: where Θ = (Λ 1 , θ 1 , Λ 2 , θ 2 ).Label 1 is referred to the frequent events, while label 2 is associated to the outlier events.Λ i (i = 1, 2) represents the mean annual number of events (with clearly Λ 1 > Λ 2 ) and θ i (i = 1, 2) is the mean value of magnitudes (with clearly θ 1 < θ 2 ).By fixing can be rewritten as a product of two EV1 CDFs: with Θ = (ε 1 , θ 1 , ε 2 , θ 2 ).
As TCEV distribution is characterized by four parameters, their estimation from an at-site data series is robust enough if the sample has an adequate size N (at least 90-100 data).On the contrary, if N is less than 90 data, a regionalization approach is necessary, that means TCEV distribution is hierachically developed into three regional levels, explained below.
In the 2nd level (10 < N < 40 data), homogeneous sub-regions are also considered, within each region identified in the previous level.In each sub-region, besides the constancy of Λ * , θ * (and then γ 1 ) of the related region, the parameter Λ 1 is also assumed theoretically invariant from site to site, and thus CV is constant too.In this case, only one parameter (θ 1 or µ) is at-site estimated.
In the 3rd level (N < 10 data), θ 1 (or µ) is also estimated in a regional way, by adopting defined empirical formulas, depending (in the case of rainfall) on the altitude of the point or other topographic covariates [31].
For the selected case study, concerning annual maxima of daily and sub-daily rainfall accumulation, from application of TCEV-ML, the whole Calabria region is a unique homogenous region, and it is subdivided into three homogeneous sub-regions, named Tyrrhenian (T), Central (C), and Ionian (I), respectively, as shown in Figure 1.The regional values of parameters Λ * , θ * and Λ 1 are reported in Table 2 [34].
Table 2. Two-Component Extreme Value (TCEV) distribution of rainfall annual maxima with daily and sub-daily durations: regional values of Λ * , θ * , and Λ 1 for the Calabria region [34].

Stationary Framework
For all the time series listed in Section 2.2, the authors considered two models belonging to the stationary framework.
The first is indicated as Model 1.It is the simplest version of Equation ( 1), in which g(.) = δ(.), that means Θ is deterministic and then F X (x) = F X (x|Θ).In this case, the authors assumed Θ = (Λ * , θ * , Λ 1 , θ 1 ) and the first level of regionalization for parameter estimation was used (as only time series with N > 40 data were analyzed).
The second is denoted as Model 2, and it represents the complete version of Equation ( 1), in which g(.) is computed as described below.
As regards Λ * and θ * , let a standardized random variable Y be defined as: If X is a TCEV random variable, then Y is TCEV distributed too, with the following expression [2]: For each considered value of N (40, 50, . . ., 500), 1000 synthetic samples of Y were generated with a Monte Carlo approach by using Equation ( 13) and the regional values of Λ * and θ * , as reported in Table 2.For each synthetic sample, estimations of Λ * and θ * , indicated as Λ * and θ * , were carried out by applying a ML procedure [22,32], and then empirical distributions for Λ * and θ * were obtained, depending on N, and assumed as valid for the whole homogenous region.
Regardless of N, the normal and log-normal distributions showed the best fitting for Λ * and θ * , respectively, as shown in Figures 2 and 3, and the median values were, in any case, equal to the assumed theoretical values Λ * and ln θ * , respectively (for the latter correspondence, please refer to [35], p. 216).
Concerning the standard deviations σ− , a clear and expected dependence on N emerged, as shown in Figure 4.
Analyzing the dispersion plot in Figure 5, obtained with all the values Λ * and θ * estimated from all the synthetic series, Λ * and θ * were assumed as independent random variables.
As regards g Λ 1 (Λ 1 ), it was theoretically derived ( [35], p. 237): as Λ 1 corresponds to the mean value of a distribution (Poisson in this case), its distribution is close to normal if N > 30 data even when the original variable is not normally distributed.Further theoretical results are: where √ Λ 1 is the theoretical standard deviation of the original (Poisson) distribution.Consequently, for a time series with a sample size N, the authors considered the sub-region where the rain gauge is placed, and then they used the associated regional value Λ 1 , as shown in Table 2, in order to compute the normal distribution g Λ 1 (Λ 1 ) with parameters expressed by Equations ( 14) and (15).Moreover, Λ 1 was also assumed as independent on both Λ * and θ * , and thus g(Λ * , θ * , Λ 1 ) was set equal to the product , where g Λ * (Λ * ) and g θ * (θ * ) are the normal and log-normal, respectively, distributions explained above and assumed as valid for the whole Calabria region.Water 2018, 10, x FOR PEER REVIEW 9 of 27       Moreover, 1 Λ was also assumed as independent on both * Λ and * θ , and thus ( ) , where

Non-Stationary Framework
The whole procedure, described below, is indicated by the authors as Model 3.Because of violation of ergodicity when a process is non-stationary [2], the authors used information external from the data series of rain gauges for non-stationary analysis.The report of the European Environmental Agency (EEA) [36] and the publication of the Italian Institute for Environmental Protection and Research (ISPRA) [37] were examined, concerning analysis of future climate in Europe and Italy, respectively, by using projections from Regional Climate Models (RCMs).Four RCMs were adopted (named as ALADIN, GUF, CMCC, LMD) and, for each one, four scenarios of Representative Concentration Pathways (RCPs) of total radiative forcing (i.e., cumulative measure of human emissions of greenhouse gasses from all sources expressed in W/m 2 ) were used as input: RCP 2.6, RCP 4.5, RCP 6, and RCP 8.5.Moreover, an ensemble mean projection from all the RCMs is also derived for each RCP.
Focusing the attention on RCP 4.5 (intermediate emissions) and RCP 8.5 (high emissions), with respect to the reference period 1971-2000, the results of the simulations, related to a future period until 2090 for the Calabria region, can be summarized as indicated below.
A modest increase is predicted for the annual maximum daily rainfall.The ensemble mean is, for both RCP 4.5 and RCP 8.5, up to 10 mm, with a maximum increase of 35-40 mm from the CMCC model.
On the contrary, a significant increase for the waiting time between two consecutive rainfall events is forecasted.The ensemble mean is: up to 15 days for RCP 4.5, with a maximum increase of 30 days from the CMCC model; up to 30 days for RCP 8.5, with a maximum increase of 40 days from the CMCC model.
It has to be highlighted that these results are, obviously, affected by uncertainties, derived from the structure of the adopted numerical schemes of RCMs, and from their spatial and temporal resolutions.In fact, RCMs generally underestimate extreme precipitation intensity, and they are better at locating an extreme rainfall event than estimating its intensity.Model accuracy clearly improves with resolution, but RCMs with spatial resolution between 10 and 30 km, typically used in climate change studies, are still too coarse to explicitly represent sub-daily localized heavy precipitation events [38,39].A very high-resolution model (typically 1-5 km), used for weather forecasts with explicit convection, has recently used for a climate change experiment for a region in the UK: this study projects an intensification of short duration heavy rain in summer, with significantly more events exceeding the high thresholds indicative for serious flash flooding [38,40,41].
However, on the basis on these results, the non-stationary framework was organized by the authors by considering the features associated to the outliers (i.e., the parameters with label 2 in Equation ( 4)).In detail, the authors established a positive trend for θ 2 and a negative trend for Λ 2 .
Clearly, the specific mathematical expression for the parameters' trend constitutes a not negligible source of uncertainty in a non-stationary framework.In this context, several expressions could be assumed [2]: linear, quadratic, regime shifts, sequence of different forms, etc.
In this work, the authors assumed the following expressions: where t 0 is the year when the trend is supposed to begin, α θ 2 and α Λ 2 are the rates of increase/decrease in one year for θ 2 and Λ 2 , respectively, while θ st 2 and Λ st 2 are the steady values of parameters before the beginning of the trend.
The authors set t 0 equal to 1980 (i.e., a year representative of the period when the scientific community began to discuss climate changes and the first conferences were organized) and to 2000 (i.e., the last year of the reference period for RCMs).
As regards α θ 2 and α Λ 2 , their values were fixed in order to reach an increase/decrease in 100 years with the same order of magnitude of RCMs scenarios forecasted for the Calabria region; consequently, the authors set α θ 2 = 0.005 (i.e., an increase of 50% in 100 years) and α Λ 2 = −0.002(i.e., a decrease of 20% in 100 years).Finally, θ st 2 and Λ st 2 were assumed as the estimations obtained by considering the time series as stationary until t 0 .
Overall, two versions were considered, depending on the assumed value for t 0 : Model 3_1980 and Model 3_2000.In detail, for a time series with N > 40, all the data which occurred when t < t 0 were used to estimate the stationary vector Θ = (Λ * , θ * , Λ 1 , θ 1 ), by using a suitable level of regionalization, and the authors calculated θ 2 = θ st 2 and Λ 2 = Λ st 2 with Equations ( 7) and (8).Then, for each year t ≥ t 0 , θ 2 (t) and Λ 2 (t) holds from Equations ( 16) and ( 17), and the parameters vector for TCEV is For large values of sample size N, the stationary parameter vector of Model 3 is usually very close to be identical to that associated to Model 1, when a user adopts a maximum likelihood (ML) procedure for parameter estimation (like in this work); in fact, ML application provides parameter values which are not so influenced by the presence of outliers.

Procedure for Models Application
For each time series, the authors tested all the previously described models.As regards to Model 1, starting from Equation ( 9) with Θ = (Λ * , θ * , Λ 1 , θ 1 ) estimated at first level of regionalization, 1000 synthetic samples (with the same N of the observed series) were generated by using the Monte Carlo method.Then, for each synthetic sample, all the generated data were sorted in ascending order and assigned to specific Plotting Position (PP) values, calculated with Weibull formula [42].It is clear that synthetic samples with the same N present the same set of PP values.Moreover, for each PP value, the 95% Confidence Intervals (CIs) of rainfall accumulation was derived and then the authors checked if the whole observed sample fell inside 95% CIs on an EV1 probabilistic plot; if the check is positive, the hypothesis of a stationary TCEV with Θ as deterministic cannot be rejected.
Concerning Model 2, the authors considered Equation (10), with µ equal to the sample mean, while Λ * , θ * , Λ 1 as random variables with probability distributions which are described in Section 2.3.2 (normal for Λ * , Λ 1 and log-normal for θ * ) and clearly depend on the value of sample size N.One thousand synthetic triads (Λ * , θ * , Λ 1 ) were derived with the Monte Carlo method, and each triad was used (togheter with µ) in Equation (10) to generate a specific sample with the same N of the observed series.Also in this case (with the same procedure of the previous Model 1), the observed data were compared with the estimated 95% CIs on an EV1 probabilistic plot, in order to evaluate the possibility of rejection or not for Model 2.
Referring to Model 3, as mentioned in Section 2.3.3, the authors assumed Equations ( 16) and ( 17) for trends regarding (Λ 2 , θ 2 ).Also in this case, the authors carried out a generation of 1000 synthetic samples (with the same N of the observed series) by using the Monte Carlo method, but considering (unlike Model 1) parameters vector as stationary until t 0 , and variable for t ≥ t 0 .The check of rejection or not for Model 3 comprised the following two conditions: (i) all the occurred events when t ≥ t 0 are inside the 95% CIs of Model 3 (which are estimated in a similar way to those of Models 1 and 2); (ii) all the occurred events when t < t 0 are inside the 95% CIs of the stationary component of Model 3 (which are very close to those of Model 1, for the reasons detailed at the end of Section 2.3.3).Only if both conditions are satisfied, then the validity of Model 3 cannot be rejected.

Results and Discussion
For all the time series, the obtained results on EV1 probabilistic plots are reported in Appendix A. In this section, as examples, the outcomes related to Vibo Valentia (1-h and daily time series) and Ardore Superiore (1-h and 24-h time series) are shown in Figures 6-9.
Water 2018, 10, x FOR PEER REVIEW 13 of 27 component of Model 3 (which are very close to those of Model 1, for the reasons detailed at the end of Section 2.3.3).Only if both conditions are satisfied, then the validity of Model 3 cannot be rejected.

Results and Discussion
For all the time series, the obtained results on EV1 probabilistic plots are reported in Appendix A. In this section, as examples, the outcomes related to Vibo Valentia (1-h and daily time series) and Ardore Superiore (1-h and 24-h time series) are shown in Figures 6-9.Concerning the Vibo Valentia rain gauge, Model 1 is not able to reproduce, for both 1-h and daily durations, the outliers are inside the 95% CIs.On the contrary, taking into account parameter uncertainty (Model 2) allows to do not reject the assumption of stationary process for both time series.
The comparison between the theoretical curve of Model 1, from which 1000 synthetic Concerning the Vibo Valentia rain gauge, Model 1 is not able to reproduce, for both 1-h and daily durations, the outliers are inside the 95% CIs.On the contrary, taking into account parameter uncertainty (Model 2) allows to do not reject the assumption of stationary process for both time series.
The comparison between the theoretical curve F X (x|Θ) of Model 1, from which 1000 synthetic samples were generated, and the theoretical curve F X (x), obtained from Equation (1) for Model 2 shows, for each duration, the similarity of these two parent distributions and, consequently, remarks on the crucial role played by taking into account g Λ * (Λ * ), g θ * (θ * ) and g Λ 1 (Λ 1 ) if a stationary framework is assumed.
As regards the non-stationary framework: for 1-h time series, only Model 3_1980 is able to reconstruct sample data inside the 95% CIs as previously explained (i.e., satisfaction of both conditions), while both models have to be rejected for 1-day time series.
Water 2018, 10, x FOR PEER REVIEW 14 of 27 As regards the non-stationary framework: for 1-h time series, only Model 3_1980 is able to reconstruct sample data inside the 95% CIs as previously explained (i.e., satisfaction of both conditions), while both models have to be rejected for 1-day time series.For 1-h time series of Ardore Superiore, Model 1 is sufficient for a suitable reconstruction of sample data, and application of the other models did not significantly improve the performances.On the contrary, concerning 24-h series recorded in the Ardore rain gauge, rejection is necessary for Model 1 and for all the versions of Model 3, while Model 2 allowed for the reproduction of the whole data sample inside the assigned CIs.Focusing on the theoretical curves    For 1-h time series of Ardore Superiore, Model 1 is sufficient for a suitable reconstruction of sample data, and application of the other models did not significantly improve the performances.On the contrary, concerning 24-h series recorded in the Ardore rain gauge, rejection is necessary for Model 1 and for all the versions of Model 3, while Model 2 allowed for the reproduction of the whole data sample inside the assigned CIs.Focusing on the theoretical curves F X (x|Θ) of Model 1 and F X (x) of Model 2, the same considerations for Vibo Valentia time series hold.Overall, the histogram in Figure 10 summarizes the performances of all the adopted models in terms of No-Rejection Rate (NRR), defined as the ratio between the number of time series for which the hypothesis of model validity cannot be rejected, and the total number of analyzed time series.In many cases, it is sufficient to use Model 1 (NRR = 0.76, i.e., 11 out of 14 time series): this is mainly true if a rain gauge recorded a significant number of outlier events also in the last century.NRR (-) Overall, the histogram in Figure 10 summarizes the performances of all the adopted models in terms of No-Rejection Rate (NRR), defined as the ratio between the number of time series for which the hypothesis of model validity cannot be rejected, and the total number of analyzed time series.In many cases, it is sufficient to use Model 1 (NRR = 0.76, i.e., 11 out of 14 time series): this is mainly true if a rain gauge recorded a significant number of outlier events also in the last century.Overall, the histogram in Figure 10 summarizes the performances of all the adopted models in terms of No-Rejection Rate (NRR), defined as the ratio between the number of time series for which the hypothesis of model validity cannot be rejected, and the total number of analyzed time series.In many cases, it is sufficient to use Model 1 (NRR = 0.76, i.e., 11 out of 14 time series): this is mainly true if a rain gauge recorded a significant number of outlier events also in the last century.NRR (-) Model 2 is able to reproduce all the time series inside the 95% CIs (NRR = 1), thus appearing as the preferable approach for a user, as it allows for only using information of rain gauge time series in the study area, without considering external sources.
with respect to a stationary model with deterministic parameters.These results are justified by two reasons: the stationary component of this framework does not correctly model some outlier events that occurred in the last century; the assumed period for trends could have an insufficient duration in order to reproduce the more recent events.
Consequently, if the hypothesis of a stationary process with parameters as random variables cannot be rejected, a user could surely prefer stationary frameworks.The obvious advantages are: the respect of ergodicity of the process (and then only the observed data series can be used); no other source of uncertainty is introduced, mainly related to the mathematical expression of the trends, the beginning t 0 , and the reliability of the external information (in particular, outputs from RCMs, characterized by a spatial resolution not compatible with rainfall phenomena at point-scale).
Overall, the obtained results in this work showed that, for annual maximum rainfall series in Calabria, it is not essential to remove the hypothesis of stationary parameters: significant trends do not appear from the observed data, as a relevant rate of heavy events also occurred in the central part of the last century.However, further analyses are needed, based on broader samples than annual maxima (for example POT-Peak Over Threshold series), in order to establish the more general validity of the results presented here.Obviously, these results do not imply any statement about climate changes (as discussed in Section 3), but they only intend to offer operating elements for an estimate of hydrological design extremes referred to a very close time window in the future.
The methodology, here applied, further suggests that a good modelling strategy should firstly consider the simplest approach, according to the principle of parsimony.In fact, using directly a more complex model instead of a simpler one could not be justified, without clear evidence from observed data, or understanding of physical dynamics.
Therefore, the stationary framework should represent the default assumption [49] and, even if there is evidence for non-stationarity, stationary models should be used as a benchmark for more complex approaches, in order to evaluate the magnitude of the improvements.
Lastly, uncertainty plays a crucial role in this context, related to model parameters (should a user assume them as deterministic or as random variables?), and to the other previously mentioned aspects referred to the use of a non-stationary model.Thus, adoption of few sources of uncertainty, or at most one (associated to the parameters assumed as random variables in a stationary framework), appears more competitive when changes of some features are not so clear from the data.(d) (c)    (d) (c)           (d) (c)    (d) (c)

Water 2018 ,
10, x FOR PEER REVIEW 7 of 27 and Ionian (I), respectively, as shown in Figure 1.The regional values of parameters * Λ , * θ and1

Figure 1 .
Figure 1.Study area of the Calabria region (southern Italy): rain gauge network (blue points), rain gauges considered in this work (red points with ID code labelled), sub-regions Tyrrhenian (green polygons), Central (yellow polygon), and Ionian (pink polygons).

Figure 1 .
Figure 1.Study area of the Calabria region (southern Italy): rain gauge network (blue points), rain gauges considered in this work (red points with ID code labelled), sub-regions Tyrrhenian (green polygons), Central (yellow polygon), and Ionian (pink polygons).

Figure 6 .
Figure 6.Vibo Valentia rain gauge-annual maximum 1-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

Figure 6 .
Figure 6.Vibo Valentia rain gauge-annual maximum 1-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

Figure 7 .
Figure 7. Vibo Valentia rain gauge-annual maximum 1-day time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

of
Model 2, the same considerations for Vibo Valentia time series hold.

Figure 7 .
Figure 7. Vibo Valentia rain gauge-annual maximum 1-day time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

Figure 8 .
Figure 8. Ardore Superiore rain gauge-annual maximum 1-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

Figure 8 .
Figure 8. Ardore Superiore rain gauge-annual maximum 1-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

Figure 9 .
Figure 9. Ardore Superiore rain gauge-annual maximum 24-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

Figure 10 .
Figure 10.Histogram of No-Rejection Rate (NRR) for each adopted model.

Figure 9 .
Figure 9. Ardore Superiore rain gauge-annual maximum 24-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

Water 2018 , 27 Figure 9 .
Figure 9. Ardore Superiore rain gauge-annual maximum 24-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

Figure 10 .
Figure 10.Histogram of No-Rejection Rate (NRR) for each adopted model.

Figure 10 .
Figure 10.Histogram of No-Rejection Rate (NRR) for each adopted model.

Figure A1 .
Figure A1.Chiaravalle Centrale rain gauge-annual maximum 1-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

FXFigure A1 .
Figure A1.Chiaravalle Centrale rain gauge-annual maximum 1-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

Figure A1 .
Figure A1.Chiaravalle Centrale rain gauge-annual maximum 1-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

Figure A2 .
Figure A2.Chiaravalle Centrale rain gauge-annual maximum 24-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

Figure A3 .
Figure A3.Santa Cristina d'Aspromonte rain gauge-annual maximum 24-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

FFigure A2 .
Figure A2.Chiaravalle Centrale rain gauge-annual maximum 24-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

Figure A2 .
Figure A2.Chiaravalle Centrale rain gauge-annual maximum 24-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

Figure A3 .
Figure A3.Santa Cristina d'Aspromonte rain gauge-annual maximum 24-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

FFigure A3 .
Figure A3.Santa Cristina d'Aspromonte rain gauge-annual maximum 24-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

Figure A4 .
Figure A4.Montalto Uffugo rain gauge-annual maximum 24-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

FFigure A4 .
Figure A4.Montalto Uffugo rain gauge-annual maximum 24-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

Figure A4 .
Figure A4.Montalto Uffugo rain gauge-annual maximum 24-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

Figure A5 .
Figure A5.Sant'Agata del Bianco rain gauge-annual maximum 1-day time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

Figure A6 .
Figure A6.Gioiosa Ionica rain gauge-annual maximum 3-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

FFigure A5 .
Figure A5.Sant'Agata del Bianco rain gauge-annual maximum 1-day time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

Figure A5 .
Figure A5.Sant'Agata del Bianco rain gauge-annual maximum 1-day time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

Figure A6 .
Figure A6.Gioiosa Ionica rain gauge-annual maximum 3-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

FFigure A6 .
Figure A6.Gioiosa Ionica rain gauge-annual maximum 3-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

Figure A7 .
Figure A7.Platì rain gauge-annual maximum 24-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

FFigure A7 .
Figure A7.Platì rain gauge-annual maximum 24-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

Figure A7 .
Figure A7.Platì rain gauge-annual maximum 24-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

Figure A8 .
Figure A8.Serra San Bruno rain gauge-annual maximum 1-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

Figure A9 .
Figure A9.Albi rain gauge-annual maximum 3-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

Figure A8 .
Figure A8.Serra San Bruno rain gauge-annual maximum 1-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

Figure A8 .
Figure A8.Serra San Bruno rain gauge-annual maximum 1-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

Figure A9 .
Figure A9.Albi rain gauge-annual maximum 3-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

Figure A9 .
Figure A9.Albi rain gauge-annual maximum 3-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

Figure A10 .
Figure A10.Cittanova rain gauge-nnual maximum 6-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

Figure A10 .
Figure A10.Cittanova rain gauge-nnual maximum 6-h time series: application of (a) Model 1, (b) Model 2, (c) Model 3_1980 (the CIs of the stationary component are represented by dotted lines), (d) Model 3_2000 (the CIs of the stationary component are represented by dotted lines).

Table 1 .
List and characteristics of analyzed time series in this work.