A Multi-Modal Machine Learning Approach to Detect Extreme Rainfall Events in Sicily

In 2021 300 mm of rain, nearly half the average annual rainfall, fell near Catania (Sicily island, Italy). Such events took place in just a few hours, with dramatic consequences on the environmental, social, economic, and health systems of the region. This is the reason why, detecting extreme rainfall events is a crucial prerequisite for planning actions able to reverse possibly intensified dramatic future scenarios. In this paper, the Affinity Propagation algorithm, a clustering algorithm grounded on machine learning, was applied, to the best of our knowledge, for the first time, to identify excess rain events in Sicily. This was possible by using a high-frequency, large dataset we collected, ranging from 2009 to 2021 which we named RSE (the Rainfall Sicily Extreme dataset). Weather indicators were then been employed to validate the results, thus confirming the presence of recent anomalous rainfall events in eastern Sicily. We believe that easy-to-use and multi-modal data science techniques, such as the one proposed in this study, could give rise to significant improvements in policy-making for successfully contrasting climate changes.


Introduction
Is it possible to detect extreme rainfall events areas by clustering spatio-temporal data?The intensification of weather extremes, which is dramatically changing the climate scenario worldwide, is currently thought to be as one of the most important factors related to green-house effect and climate change [1][2][3][4][5][6][7][8] .The increase in the frequency and intensity of daily temperatures has contributed to a widespread escalation of daily precipitation 9,10 .Moreover, severe weather and climate events, interacting with exposed and vulnerable human and natural systems, can lead to disasters which require an extraordinary adaptation ability 4 .It is therefore mainly for this reason that, nowadays, the study of climate change is not only about temperature increase, but it also focuses on catastrophic rainfall extreme events and drought 8,11 .The concept of extreme precipitation and its changes in response to warming are well described in 12 .For this reason, the scientific community faces an increasing demand for regularly updated estimations of evolving climate conditions and extreme weather events 5,11,13 .Moreover, a correlation between changes in heavy precipitation and landslides in several regions has been found in 4 .More specifically, it is possible to identify 3 examples of extreme weather events, that have raised the question of a potential link to climate change: more intense precipitation events, increased summer drying over most mid-latitude areas and increase in tropical cyclone peak winds intensities 14 .These results show that rainfall extreme events are related to climate change and represent the triggers of a chain of reactions involving several human activities.The change in temperatures will, in fact, have serious long-term effects [15][16][17] , although extreme rainfall events will also cause a short-term danger to the environment and the population.In more recent years several extreme events all over the world caused large losses of lives, as well as a tremendous increase in economic losses from weather hazards 18 .Such disasters have forced public opinion to consider climate change as the main cause of these events 19 and to deeply analyse the economic consequences of climate change in terms of investments and productivity 20,21 .A relevant example of this regards wine industry.For instance, in the past two decades, Sicilian winemakers have enhanced the biological production of wine all around the island, especially on the slopes of Mount Etna.Although wine is not essential to human survival, it is an important product of human ingenuity and its economy is rapidly growing 22 .Agricultural activities depend on climate and are interconnected to weather changes.Any shift in climate and weather patterns may potentially affect the entire local wine industry 23 and the stability of many crops, thus undermining the related economies 20 .Any shift in climate and weather patterns may potentially affect the entire local wine industry 23 .Abnormal climate changes might also undermine the stability of crops and might be critical for the related economy 20 .Considering all of these aspects, in Mediterranean areas, rainfall is probably the most important climatic variable due to its manifestation as a deficient resource (dryness) or a catastrophic agent, such as water bombs 24 .Therefore, many challenges arise during the measurement of the precipitation.For instance, in situ measurements are especially affected by wind effects on the gauge catch, particularly for snow but also for light rain 15 .Moreover, to reduce this uncertainty, it is crucial to analyze spatio-temporal data in the most efficient way 25,26 .
In this regard, over the last decades scientists conducted several studies on rainfall time series.These studies investigated potential trends in different rainfall indicators, such as total and maximum annual precipitation and mean daily intensity 27 .A tendency toward higher frequencies of heavy and extreme rainfalls emerged for some areas 28 .In most of these areas, an increase in total precipitation has also been observed, for instance in 24 , thanks to the analysis of 247 stations over the 1921-2000 period.However, the correlation between the increase of total precipitation and extreme events is not always clear, as in other areas (i.e.Italy) several authors have observed an increase in heavy precipitation, together with a tendency towards a decrease in the total amount of precipitations 29 .Among the studies mentioned, a few of them were specifically focused on the Mediterranean areas, given their peculiar climate, which is affected by interactions between mid-latitude and tropical processes, lying between the arid climate of North Africa and the temperate and rainy climate of central Europe.For these reasons, even relatively minor modifications of the general circulation can lead to substantial changes in the Mediterranean climate 27,30 , including rainfall frequency 31 , thus making these areas vulnerable to climatic changes and in particular to catastrophic precipitations.In this setting, scientists analysed the region of Sicily to identify climate change signals, as for instance in 32 .In most of those studies, the authors analysed annual, seasonal and monthly rainfall data in the entire Sicilian region, showing a global reduction of total amount of annual rainfall 32 .For example, in 27 the annual maximum rainfall for fixed time duration of 1, 3, 6, 12 and 24 h, and the daily rainfall series recorded from 1956 to 2005 in approximately 60 stations were analyzed using the non-parametric Mann-Kendall test 33,34 .
Results of this study, confirmed an increasing trend for rainfall of short duration, in particular for the 1 hour rainfall length.On the other hand, time-persistent rainfalls exhibited a decreasing course 33,34 .In particular, heavy-torrential precipitation have been reported to be more frequent at a regional scale, while light rainfall have shown negative trends at some sites.In 35 the presence of linear and non-linear trends in 16 series from rain gauge stations, mostly placed in the eastern Sicily, was studied.The results indicated a different behaviour according to the time scale: for short duration, historical series generally presented increasing trends, that switched to decreasing for longer time courses.
A total of 67 sites of daily precipitation records over the 1951-1996 period in Italy were also analyzed in 29 considering seasonal and yearly total precipitation, number of wet days and precipitation intensity with the aim of evaluating the trends both from the single-station records, and for larger areas by using averaged series.Results showed that the trend for the number of wet days in the year was significantly negative throughout Italy, particularly stronger in the north than in the south, especially in winter.A tendency towards an increase in precipitation intensity, which was globally less strong and significant than the decrease in the number of wet days was also found.
In 36 the authors identified the presence of homogeneous areas over Sicily using the Regional Frequency Analysis (RFA), which is a procedure estimating the frequency of rare events at one site by using data from several sites 37 , used frequently in the analysis of environmental data 38 .They also developed Principal Component Analysis (PCA) followed by a clustering analysis, performed by applying the K-Means method, to identify regional groups, starting from annual maximum series for rainfall duration of 1, 3, 6, 12 and 24 h over about 130 rain gauges.
One of the most interesting papers studying different rainfall time series in Sicily is 28 , where the authors investigated temporal changes in extreme rainfall by performing a regional study.In particular, a regional frequency analysis based on L-moments approach 39 was applied to 1, 3, 6, 12 and 24 h annual maxima rainfall (AMR) series grouped per homogeneous regions, identified through a hierarchical cluster analysis 40 .Changes were investigated in a long-term dynamic (from 1928 to 2009) with special reference to the last forty years.The study 40 detected an increasing trend on rainfall extreme events between 2003 and 2009 with several heavy localized storms all over Sicily and a remarkable tendency towards more intense storm events during the 2000's affecting mainly the outer western part of the region.On the contrary, the increasing trend in extreme rainfall detected in eastern Sicily, have been considered only apparent, as related to a few severe local storms.
In our work we present a spacial and temporal clustering analysis on rainfall data over Sicily, performed by the Affinity Propagation clustering algorithm.The research interest towards this topic lies in new alarming violent rainfall events occurred in East Sicily at the end of 2021 41,42 .The main objective of the research is to introduce a machine learning methodology to detect extreme rainfall events.In particular, this work analyzes where extreme events took place in the region of Sicily, highlighting the differences and the similarities both between stations and years, i.e. in a spatio-temporal analysis.Differently from 28 , in our study clustering is not only used for identifying homogeneous sub-regions, but also to detect critical rainfall sites.Moreover, while in 28 the authors focused on finding long-term trends, we concentrated our attention on short-term changes between 2009 and 2021, analyzing high-frequency data, so as to obtain clusters specifically related to extreme events.
The present work faces several steps: • Consider the RSE (Rainfall Sicily Extreme) dataset, using originally recorded data, which allowed us to highlight different climate behaviours respect to previous studies.
• Apply clustering to the Sicilian data in order to cluster regions and detect extreme sites according to rainfall data observations.To further validate the clusters obtained we defined rainfall indicators, which allowed us to validate and understand further the meaning of the obtained clusters.
• Use a multi-modal approach to merge both geographical and temporal information gathered from the dataset.
• Highlight significant characteristics of rainfall distribution over Sicily, detecting an increasing trend on extreme events in East Sicily, that reinforces the results of the state of the art in 28 .
According to this scheme, the paper is structured as follows: in Section 1 the regional dataset used in the analysis, including the data pre-processing, is presented.In Section 2 the methods applied in the study, in particular the adopted clustering algorithm, and the statistical validation methods, are introduced.In Section 3, we describe the application of methods to the considered rainfall datasets.In Section 4 we report the discussion of the results, concerning each analyzed variable, and the most relevant conclusions drawn.Furthermore, we report in the supplementary material the analysis concerning the annual histograms of specific rain gauges and local data plots at different levels, as well as the complete annual clustering results.

