Multi-Secular Trend of Drought Indices in Padua, Italy

Abstract

The aim of this work is to investigate drought variability in Padua, northern Italy, over a nearly 300-year period, from 1725 to 2023. Two well-established and widely used indices are calculated, the standardized precipitation index (SPI) and the standardized precipitation evapotranspiration index (SPEI). They are compatible with a data series starting in the early instrumental period, as both can be estimated using only temperature and precipitation data. The Padua daily precipitation and temperature series from the early 18th century, which were recovered and homogenized with current observations, are used as datasets. The standard approach to estimate SPI and SPEI based on gamma and log-logistic probability distribution functions, respectively, is questioned, assessing the fitting performance of different distributions applied to monthly precipitation data. The best-performing distributions are identified for each index and accumulation period at annual and monthly scales, and their normality is evaluated. In general, they detect more extreme drought events than the standard functions. Moreover, the main statistical values of SPI are very similar, regardless of the approach type, as opposed to SPEI. The difference between SPI and SPEI time series calculated with the best-fit approach has increased since the mid-20th century, in particular in spring and summer, and can be related to ongoing global warming, which SPEI takes into account. The innovative trend analysis applied to SPEI12 indicates a general increasing trend in droughts, while for SPI12, it is significant only for severe events. Summer and fall are the most affected seasons. The critical drought intensity–duration–frequency curves provide an easily understandable relationship between the intensity, duration and frequency of the most severe droughts and allow for the calculation of return periods for the critical events of a certain duration. Moreover, the longest and most severe droughts over the 1725–2023 period are identified.

1. Introduction

In the wake of what has happened in the past two years, spring and summer of 2024 were again warm and dry in many European regions [1]. Global warming and increased rainfall extremes are two of the main hallmarks of climate change. At the same time, future projections indicate that precipitation will likely decrease in many areas, enhancing the risk of drought [2].

Drought exhibits a high variability at both temporal and spatial scales; therefore, there is a great uncertainty in capturing trends compared to other natural hazards. Results of recent studies are quite contradictory at the global scale, while at the pan-European scale, a slightly increasing tendency has been identified in the last decades [3]. At the continental scale, trend studies have shown a decreasing trend in Northern Europe and an increasing trend in Southern Europe. In particular, in the Mediterranean coastal area, the synergy between reduced precipitation and increased evapotranspiration fostered by higher temperature will affect the hydrological balance, leading to decline in water availability [2]. Droughts are projected to become more severe, more frequent and longer under moderate emission scenarios and strongly enhanced under severe emission scenarios ([2] and references therein). Therefore, as spatial scale is one of the main sources of uncertainty, the reconstruction of climate evolution at high spatial resolution is crucial for drought variability assessment and consequently for risk mitigation, preparedness and response [2,4].

In recent research, droughts are usually subdivided into four types according to their duration and impact on different systems: meteorological, agricultural, hydrological and socio-economic drought [5]. This study deals mainly with hydrological type because only meteorological variables are used as input (precipitation and temperature), and data analysis is focused on a 12-month accumulation period, a compromise independent of individual seasons and regimes, which allows for the representation of water shortages caused by a lack of precipitation and/or high temperatures over an entire year. Besides the spatial scale, the selected time period is another important variable, together with the drought index formulation.

In the past, several indices were proposed to objectively quantify and compare droughts in different climatic and hydrological regimes. A detailed review can be found in Mishra and Singh [6]. Among the existing indices, the standard precipitation index (SPI) [7] and the standardized precipitation evapotranspiration index (SPEI) [8] are the most widely used, in particular in the case of datasets from the early instrumental period when precipitation and temperature were the only variables available. While SPI considers precipitation only, SPEI is based on water balance (BAL), i.e., the difference between precipitation and potential evapotranspiration (PET). Precipitation is certainly the main driver of drought, but the role of temperature cannot be neglected, in particular in the Mediterranean region, which can be considered a global warming hotspot [9].

Both SPI and SPEI are probability-based indices that consider the distribution of precipitation and water balance, respectively. The choice of the most appropriate probability distribution function for precipitation and water balance is crucial, as an improper distribution may impart bias on the index values, exaggerating or minimizing the severity of drought events. Several probability distribution functions (PDFs) have been tested in the literature, but a general consensus has not been reached, as the best-performing distribution depends on location, spatial scale and accumulation period. When working in an area not previously explored or with a new dataset, extensive and rigorous statistical testing is recommended.

Stagge et al. [10] compared a suite of candidate PDFs and used a new method based on the Shapiro–Wilk (SW) test to evaluate the goodness of fit. Relative distribution ranking was determined by the Akaike information criterion (AIC). For regional studies in Europe, they recommended the two-parameter Gamma distribution when calculating SPI, in agreement with previous studies [7,11]. Instead, when computing SPEI, the use of the generalized extreme value (GEV) was a better alternative specifically for the European domain than the three-parameter log-logistic (Log-log3) distribution mostly adopted in the literature [8,12]. Angelidis et al. [13] analyzed a case study in Portugal and found different PDFs depending on the accumulation period, i.e., gamma for time scales up to 6 months and log-normal for 12 months. According to Pieper et al. [14], the exponentiated Weibull distribution for SPI performs proficiently in simulations and observations for every common accumulation period and virtually everywhere around the globe. The AIC was used to penalize unnecessary complexity of the candidate PDFs. Concerning SPEI, Monish et al. [15] found that Pearson type 3 (PE3) outperforms the other distributions, especially for shorter timescales, while GEV performs well for longer timescales for most parts of India. The findings of Wang et al. [16] confirm the choice of the Log-log3 distribution for around 500 stations in China. Ramezani et al. [17] examined 65 PDFs for the estimation of SPEI in eastern Iran, and the results of Anderson–Darling (AD), Kolmogorov–Smirnov (KS) and chi-square tests indicated the four-parameter Burr distribution as better performing than the Log-log3 distribution. The generalized logistic distribution (Genlog) was the most appropriate for a case study in Ethiopia that included 125 stations and a 30-year dataset [18].

