Detecting precipitation trend using a multiscale approach based on quantile regression over a Mediterranean area

One of the most relevant and debated topics related to the effects of the climate change is whether intense rainfall events have become more frequent over the last decades. It is a crucial aspect, since an increase in the magnitude and frequency of occurrence of heavy rainfall events could result in a dramatic growth of floods and, in turn, human lives losses and economic damages. Because of its central position in the Mediterranean area, Sicily has been often screened with the aim to capture some trends in precipitation, potentially related to climate change. While Mann‐Kendall test has been largely used for the rainfall trend detection, in this work a different procedure is considered. Precipitation trends are here investigated by processing the whole rainfall time‐series, provided by the regional agency SIAS at a 10‐min resolution, through the quantile regression method by aggregating precipitation across a wide spectrum of durations and considering different quantiles. Results show that many rain gauges are characterized by an increasing trend in sub‐hourly precipitation intensity, especially at the highest quantiles, thus suggesting that, from 2002 to 2019, sub‐hourly events have become more intense in most of the island. Moreover, by analysing some spatial patterns, it has been revealed that the south and the east of Sicily are more interested in significant increasing rainfall trends, especially at the 10‐min duration. Finally, the comparison between the two procedures revealed a stronger reliability of the quantile regression in the trend analysis detection, mainly due to the possibility of investigating the temporal variation of the tails of precipitation distribution.


| INTRODUCTION
For about 30 years climate change has been, and still today is, one of the most relevant and debated topics for the scientific community. Although the signals, impacts, and effects of climate change appear to be obvious and we have been hearing more and more about matters such as global warming, rising of oceans and seas temperature, rising of global sea level, decreasing of arctic sea ice extent, melting of glaciers, increasing of extreme events (e.g., heavy rainfall events and/or droughts), and so forth, still today a part of the scientific community is sceptical about climate change (Koutsoyiannis, 2020). What is certain is that if no action will be taken to face the climate change, some climate tipping points might be reached (Lenton, 2011), resulting in strong consequences for the Earth difficult, or even impossible, to revert (e.g., the ice melt in the poles). As an example, an increasingly warm climate could result in an increase of the occurrence of longer and more intense droughts in the 21st century (Field et al., 2018) with relevant implications in water availability for human, vegetation, and agriculture (Thuiller et al., 2011). At the other side of the spectrum, an increasing trend in the magnitude and frequency of heavy rainfalls could lead to floods and, consequentially, to fatalities and economic damages (Tabari, 2020). All these aspects, in the long-term, could produce social and economic frictions between populations, leading to a growth of future risk of conflicts (Mach et al., 2019).
In this respect, already in the past, the Mediterranean region has been referenced as one of the most responsive regions to climate change, so much to be defined as a primary hotspot of climate change (Giorgi, 2006). Also, the last report from the Intergovernmental Panel on Climate Change has highlighted the Mediterranean as one of the most vulnerable regions in the world to the impacts of global warming (IPCC, 2019). One of the most debated points, still today, is whether these changes have led to an increase in frequency and magnitude of heavy rainfall events over the Mediterranean area in the last years. The detection of an increasing trend in such a kind of events is very important since, other than being an indicator of a climate alteration, it could support political decisions to mitigate their effects (Ingold and Fischer, 2014). For this reason, over the last decades, a great number of studies have been aimed to detect the presence of trends in precipitation time-series over the Mediterranean region.
In general, to study the changes in rainfall characteristics, there is no doubt that it is necessary to perform a robust statistical analysis of historical data. Over the years, one of the most used methods to detect trends in precipitation has been the non-parametric Mann-Kendall (hereinafter referred to as MK) test. It has largely been used also in trend detection for extreme events, often coupled with different methods used for extracting extremes, such as the Block-Maxima (Pujol et al., 2007;Villarini, 2012;Westra et al., 2013;Wi et al., 2016), the Peak Over Threshold (POT) (Villarini et al., 2011b;Tramblay et al., 2012;Wi et al., 2016) and the indices developed by the Expert Team on Climate Change Detection and Indices (ETCCDI) (Song et al., 2015;Panda et al., 2016;Gentilucci et al., 2020). A great number of studies involving the MK test have been applied over the Mediterranean area (Norrant and Douguédroit, 2006;Pujol et al., 2007;Chaouche et al., 2010;Valdes-Abellan et al., 2017) and particularly to the Italian peninsula (Buffoni et al., 1999;Crisci et al., 2002), its southern part (Longobardi and Villani, 2010;Caloiero et al., 2011), and its main islands, that is, Sardinia and Sicily (Bonaccorso et al., 2005;Cannarozzo et al., 2006;Arnone et al., 2013;Caloiero et al., 2019), since their central geographic location in the Mediterranean basin. Focusing on Sicily (Italy), many studies have involved the use of the MK test for detecting precipitation trends at several durations and evaluating their significance. As an example, Bonaccorso et al. (2005) applied the MK test to Sicilian rainfall annual maxima at the canonical durations (i.e., 1, 3, 6, 12, and 24 hr) for those stations having more than 50 years of data. The authors found a relationship between the trend direction and duration, highlighting an increasing trend at the shortest time scale and an opposite behaviour at the longest ones. Arnone et al. (2013) used the MK test coupled with a non-parametric estimate, presented by Hirsch et al. (1982), to determinate the magnitude of trends and thus identify changes in rainfall characteristics at the canonical durations for Sicily. The authors found out that the percentage of gauges showing a positive trend tends to decrease when duration increases, finding the maximum percentage at the 1-hr duration. For this reason, they guessed that a further increase in the percentage of stations showing a positive trend and, consequently, of extreme events could be extrapolated towards sub-hourly durations. In addition, the number of stations characterized by a non-significant trend increases with duration and consists of the greatest part of the gauges under study.
Differently from the MK test, the quantile regression (hereinafter referred to as QR) analysis (Koenker and Bassett, 1978;Koenker, 2005) allows one to perform a linear regression on the whole data time-series, taking into account those values greater than a selected quantile. This means that, if very high (low) quantiles are considered, QR allows exploring the upper (lower) tail of the probability distribution function of the data. QR method has been applied in the past to study trends at different temporal and spatial scales for rainfall (Villarini et al., 2011a;Bartolini et al., 2014;Lausier and Jain, 2018) and other climatic variables, such as temperature (Barbosa et al., 2011). Lausier and Jain (2018) applied the QR method to the annual total precipitation at a global scale, comparing results with those provided by a linear regression model. The authors, who found different precipitation trends related to the mean (linear regression model) and the median (QR model at 0.5 quantile) of time-series, asserted that a wrong trend interpretation, deriving by using an easier method as simple linear regression, could have implications for some environmental systems. Hence, they suggested to use a more robust method, such as the QR. Lastly, by detecting trends for the lower and upper tails of precipitation probability density function, they classified the whole planet into three risk classes, with the aim to identify some strategies to deal with them. Bartolini et al. (2014) started from two hourly rainfall data sets to look for trends in precipitation amount, frequency, and intensity in two locations of Tuscany (Italy) by means of a QR method. Results showed a tendency to a decrease of total rainfall and wet hours, occurring in winter and spring, and to an increase of hourly average precipitation during wet hours. With reference to rainfall at sub-daily timescale, some studies used the QR with the aim to understand changes of sub-daily precipitation with temperature. For instance, Van de Vyver et al. (2019) detected the scaling rates of sub-daily precipitation with dew point temperature at various quantiles, highlighting a general increase of rainfall extremes depicted by quantiles higher than 0.9, with dew point temperature in various cities of Europe. However, to our knowledge, there seems to be a lack of studies involving applications of QR in detecting trends of shortduration (i.e., hourly and sub-hourly) precipitation, which is one of the main purposes of the present study.
Starting from the above-mentioned conjecture of Arnone et al. (2013), this work analyses a 10-min resolution precipitation data set, different from that used by Arnone et al. (2013), in order to identify eventual statistically significant trends for Sicily at the sub-hourly durations and thus verify if the hourly and sub-hourly precipitation are characterized by the same behaviour, both in terms of trend direction and significance. To determinate an eventual existence of trends and study their magnitude, a different method than that used in Arnone et al. (2013), namely the QR method, is here used. The Student's t test is used to assess the significance of trends. Finally, the presence of any spatial patterns of trends in magnitude is globally and locally verified through a spatial autocorrelation analysis based on the Global Moran's I Index (Moran, 1950) and the Local Moran (Anselin, 1995), respectively.
Despite the observation period here investigated is not enough long to infer properly about climate change effects, the results of the study may still be a further potential signal that something is probably changing in Sicilian and, more in general, in Mediterranean climate, as several studied have already pointed out over the years (Giorgi, 2006;Giorgi and Lionello, 2007;Arnone et al., 2013;Forestieri et al., 2018). In this sense, the work attempts to contribute to the research in the field of climate change and its impact on heavy rainfall occurrence trying to answer to the following research questions: is it possible to identify statistically significant trends at the sub-hourly durations, thus verifying the conjecture of Arnone et al. (2013)? Is there an increase in the significance of rainfall trends with the duration? And if so, is it possible to assume that this aspect mainly concerns extreme rather than ordinary events? Finally, is it possible to identify any spatial pattern of trends magnitude?
The article is organized as follows. Section 2 introduces the methodology and data used in this study. Here just a brief description of the detection trends methods is provided. Section 3 is comprehensive of the achieved results and the related discussions. Lastly, Section 4 reports the conclusions concerning the study.