RSE: the Rainfall Sicily Extreme dataset
The dataset used in this analysis consists in geographical rainfall records with a 10 minutes periodicity from 2009 to 2021, provided by SIAS, the Servizio Informativo Agrometeorologico Siciliano 43 .The dataset together with the code is available at the following GitHub Repository 44 .
The most common rainfall measurement gathered from the database is the number of millimeters (mm) of rain in a given period.Accordingly, six collections were considered, as described in Table 1.C.A and C.B contain 13 datasets per stationone per year -with the original data and the weekly mean data, respectively.C.C and C.D include one full dataset per stationinvolving all the records from 2009 to 2021 -with the original data and the weekly mean data, respectively.C.A s and C.B s are subsets of C.A and C.B, respectively, since one station per time is considered, so that each of them includes 13 datasets.

Data preprocessing
We will now describe the initial data selection process, obtained through the analysis of annual data.On the basis of an initial graphical analysis reported in the SI document, we decided to select the most extreme stations.A station is considered extreme if it is possible to observe a high amount of rain in a relatively short time interval.We implemented this concept of "extremeness" using the following strategy.First, we considered the following data for all the 96 available stations in Sicily and for all the years: • The total annual precipitation in mm (tot).
• The percentage of rainy days over the year (rd).This measure indicates the percentage of days with more than 1 mm of rain.
• The mm of rain during the rainiest day in the year (dmax).