Passing to the national scale, in Italy, recent studies have focused on the analysis of drought frequency and intensity variations without addressing the issue of the distribution function of drought indices and instead simply adopting the standard approach (SA) in the literature, i.e., gamma and Log-log3 distributions for SPI and SPEI, respectively. Moreover, these studies deal mainly with drought characterization for central and southern regions [19,20,21] and less so for northern ones [22,23]. Moccia et al. [24] tested different PDFs on a 50-year dataset containing records from all over Italy. The log-normal distribution was the best-fitting model to describe almost all the monthly precipitation samples, followed by Weibull (for 1-month scale) and gamma. The fitting performance was evaluated with a modified mean square error normalized (MSEN) and the normality with the SW test.

Concerning drought trends, the two studies focused on northern Italy [22,23] are in accordance with the Europe-wide analysis, which shows increasing drought frequencies in Southern Europe [3]. Nevertheless, drought characteristics, mainly intensity, strongly depend on the geographic area, as well as the spanned period and temporal resolution. Both studies cover the period from around the mid-20th century to the early decades of the 21st century and are based on both SPI and SPEI. Crespi et al. [23] analyzed a portion of the Po Plain in the Lombardy region and found drying tendencies, in particular in summer, in the southern and western parts of the domain. The more negative trends of SPEI than of SPI were ascribed to the increasing role of evapotranspiration over recent decades triggered by rising temperature. According to Baronetti et al. [22], since 2001, drought episodes in the Po Valley have increased in terms of frequency and duration, and they seem to be mostly related to changes in the intra-annual precipitation distribution.

The present study fits into this context with the general aim to improve the understanding of drought variability in Padua, Italy, over three centuries. The Padua precipitation and temperature series are the longest ones available in Italy and are used for the first time for drought assessment.

The specific aims are as follows:

• To assess the best-performing PDFs for use in the SPI and SPEI normalization of the multi-secular series of temperature and precipitation in Padua;
• To investigate the impact of the best-fit approach (BFA) compared to the standard approach (SA) on identifying drought characteristics, i.e., duration, severity and intensity;
• To analyze the trends in drought patterns, extreme events and periodicity.

The paper is organized as follows: Section 2 is devoted to a description of the study site, the datasets and the methodology used for data analysis, divided into several steps, each with a specific purpose; results are discussed in Section 3, following the general lines of the methodology. Finally, conclusions are drawn in Section 4.

2. Materials and Methods

2.1. Study Site and Climatic Data

Padua is a city located on the eastern side of the Po Valley, about 40 km west of Venice and the Adriatic Sea (Figure 1). The municipal territory covers 93 km² entirely, with a mean altitude of 12 m a.s.l.

The climate in Padua is classified as humid subtropical (Köppen climate classification Cfa) with cold winters and hot summers, frequently associated with air stagnation [25].

Although Padua is an urban center, agriculture is intensively practiced in the suburbs that surround the historic center. In terms of production, a prevalence of arable farming is observed, and most of the orchards and vineyards have a family character. As in other parts of the Po Valley, agriculture in Padua is largely dependent on irrigation and, therefore, is strongly influenced by droughts.

The datasets used to assess drought variability in Padua consist of the unbroken series of daily temperature and precipitation, the longest in Italy, covering nearly 300 years, from 1725 to the present day. The original data were recovered, corrected and homogenized with the help of a careful interpretation of the metadata during a demanding study that started around 40 years ago and has been continued up to now. In particular, the bias due to changes in instruments, exposure location, height above the ground and observation protocols were corrected. Details can be found in the rich literary production [26,27,28,29,30,31,32,33] and references therein. The final composition of the temperature series was presented by Stefanini et al. [33]. The comparison with independent datasets, such as data collected at other Italian meteorological stations (Milan, Bologna, Turin, from 100 to 300 km away from Padua) and a reanalysis dataset (ModE-RA), confirmed the general trend and anomalies, assessing the reliability of the Padua temperature series. The completion of the precipitation series involved further work, described here for the first time (Section 3.1). The assessment of the accuracy of the precipitation series is much more difficult than that for temperature. Given the high variability in precipitation, only stations very close to Padua can be considered, and unfortunately, no station near Padua provides a dataset covering nearly three centuries.

2.2. Methodology for Calculation of Drought Indices

As monthly data are required for the calculation of SPI and SPEI, daily precipitation amounts were summed, and monthly and daily temperature values were averaged month by month. The potential evapotranspiration (PET) was calculated using Thornthwaite’s equation [34]. It was the only possible solution considering the available data, since it requires only mean temperature values. Other methods require observations of other variables, such as daily minimum and maximum temperatures, vapor pressure and wind speed, which are not available for the early instrumental period.

