Analyzing Space–Time Coherence in Precipitation Seasonality across Different European Climates

Seasonality is a fundamental feature of environmental systems which critically depend on the climate annual cycle. The regularity of the precipitation regime, in particular, is a basic factor to sustain equilibrium conditions. An incomplete or biased understanding of precipitation seasonality, in terms of temporal and spatial properties, could severely limit our ability to respond to climate risk, especially in areas with limited water resources or fragile ecosystems. Here, we analyze precipitation data from the Climate Hazards Group Infrared Precipitation with Stations (CHIRPS) at 0.050 resolution to study the spatial features of the precipitation seasonality across different climate zones in Central-Southern Europe during the period 1981–2018. A cluster analysis of the average annual precipitation cycle shows that seasonality under the current climate can be synthesized in the form of a progressive deformation process of the annual cycle, which starts from the northernmost areas with maximum values in summer and ends in the south, where maximum values are recorded in winter. Our analysis is useful to detect local season-dependent changes, enhancing our understanding of the geography of climate change. As an example of application to this issue, we discuss the seasonality analysis in a simulated scenario based on IPCC projections.


Introduction
Climate spatial variability around the world is determined by many factors, such as latitude, orography, influence of oceans, and direction of prevailing winds.Differences between climate zones [1,2] not only concern average values of temperature and precipitation, which are absolute key variables, but also their temporal variability.Climate seasonality, in particular, is a basic temporally coherent pattern for many environmental processes, as it influences biodiversity, ecological network structures [3], and evolutionary adaptation [4].Its effects on human functions, activities, and socio-economies [5] are crucial, too, especially in areas whose economy is strictly dependent on agriculture and exploitation of natural resources.
Many ecological studies have documented the negative effects of asynchronous climate changes, i.e., changes that are unevenly distributed throughout the year and therefore do not affect the annual regime uniformly.These season-dependent alterations, especially if they are also contrasting (with a decrease of the value of climate variables in some parts of the year and an increase in others), can dramatically limit the reactive capabilities of species and ecosystems.These studies suggest that climate-related predictive frameworks should also account for spatial and temporal irregularity to fully capture the spectrum of potential future climate change effects (see [6] for a review).In such a context, the assessment of the spatial and temporal scales involved in the many climate processes which condition environmental sustainability are of growing concern to scientists and policy makers interested in mitigation and adaptation strategies.Continuous monitoring and analysis of actual space-time climatic trajectories, therefore, appear ever more indispensable tools to obtain realistic information useful for refining models and improving decision support systems.
Although the most representative, synoptic climate variable is temperature, a specific focus on variability and change in precipitation regimes is needed to evaluate the risk of immediate or long-term detrimental effects on the environment (drought, flood, water scarcity, loss in land productivity, etc.) [7].The contrast between wet and dry regions and between wet and dry seasons, as well as the intensity and frequency of extreme rainfall events over most of the mid-latitude lands, are expected to increase [8].This will exacerbate land degradation processes, compromising extensive agricultural areas and increasing the vulnerability of local communities.
Satellite remote sensing data represent the main source to get climate information on the present climate [9].Earth observation satellite sensors are providing collective, high-quality data on climate variables allowing the observation of states and processes of the atmosphere, land, and ocean at several spatial and temporal scales [10][11][12].
The Global Climate Observing System (GCOS) has listed 26 out of 50 essential climate variables (ECVs) as significantly dependent on satellite observations [12,13].
Over the past fifty years, a range of satellite platforms carrying many different sensors has been constructed to collect atmospheric parameters used in meteorological and climatological studies.Several countries, including the US and member countries of the EU, are scheduling new satellite missions for climate observation, such as the Copernicus Sentinel-5P mission, which has been successfully launched with the main objective of performing atmospheric measurements with high spatial and temporal resolution, to be used for air quality, and climate monitoring and forecasting (https://sentinel.esa.int/web/sentinel/missions/sentinel-5p).Very interesting products, which merge satellite data with in-situ observations, have been developed to provide temporally consistent and near real-time datasets.Among these, the precipitation data from the Climate Hazards Group Infrared Precipitation with Stations (CHIRPS) [14] offer estimates over land with 0.05 0 resolution and quasi-global coverage and have proven to be useful to study environmental change, supporting, in particular, seasonal drought predictions (see [15] and reference therein).
The availability of such spatially continuous, long-term observations makes it possible to strengthen knowledge about precipitation variability within and across different climate zones, thereby enabling the identification of space-time trajectories of regional climate.Our study area lies in the central-southern part of Europe, which includes nine main climate classes ranging from ET (Polar, tundra) to Csa (Temperate, dry summer, hot summer) and Bsk (Arid, steppe, cold), according to the recently revised Köppen-Geiger classification [16].Here, we refer to this well-known bioclimatic classification scheme, which has been largely used to synthesize climates in relation to the ecosystem distribution and to estimate climate change effects on biodiversity (see, e.g., [17]).
We analyze CHIRPS time series from 1981 to 2018 focusing on precipitation seasonality, which is the temporally coherent pattern that usually characterizes the annual precipitation regime of an area.Studies on precipitation seasonality in Europe generally investigate the leading modes of intra-and inter-annual fluctuations to identify the main drivers of variability in the region (e.g., [18]).Our point of view is quite different.We want to study coherence between seasonality features observed across the study area by identifying homogeneous (synchronous) regions and their characteristic annual regimes.This analysis is particularly interesting for environmental researchers who need a simple and effective strategy to integrate precipitation regimes in their environmental research and to evaluate sustainability in a changing climate.Our study therefore focuses on normalized monthly precipitation values (annual averages and amplitudes do not affect coherence) and requires the use of a methodology to cluster together areas with similar timing.Here, we use the K-means clustering algorithm to regionalize seasonality and to define a discrete set of seasonal curves which describe precipitation regimes in the area.
We propose a simple procedural scheme for the continuous analysis of precipitation regimes in climate mosaics under a changing climate, able to detect peculiar transformations occurring within single climate patches and to provide a collective picture of climate change in a holistic approach to the problem.The CHIRPS dataset is still too short to demonstrate the performance of our approach in this context.Thus, we resort to a simulated scenario based on IPCC projections for explaining how the methodology is able to highlight the areas affected by change and the effects of change on the spatial distribution of the seasonal timing.

Study Area
The study area is located in the central-southern part of Europe (Figure 1a) encompassing a wide range of terrestrial and marine ecosystems mainly dominated by Mediterranean environments.However, large areas are characterized by continental and tundra landscapes and there are some limited patches containing semi-arid (steppe) ecosystems.We propose a simple procedural scheme for the continuous analysis of precipitation regimes in climate mosaics under a changing climate, able to detect peculiar transformations occurring within single climate patches and to provide a collective picture of climate change in a holistic approach to the problem.The CHIRPS dataset is still too short to demonstrate the performance of our approach in this context.Thus, we resort to a simulated scenario based on IPCC projections for explaining how the methodology is able to highlight the areas affected by change and the effects of change on the spatial distribution of the seasonal timing.

Study Area
The study area is located in the central-southern part of Europe (Figure 1a) encompassing a wide range of terrestrial and marine ecosystems mainly dominated by Mediterranean environments.However, large areas are characterized by continental and tundra landscapes and there are some limited patches containing semi-arid (steppe) ecosystems.Here, we observe different ecosystems placed across an increasing gradient of anthropic influence: from quite undisturbed areas (grassland, outcropping rocks, and snow or ice-dominated habitats in alpine ecosystems) to forestlands (spread throughout the study area) and barren or sparsely vegetated lands (typical of semi-arid landscapes) to agricultural and artificial areas with a different degree of urban density [19].Orography is rather varied and stretches from the characteristic Mediterranean coasts and the huge plains to the extensive European mountain ranges (Figure 1b).
According to the Köppen-Geiger classification ( [16], see Figure 2) the study area hosts nine main climate classes ranging from ET (Polar, tundra) to Bsk (Arid, steppe) passing through the prevalent Mediterranean (temperate) Cs and Cf and the cold (Df) climates.Here, we observe different ecosystems placed across an increasing gradient of anthropic influence: from quite undisturbed areas (grassland, outcropping rocks, and snow or ice-dominated habitats in alpine ecosystems) to forestlands (spread throughout the study area) and barren or sparsely vegetated lands (typical of semi-arid landscapes) to agricultural and artificial areas with a different degree of urban density [19].Orography is rather varied and stretches from the characteristic Mediterranean coasts and the huge plains to the extensive European mountain ranges (Figure 1b).
According to the Köppen-Geiger classification ( [16], see Figure 2) the study area hosts nine main climate classes ranging from ET (Polar, tundra) to Bsk (Arid, steppe) passing through the prevalent Mediterranean (temperate) Cs and Cf and the cold (Df) climates.The mean annual temperature in the area, reported in Figure 3a, was estimated from latitude, longitude, and elevation according to a regressive model appositely devised for Europe and parameterized with the support of geostatistical tools.The model parameters were estimated from temperature data measured at 387 European meteorological stations mostly provided by the European Climate Assessment & Dataset (ECA&D) [20].As for annual precipitation averaged over the observational period (Figure 3b), orography and latitude are the main driving factors for the precipitation amount.There is a rather evident inverse relationship between the two variables, which is particularly worrying in areas where summer months are contemporaneously the warmest and the driest.This remarkable physiographic and climatic heterogeneity is associated with a variety of economic development levels, as revealed by the EU's GDP (Gross Domestic Product) per capita [21], including underdeveloped regions mostly located in Southern Italy and in the Balkan Peninsula and The mean annual temperature in the area, reported in Figure 3a, was estimated from latitude, longitude, and elevation according to a regressive model appositely devised for Europe and parameterized with the support of geostatistical tools.The model parameters were estimated from temperature data measured at 387 European meteorological stations mostly provided by the European Climate Assessment & Dataset (ECA&D) [20].As for annual precipitation averaged over the observational period (Figure 3b), orography and latitude are the main driving factors for the precipitation amount.There is a rather evident inverse relationship between the two variables, which is particularly worrying in areas where summer months are contemporaneously the warmest and the driest.The mean annual temperature in the area, reported in Figure 3a, was estimated from latitude, longitude, and elevation according to a regressive model appositely devised for Europe and parameterized with the support of geostatistical tools.The model parameters were estimated from temperature data measured at 387 European meteorological stations mostly provided by the European Climate Assessment & Dataset (ECA&D) [20].As for annual precipitation averaged over the observational period (Figure 3b), orography and latitude are the main driving factors for the precipitation amount.There is a rather evident inverse relationship between the two variables, which is particularly worrying in areas where summer months are contemporaneously the warmest and the driest.This remarkable physiographic and climatic heterogeneity is associated with a variety of economic development levels, as revealed by the EU's GDP (Gross Domestic Product) per capita [21], including underdeveloped regions mostly located in Southern Italy and in the Balkan Peninsula and This remarkable physiographic and climatic heterogeneity is associated with a variety of economic development levels, as revealed by the EU's GDP (Gross Domestic Product) per capita [21], including underdeveloped regions mostly located in Southern Italy and in the Balkan Peninsula and rich areas extending prevalently in the so-called Mitteleuropa (Northern Italy, Austria, Switzerland, Southern Germany).Some of these areas, generally falling in the Mediterranean basin, are vulnerable to land degradation because of unfavorable biophysical factors (e.g., thin soil, sparse vegetation, etc.) and bad management practices [22][23][24], which are accompanied by even more frequent extreme events such as drought or frost [25][26][27][28][29][30].The intensity and frequency of these extreme events are likely to increase in response to climate change, further impacting economic activities of underdeveloped regions and widening the social disparities within Europe.On the other hand, risk of floods and/or landslides triggered by extreme rainfall events are particularly noticeable in the central-northern part of the study area, entailing loss of human life and severe damage to communities [31][32][33][34][35][36][37][38].

CHIRPS Data
Climate Hazards Group Infrared Precipitation with Stations (CHIRPS) is a high-resolution (0.05 • ) satellite and observation-based precipitation estimate for the quasi-global coverage (50  [43,44]; • the in situ precipitation observations obtained from a variety of sources including national and regional meteorological services.
The CHIRPS procedure is detailed in [14].In brief, CHIRPS dataset is produced in two main steps: 1. rainfall from satellite data is estimated as the percentage of time during the pentad that the IR observation is lower than 235

Cluster Analysis of Average Seasonal Patterns by K-means
Here, we want to briefly recall that K-means [45] is a simple clustering tool which can be used to subdivide data into a number of non-overlapping clusters K (fixed a priori) according to given features.The algorithm works iteratively to group data according to feature similarity.It identifies K centroids (collections of feature values which define the clusters), and then allocates every data point to the nearest cluster by minimizing the intra-cluster distance (usually Euclidean distance is used) while keeping different clusters as distant as possible.The cost function to be minimized is: that is the sum of squared distances of each data point x i from the centroid µ k of its assigned cluster c k .
In our work, K-means is used to "quantize" the spatial "continuum" of the seasonal variability in a discrete set of centroids thereby reducing the number of spatial degrees of freedom.
The precipitation field P(p, t), which represents the monthly precipitation amount estimated at the pixel p and time t, can be properly represented by splitting the time index t into two time indices: m, which represents the month of the year (m = 1, ..., 12), and y, which is the year (y = 1, . . ., 38).
The first step of the analysis concerns the estimation of the monthly mean values of P(p, m, y) : for each p and m.The sequence (P(p, m), m = 1, . . ., 12) represents therefore the seasonal curve estimated for the pixel p.In order to investigate the degree of synchronism of the seasonal patterns observed in different pixels, we estimated the 38-year monthly means and normalized them to obtain seasonal curves with zero mean and unitary variance P n (p, m), m = 1, . . ., 12: We used these normalized values to define a 12-dimensional feature space where the distance between two pixels p 1 and p 2 is the Euclidean distance between the vectors P n (p 1 , 1), . . ., P n (p 1 , 12) and P n (p 2 , 1), . . ., P n (p 2 , 12) , which synthesize their seasonality.This filtering allows us to apply K-means by using the Euclidean distance to detect areas with the most correlated seasonality.In addition, normalization minimizes most of the effects of factors, such as elevation, which act on the absolute values of rainfall thus enabling us to better focus on the annual precipitation timings.We can consequently pick up spatial aggregations that are based on the seasonal contrast in precipitation more than on the degree of dryness and wetness in absolute terms.
In order to measure similarity and separation in the clustering, we estimate the Silhouette [46], which is an index whose values range in the interval (-1, 1).A value of +1 indicates that the given pixel is in the core of the cluster assigned to it.On the contrary, a value of −1 indicates that the pixel is closer to its neighboring clusters than to its cluster.Values close to 0 mean that the given pixel is located on the boundary between two clusters.The performance of K-means for different values of K was evaluated according to the mean value (SI) of the Silhouette on the whole clustering.

Cluster Analysis of Sample Time Distributions
The first step of the analysis focuses on the statistical distributions of the monthly values of precipitation estimated per pixel for the total 1981-2018 sample and is for contextualizing the analysis of precipitation seasonality.The three descriptive quartiles, [Q1(p), Median(p), Q3(p)], were estimated per pixel and exploited as similarity parameters to cluster areas with similar distributions.Of course, we did not expect well separated clusters due to the spatial continuity of precipitation values against elevation.Silhouette analysis, in fact, did not indicate a unique optimal K.A few clusters (K = 5, 6) produce almost equivalent results (SI ~0.75), while the performance drops greatly for higher values.We decided to group pixels into six clusters (Figure 4 and Table 1) because this partition better accounts for the continuity of the altitude variability at the 0.05 • resolution of the data and for the distance from the coastline.This clustering has no implication on our final results but is useful to show that synchronized areas, which are studied by analyzing normalized values, do not necessarily match areas with a similar degree of dryness/wetness.Median values progressively decrease in the passage from C1 to C6 (Figure 4b); the range (Q3-Q1) computed for the cluster centroids, which includes 50% of the precipitation values, decreases accordingly.This means that a large variability is observed where median values are high, especially in mountainous areas, whereas narrower distributions characterize drier areas.It is interesting to note that clusters are nested inside each other; an ideal line that would connect C1 to C6 should cross all the other ones.Along the east side of the Adriatic Sea, in particular, clusters boundaries follow the coastline and, going inland, we cross all the clusters from C6 up to C2 (C1 in the area of the Dinaric Alps).The inverse succession from C2 to C6 appears if we further penetrate the continent until the Pannonian Plain.C1, C2, and C3, which include more rainy areas, mostly collect mountain climates (ET, Dfa, and Dfc); on the other hand, in C5 and C6, where precipitation are more scarce, temperate and arid climates dominate (Csa, Csb, Cfa, Cfb, and Bsk).C4 appears to be an intermediate class and is the most heterogeneous, as it is composed by pixels from all the climate zones with the exception of the extreme ones.

Cluster Analysis of Average Seasonal Patterns
The analysis of the spatial coherence of the seasonal timing was performed on normalized data (zero mean and unitary variance).The optimal number of clusters was obtained estimating the average value of the mean Silhouette SI for different values of K and selecting the K value which  Median values progressively decrease in the passage from C1 to C6 (Figure 4b); the range (Q3-Q1) computed for the cluster centroids, which includes 50% of the precipitation values, decreases accordingly.This means that a large variability is observed where median values are high, especially in mountainous areas, whereas narrower distributions characterize drier areas.It is interesting to note that clusters are nested inside each other; an ideal line that would connect C1 to C6 should cross all the other ones.Along the east side of the Adriatic Sea, in particular, clusters boundaries follow the coastline and, going inland, we cross all the clusters from C6 up to C2 (C1 in the area of the Dinaric Alps).The inverse succession from C2 to C6 appears if we further penetrate the continent until the Pannonian Plain.C1, C2, and C3, which include more rainy areas, mostly collect mountain climates (ET, Dfa, and Dfc); on the other hand, in C5 and C6, where precipitation are more scarce, temperate and arid climates dominate (Csa, Csb, Cfa, Cfb, and Bsk).C4 appears to be an intermediate class and is the most heterogeneous, as it is composed by pixels from all the climate zones with the exception of the extreme ones.

Cluster Analysis of Average Seasonal Patterns
The analysis of the spatial coherence of the seasonal timing was performed on normalized data (zero mean and unitary variance).The optimal number of clusters was obtained estimating the average value of the mean Silhouette SI for different values of K and selecting the K value which produces the highest SI.The best result (SI 0.8) was obtained for K = 2.The spatial partition so obtained (Figure 5a) roughly separates continental areas (curves with maximum values in summer) from southern areas more directly affected by the influence of the Mediterranean Sea (curves with minimum values in that period), as shown in Figure 5b.If we exclude the effect of orography, we find a reasonable agreement with the Köppen-Geiger Df and Cf macro-groups.produces the highest SI.The best result (SI ≅ 0.8) was obtained for K = 2.The spatial partition so obtained (Figure 5a) roughly separates continental areas (curves with maximum values in summer) from southern areas more directly affected by the influence of the Mediterranean Sea (curves with minimum values in that period), as shown in Figure 5b.If we exclude the effect of orography, we find a reasonable agreement with the Köppen-Geiger Df and Cf macro-groups.Starting from this hard partition, we wanted to identify the spatial process that is associated with the transformation of the curves of Figure 5b in each other.The strategy adopted here was to evalute the distance d between seasonality at the pixel level within a cluster and the centroid of the other one.High values are expected in areas that have a high degree of membership in their cluster whereas ever lower values should identify areas that are ever more similar to those of the other cluster.
The results of these estimations are reported in Figure 6a.Seasonal variability within clusters is not randomly arranged but is organized in swaths that gradually change in the North-South direction.The core areas of the two clusters (the reddest ones in Figure 6a) are linked by a sequence of regions with in-between features.Two complex regions in the west and east part of the continental area account for more complex zones.In order to segment the investigated area in the regions qualitatively identified in Figure 6a, we defined the boundaries of nine areas with the help of Kmeans applied with K = 9.
Looking at the similarity of seasonal timings, we discover spatial structures that mostly vary with latitude (over the entire study area) and with longitude but are limited to the continental interior (Figure 6b and Table 2).Starting from this hard partition, we wanted to identify the spatial process that is associated with transformation of the curves of Figure 5b in each other.The strategy adopted here was to evalute the distance d between seasonality at the pixel level within a cluster and the centroid of the other one.High values are expected in areas that have a high degree of membership in their cluster whereas ever lower values should identify areas that are ever more similar to those of the other cluster.
The results of these estimations are reported in Figure 6a.Seasonal variability within clusters is not randomly arranged but is organized in swaths that gradually change in the North-South direction.The core areas of the two clusters (the reddest ones in Figure 6a) are linked by a sequence of regions with in-between features.Two complex regions in the west and east part of the continental area account for more complex zones.In order to segment the investigated area in the regions qualitatively identified in Figure 6a, we defined the boundaries of nine areas with the help of K-means applied with K = 9.
Looking at the similarity of seasonal timings, we discover spatial structures that mostly vary with latitude (over the entire study area) and with longitude but are limited to the continental interior (Figure 6b and Table 2).produces the highest SI.The best result (SI ≅ 0.8) was obtained for K = 2.The spatial partition so obtained (Figure 5a) roughly separates continental areas (curves with maximum values in summer) from southern areas more directly affected by the influence of the Mediterranean Sea (curves with minimum values in that period), as shown in Figure 5b.If we exclude the effect of orography, we find a reasonable agreement with the Köppen-Geiger Df and Cf macro-groups.Starting from this hard partition, we wanted to identify the spatial process that is associated with the transformation of the curves of Figure 5b in each other.The strategy adopted here was to evalute the distance d between seasonality at the pixel level within a cluster and the centroid of the other one.High values are expected in areas that have a high degree of membership in their cluster whereas ever lower values should identify areas that are ever more similar to those of the other cluster.
The results of these estimations are reported in Figure 6a.Seasonal variability within clusters is not randomly arranged but is organized in swaths that gradually change in the North-South direction.The core areas of the two clusters (the reddest ones in Figure 6a) are linked by a sequence of regions with in-between features.Two complex regions in the west and east part of the continental area account for more complex zones.In order to segment the investigated area in the regions qualitatively identified in Figure 6a, we defined the boundaries of nine areas with the help of Kmeans applied with K = 9.
Looking at the similarity of seasonal timings, we discover spatial structures that mostly vary with latitude (over the entire study area) and with longitude but are limited to the continental interior (Figure 6b and Table 2).The macro-groups reported in Figure 5a are substantially subdivided in two separated sets of regions, from R1 to R6 and from R7 to R9 respectively.Therefore, this segmentation is in agreement with the optimal one (Figure 5) but adds further details to it.
In Figure 6b, colors varying progressively from blue to red, green clusters excluded, allow to follow an ideal path from North to South along which the seasonal patterns vary, as illustrated in Figure 7a.It is possible to note that all the seasonal curves along this path can be obtained by small deformation from the adjacent ones.R1, which includes the coldest climates, is characterized by a net dominance of the simple annual frequency (see the power spectrum for R1) with maximum values in the summer.Almost the same is true for R2, even if the rainy phase lengthens and a small decrease is observed in July.In R3, the summer minimum becomes more evident and the second harmonic appears in the power spectrum because of the appearance of two separate maxima, in spring and autumn.The seasonal pattern of R7 is very similar but more asymmetric than R3, as the spring maximum decreases; the second harmonic dominates over the annual one.In R8, the summer minimum expands and the annual harmonic returns to be dominant.In R9, finally, we find a seasonal curve that is in an almost perfect counter-phase with the pattern observed for R1; in this case, the power spectrum is dominated by the first harmonic with winter maxima.Note that although the analysis focused on normalized values, in Figure 7 real estimated values are reported.Therefore, the centroids are the monthly mean values characterizing the different clusters.This figure shows that our clustering provides a rich library of curves, which describe seasonality within the current climate and can be used as a baseline to evaluate future climate change in terms of passages from a cluster to another one.In this figure, error bands are also reported as a measure of the degree of representativeness of the centroids.These errors were measured by estimating the standard deviation per month of each pixel value in its region of membership.
Figure 7b illustrates the West-to-East sequence.These clusters show the most uniform profiles, especially in the westernmost part.In R5 and R6, there is no sign of seasonality properly defined and precipitation levels are high along all the year and their centroids are rather close.We did not group them in a unique region as they belong to C1 k = 2 and C2 k = 2 , respectively (see Figure 5a).This lack of a proper seasonality is picked up by the power spectrum, which shows contributions from all the frequencies.This behaviour significantly changes at the boundary of R6 where this nearly uniform distribution seems to be split in the three R1, R2, and R3 seasonal patterns.In Figure 7b we have reported the already discussed two-season pattern of R3 which exhibits the higher analogy with R4, in the easternmost part of the region, which again shows an increase in irregularity with only a slight maximum in the late spring.Within clusters R7, R8, and R9, summer is the driest season in the warmest period of the year.In practice, all the Mediterranean areas show this feature.About 90% of the minimum values fall in July and 5% in August.If temperature would further increase concomitantly with drought frequency and persistence, these areas would be very vulnerable to climate risk.Already today the mean summer values in the southernmost and insular areas are very low (Figure 8), especially along the coastlines, where often minimum values are lower than 10 mm.

Detection of Seasonally Asynchronous Climate Change in a Simulation-Based Scenario
Results discussed in Section 3.2 allowed us to synthesize the spatial variability of the annual precipitation regime in a set of temporal curves characterizing regions with synchronous seasonality under the current climate.Any future change, affecting just some areas and just some parts of the year, would destroy synchronism between the changed areas and their class of membership, thereby determining a new partition.In order to test the potential of our approach to detect such asynchronous change, we used simulated data, as satellite observations assuring spatially continuous information are still too short for climate change analyses.
According to the IPCC climate projections for the study area [48], the precipitation amount should decrease in some areas of Southern Europe in the three-decade period 2071-2100 with respect to 1961-1990.Different Regional Climate Models (RCM) predict significant changes in the precipitation regimes characterizing the areas of Calabria and Sicily (Southern Italy, cluster R9).In particular, the RCA3-BCM (Bergen Climate Model), an RCM including a full description of the atmosphere and its interaction with the land surface, predicts a statistically significant precipitation decrease (20-30%) in winter values for both the regions and in summer values for Calabria [49].This is an example of geographically confined and seasonally asynchronous change, as different months are differently affected by change.Thus, we modified the CHIRPS data by decreasing the current winter values in Calabria and Sicily as well as the summer values in Calabria by 25%.Then we repeated the clustering described in Section 4.2, selecting the current precipitation curves as initial centroids for the K-means iterative algorithm.The first interesting result is that the algorithm did not modify these initial centroids, which still represent characteristic seasonal curves in the area.In the

Detection of Seasonally Asynchronous Climate Change in a Simulation-Based Scenario
Results discussed in Section 3.2 allowed us to synthesize the spatial variability of the annual precipitation regime in a set of temporal curves characterizing regions with synchronous seasonality under the current climate.Any future change, affecting just some areas and just some parts of the year, would destroy synchronism between the changed areas and their class of membership, thereby determining a new partition.In order to test the potential of our approach to detect such asynchronous change, we used simulated data, as satellite observations assuring spatially continuous information are still too short for climate change analyses.
According to the IPCC climate projections for the study area [48], the precipitation amount should decrease in some areas of Southern Europe in the three-decade period 2071-2100 with respect to 1961-1990.Different Regional Climate Models (RCM) predict significant changes in the precipitation regimes characterizing the areas of Calabria and Sicily (Southern Italy, cluster R9).In particular, the RCA3-BCM (Bergen Climate Model), an RCM including a full description of the atmosphere and its interaction with the land surface, predicts a statistically significant precipitation decrease (20-30%) in winter values for both the regions and in summer values for Calabria [49].This is an example of geographically confined and seasonally asynchronous change, as different months are differently affected by change.Thus, we modified the CHIRPS data by decreasing the current winter values in Calabria and Sicily as well as the summer values in Calabria by 25%.Then we repeated the clustering described in Section 4.2, selecting the current precipitation curves as initial centroids for the K-means iterative algorithm.The first interesting result is that the algorithm did not modify these initial centroids, which still represent characteristic seasonal curves in the area.In the light of this result, we did not increase the number of clusters.Figure 9 shows that both the regions affected by change move from R9 to R8 as seasonality in these areas becomes more asymmetric due to the decrease in winter values.Such a kind of change breaks the continuity of cluster R9 slightly redesigning the geography of the precipitation seasonal timing in Southern Italy, where R9 appears to be confined in isolated aggregations surrounded by R8.Changes in summer values predicted for Calabria do not produce effects as this decrease does not modify the seasonal timing but only the magnitude of minimum values.This change in summer values, however, would trivially modify the map of the most climatically vulnerable areas illustrated in Figure 8, showing an increase of the climatic risk in Calabria due to the decrease of precipitation values in the warmest period in the year.
light of this result, we did not increase the number of clusters.Figure 9 shows that both the regions affected by change move from R9 to R8 as seasonality in these areas becomes more asymmetric due to the decrease in winter values.Such a kind of change breaks the continuity of cluster R9 slightly redesigning the geography of the precipitation seasonal timing in Southern Italy, where R9 appears to be confined in isolated aggregations surrounded by R8.Changes in summer values predicted for Calabria do not produce effects as this decrease does not modify the seasonal timing but only the magnitude of minimum values.This change in summer values, however, would trivially modify the map of the most climatically vulnerable areas illustrated in Figure 8, showing an increase of the climatic risk in Calabria due to the decrease of precipitation values in the warmest period in the year.In summary, seasonally asynchronous change is detected by our procedure in terms of pixel transfer between clusters, as the seasonality library described in Figure 7 includes rather fine details that are able to account for most of the possible distortions introduced by climate change, if sufficiently significant.

Conclusions
This study examined the seasonality and spatial coherence of precipitation by analyzing CHIRPS multi-source data from 1981 to 2018 at 0.05 0 resolution in central-southern Europe.Such a resolution appeared to be suited to studying climate at the limit between regional and local scale and is potentially useful for end-users involved in the land management of fragmented areas, which cannot straightforwardly benefit from the results from global and regional climate models.
As far as we know, this is the first work which explicitly synthesizes the spatial variability of the precipitation annual regime in a collection of seasonal curves characterizing seasonally synchronous regions.The overall picture provided by our investigations accounts for a high degree of coherence in precipitation seasonality when we move across the study area.At the 0.05 0 climatological resolution and monthly sampling, the mean annual precipitation regime estimated over 30+ years shows rather homogeneous features over large areas.If aggregated in few homogeneous clusters, the annual curves display very interesting features.Continental inner clusters are the rainiest but only in In summary, seasonally asynchronous change is detected by our procedure in terms of pixel transfer between clusters, as the seasonality library described in Figure 7 includes rather fine details that are able to account for most of the possible distortions introduced by climate change, if sufficiently significant.

Conclusions
This study examined the seasonality and spatial coherence of precipitation by analyzing CHIRPS multi-source data from 1981 to 2018 at 0.05 0 resolution in central-southern Europe.Such a resolution appeared to be suited to studying climate at the limit between regional and local scale and is potentially useful for end-users involved in the land management of fragmented areas, which cannot straightforwardly benefit from the results from global and regional climate models.
As far as we know, this is the first work which explicitly synthesizes the spatial variability of the precipitation annual regime in a collection of seasonal curves characterizing seasonally synchronous regions.The overall picture provided by our investigations accounts for a high degree of coherence in precipitation seasonality when we move across the study area.At the 0.05 0 climatological resolution and monthly sampling, the mean annual precipitation regime estimated over 30+ years shows rather homogeneous features over large areas.If aggregated in few homogeneous clusters, the annual curves display very interesting features.Continental inner clusters are the rainiest but only in the west areas is the month-to-month variability the minimum.In all the other cases, well-defined cycles are distinguishable.Along a North-to-South path starting from the central Alps, these cycles can be obtained by means of a deformation process which progressively transforms the alpine seasonality in the counter-phase cycle characterizing the southernmost Mediterranean areas.In such a picture, orography appears to have a minor role, acting more on mean precipitation levels and ranges than on timing.

Figure 2 .
Figure 2. Climate map according to the Köppen-Geiger classification.

Figure 3 .
Figure 3. Annual values for: (a) simulated mean temperature and (b) total precipitation averaged over the period 1981-2018.Grey pixels indicate areas where CHIRPS data were not analyzed because of the presence of artefacts.

Figure 2 .
Figure 2. Climate map according to the Köppen-Geiger classification.

Figure 2 .
Figure 2. Climate map according to the Köppen-Geiger classification.

Figure 3 .
Figure 3. Annual values for: (a) simulated mean temperature and (b) total precipitation averaged over the period 1981-2018.Grey pixels indicate areas where CHIRPS data were not analyzed because of the presence of artefacts.

Figure 3 .
Figure 3. Annual values for: (a) simulated mean temperature and (b) total precipitation averaged over the period 1981-2018.Grey pixels indicate areas where CHIRPS data were not analyzed because of the presence of artefacts.

Figure 4 .
Figure 4. Distribution-based clustering: (a) cluster map; (b) descriptive triplets [Q1, Median, Q3] (centroids) estimated on the total observational period; coloured points indicate the median values whereas the error bars indicate Q1 and Q3 respectively.(a) has to be read according to the colour/cluster association rule in (a); grey pixels mask CHIRPS artefacts, which will be corrected in the next v2.1 version (Pete Peterson, private communication, 2019).

Table 1 .
Clusters based on the values of triplets [Q1, Median, Q3] estimated over the period 1981-2018 (Figure 4) with the respective number of pixels.

Figure 4 .
Figure 4. Distribution-based clustering: (a) cluster map; (b) descriptive triplets [Q1, Median, Q3] (centroids) estimated on the total observational period; coloured points indicate the median values whereas the error bars indicate Q1 and Q3 respectively.(a) has to be read according to the colour/cluster association rule in (a); grey pixels mask CHIRPS artefacts, which will be corrected in the next v2.1 version (Pete Peterson, private communication, 2019).

Figure 5 .
Figure 5. Clusters (a) and centroids (b) obtained for the optimal clustering with K = 2.Note that real values are reported in (b).

Figure 5 .
Figure 5. Clusters (a) and centroids (b) obtained for the optimal clustering with K = 2.Note that real values are reported in (b).

Figure 5 .
Figure 5. Clusters (a) and centroids (b) obtained for the optimal clustering with K = 2.Note that real values are reported in (b).

Figure 6 .
Figure 6.Seasonal clusters: (a) distances d between each pixel in C1 k = 2 from the centroid of C2 k = 2 and vice versa; (b) results from the K-means algorithm applied with K = 9.
Remote Sens. 2019, 11, x FOR PEER REVIEW 10 of 15 climate risk.Already today the mean summer values in the southernmost and insular areas are very low (Figure 8), especially along the coastlines, where often minimum values are lower than 10 mm.

Figure 7 .
Figure 7. Seasonal curves (centroids) for the regions in Figure 6b: (a) North-to-South sequence; (b) West-to-East sequence.Dotted areas sign the ±σ(m) band around the centroids, where ±σ(m) is the standard deviation of the monthly means within the given cluster.Plots include the variance power spectra S(f ) [47] estimated for the seasonal curves and normalized to S(1).The height of the bars represents the variance explained by each frequency (cycles per year) relative to the variance explained by annual harmonic.

Figure 7 .
Figure 7. Seasonal curves (centroids) for the regions in Figure 6b: (a) North-to-South sequence; (b) West-to-East sequence.Dotted areas sign the () band around the centroids, where () is the standard deviation of the monthly means within the given cluster.Plots include the variance power spectra S(f) [47] estimated for the seasonal curves and normalized to S(1).The height of the bars represents the variance explained by each frequency (cycles per year) relative to the variance explained by annual harmonic.

Figure 8 .
Figure 8. Summer (July and August) precipitation mean values for the only pixels in the clusters R7, R8, and R9, which exhibit the annual minimum in summer.

Figure 8 .
Figure 8. Summer (July and August) precipitation mean values for the only pixels in the clusters R7, R8, and R9, which exhibit the annual minimum in summer.

Figure 9 .
Figure 9. Analysis of precipitation projections: seasonal clusters, update of Figure 6.Hatched areas identify regions moving from R9 to R8.

Figure 9 .
Figure 9. Analysis of precipitation projections: seasonal clusters, update of Figure 6.Hatched areas identify regions moving from R9 to R8.

Table 1 .
Clusters based on the values of triplets [Q1, Median, Q3] estimated over the period 1981-2018 (Figure 4) with the respective number of pixels.

Table 2 .
Seasonally homogeneous regions with the respective number of pixels.