3/14
Afterwards, a selection strategy has been applied.Extreme rainfall events are generally characterized by the increasing of either drought and/or excessive wetness 24 .The logical rule below highlights precisely such characteristics: 1) Fix a station.
2) Compute µ 1 : the mean over years of the rd annual indicator.
3) Compute µ 2 : the mean over years of the dmax annual indicator.4) Fix a year y.

5)
If the rd value in the year y is less than µ 1 and the dmax value in the year y is grater than µ 2 , then the year y is considered as extreme.Otherwise no.Since the procedure works year by year, we selected the stations satisfying the extreme events detection rule for at least 3 years (the stations respecting this condition for at least one year were 85 out of 96, almost all).In this way, we obtained 32 stations out of the 96 rain gauges.Furthermore, we decided to include all of the provincial capitals in the region, thus obtaining the 34 stations shown in Fig. 1.
After the selection, we observed rainfall data time series, by fixing a station and using full, annual, and monthly data plots, as well as mean data graphics (all details regarding these initial observations are reported in the SI document).This preliminary analysis lead to different reasoning.The full plots proved the necessity of quantifying and understanding variation in the stations time series behavior.In contrast, the annual plots showed a typical seasonality pattern.Moreover, the graphics observation led to the idea of comparing annual time series.Finally, a similar reasoning has been done with regard to the monthly view.
All of the above considerations suggested us to highlight the differences and the similarities both among stations and years, in order to identify multi-modal (geographical and historical) rainfall changes.Instead of performing classical time series analysis, we proceeded by applying the suitable clustering algorithms described in the following sections.