The evaluation of the best-performing distribution function is the key point in the methodology followed to calculate the drought indices, as the subsequent climatic analysis is based on the index values. The overall methodology is visualized in Figure 2 and includes the main steps described hereunder. The input datasets of the process are the monthly precipitation totals for SPI and the monthly water balance (BAL) for SPEI.

• Application of a fitting distribution to the input dataset;
• Calculation of the drought index;
• If the index results are normally distributed, go to the next step; otherwise, try another distribution function;
• Application of goodness-of-fit tests to the normally distributed index values.

SPI and SPEI are calculated at a monthly temporal resolution, and the accumulation periods considered are 3, 6 and 12 months. The same nomenclature is used as in previous studies, i.e., the index acronym is followed by a number indicating the accumulation period, e.g., SPI3 of a certain month is calculated cumulating the precipitation over the 3-month period ending in that month. Similarly, SPI12 of a certain month is calculated cumulating the precipitation over the 1-year period ending in that month.

The end of each run of the process is the dataset of the values of one index, calculated with the PDF that (i) provides a normally distributed index and (ii) exhibits the best-fitting parameters. The process was applied separately to each index (i.e., SPI and SPEI) for every accumulation period k, first considering the whole dataset (i.e., SPI_k and SPEI_k, where k = 3, 6, 12) and then separating each month (i.e., SPI_j,k, SPEI_j,k, where k = 3, 6, 12 and j = 1,…,12).

2.3. Goodness-of-Fit and Normality Tests

In addition to the PDFs recommended in previous SPI/SPEI studies (see Section 1), new distributions were also considered, limiting the number of parameters to two for SPI and three for SPEI. Some functions were directly available in the R extraDist package [35], while new scripts were created to calculate the others. Due to the length of the input series, the whole period from 1725 to 2023 was chosen as the reference period for both precipitation and temperature.

The distribution that fit the data best was assessed using the modified mean square error normalized function (MSEN) [24,36,37,38]. This method has the main advantage that it allows for the simultaneous estimation of the unknown parameters and the identification of the best-performing distribution between the candidates. The AIC was also applied to provide robustness to the results. Once the best PDF was identified for each input series, the SW normality test [39] was applied to the calculated values of SPI and SPEI. The distributions that yielded values of SPI and SPEI meeting the requirements described by Wu et al. [40] were considered non-normal and consequently rejected.

2.4. Standard and Best-Fit Approaches

At this point, it is possible to distinguish between two types of approach for calculating the drought indices, i.e., the standard approach (SA) and the best-fit approach (BFA). The SA foresees the use the same PDF for the whole input series, which is gamma for SPI and log-logistic 3 for SPEI. The BFA instead considers the best-fit distribution. This approach can be further split into two: one case considers the same function for all the months, called “general” in the text and indicated as BFAg; the other uses a different PDF for each month, named “monthly” and indicated as BFAm.

Drought and wet events are classified based on SPI/SPEI values [7] as reported in

Table 1. Events classification based on SPI/SPEI values.

SPI/SPEI Values	SPI/SPEI Classes
≥2.00	Extremely wet
1.50 to 1.99	Severely wet
1.00 to 1.49	Moderately wet
−0.99 to 0.99	Near normal
−1.49 to −1.00	Moderate drought
−1.99 to −1.50	Severe drought
≤−2.00	Extreme drought

.

The impact of the different approaches on the main drought characteristics, i.e., duration, severity, intensity and frequency, was then evaluated. Drought events were identified following the definition of McKee et al. [7], i.e., a drought starts in the month (included) when the index value falls below −1 and ends in the month (not included) when its value becomes positive for at least two consecutive months. Duration (D) refers to the period from the beginning to the end of a drought event and is expressed in months. Drought severity (S) is the sum of SPI/SPEI values within the period D, and intensity (I) is the average of SPI/SPEI values, i.e., the ratio of severity and duration. Finally, the number of drought events in a period is defined as the drought frequency (F).

2.5. Trend and Periodicity

One of the most commonly used trend tests is the Mann–Kendall (MK) test [41,42], the power of which depends on the preassigned significance level, magnitude of the trend, sample size and amount of variation within a time series [43]. Moreover, its application is limited by a set of restrictive assumptions, e.g., independent structure of the time series, normal distribution and data length. Sen [44] introduced a new simple and immediate method to investigate trends in time series, free from the abovementioned assumptions. This innovative trend analysis (ITA) has been already applied to drought trend assessment [45,46]. The procedure is quite simple: the time series is divided into two equal parts, and each one sorted in ascending order separately. Then, a scatter plot of the two parts is created, with the first half on the X-axis and the second half on the Y-axis. Data points plotted on the 1:1 (45°) ideal line indicate no trend in the time series. There is an increasing or decreasing trend if data are plotted above or below the ideal line, respectively. In addition, the distance from the 1:1 line indicates trend acceleration or strength. Both the MK and ITA tests were applied to assess the presence of trends in the Padua drought series, and the results were compared.

Before performing ITA, the presence of change points in the SPI and SPEI series was investigated using the standard normal homogeneity test (SNH) [47], a test sensitive to the beginning and end of the series. The SNH test was successfully applied to the temperature and precipitation time series for Padua [26,32].