| METHODOLOGY AND DATA SET
In this paper, two different procedures are used to identify trends in rainfall time-series: the QR method, originally developed for econometric and statistical applications and subsequently applied also in the environmental field, associated with the Student's t test, and the MK test, largely used in extreme events analysis, coupled with the Sen's slope estimator. The two procedures are briefly described in the following.
2.1 | Quantile regression procedure QR method, as introduced by Koenker and Bassett (1978), can be considered as a natural extension of the standard linear regression models, due to the possibility to perform a regression on quantiles rather than just on the mean. The capability to investigate, at any quantile level, the linear relationship between two or more variables provides a more complete view of the statistical properties of a sample, also inspecting the tails of its distribution. Furthermore, standard regression models are strongly influenced by the presence of outliers, aspect that could be quite annoying especially in detecting trends along time.
The main difference between a simple linear regression and the QR method is on the evaluation of coefficients. While in a classical bi-dimensional simple linear regression the intercept and the slope of the regression line are evaluated through the least square minimization problem, the QR model is based on a minimization of the sum of the weighted absolute value of a difference between the ith observation (y i ) and the τth quantile line (β 0 (τ) + β 1 (τ)x i ) at x i : From Equation (1) it is possible to evince the dependency of the intercept, β 0 , and the slope, β 1 , of the regression line on the quantile level τ, for any 0 < τ < 1. The role of τ and (1 − τ) is to weight the vertical distances, that depend on the position of the observations with respect to the τth quantile line. In particular, points above the quantile line are weighted by τ, while those below the quantile line are weighted by (1 − τ), thus meaning that the greater the considered quantile, the more relevant are points above the quantile line in the evaluation of the slope and intercept of the regression line. To evaluate the trends significance, the Student's t test has been applied to the QR results. The test is here used to reject the null hypothesis, with a significance level of .05 and .1, that the slope of the quantile line is equal to zero. In order to have a measure of the accuracy in estimating the slope and the intercept of the QR line, the standard error is also computed with a sparsity method, known as 'nid' ('not independently and identically distributed error') (Koenker, 2004). This method, as considered by Koenker and Machado (1999), allows one to estimate the sparsity function for data that are not independently and identically distributed, assuming local linearity of the conditional quantile function Q(τjx) in x, where Q(Á) indicates the conditional probability of the quantile τ given the observed variable x. The 'nid' method is sensitive to the presence of many equal values in the analysed data set; indeed, such a condition could generate a singular matrix, from which the algorithm cannot compute the SE and, consequently, the confidence interval.
For further details on the QR method the reader is referred to Koenker (2005) and Hao and Naiman (2007).