Methods
Clustering is an unsupervised machine learning methodology 45,46 .Its goal is to detect groups of observations sharing similar characteristics.More precisely, it consists in the partitioning of a dataset into subsets, so that the data in each subset are characterized by a higher similarity than elements in different sets, according to some defined distance measure.
Two main types of clustering techniques can be defined: methods in which the number of clusters needs to be established a priori, as, for instance, the K-Means algorithm 47 , and algorithms in which, instead, the number of clusters is inherently estimated during the optimization phase, such as the Affinity Propagation 48 .The last one has been used in this work, since no prior knowledge on the number of expected clusters was available.

Affinity Propagation algorithm
Affinity Propagation (AP), introduced by Frey and Dueck in 2007 48 , and its extension to Hierarchical Affinity Propagation 49 , are nowadays becoming extremely popular due to their simplicity, general applicability, and performance.
AP takes as input the measures of similarity between pairs of data points, and simultaneously considers all of them as potential exemplars.The number of clusters does not need to be defined in advance, indeed the algorithm is based on the hypothesis that the so called "real-valued messages" are exchanged between data points until a high-quality set of exemplars, together with the corresponding clusters, gradually emerges.Given that no assumption on the number of clusters was requested in our case, AP has been a natural choice.
The algorithm requires two inputs parameters 48 : • Similarities s(i, k) between data points, representing how similar a point is to be another one's exemplar.If there is no similarity between two points, as in this case they cannot belong to the same cluster, this similarity can be omitted or set to -∞ depending on the implementation.
• Preferences s(k, k), indicating each data point's suitability to be an exemplar.Since some prior information which points could be favored for being an exemplar can be available, it can be represented through preferences.
Similarity is usually defined starting from the negative Euclidean distance or the Pearson correlation coefficient, depending on the considered situation.
If all data points are supposed to be equally suitable as exemplars, the preferences should be set to a common value, such as for example the median of the input similarities, thus resulting in a moderate number of clusters, or their minimum, thus resulting in a small number of clusters 48 .In this work we initialized the preferences to the median and the availabilities to zero, a(i, k) = 0.Each iteration step of the optimization performance is composed by 2 main message-passing steps: 1. Computing responsibilities: where s(i, k) and s(i, k ) are similarities, while a(i, k ) are availabilities.

Computing availabilities
where r(k, k) are the self-responsibilities, while r(i , k) are general responsibilities.To limit the influence of strong incoming positive responsibilities, the total sum is lower bounded, so that it cannot be negative.
The "self-availability", a(k, k) is updated differently, as follows: The way for calculating how suitable a point is for being an exemplar is that it is favored more if the initial preference was higher, but the responsibility gets lower when there is a similar point that considers it as a good candidate, so there is a 'competition' between the two, until one of the two options is chosen in some iteration.The above procedure may be terminated after a fixed number of iterations, after changes in the messages fall below a threshold, or after the local decisions stay constant for a given number of iterations 48 .

Statistical validation
To assess the presence of statistical differences between 2 communities we made use of the well known Kruskal-Wallis test [50][51][52] .