Finally, a continuous wavelet transform (CWT) analysis was performed to examine the periodicities in drought index time series [48,49]. The wavelet spectrum was evaluated using the R package biwavelet [50], while the red-noise spectrum was analyzed using the R package dplR [51].

2.6. Critical Drought Intensity–Duration–Frequency (IDF) Curves

Drought IDF curves are helpful in evaluating the relation between intensity and frequency for several drought durations, providing useful information about drought events using a single graph that is easily understandable by end users for practical purposes. Aksoy et al. [52] developed a procedure for estimating IDF curves for critical drought and applied it to the SPI. The same method was recently applied to an SPI-based Italian case study by Moccia et al. [24] and extended to the SPEI by Arra et al. [53]. According to Aksoy et al. [52], a drought of duration D is defined as critical for a certain year when it is the most severe drought in terms of the total deficit over that duration. It is likely that no drought was observed in some years; therefore, a zero value is taken for the critical drought severity in a no-drought year. Critical drought IDF curves were implemented following the procedure described in detail by Aksoy et al. [52], which was developed through the following steps:

• Identify the critical droughts of each duration D for each year;
• Compose the time series of the intensities of the critical droughts for each duration D; if a drought of a certain duration was not critical for a given year, the intensity for that year is set equal to zero;
• Identify the best-fit probability distribution for each time series using the MSEN method;
• Calculate the intensity of critical droughts of a given duration for a fixed return period using the probability distribution identified in (3);
• Perform a linear regression for each return period between the intensity and duration: $I=aD+b$, where D is the duration in months (it varies between 1 and the maximum observed), and $a$ and $b$ are the regression parameters.

Then, from the critical drought IDF curves, it was possible to find a relation between drought intensity and return period (RP) for each drought duration and thus calculate RP.

3. Results

3.1. Precipitation Series

The precipitation series from 1725 to 1950 has been already published in recent works [26,29,30,31]. The missing months of September and October 1825 were filled using the contemporary monthly amounts recorded in Venice. The procedure followed to build the series from 1951 to 2023 is described hereunder and visualized in Figure 3. The locations of the meteorological stations mentioned in the list are shown in Figure S1.

• 1951–1954 Meteorological Observatory of the Water Magistrate (PD_WMa)
  The climatological day starts at 00 local time (LT).
• 1955–1978 Meteorological Observatory of the Water Magistrate (PD_WMb)
  During this period, the climatological day starts at 09 LT; therefore, 9-9 daily amounts were adjusted to 0–24 following the methodology described in [26].
• 1979 Meteorological Service of the Italian Air Force at Padua Airport (PD_AM)
  The climatological day starts at 00 UTC, i.e., 01 LT; therefore, daily amounts were adjusted to central European time (CET) [26].
• 1980–2023 Botanical Garden Station of the Regional Agency for the Prevention and Protection of the Environment in the Veneto Region (ARPAV) (PD_OB)
  Gaps were filled with the PD_WM data adjusted to 0–24 as in (1), when available; otherwise, they were filled with the daily amount recorded at Legnaro station, the ARPAV meteorological station closest to PD_OB in the Padua outskirts.

PD_AM data, except those for 1979, were excluded because this period is documented by a number of other contemporary series, which proved to be more reliable, such as PD_WMa, PD_WMb and PD_OB. As an example, a cloudburst that occurred in September 1967, which caused an accumulation of over 100 mm in a single day, was completely absent in the PD_AM series, while it was widely documented in metadata, and it was present in the observations from other meteorological stations.

Before proceeding with the calculation of drought indices, the most widely used homogeneity tests were applied to the whole series of the yearly precipitation amounts from 1725 to 2023. The same tests had already been applied to the last 30 years of the series, i.e., 1993–2022 [26]. These tests are the standard normal homogeneity test [47], the Buishand U and range tests [54,55], the Pettitt test [56,57] and the Von Neumann ratio test [58]. No change points were found; thus, the overall series can be considered homogeneous.

3.2. Calibration Dataset

Estimating a drought index requires a calibration step. The choice of the calibration dataset may influence the assessment of drought characteristics, in particular in regions with dominant drying trends [59]. It is quite common to use self-calibrated indices (i.e., using the same dataset for calibration and index estimation); nevertheless, in some studies, reference climate data are proposed to allow for intercomparisons among stations or periods [60].

Different calibration datasets were tested for the calculation of SPI12 and SPEI12 for Padua. To make the interpretation of the results easier, the index values calculated with 1961–1990 and 1994–2023 calibration datasets were subtracted from the values calculated with the whole 1725–2023 series. The indices were calculated with the R package SPEI [61], which uses the gamma and log-logistic distribution functions, respectively, for SPI and SPEI [8]. The most significant differences were for the “normal” class, i.e., when the index was in the range (−1;1] (Figure S2), which was the least significant in drought assessment. Concerning the extreme classes, the 1994–2023 calibration dataset provided negative differences, while the 1961–1990 one resulted in positive differences with respect to the whole dataset, both for SPI and SPEI. This means that according to the 1961–1990 calibration dataset, there were more extreme periods, both wet and dry, than observed using 1725–2023 dataset; the opposite was true for the 1994–2023 calibration dataset. However, the differences were within 5%.

Since the series is very long and characterized by different levels of accuracy over the various subperiods, while still being homogeneous, the most reasonable choice is to use the entire series as calibration dataset. Moreover, trends toward more extreme conditions are amplified when the calibration period does not include recent data, including the recent effects of climate change [62].