| Mann-Kendall procedure
The MK test (Mann, 1945;Kendall, 1948) has been largely used in the history of hydrological trends, mainly because it is a non-parametric test. Due to this aspect, trends related to a generic variable can be obtained without any assumption about the properties of its distribution. Furthermore, this method is not so much affected by the presence of outliers.
Unlike the QR method, the MK test is usually not applied to raw precipitation data (i.e., time series of rainfall depth originally recorded by the gauges), but it is generally used with specific data sets, such as annual maxima (Bonaccorso et al., 2005;Arnone et al., 2013). The null hypothesis in the test indicates that the population from which the sample is extracted has no trend, while the alternative hypothesis is that a trend exists. To accept or decline the null hypothesis at a fixed significance level (i.e., α sig ) a comparison between α sig and a local significance level (i.e., p value ) is required. This last term is obtained as follows: where Φ(Á) is the CDF (Cumulative Distribution Function) of a standard normal variate. The standardized test statistic, Z S , follows a standard normal distribution and can be computed as reported below: In Equation (3), σ is the variance of the standardized normal distribution function followed by the Kendall's S statistic, under the null hypothesis. The S statistic is computed as the sign function of the difference between two consecutive observations, namely x i and x j : In the case of an autocorrelated series, the MK test could detect a trend even if it is not real, as demonstrated by Von Storch (1999). For this reason, a pre-whitening procedure is usually suggested, especially when the observed data set is shorter than 50 elements (Yue and Wang, 2002). It consists in removing from each observation, x i , a component given by the product of the previous one and the lag 1 serial correlation coefficient.
In literature, MK is quite often coupled with the Sen's slope method (Almeida et al., 2017;Güçlü, 2018), as an estimator of the trend magnitude. This latter assumes that the slope of the regression line is estimated as the median of the ensemble of slopes derived by linking the pairs of consecutive observed data (Sen, 1968).

| Spatial autocorrelation analysis
The Global Moran's I statistic (Moran, 1950) can be interpreted as an extension of the autocorrelation coefficient, by means of a symmetric spatial weights' matrix filled by the inverse of the geographic distance between the pairs of points, where data are recorded. The value of Moran's I ranges from −1 to 1 (same as the autocorrelation). The Global Moran's I has been here used to perform a spatial autocorrelation (hereafter SAC) analysis in order to verify the presence of specific patterns in the spatial distributions of the trend magnitude.
In order to establish whether data are randomly distributed or not, the Moran's I must be compared with its expected value, E[I]. In general, if I is less than E[I] data are dispersed, while I values greater than E[I] indicate a clustered pattern. When I is close to its expected value, instead, data tend to be randomly distributed across the points. However, since SAC analysis is an inferential statistic, results must be interpreted on the base of null hypothesis, stating that data are randomly distributed across points. Therefore, the p-value and the Z-score of the analysis must be assessed. The null hypothesis cannot be rejected if the pvalue is not statistically significant, thus meaning that data are randomly distributed across the points. On the other hand, for those cases in which the p-value is statistically significant, that is, the spatial distribution of the processed variable is not the result of a random spatial process, a positive sign of the Z-score statistic reveals the presence of a clustered pattern of the analysed feature, while a negative one means that the spatial pattern is dispersed.
The Global Moran's I statistic is useful to conclude if the spatial distribution of a given variable is globally clustered or not, but a local statistical analysis is necessary to derive the geographical position of a cluster. To this regard, Anselin (1995) introduced a class of local indicators of spatial association, known as LISA. Among all the possible LISA, in this study the Local Moran has been applied. The link between the two global and local indicators is that the global I can be seen as an average value (up to a factor of proportionality) of local I i , where the subscript i stands for the location of the measurement. It is worth to mention that, even if the Local Moran reflects the presence of significant local clusters, it does not mean that the global spatial distribution of the feature needs to be necessarily clustered. The local clusters are identified by means of the Moran scatterplot, in which points have x-y coordinates given by the original and the spatially lagged variable, respectively. This scatterplot is divided into four quadrants, in which the axes intersect in the centroid of the point cloud. The upper-right and thelower-left quadrants refer to a positive spatial autocorrelation, representative of similar values at geographically near location. On the opposite, the other two quadrants represent a negative spatial autocorrelation, meaning that there are dissimilar values at neighbouring locations. Combining the information provided by the Moran scatterplot with an indication of significance, it is possible to make a classification into four classes. In particular, significant clusters, identified in the upper-right and the lower-left quadrants, are denoted as high-high (HH) and low-low (LL), while significant outliers, identified in the lower-right and the upper-left quadrants, are denoted as high-low (HL) and low-high (LH), respectively.