Kruskal-Wallis test
It is a non parametric statistical test that assesses the differences among three or more independently sampled groups 53 .Kruskal-Wallis test is used to determine whether or not there is a statistically significant difference between the medians of three or more independent groups.It does not assume normality in the data and is much less sensitive to outliers than the standard analysis of variance (ANOVA) 54 .The test is based on the null hypothesis H 0 55 , which allows one to state whether the considered samples are realizations of identical populations.The application of the test returns a p-value which confirms or rejects the null hypothesis.If p < 0.05, then the null hypothesis is rejected, on the contrary, if p ≥ 0.05, then the null hypothesis is confirmed 54 .The related p-value for the test is computed using the assumption that H has a χ 2 distribution.

Experiments
This section explains the experimental procedure followed to design and apply the clustering algorithm.Global and local clustering analysis have been preformed by means of high frequency (measurements collected every 10 minutes) and weekly averaged data.The reason of this choice lied in the need to reduce the dataset dimension for minimizing sensitivity to outliers and oscillations in the original time series of observations.Two main streams of experiments were performed: 1. Geographical (or spacial) clustering: it consists of grouping similar geographical stations together along different time horizons.

2.
Local (or temporal) clustering: it consists of grouping similar years together on each single location.
For the first category, we ran the algorithm four times, according to the four Collections of datasets C.A, C.B, C.C and C.D described in Table 1.In contrast, the second category involves C.A s and C.B s of Table 1.The Affinity Propagation algorithm has been implemented in Python programming language (version 3), making use of with the Scikit-learn library (V.1.0.2),which is a free software machine learning library for Python 56 , designed to inter-operate with the Python numerical and scientific libraries NumPy (V.1.21.4) 57 , SciPy (V.1.8.0) 58 and Pandas (V.1.3.5) 59.The algorithm was implemented using default hyper-parameters.For instance, convergence_iter was set to 15, that is the number of iterations with no change in the number of estimated clusters, that stops the convergence.Moreover, the preference value was set to the median of the input similarities.

Metrics for Affinity Propagation
Two different metrics of similarities were used in the Affinity Propagation algorithm: the Euclidean Metrics and the Correlation Metrics.We present in the following subsections the results obtained for both.
Euclidean Metrics We first conducted clustering using the Euclidean affinity metric, which resulted in a principal large cluster and few smaller communities consisting in one element each.For this reason, we decided to apply an iterated version of the AP algorithm in order to detect new geographical clusters, at first glance hidden by the anomalies.To this aim, the AP algorithm based on a particular multi-step structure was implemented as follows: 1) AP is applied to the whole considered collection of datasets.
2) The exceptions found at level one from the data are removed, and the AP algorithm reiterated over the remaining datasets.
3) The process is repeated from Step 1.
Correlation Metrics On the basis of the theoretical arguments reported in the Affinity Propagation algorithm subsection, the Correlation distance was also chosen as a affinity metric for a second exploration analysis.In this case no multi-step procedure was needed.

Clusters validation procedure
In order to understand the rainfall phenomena that mostly characterize the clusters, several rainfall indicators over the time series were introduced, as reported in Table 2.
We assembled the original 10 minutes records according to specific needs: naturally an hour data includes six consecutive records summed up together, whereas a day consists of the sum of 144 consecutive data.We also computed the total number of rainy hours, where one hour is considered "rainy" if its amount of rain is higher than zero.Hence, some of the introduced indicators are the percentages of light (l), moderate (m), and heavy (h) rainy hours over the total.Moreover, we considered the absolute number of violent rainy hours v, which is not expressed in percentage since it represents very rare events.The Kruskal-Wallis test was then applied to the indicators, in order to understand which of them better characterize clusters.To this aim, the SciPy scientific library has been used.Indeed, it provides algorithms for many classes of problems, extends standard tools of array computing, wraps up highly-optimized implementations, is easy to use, and enlarges NumPy 58 .
The Kruskal Wallis test was applied to each indicator and to all the experiments described above, according to the following logical evaluation steps: 1) Fix an indicator i.
2) Run the clustering algorithm.
3) Create an array k with one element for cluster.Every element of k is in turn an array a c , containing the indicator values of the stations belonging to that cluster.
4) Run the Kruskal-Wallis test on k.
5) If the p-value is less than 0.05: i is considered as characterizing for the clusters.Otherwise no.