3.3. SPI and SPEI Best-Fitting Distributions

According to the methodology described in Section 2.2, only the distributions that satisfied the SW test (p-value > 0.05) were considered in the next step, i.e., the application of goodness-of-fit tests. The fitting distributions that satisfied the SW test and showed the best performance based on MSEN and AIC indicators (cf., Section 2.3) are reported in

Table 2. Best-fit distributions for each drought index and time scale according to the standard approach (SA) and “general” best-fit approach (BFAg).

Method	SPI3	SPI6	SPI12	SPEI3	SPEI6	SPEI12
SA	Gamma	Gamma	Gamma	Log-log 3	Log-log 3	Log-log 3
BFAg	Gamma	Gamma	Gamma	GEV	GEV	Skew-normal

and

Table 3. Best-fit distributions for each month, drought index and time scale according to the “monthly” best-fit approach (BFAm).

	Method: BFAm
Month	SPI3	SPI6	SPI12	SPEI3	SPEI6	SPEI12
January	Gamma	Gamma	Gamma	Pearson type 3	GEV	Skew-normal
February	Log-normal	Gamma	Weibull	Skew-normal	Skew-normal	Skew-normal
March	Log-normal	Gamma	Weibull	Pearson type 3	GEV	Skew-normal
April	Gamma	Gamma	Weibull	Skew-normal	Skew-normal	Skew-normal
May	Gamma	Gamma	Weibull	Pearson 3	Skew-normal	Skew-normal
June	Gamma	Gamma	Gamma	Skew-normal	Skew-normal	Skew-normal
July	Gamma	Gamma	Weibull	Skew-normal	Skew-normal	Skew-normal
August	Log-normal	Gamma	Normal	Log-log 3	Skew-normal	Skew-normal
September	Log-normal	Gamma	Gamma	Log-log 3	Skew-normal	Skew-normal
October	Log-normal	Gamma	Weibull	Pearson type 3	Skew-normal	Skew-normal
November	Gamma	Gamma	Weibull	GEV	Skew-normal	Skew-normal
December	Gamma	Gamma	Gamma	GEV	Skew-normal	Skew-normal

for the “general” and “monthly” best-fit approaches, respectively. Table 2 includes the distributions calculated by the standard approach.

Concerning SPI, the best-fitting distributions according BFAg and SA are the same, i.e., gamma, for each accumulation period. For SPEI 3 and 6, BFAg indicates GEV, while SPEI12 indicates skew-normal. Passing to BFAm, the results differ from month to month. In the case of SPI, in addition to gamma, log-normal and Weibull distributions appear, respectively, for SPI3 and SPI12. In the case of SPEI, for the 6- and 12-month accumulation periods, skew-normal is the best-fitting distribution with few exceptions, while results for SPEI3 are more heterogeneous and also include Pearson type 3 and Log-log 3.

3.4. Comparison of SPI and SPEI Values Estimated with Different Approaches

The ability of the different approaches to detect wet and dry conditions can be tested by analyzing several aspects. As a first comparison, the Pearson correlation coefficient ρ and the mean bias error (MBE) = $-I_{BFAg,m})/N$ were calculated, where I is the index value (SPI or SPEI) calculated with the different approaches, i.e., SA, BFAg and BFAm, and N is the number of months. In particular, ρ estimates the linear correlation between two variables, while MBE estimates the degree to which the mean of a variable differs from that of another. All the ρ values were higher than 0.99, indicating a high correlation between all three approaches, both for SPI and SPEI. MBE values are reported in Figure 4a, excluding the ones for SPI, as they were not significant, being all lower than 0.005. Concerning SPEI3 and SPEI6, BFAg tended to estimate lower values than SA; the opposite was true for BFAm, with the difference increasing with the accumulation period. For SPEI12, BFAg matched with BFAm.

To further investigate the differences between the approaches, each SPI and SPEI values was assigned to a specific class according to Table 1. Then, for each class, the relative error was estimated: $RE\left(c\right)=({ne}_{BFAg,m}\left(c\right)-{ne}_{SA}\left(c\right))/{ne}_{BFAg,m}\left(c\right)$, where c is the class, and ne (c) is the number of events detected by SA and BFA for each class. The results concerning BFAg and BFAm were averaged, as they were very similar, and considered a single approach, named BFA for simplicity (Figure 4b). The most significant errors were related to the extreme classes: BFA detected more extreme drought events than SA; the opposite was true for the extremely wet events (Figure 4b).

The comparison between the different approaches was completed with the analysis of the drought characteristics defined in Section 2.4, i.e., duration, severity, intensity and frequency. The results are summarized in

Table 4. Main drought characteristics, i.e., duration (D), severity (S), intensity (I) and frequency (F), according to SPI calculated with the standard (SA), “general” (BFAg) and “monthly” (BFAm) best-fit approaches for k = 3, 6 and 12 months. The results for SA and BFAg are the same.

		SA = BFAg			BFAm
k		3	6	12	3	6	12
D (month)	Min	1	1	1	1	1	1
	Mean	5.0	8.3	12.5	5.4	8.3	12.3
	Max	27	36	68	28	36	68
	sd	3.1	5.9	10.2	3.3	5.7	9.7