| Study area and rainfall data set
The trend analysis detection here proposed is applied for Sicily, which is the largest island of the Mediterranean Sea. The island has an extension of about 25,700 km 2 and a morphology that is characterized by a mountain range along the longitudinal direction on the northern side and the Etna volcano on the eastern side. Elevation ranges from 0 to about 3,300 m a.s.l. in correspondence of the volcano Etna. Precipitation across the island has a significant spatiotemporal variability, with a mean annual precipitation (MAP) ranging between about 360 mm in the southeast part of the island and about 1,900 mm at the volcano Etna, situated in the northern east side (Caracciolo et al., 2018;Cipolla et al., 2020). Concerning temporal variability, the greatest part of MAP occurs during the winter seasons, while the summers are generally drier. Being Sicily the widest island of the Mediterranean area and lying in the heart of the Mediterranean Sea, it has a climate that can be considered as representative of an area quite extensive of the Mediterranean region.
The precipitation data set used for the study was provided by the regional agency SIAS (Servizio Informativo Agrometeorologico Siciliano; i.e., Agro-Meteorological Information Service of Sicily) and covers the period 2002-2019. SIAS rain gauge network consists of 107 tipping bucket rain gauge stations rather homogeneously distributed across the island with an average density equal to about 250 Km 2 Ágauge -1 . Data are retrieved with high temporal resolution (10 min) allowing time aggregation when necessary. The time-series are characterized by a high level of homogeneity, as declared by the SIAS Agency. The survey sites were chosen according to the World Meteorological Organization (WMO) standards and the strict criteria used in data detection, management and validation provide a high degree of uniformity over the whole island.
In order to avoid any misinterpretation in calculating trends, the original data set from SIAS was preprocessed by removing all those gauges with at least 1 year of missing data. This operation has led to take into account only 72 stations for further analyses (Figure 1). The metadata (i.e., name, ID, location, and altitude) of the selected stations are reported in Table S1 in the Supporting Information.
Before proceeding with the QR analysis, it was necessary to carry out a further two more operations. Firstly, in order to prevent a great number of null values from being weighted in Equation (1), all the zero precipitation values were removed from the data set. Moreover, to guarantee the correct sequence of rainfall in time, each record has been previously associated to a 'timestamp' that fixes its position in the timeline. Secondly, since the 'nid' method is sensitive to the presence of many equal values, as stated in Section 2.1, a Gaussian white noise (i.e., zero mean with a negligible standard deviation, here fixed equal to 10 −5 mm) was added to the original timeseries. Indeed, the raw 10-min data set includes plenty of values equal to the rain gauge resolution (i.e., 0.2 mm) and its multiples that invalidate such a method.
Starting from the modified data set, data have been aggregated to coarser time resolutions (20, 30, 40 min and 1, 3, 6, 12, 24 hr), for the further trend analysis detection with the QR. In order to leave unchanged the total precipitation amount at the different durations, a moving window with the size of the chosen duration and that moves with a time step equal to its size has been considered; at each step, all the data within the window have been summed up to return the value of the aggregated precipitation. The use of the above-mentioned 'timestamp' to fix the position of the 10-min record in the timeline, guarantees the correct sequence of rainfall in time also for the coarser time resolutions. Moreover, in order to compare results with those obtained with the MK analysis, the rainfall annual maxima for the previous nine durations have been extracted as well.

| Precipitation trends through the QR method at the gauge level
In order to derivate the slope (i.e., the trend magnitude) and the intercept for a various range of quantiles, the Rpackage 'quantreg' (Koenker, 2004) was used to apply the QR method to the rainfall time-series generated as in Section 2.4. Although the QR analysis was carried out for all the gauges displayed in Figure 1, for the sake of length, only the results related to the station of 'Palermo' are reported and discussed in this section. The significance of trends has been assessed by means of the Student's t test with reference to all the stations under study.
The results of such an analysis are reported in Section 3.2, where the procedure to identify the presence of spatial patterns has been carried out as well.
In order to provide a full view of the rainfall behaviour at different timescales aggregations, Figure 2 a, b, and c shows the results only for the shortest (10 min), F I G U R E 1 Spatial distribution of SIAS rain gauges overlaid on the digital elevation model (DEM) of Sicily. The gauge named 'Palermo' is highlighted with a bigger red point as compared with the remaining points [Colour figure can be viewed at wileyonlinelibrary.com] intermediate (1 hr), and the longest (24 hr) durations, respectively. Grey points represent the hourly rainfall intensities, i h , obtained from aggregated rainfall depths, while QR lines (coloured-solid lines) for the quantile levels 0.25, 0.5, 0.9, 0.95, and 0.99 are shown contextually to the ordinary least square (hereafter OLS) method line (black-dashed line).
Results confirmed that, as the considered quantile level grows, the intercept value grows as well. Indeed, for the highest quantiles, the highest rainfall intensity values have more weight than other values in the QR procedure, as expressed by Equation (1). This behaviour is not always valid for the slope since, depending on the data, there could be a trend inversion at any quantile level.
From Figure 2, it is noteworthy to highlight that the OLS regression line is not suitable to describe the behaviour of extreme events, both in terms of high and low intensity precipitation, since the sample includes a great number of ordinary low rainfall events that affect the intercept of the regression line, thus pushing the regression line towards low values of intensity.
By observing the slope of the regression lines, it is possible to make some inferences about the trend magnitudes for the durations considered. In addition, to make possible the comparison of the trend magnitudes provided by the slopes of the regression lines at different durations, it was decided to refer the slope units to the annual increment of hourly intensity (Δi h Áyear -1 ). A graphical representation of this transformation can be observed in the lower right panel of Figure 2. With this regard, Figure 3 represents the variation of slope with quantiles for the same durations shown in Figure 2. Black points are representative of the slope of the regression lines for different quantile levels, while the grey bands delimit the 90% confidence interval for the estimated slopes. The QR has been applied to 0.02 equal step quantiles ranging from 0.02 to 0.98, together with the quantiles 0.95 and 0.99 in order to extensively explore the relationship between slope and quantile. The slope of the OLS regression line, which is constant at a fixed duration, is indicated using a red-solid line whereas the related confidence intervals are displayed by means of some red-dashed lines. In general, if one observes the results of QR analysis, an increasing (decreasing) trend at the higher quantiles corresponds to the probability (with a certain level of significance) to have more (less) severe events. With this in mind, from the three panels of Figure 3, it is possible to notice that for the lower quantiles the slope is always close to zero, while it tends to generally increase for the higher quantiles. Such a behaviour is mainly due to the high number of low-intensity rainfall events which, for the lower quantiles, have more weight than the higher intensities. Moreover, with reference to the slope of the ordinary least square method (red-solid line in Figure 3), this is characterized by a slope significantly different from those relative to the highest quantiles. Such a behaviour is not observable for the lower quantiles, which have points mostly within the confidence interval of the mean (red-dashed lines). This certifies the fact that the ordinary least square method is only representative of the average behaviour of the sample and is not effective to characterize any trend at higher quantiles, which correspond to heavy and very heavy rainfall.
From the analysis of the confidence interval (grey bands in Figure 3), for each duration and rain gauge analysed, it can be observed that the standard error grows with the quantile; this is mainly due to the reduced extent of the sample used in the QR procedure for the higher quantiles.
For completeness, it is noteworthy to highlight the presence of some peaks in the slope at the 10-min duration for the higher quantiles (Figure 3a). This aspect can be explained by comparing the shape of the empirical pdf of rainfall intensity data at 10-min, 1-and 24-hr durations, as shown in Figure S1. In the case of the 10-min duration ( Figure S1a), the pdf shows the highest peak in correspondence of the gauge resolution (i.e., 0.2 mm) and several other peaks for the multiples of the gauge resolution that become smaller and smaller as the multiple of the gauge resolution increases. This reflects in weighting data in the QR procedure, thus generating the spikes shown in Figure 3a. This aspect is lost at the higher durations (i.e., 1-and 24-hr durations in Figure S1b and c, respectively) since data are aggregated and the empirical pdf of rainfall intensity assumes a smoother shape.