Results and discussion
In this section we report the main results obtained from the study for both geographical and single station investigations.Additional details on the whole dataset results are reported in the supplementary information document.

Geographical investigation
Results of this set of experiments are visualized in the Sicily map of Fig. 2, where the 34 stations with names or symbols coloured according to their relative clusters are drawn.When the multi-step version of the algorithm is applied, different shapes for the points are used.Specifically, circle, squares and diamond markers represent clusters resulting from the first, second, and third iterations, respectively.

Annual clustering
The results of the geographical clustering year by year for C.A and C.B (Table 1), both with the Euclidean and Correlation similarities are included in the SI document.Additionally, a video showing the clustering results proceeding in years is available for each collection.We hereby report the main results drawn from the several performed experiments: • Euclidean metrics -C.A (Video_1): in this case the annual results consist mostly of a principal cluster (at most two) and some anomalous stations.Proceeding in the years, a flow in anomalies that goes from western to eastern Sicily is detectable.We believe that anomalous clusters seem to be more susceptible to extreme events.In fact, since them change drastically in intensity according to the specific location, this legitimizes the algorithm in finding exceptions, namely only one station in one cluster.In order to validate this observation, a case by case analysis has been carried out and reported in the SI document for the year 2021.
• Euclidean metrics -C.B (Video_2): the reduction in the dataset size led the clusters to be more uniform and referred to geographical divisions.However, there are some exceptions, mainly in the South-East Sicily, and in the neighbourhood of Palermo.As in the previous case, this represents a trend on extreme events, more diffused in the East side of the island.
• Correlation metrics -C.A (Video_3): in this case a geographical clustering pattern was obtained, identifying eastern and western Sicily.This is coherent with the fact that the Correlation metrics finds shape similarities and it is less sensitive to the micro-climatic differences.
• Correlation metrics -C.B (Video_4): here the combination between dimensionality reduction and correlation metrics brings to a rough splitting of the island.The number of clusters does not exceed 3 and often very far away stations are grouped together in the same cluster.
The results of the 4 settings for the year 2021 are in line with our initial research hypothesis, for which the anomalies correspond to extreme stations.This behaviour was further confirmed by the Kruskal-Wallis test, as reported in the SI document.
Moreover, East Sicily emerges as the most extreme zone of the island, confirming the real occurred events reported in 42 and discussed in 41 .
Finally, Video_2 shows that the use of weekly averaged data (C.B) reduces noise, thus giving rise to balanced presence of both anomalies and territorial clusters in the Euclidean case.On the other hand, Video_4 shows that CB, differently from C.A, does not provide clear geographical splitting in the Correlation case.
In conclusion, in the annual case the use of Euclidean metrics led to detect the anomalies, while in contrast, the use of the Correlation metrics as a similarity measure allowed us to identify more uniform clusters.

