Daily Precipitation and Temperature Extremes in Southern Italy (Calabria Region)

Abstract

We apply extreme value theory (EVT) to study the daily precipitation and temperature extremes in the Calabria region (southern Italy) mainly considering a long-term observational dataset (1990–2020) and also investigating the possible use of the ERA5 (ECMWF Reanalysis v5) fields. The efficiency of the EVT applied on the available observational dataset is first assessed—both through a punctual statistical analysis and return-level maps. Two different EVT methods are adopted, namely the peak-over-threshold (POT) approach for the precipitation and the block-maxima (BM) approach for the temperature. The proposed methodologies appear to be suitable for describing daily extremes both in quantitative terms, considering the punctual analysis in specific points, and in terms of the most affected areas by extreme values, considering the return-level maps. Conversely, the analysis conducted using the reanalysis fields for the same time period highlights the limitations of using these fields for a correct quantitative reconstruction of the extremes while showing a certain consistency regarding the areas most affected by extreme events. By applying the methodology on the observed dataset but focusing on return periods of 50 and 100 years, an increasing trend of daily extreme rainfall and temperature over the whole region emerges, with specific areas more affected by these events; in particular, rainfall values up to 500 mm/day are predicted in the southeastern part of Calabria for the 50-year-return period, and maximum daily temperatures up to 40 °C are expected in the next 100 years, mainly in the western and southern parts of the region. These results offer a useful perspective for evaluating the exacerbation of future extreme weather events possibly linked to climate change effects.

1. Introduction

The latest IPCC (Intergovernmental Panel on Climatic Change) report published in 2022 (IPCC Sixth Assessment Report (AR6)) [1] describes an unequivocal increase of the Earth’s surface temperature over the next decades. Particular attention in the report was given to extreme weather and climate phenomena, whose increased frequency and intensity have the potential to cause widespread and pervasive effects on ecosystems and infrastructure [1]. Furthermore, recent works show that significant information about multi-decade climate change can be obtained by quantifying shifts in the overall probability distribution of daily weather conditions—that is, by using the study of extreme events as a magnifying lens for understanding climate change [2,3,4,5,6,7].

The variation in both frequency and intensity of extreme events has a strong spatial dependence, [8,9]. The intensification of the hydrological cycle due to global warming in the Mediterranean area [10] has raised the need to conduct several studies based on future climate projections. These studies mainly reveal that the occurrence of extreme precipitation in a (future) warmer climate can possibly increase [11], even if the mean (light/medium) precipitation amount is predicted to decrease in various zones [12,13]. Similar considerations can be applied for temperature extremes, clearly more linked to the expected global warming. Heatwaves in the Mediterranean Basin have been assessed to occur more frequently in recent decades [14,15].

Based on regional climate model simulations under the more pessimistic greenhouse gas emission scenario (RCP8.5), the authors found an increase of heatwaves over the Mediterranean for the future and that the intensity of these events will be strongest in the southern and eastern part of the Basin. The spatial dependence of extreme events reveals the main limitation of a statistical approach to extreme climate events, namely the availability and quality of meteorological data. In this work, temperature and precipitation data for a period longer than 30 years were used for the first time in southern Italy—the Calabria region in particular.

The present study focuses on the Calabria region (southern Italy), in the central Mediterranean, an area particularly prone to heavy rain events and heatwaves. The region has a complex orography with several mountain chains and a marked land–sea contrast that leads to a coexistence of different local atmospheric circulation regimes. Several papers have dealt with the topic of severe weather events in the Calabria region—in particular, heavy rains—all adopting statistical/climatological approaches [16,17,18] and through in-depth analyses of remarkable extreme episodes [19,20,21,22] in order to study their characteristics and dynamics.

The study of [16], in particular, showed a first exploratory result on a 30-year (1978–2007) homogeneous precipitation database for the Calabria region (daily precipitation/88 rain gauges), to assess the key roles of the orography and the sea as well as the seasonal dependence of rainfall unequivocally linked to the synoptic scale conditions. As a general deduction, despite yearly precipitation being larger on the west side of the region, the most intense rainstorms are more frequent on the east side.

In continuation of a work [23] that provided a first classification of atmospheric patterns for the Calabria region, a recent study [18] classified the main precipitation systems through the analysis of selected heavy rainfall events, considering a high-resolution rain gauge network. This work also assessed the relationships between the selected events and the main synoptic atmospheric patterns derived by the ERA5 Reanalysis dataset.

The works that have dealt with (extreme) temperatures in Calabria are few. Among them, considering a long-time period of monthly mean values and extreme daily temperatures in southern Italy, [24,25] revealed a positive trend in spring and summer and a negative trend in the autumn–winter period; specifically, regarding the extreme temperatures, the authors observed an increase in the frequency and intensity of the highest temperatures and some negative trends for the lowest ones, i.e., a major (minor) probability of heatwaves (cold extremes) throughout the years.