| Variability of precipitation trends
The QR analysis was used to detect eventual relationships between the magnitude and significance of the trends with their durations and quantiles. In particular, Section 3.2.1 will regard the detection of the relation between significant trend magnitudes with duration/quantile, grouping all the gauges. Section 3.2.2, instead, will be focused on analysing the previous relationships from a spatial point of view.

| Effects of duration and quantile on precipitation trends
Following the conjecture provided by Arnone et al. (2013) about the existence of trends in sub-hourly precipitation, a particular attention was paid to identify eventual positive trends for the events at the sub-hourly durations and for the higher quantiles. Moreover, a further and most crucial aspect is linked to the fact that eventual positive trends for these events would correspond to an increase in rainfall events with a shortduration and high-intensity that, in certain cases (e.g., small catchments with low times of concentration), could cause flash flood events and, as a consequence, higher risks of economic damages and fatalities. Figure 4 summarizes the QR results for all the considered rain gauges. It provides the percentage of gauges showing a significant positive (red), negative (green), or non-significant (grey) trend for the whole ensemble of F I G U R E 3 Slope of the regression lines versus quantile level for the station of 'Palermo' at (a) 10-min, (b) 1-hr, (c) 24-hr durations. Black points are representative of the slopes for various quantiles, while the grey bands stand for the 90% confidence intervals. The figure also displays the slope of the OLS regression line (red-solid line) and the related confidence intervals (reddashed lines) at the three considered durations [Colour figure can be viewed at wileyonlinelibrary.com] considered durations and 0.2, 0.5, 0.9, 0.95, and 0.99 quantiles. The two columns of Figure 4 are representative of two different levels of significance (i.e., α sig = .05 and α sig = .1, respectively).
Focusing on the lowest and intermediate quantiles, namely 0.2 and 0.5, respectively, it is possible to notice that a high percentage of rain gauges show a nonsignificant trend for all the considered durations for both the levels of significance. However, it is noteworthy that a noticeable percentage of stations has a positive trend at the sub-hourly durations, reaching the maximum value at the 10-min duration (i.e., 26.4 and 33.3% for α sig = .05 and α sig = .1, respectively). As the duration increases, there is an increase of the percentage of stations showing a negative trend and, at the same time, a corresponding decrease of the percentage of gauges manifesting a F I G U R E 4 Percentage of rain gauges characterized by a positive (red), negative (green), and non-significant (grey) trend coming out from QR procedure at 0.2, 0.5, 0.9, 0.95, and 0.99 quantiles at all durations [Colour figure can be viewed at wileyonlinelibrary.com] positive trend. This is more evident for the α sig = .1. When the percentages related to the shortest (10 min) and longest (24 hr) duration (Figure 4, top right panel) are taken into account, for instance, one can observe that the first bin (10 min) provides about 33% of gauges characterized by a significant positive trend and about 6% of gauges with a significant negative trend, while, on the other hand, the last bin of the panel (24 hr) displays an opposite behaviour (about 4% with a significant positive trend and about 20% with a significant negative trend). Furthermore, it seems that as the duration increases there is a reduction of the very ordinary (i.e., 0.2 quantile) and ordinary (i.e., 0.5 quantile) long-duration events. If these negative trends for the lower quantiles (e.g., 0.2) will persist especially at the longer durations, they could result in an increased risk of dry conditions (Lausier and Jain, 2018) with a reduction of available water resources, impacting, for instance, the agricultural sector (Field et al., 2018), hence causing noticeable economic damages.
Regarding the higher quantiles (i.e., 0.9, 0.95, 0.99), as quickly as the duration decreases, it is possible to observe a clear increasing pattern of the percentage of rain gauges with a positive trend. In particular, for the two highest quantiles (i.e., 0.95 and 0.99, respectively) and both the significance levels (i.e., α sig = .05 and α sig = .1, respectively), at least the 50% of the stations shows an increasing trend at the 10-min duration. Regarding the quantile level of 0.9, unless for the 10-min duration, which shows a percentage slightly lower than the other sub-hourly durations, the increasing pattern of positive trends with the decreasing of durations is still maintained. On the opposite, for the 24-hr duration, the stations with notrend strongly prevail, especially for 0.99 quantile level, where about 98% of stations reveal a non-significant trend at α sig = .05. All of the previous considerations suggest an increase in rainfall intensity provided by subhourly extreme events in Sicily, thus confirming what Arnone et al. (2013) had already guessed for Sicily region and other studies had found for other parts of the world (De Toffol et al., 2009;Adamowski et al., 2010). As already said at the beginning of this section, such a condition could lead to an increase of flash floods, with all of its consequences, but also to other consequences such as an increase in the soil loss, due to its erosion, and a consequent decreasing in production of these soils as found by Wei et al. (2009) in the North-West of China.
Differently from the case of the lower quantiles, in the last cases (i.e., 0.90, 0.95, and 0.99) less than the 10% of stations show a negative trend at all the durations under study; this percentage becomes smaller and smaller as quickly as the quantile increases. Moreover, differently than for the positive trend, it is not possible to recognize any pattern with duration. Furthermore, for the 0.99 quantile and α sig = .05, almost no station reveals a significant negative trend at all durations.
In order to describe the variation of trend magnitudes with quantiles and durations, Figure 5 reports, the empirical cumulative distribution function (hereafter ECDF) of the trend magnitude for the gauges showing a statistically significant (positive and negative) trend for α sig = .1 and the durations of 10, 30-min and 1, 6, and 24-hr and the 0.2, 0.95, and 0.99 quantiles.
Firstly, it is worth to focus on the entity of the trend magnitude, which considerably varies with quantiles. Indeed, at the 0.2 quantile (Figure 5a), it is possible to notice that, for all the durations, the trend magnitude values are close to zero, so that hourly and sub-hourly ECDFs, apart from those at 6 and 24-hr durations, cannot be distinguished. An opposite behaviour is illustrated in Figure 5b,c, related to the 0.95 and 0.99 quantiles, respectively. In particular, as the duration increases, the sample size becomes even smaller, due to a general loss of statistical significance (see also Figure 4) and, at the same time, the trend magnitude grows. Furthermore, as the quantile increases, the trend magnitude increases as well. As an example, focusing on the quantile 0.99 and sub-hourly durations, it is possible to notice that the 75% of stations present a trend magnitude greater than about 0.3 mmÁh -1 Áyear -1 , which is approximately the maximum magnitude value which characterizes the analyses carried out for 0.95 quantile. These findings confirm that the short-duration and high-intensity rainfall events are occurring more frequently in Sicily, at least with reference to the considered period, thus confirming that in the last years we are experiencing an increase in shortduration high intensity rainfall events that could be probably a consequence of climate change, as affirmed by Field et al. (2018).
The remaining ECDFs, related to all the durations and quantiles considered, are reported in Section S3.