Full clustering
We report here the full clustering results, obtained using C.C and C.D of Table 1.Similarly to the annual case, the use of Euclidean metrics brings to exceptions detection, by highlighting the presence of clusters composed by a unique site.We believe that the reason why anomalous clusters are independent lies on the fact that extreme events intensities are very different among sites.Fig. 2a shows the presence of one principal cluster and many anomalies, such as Pedara, Augusta and Siracusa.In contrast, Fig. 2b reports different principal clusters -geographically distributed -and only one exception: Pedara.
In order to validate results, we carried out a case by case analysis.First of all, the full and the annual clustering results in the case of Euclidean metrics (C.A and C.C, respectively) are compared in Fig. 2a by counting how many times the stations has been clustered as anomalous in the annual case.It turns out that the stations with an higher counter are the ones clustered as anomalies in the full case as well, except for Catania.
In any case, Fig. 2a confirms the results consistency, since East Sicily emerges as the most extreme side of the island.The Kruskal-Wallis test was applied to the full case, providing similar results to the annual case.Among the characterizing indicators (subsection Clusters validation procedure), md (Maximum per day) and h (Heavy rain (%)) are particularly relevant in the full case.
Fig. 2c and 2d show the md and the h heat-maps, in the full case.In contrast, Fig. 2a and 2b show the clustering results.Several similarities among the maximum values of the indicators and the anomalies can be observed.Therefore, also in the full case, the extreme stations coincide with the anomalous clusters.Moreover, the red cluster in Fig. 2b represents a cluster of extreme stations, confirmed by Fig. 2d.In fact, apart for the anomaly of Pedara, these stations retain the highest values of the h indicator.
In conclusion, the different implemented experimental settings allowed us to discover several different aspects of extreme events.Certainly, the presence of these phenomena in eastern Sicily emerges both from the annual and the full clustering, especially when the Euclidean metrics is used as a similarity measure in the AP algorithm.On the other hand, the use of Correlation metrics brings to consider Sicily composed of two different climatic areas: West side and East side, as shown in Fig.  3, where there are only 2 large clusters.Moreover, in this case no similarities between characterizing indicators and clusters are found (compare Fig. 3a and Fig. 2c, Fig. 3b and Fig. 2d).
Eventually, the clustering involving C.D of Table 1 with the use of the Euclidean metrics seems to be the most suitable setting among those tested.

Local investigation
In the local case we investigated the temporal evolution of rainfall events.In particular, similar years in the entire observed period were detected.To this aim, the AP algorithm was applied only to C.A s and C.B s , analysing one station per time.As in the geographical investigation, we chose to use both the Euclidean and the Correlation metrics.
In order to understand the most anomalous years, we counted (over stations) how many times one year appears as exception when using the Euclidean distance.cases the most anomalous years were 2015 and 2018.This means for instance that 2018 is clustered as anomalous in about 20 over 34 stations for Euclidean distance and C.A s .We see that also the 2021 counter increases after the years 2019 and 2020.Additionally, Fig. 4a and 4b also report the heavy rain (%) mean values (in blue).In this case, we fix a year y and we compute the mean of the heavy rain (%) values for all the stations that cluster the year y as anomalous, thus obtaining, for instance, that 2020 and 2021 have the highest mean values.Summarizing, in the case of C.A s and Euclidean distance, an increasing trend on anomalous years was found concerning the heavy rain mean indicator (see Fig. 4a in blue).On the other hand, in C.B s and Euclidean distance, the trend is less detectable and the highest value of heavy rain mean is measured in 2012 (see Fig. 4b in blue).

Conclusions
The main goal of this work was to introduce a clustering approach detecting extreme rainfall events occurred in Sicily, from 2009 to 2021 and to identify communities of sites with similar behaviors.
To the best of our knowledge, we are presenting for the first time in the literature, the use of multi-modal clustering analysis to detect extreme rainfall events in Sicily.With this approach, we were able to confirm and expand some preliminary observations presented in 28 , where a statistical approach has been applied to analyse rainfall trends.Specifically, in our work a clustering technique, the Affinity Propagation algorithm, was employed to confirm and discover geographical and historical rainfall changes.
In order to understand the rainfall phenomena mostly characterizing the geographical communities identified, several rainfall indicators were introduced and evaluated over the available time series.In addition, the obtained results were validated by means of the Kruskal-Wallis statistical test.
Three types of clustering analysis were conducted -full, annual and local: firstly applied to the entire high frequency time series and then applied to the weekly averaged data.The reason of this choice lied in the need of reducing the dataset size in order to minimize sensitivity to outliers.Eventually, we investigated both the Euclidean and the Correlation metrics as distance measures for the AP algorithm.