S	Min	−26.9	−43.8	−108.7	−27.4	−45.1	−100.4
	Mean	−5.6	−9.1	−14.6	−6.0	−9.2	−13.5
	Max	−1.7	−1.3	−1.6	−1.2	−2.0	−2.4
	sd	3.3	7.2	16.7	3.5	7.3	14.9
I (month−1)	Min	−2.3	−2.2	−1.9	−2.2	−2.3	−1.8
	Mean	−1.1	−1.1	−1.0	−1.1	−1.1	−1.0
	Max	−0.6	−0.3	−0.5	−0.4	−0.5	−0.4
	sd	0.3	0.3	0.3	0.3	0.3	0.3
F		239	136	78	226	136	76

and

Table 5. Main drought characteristics, i.e., duration (D), severity (S), intensity (I) and frequency (F), according to SPEI calculated with the standard (SA) and best-fit approaches (BFAg, BFAm) for k = 3, 6 and 12 months.

		SA			BFAg			BFAm
k		3	6	12	3	6	12	3	6	12
D (month)	Min	1	1	1	1	1	1	1	1	1
	Mean	3.9	4.9	12.7	3.9	4.9	13.3	5.4	8.0	13.3
	Max	8.0	8.0	41.0	7.0	8.0	41.0	27.0	27.0	41.0
	sd	0.9	1.2	7.8	0.9	1.1	7.9	3.2	5.1	7.9
S	Min	−11.4	−16.9	−56.0	−10.4	−15.9	−50.1	−27.3	−37.8	−50.2
	Mean	−5.3	−6.1	−14.7	−4.8	−5.7	−14.3	−6.0	−8.6	−14.3
	Max	−1.8	−1.5	−1.8	−1.5	−1.8	−2.1	−1.2	−1.8	−2.1
	sd	2.0	2.9	12.3	1.8	2.9	11.1	3.9	6.6	11.2
I (month−1)	Min	−2.6	−2.3	−2.1	−2.4	−2.1	−1.9	−2.3	−1.9	−1.9
	Mean	−1.4	−1.2	−1.1	−1.2	−1.1	−1.0	−1.1	−1.1	−1.0
	Max	−0.6	−0.5	−0.4	−0.5	−0.6	−0.5	−0.4	−0.4	−0.5
	sd	0.4	0.4	0.4	0.3	0.3	0.3	0.4	0.3	0.3
F		268	214	81	250	192	75	221	139	75

for SPI and SPEI, respectively.

Coherently with the previous findings, the main statistics for SPI are very similar, regardless of the approach used (Table 4). In contrast, passing to SPEI (Table 5), durations and severity are higher for BFAm than for the other approaches for a 3-month accumulation period. Moreover, the number of droughts according to SPEI3 and SPEI6 follows the scale SA > BFAg > BFAm, while for SPEI12, the frequency is similar (Table 5).

3.5. Comparison of SPI and SPEI Values Estimated via “Monthly” Best-Fit Approach

The comparison between the two drought indices was performed considering the best-fit approach applied month by month, i.e., BFAm. Results can be better interpreted by looking at the time series of the differences between SPI and SPEI (Figure 5). The SPI and SPEI time series, as well as their 30-year moving averages, are reported in Figures S3 and S4 of the Supplementary Material, respectively. SPI and SPEI values were limited to vary between −3 and 3, as values outside this range are characterized by high uncertainty [10]. Values below −3 and above 3 were assigned −3 and 3, respectively. The 3- and 12-month accumulation period records were selected, as they represent seasonal and annual series, respectively. Specifically, seasonal series were defined by considering SPI3 (SPEI3) values in February for winter (December–January–February), May for spring (March–April–May), August for summer (June–July–August) and November for autumn (September–October–November), while the annual ones were set by selecting the SPI12 (SPEI12) values in December of every year.

Looking at Figure 5a, there is a remarkable negative peak in 1740 and a clear increasing trend starting in the mid-20th century, both due to the temperature, which factors into the calculation of SPEI through water balance as the input variable. In fact, 1740 was the coldest year in Central Europe in 600 years [63], and it was the coldest year of the 1725–2023 Padua series, −3.7 °C lower than the 1901–2000 climatological mean. The coldest winter (Figure 5d) and coolest spring (Figure 5c) also occurred in 1740 [33]. Moreover, in the 1970s, a significant warming trend started in Padua [33]. The increase is part of ongoing climate change, of which global warming is a well-known aspect. To provide a few numbers, the annual mean temperature of the air in the Mediterranean region is +1.54 °C higher than the 1860–1890 preindustrial level for land and sea areas, i.e., +0.4 °C more than the global average [2]. Over the 20th century, climate reconstructions, ground-based observations, reanalyses and remote-sensing datasets all corroborated the transition to warmer conditions and showed that warming has accelerated during the last decades, with significant positive trends of the order of +0.1 °C to +0.5 °C per decade [2].

According to BFAm, the main differences between SPI and SPEI statistics are related to the longest and most severe drought for the 12-month accumulation period, the duration and severity of which are reported in

Table 6. Duration (D), severity (S) and covered period of the most severe droughts from 1725 to 2023 according to SPI and SPEI, calculated with the “monthly” best-fit approach for k = 3 and 12 months.

Most Severe Drought over the 1725–2023 Period
	SPI3	SPEI3	SPI12	SPEI12