| Spatial analysis of precipitation trends
In order to assess an eventual spatial distribution of the trends in Sicily, for each rain gauge station, Figure 6 represents the magnitude, direction, and significance of the trends detected with the QR approach. For the sake of clarity, the trend magnitude, which is expressed as the annual increment of hourly intensity, is symbolized by a graduated colour, the trend direction is represented by the symbol orientation (i.e., positive or negative), while the dimension of the triangle is related to three different classes of significance level (i.e., ≤.05; .05 Ä .1; > .1). The panels A, B, C, D and E in Figure 6 represent the spatial distribution of detected trends for the quantiles 0.2, 0.5, 0.9, 0.95, 0.99, respectively, while the columns 1, 2, and 3 are relative to the durations of 10-min, 1-, and 24-hr, respectively. For the sake of the length of the paper, the plots related to all the other durations considered for the previous analyses are reported in Section S4.
By observing this figure, it is possible to notice that the number of stations showing a significant trend (positive or negative) increases with the quantile level and decreases as the duration increases, thus confirming what has been already shown in Figure 4. Moreover, it is possible to observe that the trend magnitude values decrease moving from the higher to the lower quantile for each duration. The low trend magnitudes, which have been obtained at the lower quantile levels, are strictly connected to the high number of similar low rainfall intensity values weighted in the QR process. For this reason, these trends are expected to be low but, at the same time, not negligible, since they refer to very ordinary rather than extreme events. Indeed, even small changes in events that frequently occur throughout the year could lead to an important alteration of the local hydrological cycle.
In order to objectively identify the potential spatial clustered situations in the trend magnitude, the SAC analysis was carried out on those rain gauges showing a significant trend with α sig = .1 (i.e., all the medium and large colourful triangles in Figure 6) using the GeoDa software (Anselin et al., 2006). Table 1 shows the SAC analysis results in terms of Global Moran's I estimated value, Z-score, p-value, and sample size, related to all the possible combinations between durations and quantiles shown in Figure 6. It is important to remark that this analysis aims to determine whether the trend magnitude is globally clustered or not over the island. In such a kind of analysis the sample size plays a very relevant role;  Legendre and Fortin (1989) suggest to apply the SAC analysis to samples composed of about 30 values, in order to have a sufficient amount of data to reliably identify potential spatial patterns. Nevertheless, this threshold is not a strict limit and, in fact, Griffith (2010) suggests a minimum sample size of 25 elements.
The results of such an analysis are reported in Table 1, where the cases with a p-value lower than .1 but a sample size smaller than 25 are indicated with an italic font, while the cases with a p-value lower than .1 and a sample size greater than 25 data are reported with a bolditalic font. Normal font face, instead, is related to the quantile-duration combinations in which the p-value is higher than .1 and, consequentially, in which the spatial distribution of the magnitude can be considered as the result of a random spatial process.
F I G U R E 6 Spatial distribution of the gauges under study and magnitude (colour), expressed in mmÁh -1 Áyear -1 , direction (triangles orientation) and significance (triangle size, large for α sig ≤.05, medium for .05 < α sig ≤.1 and small for α sig >.1) from QR at 0.2, 0.5, 0.9, 0.95, 0.99 quantiles for 10-min, 1-and 24-hr durations. The bold letters, A-E, stand for the quantiles, while the bold numbers 1-3, denote the durations [Colour figure can be viewed at wileyonlinelibrary.com] Focusing on the lower (i.e., 0.2 and 0.5) and highest quantiles (i.e., 0.99), the results of Table 1 indicate that no spatial patterns of the trend magnitude can be found at any duration, given the resulting p-value and Z-score statistics. Indeed, regarding the 0.99 quantile and 10-min duration, the p-value of the SAC analysis is higher than .1, probably due to a great variability of the trend magnitude. Therefore, it is not possible to assume that, globally, the distribution of the trend magnitude is significantly clustered.
On the contrary, when one observes the quantile 0.9 it is possible to highlight that for the 1-and 24-hr durations, the p-value is less than .1 and, contemporary, the Z-score is positive. This means that the spatial distribution of the trend magnitude is clustered, even though the result related to the 1-hr duration may be considered more reliable than the 24-hr one due to the greater sample size.
Regarding the 0.95 quantile, while all the durations exhibit a p-value lower than α sig = .1 only the 10-min and 1-hr durations are characterized by a sample size greater than 25. Furthermore, the Z-score is always positive, indicating that the trend magnitude is clustered.
The results of the SAC analysis for the remaining quantiles and durations are reported in Section S4. With reference to all the durations under study, the Global Moran's I generally resulted in a spatially clustered distribution of the trend magnitude especially for the higher quantiles. This behaviour is particularly enhanced at the 0.95 quantile for the whole ensemble of sub-hourly durations (i.e., 10, 20, 30, and 40 min), which are characterized by the greater sample sizes. Furthermore, for the 0.9 quantile, the same considerations are valid for the sub-hourly durations of 30 and 40 min. The reader is referred to Table S2 for more details.
In order to identify the clusters' position, a Local Moran analysis (Anselin, 1995) has been applied at the same samples and for the same duration-quantile combinations used in Figure 6. In particular, panels a to c in Figure 7 show the results for 95th percentile at 10-min, 1-, and 24-hr duration, respectively. This quantile has been chosen as representative mainly because the sample size criterion is satisfied both for 10-min and 1-hr duration (i.e., 44 and 26 locations), while the sample size is the largest among those relative to the 24-hr duration (i.e., 17 points), even if it is smaller than 25 elements.
In Figure 7, red and blue circles are relative to HH and LL clustering cases, respectively, and characterized by a p-value lower than .1. The significant outliers, instead, are marked with the diamond shape and filled with pink (HL) or light blue (LH) colours. Therefore, the crosses represent those locations in which the significance level exceeds .1.
It is worth to mention that for all the panels in Figure 7 it is possible to distinguish two significant clusters. In particular, in the south-east there is a HH cluster, while the northern-west part is characterized by a LL cluster. Looking at the dimension of these two clusters, it is possible to observe as their extension changes with the T A B L E 1 Global Moran's index estimated value, Z-score, p-value, and sample size for 0. 2, 0.5, 0.9, 0.95, and 0.95 quantiles and 10-min, 1-and 24- Note: All cases in which the p-value is less than .1 but the sample size is less than 25 are indicated with an italic font, while those cases in which the p-value is less than .1 and the sample size is greater than 25 are written in bold-italic font. Normal font face is related instead to those cases in which the spatial distribution of the trend magnitude is the result of a random spatial process (p-value greater than .1).
duration. At the 10-min duration (Figure 7a), for instance, a relevant number of stations composes both the clusters (i.e., 8 and 12 sites for the LL and HH, respectively). Moreover, by comparing this panel with the correspondent one in Figure 6 (i.e., panel D1) it is possible to notice that the HH cluster includes those stations characterized by the highest values of trend magnitude. Therefore, the LL cluster is composed by a group of gauges where the trend magnitude is moderately positive. Starting from this consideration, it is worth to highlight that the acronyms HH and LL do not necessarily refer to a cluster of gauges in which an increasing or a decreasing trend is observed. Focusing on the 1-and 24-hr durations (i.e., Figure 7b,c, respectively), both the clusters reduce their dimension. Furthermore, the LL cluster seems to be more confined in the north-west side of the island. Looking at the trend magnitude of the gauges forming the LL cluster (panels D2 and D3 in Figure 6, respectively), it is also possible to observe that, in this case, they match with a decreasing trend magnitude cluster (i.e., except for the upper-left point in LL cluster of Figure 7c).
Looking at the Figure S5, it is worth to focus that, when the HH and LL clusters are visible, their location is, more or less, the same of that highlighted in Figure 7. This consideration is not valid for 0.2 quantile at 1-hr duration, where different clusters can be detected.