10/14
The paper presents several significant findings: • East Sicily is increasingly becoming a protagonist of extreme events, both in the full period of recordings and in the annual cases.This result is more evident choosing the Euclidean metrics in the implementation of the Affinity Propagation algorithm.
• High frequency data (C.A) with the Euclidean metrics brings to the detection of an increasing trend over years of extreme events; in contrast, C.B of weakly averaged data does not provide the same evidence.
• 2021 emerges as one of the most anomalous years in the local investigation over time.Moreover, we found from the geographical analysis that it is characterized by extreme events in the East side of the island, particularly in the cities of Catania, Siracusa and Augusta.
• Using a statistical validation approach, we found out that three indicators describe the anomalous clusters finely: the maximum per day (md), the maximum daily variation (mv) and the heavy rain percentage (h).This entails that most of the time anomalous clusters are characterized by the presence of extreme events.
• The Affinity Propagation algorithm allowed to detect anomalies, namely extreme stations, considering the full dataset.Specifically, using the Euclidean metrics, the cities of Augusta, Siracusa, and Pedara were identified clearly as anomalous at the first iteration of the algorithm.In contrast, Palazzolo Acreide and Messina have been detected at the second and third run of the algorithm, respectively.Catania does not emerge as an anomalous cluster, however, similarly to the previous sites, it presents 6 out of 13 anomalous years.
• The Euclidean metrics is sensitive to micro-climatic changes, i.e. geographically close stations are clustered in different groups.This is consistent with what can be experimentally observed since there are actually rainfall events of different frequency and magnitude a few kilometers apart.
• The Correlation metrics allowed to identify uniform clusters.This is particularly evident in the full case, in which the algorithm splits Sicily in East and West parts.
• The dataset size reduction by weekly means is successful in the geographical clusterings.In particular, it merges together anomalies and territorial clusters in a balanced way, finding appropriately extreme clusters, such as the one in eastern Sicily.
We are aware that those obtained results are no more robust than in previous before-mentioned studies, given the short observational period and therefore potentially affected by the large unforced internal/natural variability.It is clear that a wider spatial and temporal range would be needed to fully validate changes in heavy rainfall trends [60][61][62][63] .However, this paper has a methodological focus aimed to introduce an elegant approach to detect extreme events.In fact, the clustering approach is promising to interpret spatial patterns of heavy rainfall.
Further research is necessary to determine which dimensionality reduction procedure is the best for having a more precise local investigation.For instance, the total per hour datasets, the daily averaged datasets or some features refinement techniques, as for instance the Principal Component Analysis, can be advantageous to reduce the features-set dimensionality.More accurate rainfall indicators could be applied and derived, in order to entirely characterize rainfall extreme events.Additionally, in order to have more robust results, a more general analysis can be conducted, merging stations by provinces, considering a temporal clustering over the entire region, or even increasing the temporal range of investigation.In this way, the methodology introduced may potentially also yield robust findings in terms of a climate change signal.
This study could contribute significantly to the development of the decision support systems based on multi-modal and easy to use data science tools for policy makers, stakeholders, and social actors.We are currently working on the possible development of a demo based on the proposed method, to be later distributed to potential investors and researchers interested in the field.

Figure 1 .
Figure 1.Location of rainfall gauging stations in Sicily.

Figure 2 .Figure 3 .
Figure 2. Full case -Euclidean metrics.In (a) and (b) different colors represent different clusters, both in the maps and in the histograms.Square and diamond points represent results from, respectively, the second and the third iteration of the algorithm.(a) C.C.The principal cluster is reported in blue.The numbers indicate how many times the stations has been clustered as anomalous in the annual case.(b) C.D. The five main clusters are reported in green, blue, dark blue, red and yellow.(c) C.C. Maximum per day (md) heatmap.(d) C.D. Heavy rain (%) (h) heatmap.

Figure 4 .
Figure 4. Anomalous years -Euclidean metrics.C.A s (a).C.B s (b).The heavy rain mean of the year y consists of the mean of the heavy rain (%) values for all the stations that cluster the year y as anomalous.

Table 1 .
Dataset Collections.The number of considered stations is 34, except for the Single stations Collections.* 52704 for leap years.

Table 2 .
Description of the Indicators.
vViolent rain number of violent (> 50 mm) rainy hours in the time series.It was not reported in percentage since it represented very rare events.