D (month)	28	27	68	41
Covered period (m/y)	11/1892–2/1895	11/1892–1/1895	2/1737–9/1742	10/2019–12/2023
Severity	−27.4	−27.3	−100.4	−50.2

. The comparison of the periods in which these extreme events occurred confirms the previous result. In fact, the driest period for SPI12, i.e., from 2/1737 to 9/1740, includes the coldest year of the whole series. At the same time, according to SPEI12, the most severe drought occurred in the last decades, characterized by rising temperature at the global scale. During the most severe drought according to SPI12, there were three drought events for SPEI12: 1) from 2/1737 to 6/1738, D = 17, S = −24.3; 2) from 7/1739 to 5/1740, D = 11, S = −8; and 3) from 10/1740 to 5/1742, D = 20, S = −17.4. In turn, around the period of the most severe drought for SPEI12, there was the second most remarkable drought for SPI12 of the whole period, from 5/2020 to the end of 2023, with a severity of −52.6. Results for k = 3 are in accordance with the two indices, in contrast (Table 6).

3.6. Yearly and Seasonal Trends Investigated with Traditional and Innovative Methods

Before investigating the presence of trends, some homogeneity tests were applied to SPI and SPEI time series calculated with the “monthly” best-fit approach to investigate the presence of climate signals. As in Section 3.5, the 3- and 12-month accumulation period records were used for seasonal and annual series, respectively. The results are summarized in

Table 7. Change points for SPI and SPEI series.

Series	Change Point * (year)
SPEI3-JJA	1997
SPEI12	2014

* p-value < 0.05.

. In 1997, a change point was observed for SPEI3 in summer (JJA), while this occurred in 2019 for SPEI12. Long-term trends and significance were evaluated by means of Theil–Sen (TS) [44,64] and MK tests over the entire period. The MK test yielded significance only for SPEI12 (p-value < 0.05), and the TS test indicated a decreasing slope of −0.2 per decade.

Among several restrictions related to the MK test (Section 2.5), the existence of positive serial correlation in the time series increases the probability of detecting a false trend [44]. Hence, innovative trend analysis was also performed to lend robustness to the results. According to ITA, SPEI12 showed a general decreasing trend (Figure 6b), as the data were below the ideal 1:1 line. For SPI12, the decrease was notable for the values belonging to the “severe” classes (Figure 6a), i.e., at the extremes of the range values (i.e., colored rectangles). A decrease in the positive values indicates a decrease in the severe wet events. In contrast, a decrease in the negative values indicates an increase in severe droughts, as the more severe a drought is, the larger the value of the index is in absolute value and the smaller in relative value. The results of the same analysis carried out at seasonal level are represented in Figure S5 and summarized in

Table 8. Trends in SPI3 and SPEI3 series for severe droughts. Symbols “+”, “−“ and “0” indicate increasing, decreasing and no trends, respectively.

Season
	Winter	Spring	Summer	Fall
SPI3	+	+	−	−
SPEI3	0	+	−	−

for “severe” droughts. A decrease in SPI3 and SPEI3 values, i.e., an increase in the severity of drought events, was found in summer and fall.

The CWT was applied to identify the dominant and significant periodic components in SPI12 and SPEI12. As the traditional approach leads to a consistent reduction in power at low periodicities, the new biwavelet version was used. It implements the bias-correction method developed by Veleda et al. [65], which does not affect the wavelet coherence. Results for SPI12 and SPEI12 were similar; therefore, Figure 7 shows only CWT applied to SPEI12. The most powerful periodic components that emerge from the scalogram in Figure 7a (red color) are more clearly identified by the spectrum in Figure 7b. Four periodic components were above the 99^th confidence limit, and they corresponded to the return periods of 1.7, 2.3, 3.3 and 4.6 years. In contrast, the 27-year RP component was not significant. Therefore, only short-term periodic components were found as significant in the SPEI12 series.

3.7. SPI and SPEI Critical IDF Curves

Critical drought IDF curves were produced for SPI and SPEI obtained with BFAm when k = 3 and 12 months and six fixed return periods (i.e., 5, 10, 25, 50, 100 and 300 years), following the procedure described in Section 2.6 (Figure 8).

For each drought duration D, a relationship was established between critical drought intensity and return period. The best approximation resulted in an exponential equation:

$$RP=e^{{\displaystyle \frac{I-b}{a}}}$$

As an example, regression equations at the time scales k = 3 and 12 months are visualized in Figure S6 for D = 12.

In addition, for D from 1 to 12 months, RPs of the most severe droughts were calculated for k = 12 and are reported in

Table 9. Return periods (RPs) and related year of occurrence of critical droughts of duration D for k = 12.

Index	SPI12		SPEI12
D (months)	RP (years)	Year	RP (years)	Year
1	231	1737	306	2022
2	421	1922	330	2022
3	336	1922	342	2022
4	331	1922	348	2022
5	272	1893	360	2022
6	265	2022	330	2022
7	322	1922	274	2022
8	280	1922	263	2022
9	251	1922	260	2023
10	224	1922	318	2023
11	304	1922	337	2023
12	304	1922	346	2023

together with the year they were recorded.