| Mann-Kendall test for rainfall annual maxima trends
This section shows the results of the MK method, applied to the pre-whitened annual maxima extracted by the SIAS data set, as well as a qualitative comparison with the QR method results, with the aim to explore the advantages and drawbacks of both methodologies in trends detection.
As an example, the results of such an analysis for the rain gauge of Palermo are shown in Figure 8, where the subplot a) shows the annual maxima extracted for the durations of 10 min, 1, and 24 hr, while the remaining plots show the results regarding the trend-duration dependency ( Figure 8c) and the spatial variability of trends magnitude, direction, and significance ( Figure 8b). The magnitude, in this case, is obtained through the Sen's slope method and represents the variation of rainfall annual maxima per year (mmÁyear -1 ).
Focusing on canonical durations, it is possible to observe two different behaviours between positive and negative trends with duration. Indeed, as the duration increases, the percentage of stations characterized by a negative trend becomes greater, while the one featured with positive trend decreases. In particular, at the 24-hr duration, these percentages reach 18 and 0% for negative and positive trend, respectively. Regarding the sub-hourly durations, instead, no patterns with duration are noticeable. In any case, it is worth to observe that the majority of the rain gauges is characterized by a non-significant trend.
Although it is not possible to make a direct comparison between the results obtained with the QR and MK methods because of the different data sets they work with, it may be useful to highlight some differences, strengths and weaknesses of the two methods. First of all, focusing on the data, since the MK test works with the annual maxima (i.e., only a value of rainfall per year), it might be unsuitable to work with short data sets (i.e., few years of observed data) since it would return results statistically not significant. Moreover, working with annual maxima implies that all the rainfall depths slightly lower than the annual maxima, despite the fact that could be F I G U R E 7 Local Moran analysis for the 0.95 quantile at (a) 10-min, (b) 1-hr, and (c) 24-hr duration. The LISA is applied to both the positive and negative trend magnitudes with a significance level of .1. Red and blue circles are relative to high-high (HH) and low-low (LL) clustering cases, respectively, and characterized by a p-value lower than .1. The significant outliers are marked with light blue and pink diamonds for low-high (LH) and high-low (HL) clustering cases, respectively. The crosses represent those locations in which the significance exceed the level of .1 [Colour figure can be viewed at wileyonlinelibrary.com] even higher than the annual maxima of other years, are discarded from the analysis.
In the QR method, instead, all data are processed to extract information about trends, even if only those exceeding the threshold related to the examined quantile are weighted in a more significant way. Nevertheless, by looking at the ECDFs of the aggregated time-series, related to 10-min, 1-and 24-hr duration, for 'Palermo' rain gauge, it is possible to observe that about 1,200, 470, and 100 values are above the 95th percentile, respectively. Moving to the 99th percentile, the number of observations above this percentile drops to about 240, 90, and 20, respectively, but, in any case, higher than those used in the MK procedure that are 18 values (i.e., the annual maxima) for each of the considered durations. This could be the main reason why the QR approach is able to detect a consistent number of stations featured by a positive trend for the 0.95 and 0.99 quantile (i.e., about 57 and 56% at α sig = .1, respectively), while this information is not captured by the MK procedure.
With reference to the spatial patterns (Figure 8b), most of the stations shows a non-significant trend at all durations. Since the maximum number of rain gauges showing a significant trend is equal to 13 (at 24-hr duration), the SAC analysis has not been carried out to avoid inconsistent and/or unreliable results.
For all of these reasons, the QR could be a valid alternative to the MK procedure to detect trend in extreme rainfall, especially when the period under study is short. It is further worth to highlight that, also in the QR procedure, the loss of significance recorded especially at the longest duration could be attributable to a too small sample size.