Although climatological studies are mainly performed using global or regional climate models, a common approach to study extreme weather events is the adopting of statistical procedures, such as the statistical extreme value theory (EVT), on different datasets. The EVT has emerged as one of the most important statistical theories for climate observations and for understanding the outputs of numerical models. A challenge already anticipated by Wigley [26] and later taken up in a commentary by Coley [27] indicates the great potential for applications of extreme value theory on climate change.

One of the most important results of EVT, which differentiates it from other statistical approaches, is its ability to estimate the distribution of extreme values using the asymptotic argument. The origins of the asymptotic characterisation of the maximum sample go back to Fisher and Tippett [28]. The first formalisation of a model of extreme events was proposed by Jenkinson [29,30], showing that there are only three families of possible limit laws for the distribution of extremes, namely the Gumbel, Frechet and Weibull distributions. These three distributions can be expressed in a single distribution function known as the “generalized extreme value” distribution (GEV) [31].

EVT is currently used in several fields and particularly in the analysis of environmental extreme events, e.g., extreme value theory provides a solid theoretical basis for the statistical modelling of extreme hydrological events due to its many applications on floods [32,33,34]. Several works have used EVT to investigate atmospheric phenomena in order to understand climate extremes [35,36,37], hydrology [38,39] and oceanography [40,41]. There are also several works using EVT on rainfall datasets [42,43,44,45] and on temperature extremes [46,47,48,49].

Observational datasets represent the best data source for performing reliable statistical studies; at the same time, the availability of long-term data series is often a problem in several zones, as there are still large areas not sufficiently covered by measurements. In this context, the use of regularly gridded data, such as those derived from modelling products (reanalysis), would help as they are not affected by problems of spatial/temporal availability. Despite this, previous works demonstrated the limited skill of the reanalysis in correctly reproducing extreme values when compared with other datasets (observations in particular).

Ref. [50] considered ERA5 hourly precipitation over Europe to study extreme values and found that reanalysis represents a good reference for general mean statistics (e.g., spatial patterns of annual precipitation and multi-year cycles of monthly precipitation); however, the coarser resolution of the dataset tends to generate a smoothing of extreme precipitation, confirming the need for adopting highly localised data and/or dynamical downscaling procedures. A similar study, based on Germany [51] and considering the precipitation field, assessed a general underestimation of precipitation of various gridded data with respect to observations; the reanalysis datasets (ERA5, in this case) produced generally worse extreme value statistics of daily precipitation—in particular, failing in reproducing the accurate timing of observed daily precipitation extremes.

Concerning extreme temperatures over Europe, [52] also considered ERA5 fields; while the reanalyses captured the mean temperatures very well, the results over certain European sub-regions (e.g., the Alps and the Mediterranean, in particular) revealed that ERA5 underestimated temperatures. The main differences with observations can be attributed, according to the authors, to the altitude differences between ERA5 grid points and stations.

This work presents the first use of the EVT (extreme value theory) to study temperature and heavy rain in the Calabria region considering a long-term observational dataset and also investigating the possible use of gridded reanalysis data. The paper is organised as follows. Information about the study area and the adopted data, both observations and reanalysis fields, are provided in Section 2. The extreme-value-theory methods are described in Section 3. Section 4 reports the results and the general discussion about the reliability of the proposed EVT methods on extreme temperature and precipitation fields using both observations and gridded data. Our conclusions are summarised in Section 5.

2. Data and Study Area

The study area was the Calabria Peninsula, in southern Italy. The region is surrounded by the Tyrrhenian Sea (west) and by the Ionian Sea (east and south). The Apennines Mountains ideally separate the region into two sectors, crossing it from north to south more or less symmetrically; a maximum elevation of about 2000 m is reached. Daily extremes were computed considering a wide (regional) network of rain-gauges and temperature sensors (see next paragraph). State-of-the-art atmospheric reanalyses were also considered to evaluate the ability of the modelled gridded data to reproduce extreme daily values over Calabria.

2.1. Observational Dataset