In general, the RPs provided by SPEI12 are higher than the ones calculated using SPI12 for the same duration D. Moreover, SPI12 indicates 1922 as the year in which critical droughts of almost all durations occurred, while SPEI12 indicates the year 2022. This result can be explained by the different nature of the related indices, as SPEI is able to capture a broader measure of the available water through the water balance as an input variable, as previously discussed. Looking at SPI12 values, a drought event started in August 1921 and lasted 22 months, with a severity of −36.6. Within this period, the most critical droughts were recorded regardless of the duration. This was not the case for SPEI12, according to which the most intense event for each duration occurred in 2022/2023. In fact, in summer 2022, several areas of the Northern Hemisphere experienced exceptional droughts, extreme heat and a dire water shortage [1]. The critical situation continued and worsened. The driest event indicated by SPEI12 (Table 9) ended in 2023 just because the series under study ended in 2023. In reality, the exceptional drought that started in 2019 continued until April 2024, with a duration of 55 months and a severity that was twice as great as that indicated in Table 9, i.e., −92.

4. Conclusions

In this paper, the drought variability in Padua over nearly 300 years was investigated through the study of two well-established indices used to evaluate meteorological drought worldwide. SPI and SPEI have the advantages of being compatible with series starting in the early instrumental period, when only daily mean temperature and precipitation amount are the available variables. As a part of this study, the daily precipitation series of Padua from 1951 to 2023 was completed.

For calculating the drought indices, two alternative approaches to the standard one were considered, both based on the identification of the best-fitting distribution: one is “general”, i.e., it considers the same probability distribution function (PDF) for all the months (BFAg); the other uses a different PDF for each month (BFAm). Concerning SPI, the best PDF according to BFAg and SA was the same, i.e., gamma, for each accumulation period. For SPEI3 and SPEI6, BFAg indicated GEV, while SPEI12 indicated skew-normal instead of the standard log-logistic 3. When passing to BFAm, the results depended on month and accumulation period. In the case of SPI, in addition to gamma, log-normal and Weibull distributions appeared, respectively, for SPI3 and SPI12. Skew-normal was the best PDF for SPEI6 and SPEI12 with few exceptions, while the results for SPEI3 were heterogeneous.

Comparing the standard to the best-fit approaches, the mean bias error was significant only for SPEI. In particular, in case of SPEI3 and SPEI6, BFAg tended to estimate lower values than SA; the opposite was seen for BFAm, with the difference increasing with the accumulation period. For SPEI12, BFAg matched with BFAm. The results for the relative error in associating an index value with a particular class were the same for BFAg and BFAm. With respect to SA, the most significant errors were related to the extremes: the BFAs detected more drought events belonging to the extreme class than the SA. Coherently with the previous findings, the main statistics for SPI were similar, regardless of the approach type. In contrast, the duration and severity of droughts estimated using SPEI3 were higher for BFAm than for the other approaches. Moreover, the drought frequency according to SPEI3 and SPEI6 followed the scale of SA > BFAg > BFAm, while for SPEI12, values were similar.

The comparison between the two drought indices was performed considering both the best-fit approach applied month by month, i.e., BFAm, and the 3- and 12-month accumulation periods as representative of seasonal and annual time scales. The difference in SPI and SPEI time series showed an increase since the mid-20th century, in particular in spring and summer, and can be related to ongoing global warming. Moreover, the longest and most severe drought over the 1725–2023 period identified by the two indices was not the same: according to SPI12, it was from 1737 to 1742; according to SPEI12, it was from 2019 to 2023. This last event continued until the first months of 2024.

The innovative trend analysis indicated a general increase in drought trend, which according to SPI12, was significant only for severe events. Summer and fall were the most affected seasons.

Drought intensity–duration-frequency (IDF) curves provide a straightforward graphical representation of the relationship between intensity, duration and frequency that stakeholders, decision-makers and water resource managers can easily understand. In particular, for each duration, it was possible to calculate the return period of the critical droughts, i.e., the most severe ones. In general, the return period provided by SPEI12 was higher than the one calculated using SPI12 for the same duration. Moreover, SPI12 indicated 1922 as the year in which critical droughts of almost all durations occurred, while SPEI12 indicated the year 2022.

As a general remark, it is recognized that the use of the best-fitting distribution approach (BFA) may add a level of complexity and complicate comparisons across space and time. Nevertheless, in the case of a specific location for which a long, high-resolution meteorological series is available, BFA is preferred to the standard approach, in particular when investigating critical events. In addition, in the context of ongoing climate change, even if precipitation is the main driver of droughts, the role of temperature cannot be neglected, especially in the Mediterranean area. The results of the present study confirm that SPEI deserves adequate research attention as a reasonable and useful drought index.

The evaluation of a drought index for such a long series, which starts in the early instrumental period, i.e., at the beginning of the 18th century, has inherent limitations due to the nature of the input data, which may affect the results. First, the accuracy of the meteorological observations is a crucial point: the further back in time one goes, the more careful reconstruction work and critical analysis must be to compose a homogeneous series and correct biases. Second, the lack of observations for several meteorological variables, in particular in the early instrumental period, when generally only mean temperature and precipitation were recorded, limits the possibility of using methods to calculate PET and thus to estimate SPEI. Last but not least, a long series such as that for Padua is an exception, so it is not possible to find series from neighboring locations to use for comparison.

Finally, this study provides interesting insight into drought characterization in northern Italy over a long-term perspective. The availability of information on drought evolution over the centuries can support risk assessment and mitigation plans, helping to cope with one of the most damaging climate events of the last decades.