| CONCLUSIONS
Climate change is still today considered as one the most critical issues to face by the scientific community. Its F I G U R E 8 Annual maxima for the rain gauge of 'Palermo' of the SIAS data set at 10-min, 1-hr, 24-hr duration (a), spatial distribution of magnitude, direction and significance of trends (b), and percentage of stations having a negative (green), positive (red), or non-significant (grey) trend obtained through the MK test (c) [Colour figure can be viewed at wileyonlinelibrary.com] manifestations could exasperate different kinds of extremes, such as water scarcity and/or heavy rainfall events.
In this perspective, this work has been conducted with the aim to identify potential changes in rainfall intensity for Sicily (Italy) over the last 20 years. The rainfall data time-series at different durations (i.e., 10, 20, 30, 40 min and 1, 3, 6, 12, and 24 hr), aggregated from the original 10-min observations provided by the regional agency SIAS, have been elaborated using the QR method, at the quantiles ranging from 0.02 to 0.98 and considering also 0.95 and 0.99. The achieved trends have been analysed under two different significance levels (i.e., α sig = .05 and α sig = .1). Moreover, in order to give a comparison of QR with one of the most employed methods in trends detection, the annual maxima at the above-cited durations have been extracted and processed with the MK test coupled with the Sen's slope method.
When considering the 0.2 quantile, results revealed the existence of a decreasing (increasing) pattern of the positive (negative) trend with duration, even though the greatest part of the rain gauges have shown a nonsignificant trend. A decreasing trend at this quantile, which can be assumed to be representative of the dry tail of precipitation distribution, could lead in the future to a higher risk of the occurrence of too dry conditions, thus to water scarcity.
Moving to the highest quantiles, the dependence of the negative trends with duration disappears and, in addition, the percentage of stations showing a negative trend tends to zero. On the opposite, a clear increase of the percentage of stations characterized by a positive significant trend has been achieved moving towards the sub-hourly durations. This is especially highlighted at the 10-min duration, where more than a half of the considered stations manifests a positive trend in rainfall intensity. In addition, a general loss of significance of trends with duration is registered. These results confirm what Arnone et al. (2013) guessed years ago about extrapolating the increase in extreme events towards sub-hourly durations.
Not only the temporal, but also the spatial variability of precipitation trends has been investigated in this work, by means of the Global and the Local Moran Indexes. According to the SAC analysis results, it is possible to assume that the spatial distribution of the magnitude is clustered at any durations for the 0.95 quantile, while only at 1-and 24-hr durations for the 0.9 quantile. For the LISA analysis at 0.95 quantile, the majority of the considered stations revealed a positive significant trend, mainly located in the southeast of the island whereas the west coast of the region mostly presents a cluster characterized by a decreasing trend.
The MK test on annual maxima did not confirm the QR results, since about the 90% of the stations revealed a non-significant trend and this is probably due to the small extension of the annual maxima data set used.
On the contrary, since regression procedure in the QR method is based on the whole precipitation data set, it is demonstrated to be less sensitive to outliers in the data (Koenker, 2005;Roth et al., 2015) and to better look at the tails of the probability distribution functions (Lausier and Jain, 2018;Villarini and Slater, 2018). This provides a more complete and reliable estimation of trends respect to the MK procedure.
In conclusion, the general increase of sub-hourly rainfall intensity, especially the more severe ones, over the considered period (from 2002 to 2019) could be considered as a typical sign of a changing climate. Particular attention is given to this result, since this kind of rainfall could cause several damages due to the high quantitative of rainfall that can occur also for a short duration event. Especially in regions like Sicily, since the presence of many small catchments with low times of concentration and sometimes highly urbanized, the increasing occurrence of short-duration and high-intensity rainfall events could lead to a higher number of flash floods, thus increasing the risk of economic damages and fatalities.