Observed data are provided by ‘Centro Funzionale Multirischi della Calabria’ of the ‘Agenzia Regionale per la Protezione dell’Ambiente della Calabria (Arpacal)’ (http://www.cfd.calabria.it, accessed on 1 March 2022). The whole dataset of daily temperature and precipitation data from 1990 to 2020 (31 years), distributed quite uniformly over the whole region, was initially considered. Although the data are quality-controlled by the centre before being granted for research activities, we also adopted further selection criteria that permitted to retain 38 points/stations, starting from a larger number (254 rain gauges and 137 temperature sensors).

In particular: (i) each station was equipped with both a rain gauge and thermometer and (ii) each station had possible missing data of less than 3 years, even if non-consecutive, considering the whole 31-year period. Using these criteria, the percentages of missing data were 3.5% for temperature and 3.7% for precipitation. Starting from the hourly data, we computed the daily accumulated (0–24 h) precipitation and the daily maximum (Tmax), minimum (Tmin) and mean (Tmean) temperatures. The spatial distribution of the selected stations is shown in Figure 1. In the figure, eight zones in which the region has been ideally divided are also visible; the spatial extensions of these zones are comparable, and they divide the Calabria region more or less uniformly into the western (zones 1, 2, 3 and 4) and eastern (zones 5, 6, 7 and 8) side as well as in the northern (zones 1–5) and southern side (zones 4–8).

We used this subdivision as well, following that proposed by the Calabrian Regional Civil Protection (https://www.protezionecivilecalabria.it, accessed on 1 March 2022), which is responsible for disseminating weather-marine warnings on the region, and in order to make our results easily accessible (and usable) to the bodies responsible for the managing of weather alerts and emergencies.

From each of these zones, one station was selected and taken as representative of the area (see

Table 1. The eight zones and the related representative stations/ERA5 points used in the analysis. Longitude, latitude and altitude (only for observational stations) are also reported. The positions of each station for each zone are identified by the blue filled circles in Figure 1.

Zone	Station Name	Lon-Lat	Lon-Lat (ERA5)	Alt. (msl)	${\mathbf{T}}_{\mathbf{avg}}{(}^{°}\mathbf{C})$	ANP (mm)
1	Castrovillari	16.25–39.77	16.15–39.80	353	16.3	646
2	Montalto U.	16.13–39.40	16.15–39.30	468	15.7	1459
3	Nicastro	16.30–38.99	16.40–39.05	200	14.5	1139
4	Reggio Calabria	15.65–38.11	15.65–38.05	15	18.7	591
5	Acri	16.39–39.48	16.40–39.55	790	13.2	901
6	Crotone	17.13–39.09	17.15–39.05	5	17.9	630
7	Palermiti	16.45–38.75	16.40–38.80	480	14.6	1230
8	Serra S. Bruno	16.32–38.57	16.40–38.55	790	11.1	1630

) considering the stations with less missing data and closest to the ERA5 gridded points. For each of these eight stations, the extreme values of daily temperature and precipitation were computed following the extreme value theory (see Section 3). In the same table, for each station, we report the average mean temperature (${\mathrm{T}}_{\mathrm{avg}}$) and average annual precipitation (ANP), considering the 31-year period, in order to provide information about the climatic conditions of the different stations/zones.

2.2. ERA5 Reanalysis

With the aim of evaluating the skill of large-scale atmospheric reanalyses in reproducing daily extreme values, we used the global climate monitoring dataset ECMWF ReAnalysis (ERA5; [53]). The ERA5 fields are available hourly on regular latitude–longitude grids at ${0.25}^{°}\times {0.25}^{°}$ resolution. The fields taken into account are the hourly total precipitation and the 2 m temperatures for the whole 31-year period. Starting from these data, the daily ones were computed. The points of the ERA5 dataset are shown in Figure 1. As of the coarse horizontal spatial resolution, the number of ERA5 points is limited (about 30 points on the whole region, clearly uniformly arranged) but comparable with the retained measuring stations.

3. Methods

3.1. Extreme Value Theory

We analysed extreme daily temperature and precipitation using the statistical extreme value theory (EVT). The theory is often applied to quantify the stochastic behaviour of unusually large (i.e., extreme) phenomena. The procedure for establishing what is meant by extreme may be different. In this paper, two different approaches were used for the two different temperature and rainfall datasets. The motivation is inherent in the type of signal for the two different cases. Whereas the temperature daily data show a strongly periodic pattern dominated by seasonality, the rainfall data are more variable with non-regular peaks.

3.1.1. Peak over Threshold (POT)

In the POT approach, an observation is treated as an extreme if an associated measurement exceeds a predetermined threshold. For independently and identically distributed random observations $\left\{X_1,X_2,\setminus ,X_n\right\}$, the distribution function of exceedances X over a threshold u is $F_u\left(x\right)=Pr(X-u\le x|X>u)$ [31]. With a high enough threshold u, the $F_u\left(x\right)$ can be fitted by

$$G\left(x\right)=1-{\left(1+\frac{\xi x}{\widehat{\sigma }}\right)}^{-\frac{1}{\xi }},$$

(1)

defined on {$x:x>0$ and $(1+\xi x/\widehat{\sigma })>0$}, where

$$\widehat{\sigma }=\sigma +\xi (u-\mu ).$$

The family of distributions defined by Equation (1) is called the generalized Pareto distribution (GPD). Moreover, the parameters of the generalized Pareto distribution of threshold excesses are uniquely determined by those of the associated GEV distribution of block maxima. In particular, the parameter $\xi$ in (1) is equal to that of the corresponding GEV distribution [31].

In this work, the threshold approach was used to analyse the daily rainfall. This choice is due to the fact that the rainfall dataset is significantly more irregular than the temperature signal, and the influence of seasonality is significantly reduced. Due to its irregularity, the use of block maxima for rainfall is an expensive approach to the analysis of extreme values as many extreme events are discarded [31,54].

The threshold approach we described has an obvious dependency on the choice of threshold itself. Ref. [31] asserted that the threshold must be “as large as possible”, indicating some methods for determining this threshold from the data itself and not a priori. The first consequence of this approach is that the threshold u is not the same for every station but depends on the individual signals.

To determine the thresholds, in this work, we used the procedure known as “graphical method” [31]. This method consists of observing the evolution of the parameters $\sigma$ and $\xi$ as the threshold u changes. The choice of the threshold will be the highest value for which the parameters remain almost constant. If a GPD is used for the excesses of a threshold $u_0$, then the excesses of a different threshold u can also be described by a GPD. The shape parameters of the two distributions are the same. If ${\sigma }_u$ is the scale parameter for a threshold $u>u_0$, then

$${\sigma }_u={\sigma }_{u_0}+\xi (u-u_0).$$

The scale parameter changes with u unless $\xi =0$. To avoid this problem, the scale parameter is reparameterized as

$${\sigma }^{*}={\sigma }_u-\xi u.$$

This implies that ${\sigma }^{*}$ is constant with respect to u. Consequently, $\xi$ should also be constant with respect to u, and indeed they should be constant above $u_0$, if $u_0$ is a valid threshold for the excesses that follows the GPD [31].

3.1.2. Block Maxima

In the method known as “block maxima” (BM), extremes are created by dividing the analysis period into non-overlapping periods of the same size and then choosing the maximum observation of each new period. The choice of the block size is crucial because a very small block could create distortions, while, from too large a block, only certain extreme values could be selected [31]. The choice of the block length is subject to several considerations. Choosing a time period of one year and then extrapolating the annual maximum eliminates the problem of seasonal periodicity that strongly influences maximum temperatures. This choice, however, greatly reduces the number of extreme events, e.g., in a 30-year dataset, there will be only 30 values.

On the other hand, selecting monthly maxima allows for many more points but does not eliminate the seasonal contribution, clearly changing the distribution of extreme events. For example, the temperature block maximum in a spring season (defined as March to May) in the upper Midwest will most likely be the same as the May maximum, and thus the three-month block maximum may not be better approximated than only the May maximum. A block length of 1 year or one season has been used in a number of studies [55,56,57]. In other works, in an attempt to use as many points as possible and overcome the effect of seasonality, temperature anomalies were used [58,59].

Daily anomalies are calculated as the difference between each daily temperature and an average value. This value varies depending on the day of the year and is calculated as the average of the daily temperatures recorded over the time interval analysed, namely $\mathrm{\Delta }{T}_i=T_i-\langle T_i\rangle$, where $\langle T_i\rangle$ represents the temperature mean value for the ith calendar day. This definition implicitly assumes the annual seasonal cycle to be constant and generated by a set of stationary processes. The validity of this assumption is often questionable due to the non-linear response of the Sun–Earth system. Irregularities in the seasonal cycle have been observed as both amplitude [60] and phase variations [61,62,63].

Therefore, in this paper, following the ideas of Vecchio and Carbone [64], we use a new definition of temperature anomalies based on the empirical mode decomposition (EMD) [65]. The EMD has been successfully used in many fields [66,67,68], including geophysical systems [69,70,71]. The EMD decomposes the temperature signal into a finite number of intrinsic mode functions (IMF) and a residual, which describes the trend, using an adaptive basis derived from each dataset [65], i.e.,

$$T\left(t\right)=\sum _{j=0}^{m}IM{F}_j\left(t\right)+r_m\left(t\right)$$

(2)

Each IMF represents a zero-mean oscillation with amplitude and frequency modulations that both depend on time [65]. For each temperature signal, therefore, it is possible to identify the IMF representing the seasonal oscillation and subtract it from the original signal, so that the new definition of temperature anomaly is: $\mathrm{\Delta }{T}_j=T_j-S_j$, where $S_j$ is the IMF that describes the seasonal oscillation, which (in our dataset) is identified for all stations by the index $j=1$ and, for simplicity, we indicate by IMF 1.

Using this definition, the monthly maximum temperatures used in the EVT analysis maintain a satisfactory number of points and, at the same time, do not lose useful information from the subtraction of a value set by a simple average. In Figure 2, we show the maximum daily temperature recorded at the Crotone station (black dots); the purple signal is the one used in the EVT analysis obtained by subtracting from the monthly maxima the IMF 1 extrapolated from the monthly mean temperature.

The block maxima signal (purple line of Figure 2) can be fitted with the GEV (generalized extreme value) distribution:

$$G\left(x\right)=exp\left\{-{\left[1+\xi \left(\frac{x-\mu }{\sigma }\right)\right]}^{-1/\xi }\right\},$$

(3)

where x are the maxima, while the parameters $\mu$, $\sigma$ and $\xi$ are free parameters of the distribution that are constrained by the relation: $1+\xi \left(\frac{x-\mu }{\sigma }\right)>0$. This method was applied to the temperature signals for each station, and the results obtained are shown in the next section.

The choice of using the BM method on temperatures and the POT method on rainfall is, as mentioned above, determined by the different types of data. The temperature signals are more clustered (i.e., an extremely warm day is likely to be followed by another similar day), and the POT method requires the exceedances to be mutually independent. The rainfall signal is definitely more irregular and not clustered, and therefore with the BM method one could ’miss’ some extreme values occurring at close but independent times.

4. Results and Discussion

4.1. EVT Reliability

The purpose of this section is to evaluate the reliability of the EVT on the available datasets in order to assess the daily extreme rainfall and temperatures over the whole region.

4.1.1. Diagnostic Plots for the Selected Stations

We present the results obtained with the application of the EVT on the observational dataset. In particular, in this section, we primarily comment on the application of the theory on the selected eight stations previously defined. As already said, we used the POT method for the precipitation and the BM method for the temperatures; in the following, we use these terms to indicate the two EVT techniques applied to the two datasets. We used two different python packages in order to apply the POT and BM methods on the available dataset.

For the POT method, we used the software developed by [72], while for the BM approach, we used the software developed by [73]. In Figure 3, the diagnostic plot obtained via the application of the POT method on the rainfall data of Crotone station is reported. This station is chosen as representative, in terms of extreme events, for zone 6; the plots for the other zones are shown in Appendix A. The diagnostic plot is composed of four sub-plots: the probability density function plot (PDF), the probability plot (P-P plot), the quantile plot (Q-Q plot) and the return level plot.

The distribution of the values above the threshold is reported in the PDF plot, while the P-P plot and the Q-Q plot are graphical methods for comparing two distributions (theoretical and empirical), and they are basically the same plots but expressed on a different scale; if the data are adequate to model the extreme daily values of precipitations through the adopted theory, they have to lie on the diagonal of these plots. Through Spearman’s correlation coefficient, we quantified this correspondence between the empirical and theoretical return level obtaining, for the station of Crotone, a value of 0.94. The final plot is the return level plot that gives a probabilistic estimation of the repetition of an extreme event (see [49] for more details).

Figure 3a shows the distribution of the extreme values, extreme rainfall in this case, above the threshold values. They are presented as $excess-u_0$, where $u_0$ is the threshold value used for the station (in this case, the threshold is set to be $u_0=10$ mm), and it differs from station to station. The distribution has the maximum values located between 10 mm and 60 mm, and it is possible to observe a long tail of values that represents extreme rainfall events (less probable but more extreme).

Figure 3b,c permit us to assess if the adopted POT method is appropriate to describe the considered event. In these two plots, most of the points lie to the diagonal except for one. This point appears to be not well described by the POT model, and this could be due to the fact that this type of event is extremely rare. A similar behaviour is also shown in the other stations (see Appendix A) in which more points do not lie on the diagonal. It is possible to define them as “abnormal” because they are out of description for the EVT. In panel (d) of Figure 3, the return level plot is reported, which gives a probabilistic forecast of extreme rainfall events based on the extreme events detected by the station.

All values (also the abnormal ones) are in the confidence interval of $2\sigma$, and the POT model (black solid line) seems to describe the behaviour of the data very well, except for a few points. For these cases, the EVT theory seems to either overestimate and underestimate some observed extreme values. In the case of the overestimation, the extreme data events appear below the theoretical line, and we can define them as “delay” because it seems they happen after the time expected from the theory. In the case of underestimation, the data are located above the theoretical line, and we can define them as “advances” because they occur before what is expected from the theory.

Figure 4 shows the diagnostic plot (obtained via the application of the BM method) for the extreme temperature and, in this case, for the Crotone station (the plots for the other stations/zones are reported in Appendix A). The plot composition is the same as Figure 3. The density plot differs from that seen for the rainfall case because we used a different approach of EVT (i.e., the BM approach). We used the GEV model, which takes into account the maximum values of the temperature for each month of the dataset. The density plot (Figure 4a) shows the distribution of the maximum temperature values.

The distribution does not appear symmetric; however, it exhibits a pronounced tail in the right side of the figure. This indicates the presence of high extreme temperatures. Figure 4b,c are similar to the plot shown for the rainfall case. Most of the points lie on the diagonal except for a few of them, indicating that the adopted GEV approach well describes the dataset of extreme temperature. Again, Spearman’s correlation coefficient between empirical and theoretical return level was calculated, obtaining, for the temperatures in Crotone, a 0.95 value. The return level plot (Figure 4d) shows a similar behaviour as the rainfall data; however, in this case, some values do not belong to the confidence interval of $2\sigma$. These points are probably the ”abnormal” data discussed before, and the BM method fails in their description.

As already said, we report in Appendix A the diagnostic-plots for all the other seven stations/zones (Section 2.1) for both precipitation and temperature extremes. Considering the other zones, the results are quite similar to zone 6. It is evident that the higher values of extreme precipitation are detected in the zones located in the eastern side of the Calabria region (along the Ionian sea). These zones are characterised by more extreme rainfall events, with respect to the zones along the Tyrrhenean sea, and this is in agreement with several previous works confirming that the east side of the region is mainly affected by higher precipitation events. An inverse behaviour was found for extreme temperatures; from the results obtained with the application of the BM method, the highest extreme values appear located on the western side of the region, i.e., along the Tyrrhenian sea.

4.1.2. Return-Level Maps for the Observational Dataset

In this section, we discuss the application of the EVT method on the whole observed dataset. The purpose of this analysis is to qualitatively compare, in terms of return periods, the areas affected by daily precipitation and temperature extremes. Clearly, the maximum time interval for which such a comparison can be made is 31 years (i.e., the time interval in which the observations are available).

The return level plots give us the probability that a certain event will be expected into a certain return period interval; this probabilistic method allows us to better identify the regions subjected to extreme phenomena. In Figure 5a,c, we report the map of the empirical (i.e., considering the observed daily precipitations) return level for 10 and 30 years, respectively, while Figure 5b,d are the theoretical (i.e., applying the POT method on the observed daily precipitations) return-level maps for 10 and 30 years, respectively. These maps, as with all those in this section, are produced using the inverse distance interpolation method, considering all the available points (Figure 1).

The maps have similar characteristics with the same features of daily extreme rainfall as also anticipated in the diagnostic plots of Figure 3 (although they refer to a single station). Comparing the 10- and the 30-year-return periods, a general increase in extreme rainfall values is visible. These probabilistic results, based on observations, can be used to assess the expected intensity of extreme precipitation events, providing a useful complementary tool to climate projections and models, that hypothesise, in the coming years, an increase of such events.

For these fixed return periods, it is easy to identify a zonal gradient of extreme precipitation over Calabria, with the south-eastern part of the region most affected by such events (the Ionian areas, in particular), as opposed to the north-western part. This occurrence is found by considering only the observational data and the application of the POT method, with values up to 300 (400) mm/day for the fixed return period of 10 (30) years.

The result is not surprising since, as said in the introduction, several studies previously assessed how the most intense rainstorms are more frequent on the east side of the region. In this case, and also confirmed by the aforementioned works (e.g., [16]), the role played by the orography of the region and by the prevailing synoptic conditions associated with extreme precipitation events in southern Italy is visible. In such situations, the presence of a cyclonic area located in the southern Ionian Sea [23] prevails, which draws more unstable and humid air on the Ionian coast of Calabria.

The contrast between these air masses with the terrain and the orography of the region, in particular the Aspromonte mountain range, causes an uplift of the air masses and a consequent increase in convective instability conditions (see, for example, the case study analysed in [19]), and this results in a significant amount of precipitation windward of the orographic reliefs. In Figure 6a,c, we report the map of the empirical (i.e., considering the observed daily maximum temperatures) return level for 10 and 30 years, respectively, while Figure 6b,d are the theoretical (i.e., applying the BM method on the observed daily maximum temperatures) return-level maps for 10 and 30 years, respectively.

Furthermore, in this case, the maps have similar characteristics with the same extreme temperature features as also anticipated in the diagnostic plots of Figure 4. Considering the fixed return period of 10 years, the areas most affected by extreme daily temperatures, and therefore possible heatwaves, are the flat areas of the region, in particular those on the Tyrrhenian side. Observing the differences between the 10- and the 30-year-return periods, we can see a general increase of the daily extreme temperature values.

As for the case of daily precipitation, it is possible to use these results to obtain information on the local effects of general global warming possibly linked to the ongoing climate change. Considering the 30-year-return period, the daily extreme temperatures, apart from being higher, seem to be more concentrated in the southern part of the region, although it is possible to note a small area with high temperature values also in the north part. The observed data (Figure 6c) suggest an intensification of extreme maximum temperature events in the south-eastern areas of Calabria, while those obtained from the application of the BM method (Figure 6d) confirm the feature seen in the 10-year-return period maps, with the Tyrrhenian part most affected by these extreme thermal features.

It is important to note that the empirical return values are not known for all the stations, especially for larger return times. This is because, being an empirical value, it is not always possible to detect it at the same return time for all stations. This implies that, for higher return times, the spatial interpolation was calculated with fewer points. This problem does not arise for the theoretical return value, since extreme events are estimated theoretically and for each return time chosen by the numerical procedure (in our case, with a variation of 0.1 years).

The empirical (based on observed data) and the EVT maps are very similar; this implies that, for both daily extreme precipitation and temperature, the EVT applied on the observations describes their behaviour well.

In considering these first results, it is important to note the significant difference in the number of available stations over the region. In particular, the eastern-northeastern part of Calabria (zone 6 and 5, above all) is much less covered by observations (see Figure 1) and, therefore, also the interpolation procedure for creating the maps is inevitably conditioned by this fact.

4.1.3. Return-Level Maps for the Reanalysis Dataset: Comparison with Observations

In this section, we evaluate the possible use of ERA5 reanalysis data in studying extreme weather events, since regular gridded data are not affected by spatial/temporal availability problems. We consider the results obtained from observed data and from ERA5, comparing the return-level map for a fixed return period of 10 years, both for temperature and precipitation.

In Figure 7, we show the return-level maps of extreme daily rainfall; Figure 7a is the same as Figure 5a (here duplicated, to allow a direct comparison), while Figure 7b is the map obtained by the EVT theory application (POT method) on the reanalysis. In order to better comment on the results, we reduced the colour scale values on the maps, mainly because the maximum extreme rainfall values for the ERA5 data are significantly lower. Figure 8 is the same as Figure 7 but for the daily extreme temperatures (BM method).

A general and unquestionable underestimation is evident for the ERA5 data with respect to observations, both for precipitation and temperature. For the extreme daily rainfall, differences between observations and reanalysis even higher than 200 mm are clearly visible in several zones (mainly in the southeast). The general behaviour of the rainfall is similar with the eastern part of the region most affected by extreme precipitation. A poor agreement is noticeable in the southeastern part of the region, the most affected by heavy rain based on observations, where the reanalysis failed to correctly reproduce the high rainfall values.

As usual, a separate comment applies to the northeastern part of Calabria where, as already said, there is a general lack of weather stations (see Figure 1); for this reason, we cannot obtain useful information using the POT method application in this area. Furthermore, for the daily extreme temperature, no useful information can be drawn on the northeastern part of the region, due to the lack of stations. The areas most affected by extreme daily temperatures (based on observations) are the Tyrrhenian ones, as already seen in Section 4.2. The ERA5 fields fail to correctly identify this behaviour, except for the southwestern side of the region, instead showing a greater predisposition to extreme temperatures in the southeastern area (not confirmed by observations).

This general underestimation is not surprising, and it is in agreement with previous works ([50,51,52]) that have attempted to use ERA5 fields for highly localised climatic studies. To further highlight this aspect, also providing a more quantitative result, we show, in Figure 9, a punctual comparison (only for the Crotone station) between the return levels obtained using both the observed data and ERA5 reanalysis. In this case, we show the return levels computed considering both the original datasets (OBS and ERA5) and the EVT results.

Figure 9a refers to extreme precipitation and shows the return level of the observed data (black dotted points) and the return level of the POT method (black solid line) applied on the observations. In the same plot, we present the results from ERA5 data (blue dotted points) and the return level of the POT method (blue solid line) applied on the reanalysis. Dotted lines indicate the error confidence intervals.

The discrepancy between the plots and the general ERA5 underestimation with respect to observations is evident, although the application of the POT method provides, also on the reanalysis, valuable results (see blue dotted points vs. blue solid line). Similar findings can be seen in Figure 9b for the temperature daily maxima. Here, the differences, in terms of temperature values, appear more pronounced, and also in this case, the application of the BM method seems to work quite well.

Regarding the ERA5 underestimation of the extreme values, the reanalysis does not seem capable of well grasping the localised dynamic characteristics influencing the temperature and rainfall fields in the region. Regular gridded modelled data represent a valid and promising tool for these types of analysis, even irreplaceable in the vast worldwide areas not sufficiently covered by observations, and the EVT technique also applied to ERA5 seems to work well. However, the fact that the reanalysis unequivocally underestimates the observed extreme values prompts us to perform our further EVT analyses only considering the observational dataset as a reference.

4.2. 50 and 100 Years Return Levels

We present in this section, the 50- and 100-year return levels by applying the EVT to the whole observational dataset. This probabilistic approach allows us to identify the Calabrian areas that are likely the most affected by daily extreme precipitation and maximum temperatures in the future. We show the return-level maps obtained from the application of the POT (for rain) and the BM (for temperature) methods for each available station. The analyses were conducted for the fixed return periods of 10, 30, 50 and 100 years. Since the return-level maps for 10 and 30 years were already presented in Section 4.1.2 (Figure 5b,d and Figure 6b,d), we report in this section only the maps referring to the 50- and 100-year-return periods.

Figure 10a,b show the return levels at 50 and 100 years, respectively, for the extreme daily rainfall. The figures clearly show two distinguished regions in which the rainfall extreme events have different features. The southeastern part of the region represents the area most affected by extreme daily precipitation, unlike the northwestern part. The behaviour, in terms of precipitation patterns, is very similar for all the return level periods. Observing Figure 1, we can see how the stations are well-distributed over the region except for the northeastern part; in this area, there are no predicted extreme precipitation values; however, we cannot truly quantify the role of the EVT theory, simply because of the lack of observational data as already mentioned in the previous sections.

On the other hand, the fact that the most intense precipitation in Calabria occurs in the southeastern part of the region is a further confirmation of previous work results ([16,18]). A sensible intensification of the extreme rainfall is evident as the return periods increase in qualitatively agreement with the increasing incidence expected by the climate projections. The Aspromonte and the areas windward the mountain range are the most affected zones for daily extreme precipitation. For the 50-year-return period, extreme values up to 500 mm/day are predicted in this zone, and the daily extreme precipitation assumes even greater values (up to 600 mm/day) if a return period of 100 years is considered.

Figure 11a,b show the return levels 50 and 100 years, respectively, for the extreme daily temperatures. As already seen for the return periods of 10 and 30 years (Figure 6b,d), these probabilistic projections suggest that the areas most affected by extreme daily temperatures will be the Tyrrhenian flat areas. The central part of the region seems not to be affected by daily extreme temperature, also because it is in the central areas that the main mountain ranges of the region are distributed, while the lack of observations in the northeastern part does not permit extracting very useful information for this area in applying the EVT. Furthermore, for temperatures, an intensification of the daily extreme is visible as the return periods increase in agreement with the ongoing global warming.

Heatwaves with maximum daily temperatures up to 40 °C are expected in the next 100 years and the southern part of Calabria is the one mainly affected; it is even more interesting to note how, in the same areas, differences of maximum daily temperatures of about 2 °C are visible, considering the 50- and 100-year-return periods. The maps obtained at 50 and 100 years (Figure 10 and Figure 11) represent a great advantage of EVT theory, namely to be able to provide a forecast of the return value exceeding the time scale imposed by the dataset (in our case of 31 years).

The values shown in the maps refer to the most probable RL; however, each has a confidence interval highlighted by the dotted lines in the return level plot (see Figure 3d and Figure 4d). The error for each station is impossible to be shown in a single map; thus, we decided to report two further figures showing the precise RL values and the related confidence intervals for all the stations listed in Table 1 (Figure 12 and Figure 13). It is important to highlight that these results are strongly conditioned by the choice of the stations.

From the box-plots, it is easy to see how the representative stations of zones 7 and 8 are those more prone to an intensification of extreme daily precipitation events over the years, i.e., the areas on the southern Ionian side, where values between 350 and 400 (450 and 500) mm/day are predicted for the return period of 50 (100) years. Considering instead temperatures, the future incidence of heat waves is more likely in the representative stations of zones 4 and 6, i.e., the southern Tyrrhenian side and the eastern part of the region, respectively, where maximum daily temperatures between 36 and 37 (36 and 38) °C are predicted for the return level of 50 (100) years.

5. Summary and Conclusions

In this work, we presented a statistical study on the daily precipitation and temperature extremes in Calabria, southern Italy, applying the extreme value theory (EVT) on a long-term series (31 years) of regional observations. For the precipitation, we adopted the peak-over-threshold method, while for the temperature, the block-maxima approach was considered.

Both datasets (observations and reanalysis) present certain limitations. The first, regarding the observations, is related to the lack of stations in a specific area (the northeastern part of the region) and only affects the procedure (inverse distance) adopted for generating the maps. A further limitation is related to the time series (or limited to 30 years) and can affect the statistical uncertainty leading to an increase of the confidence interval in the return level plot.

The main findings from this work can be summarised as follows:

• The proposed EVT methods were suitable to describe the precipitation and temperature extremes over the study area.
• The reanalysis fields showed a systematic underestimation of daily precipitation and temperature extremes with respect to the observations.
• Extreme precipitations over Calabria mainly affected the southeastern part of the region; values up to 500 mm/day were predicted for the 50-year-return period in this area.
• Extreme temperatures, instead, mainly affected the Tyrrhenian side of the region with values up to 40 °C predicted for the 100-year-return period in this area and in the southern part of the region.

In order to perform punctual statistical investigations, the region was divided in eight climatic zones, in agreement with the spatial subdivision adopted by the Calabrian Regional Civil Protection; from each of these zones, one station was selected and taken as representative of the area. To assess the reliability of the method, we first reported the diagnostic plots (i.e., the probability density function, the probability plot, the quantile plot and the return level plot) for each selected station, considering both the empirical (the directly observed values) and the theoretical (the EVT results) data.

The analysis of the diagnostic plots confirmed that the proposed methods are appropriate to describe the considered extreme weather events over the study area. In order to qualitatively compare the areas most affected by such extremes, again considering both empirical and theoretical data, we reported the regional maps of daily precipitation and temperature extremes in terms of fixed return periods (10 and 30 years). Furthermore, in this case, a good agreement was obtained between the two datasets, thus, confirming the reliability of the EVT applied on the long-term observations.

An evident zonal gradient of extreme precipitation over Calabria was found, with the south-eastern part of the region most affected by such events, in agreement with several previous studies. The flat areas of the region, in particular the Tyrrhenian ones, were instead the most affected by extreme daily temperatures and possible heatwaves.

The possible use of reanalysis in studying the daily extremes was also investigated, comparing the results obtained considering the observations with those derived by ERA5. Although the use of modelling products represents a great opportunity, mainly in areas not covered by observations, our comparisons show how the ERA5 fields are not suitable for reproducing the occurred extreme values, as they showed a systematic and substantial underestimation of both the daily precipitation and temperature extremes.

The assessed reliability of the theory applied on the observations allowed us to perform a probabilistic analysis aimed at identifying the areas most likely to be affected by daily extreme precipitation and maximum temperatures in the future; with this aim, we reported the return-level maps for the fixed periods of 50 and 100 years for both precipitation and temperature. A general intensification of extreme rainfall events is expected in future years, in agreement with climate projections, with values up to 500 mm/day for the 50-year-return period in the south-eastern part of Calabria. Considering the temperature, heatwaves with maximum daily values up to 40 °C are expected in the next 100 years, mainly in the Tyrrhenian and southern part of the region.

The unequivocal increase of the Earth’s surface temperature and the predicted exacerbation of instability and convective atmospheric conditions in the future emphasise the importance of performing studies, such as the one presented, devoted to the analysis of the environmental variables most connected to the risks deriving from extreme meteorological phenomena.