Drought Characterization in Croatia Using E-OBS Gridded Data

1. Introduction

Droughts are still among the least understood and complex extreme climate events, affecting large worldwide areas and having serious impacts on the society, the environment, and the economy [1]. The Convention to Combat Desertification from the United Nations 15th Conference of Parties (COP15) states in its 2022 report [2] that globally the information points in the same direction, which is an increase in the duration and severity of droughts and their impacts affecting both human societies and ecological systems.

Due to the recent more pronounced climate variability, drought patterns have been changing, becoming more difficult to follow their evolution, which leads to a constant need to update the understanding about their behaviour and characteristics. To better cope with drought it is important also to distinguish concepts between drought and water scarcity, the two major phenomena that challenge water security [3]: while drought is a natural climate event caused by climate variability that cannot be avoided by water management practices; water scarcity refers to the long-term unsustainable use of water resources, where consumption exceeds water availability and where water management practices can have their influence [4].

To characterize the drought phenomenon and occurrence in Croatia and to understand the historical and recent climate variability in the country, it is important to study the long-term time series of precipitation and temperature since they portrait the behaviour of two of the most important climate variables. Around the world, several authors have characterized the phenomenon of drought using well-known scientific indices, such as the Standardized Precipitation Index (SPI) [5] and the more recent Standardized Precipitation Evapotranspiration Index (SPEI) [6], and by studying the spatio-temporal patterns of droughts, their duration and frequency.

In this regard the pioneer work of [7] for the Iberian Peninsula show that using the SPI at a 12-month temporal scale in 51 monthly precipitation time series within a period from 1910 to 2000 allowed the identification of four different regions showing internal homogeneity and well-defined boundaries (north, northeast, southeast and central/western). Specifically for Portugal and using multiple SPI time-scales, [8] have analysed droughts based on monthly precipitation data from a meteorological network concluding that the country could be divided into three homogeneous regions (North, Center and South), the same spatial results confirmed by [9] using gridded data and the same drought index in 2015. More recently, also based on SPI computed on rainfall grid data for a time window of 100 years, [10] studied the regionalized droughts in three regions of Portugal concluding that the frequency of the moderate and severe droughts is increasing. Also in Southern Europe, [11] studied the spatiotemporal drought variation using clustering techniques and the SPEI6 based on gridded datasets allowing the identification of five well-distinguished regions in the south of Italy characterized by drought events different in terms of duration and severity. In China, Xinjiang Province, a study using the Mann–Kendall trend test, cluster analysis and Morlet wavelet analysis concluded that historical drought evolution is better characterized by three regions in close connection with the topography and climate [12]. Also the recent work of [13] is a good example of climate regionalization by decomposing a field into its spatial and temporal terms, allowing the study of the variability of drought in Mozambique.

For the period 1950-2012, in a study about European drought climatologies and trends, [14] found that in Central Europe and the Balkans the precipitation increase was not significant, but the temperature rise and hence of the potential evapotranspiration (PET), have driven the drought severity to increase. Especially in the last decade, the occurrences of drought have been increasing in the Balkans region so that in 2018-2019, winter precipitation deficits coupled with high record temperatures in the summer across much of the Balkan Peninsula, i.e., Croatia, Albania, Slovenia, North Macedonia, Montenegro and Hungary, as well as in Slovakia led to an uncommon hydro-climatological conditions that clearly differed from the rest of Europe [15].

Especially relevant for Croatia and according to the Center for Climate and Security [16] droughts pose the largest threat to the Balkans region’s stability and prosperity, especially through their impacts on the agriculture and hydroelectric sectors. The Western Balkans are also a hotspot of biodiversity in Europe and contain a great diversity of ecosystems, which can be explained by their extremely diverse geology, soil, climate ranges and topography, being one of the six European centres of biodiversity [17].

Some studies have been conducted in Croatia regarding drought analysis, namely the works of [18] applying principal component analysis to a drought index, or [19] which compared the SPI, SPEI, and the standardized groundwater index (SGI) in different time-scales (1, 3, 6, 12, 24, and 48 months). Drought indices (SPI and Palmer Drought Severity Index, PDSI) were also compared for one meteorological station in North Western Croatia and the agricultural impacts were assessed [20].

In this study, a comprehensive characterization of droughts over Croatia using the SPEI is achieved with a combined classification methods of Principal Component Analysis (PCA) and Cluster Analysis (K-means). The temporal characteristics of drought variability, the changes in the yearly occurrences rate, the drought areal evolution and the characterization of trends via robust and commonly used techniques within the climatology and hydrology fields is examined.

While different studies have been conducted on the Croatian continuum, a long-term analysis aimed at identifying spatial and temporal patterns of drought has never been undertaken. This limitation has largely been due to the unavailability of ground records of climate variables, as the number of meteorological stations with continuous long time series is relatively small. Therefore, the primary contribution of this research is to address this constraint by utilizing a highly dense gridded meteorological dataset and employing well-established, validated methodologies to investigate the drought phenomenon across the entirety of Croatia. Importantly, this study represents the first instance in which complete long-term data have been used for this purpose.

The results are intended to be of practical value for planning activities, water resources management, and stakeholders, providing insights that could be used to create risk-based drought management plans at a regional scale but also to understand the connection between drought risk and climate variability.

2. Materials and Methods

By using gridded climate data within a 73-year period, from September 1950 to July 2022, the drought characterization used the following research steps and complementary methods:

i)

Data validation. E-OBS precipitation and temperature datasets are first validated with records from meteorological stations at different locations of Croatia.

ii)

Drought index calculation. SPEI with a 6- and 12-month time-scale (SPEI6 and SPEI12) are calculated ascertaining the sub-annual and annual temporal variability of droughts.

iii)

Drought regional patterns. Principal Component Analysis (PCA) was applied to the previous SPEI time series aiming at identifying homogeneous regions. The K-means clustering method (K-means) was used to validate the regions identified from the PCA.

iv)

Temporal evolution of drought areas. Drought areal evolution of the SPEI6 and SPEI12 fields in each of the identified regions is achieved by assigning an area of influence to each grid cell.

v)

Yearly frequency analysis of drought occurrences. A kernel occurrence rate estimator (KORE) is used to analyse the yearly frequency of the periods under drought conditions for different drought categories according to the regionalized SPEI time series given by the factor scores previously obtained by the PCA.

vi)

Trend analysis. The Modified Mann-Kendall (MMK) trend test coupled with the Sen’s Slope estimator test is used to detect the temporal variability of drought intensities within the SPEI6 and SPEI12 fields in each of the regions.

The framework used for the application of the algorithms was developed using the R Software (RStudio Team, 2023), namely by considering the following principal packages: "SPEI" [6], “psych” [21], “modifiedmk” [22], “stats” (developed by the R Core Team and contributors worldwide) and others.

2.1. Study Area

The extensive drought analysis considered mainland Croatia (Figure 1a), with approx. 56 594 km², as the study area. The country is located in the middle of the Central and Southeast Europe transition (between latitudes 42° and 47° N and longitudes 13° and 20° E), faces the Adriatic Sea in the southern part and has a complex topography, a complex and variable climate, and diverse terrestrial ecosystems.

The major natural factors that condition the climate of the country are the latitude, the topography, ranging from approx. 0 to 1 831 meters above sea level, m.a.s.l. (the highest mountain in Croatia, the Dinaric Alps), and the proximity with the Adriatic Sea to predominant Mediterranean influence [23].

The majority of Croatia is characterized by a moderately warm and rainy climate, but it has experienced significant spatial heterogeneity in precipitation with large variability across the country, over the last decades [24]. The mean annual precipitation in Croatia is approximately 1 100 mm [25] (Figure 1b), ranging from about 3 900 mm on the summits of the southern Velebit mountain located along the northern Croatian Adriatic coast to about 300 mm on the outlying islands in mid-Adriatic, in close connection with the relief. Mean annual temperature in the lowland area of northern Croatia is 10°C–12°C, the mountain regions experience mean temperatures of 3°C–4°C, and coastal areas experiencing temperatures of 12°C–17°C [23]. In the last 20-30 years the average annual temperature has increased considerably over Croatia, where its anomalies have reached + 1.8°C (Figure 1c), a consistent warming trend in the whole country already confirmed in recent studies [26]. The precipitation and temperature anomalies in Figures 1b and c represent deviations in relation to the average of the reference period of 1950–2021.

2.2. E-OBS Data

Climate variables, namely precipitation and temperature, used to calculate the SPEI over a long base period 73-year period from September 1950 to July 2022, to correctly sample the natural variability, were obtained through the Copernicus Climate Change Service (C3S) program, operated on behalf of the European Union, by the European Centre for Weather Forecasting (ECMWF). This program combines data from observations of the climate system with recent data and provides high quality information on the past, present and future state of the climate for Europe and the World¹. The product used was the "E-OBS daily gridded meteorological data for Europe from 1950 to 2022 derived from in-situ observations"². This is a daily product for Europe obtained from observations at stations of the European National Meteorological and Hydrological Services (NMHSs) or other institutions, available on a regular grid of 0.25°x 0.25° [27], being represented in Figure 1a the location of the center points of the respective grid.

The monthly time series needed to compute the SPEI for Croatia were generated from the daily temperature and precipitation E-OBS grid datasets by averaging or aggregating daily data into a monthly time-scale, respectively. The resulted grid spatially covered Croatia with 103 grid cells all within the Croatian border. The 0.25°x 0.25° spatial resolution means a vertical distance that varies with latitude from 27 660 m in the grid cell of the most southern point to 27 814 m in the topmost northern grid cell. In order to validate the E-OBS grid dataset, the meteorological station data from Croatian Meteorological and Hydrological Service (DHMZ) was used and the monthly boxplots between each station and the nearest grid center point were obtained for both climate variables. For this purpose, a common period with available ground data, spanning from 1981-2022, was adopted (Figure 2).

The chosen meteorological stations were those that represent Croatian climate zones: the cool temperate moist in the mountainous regions, close to the Adriatic Sea (Gospić); the warm temperate moist of the northern central region (Bjelovar); and the warm temperate dry of the northeastern region (Gradište). The grid center point (E-OBS cells) that was assigned, for comparison purposes, to each of the previous stations was the nearest one, which resulted in a maximum distance of 8 000 m from the coupled station. The respective elevations from each data source location were also quite similar, reinforcing the comparison between datasets.

From Figure 2, it is clear that both datasets showed similar behaviour through all the months, namely by comparing the central tendency median but also similar monthly variability by the equivalent interquartile range. The monthly mean temperatures show the normal yearly distribution of a temperate northern hemisphere country, with the hotter months being those of summer and the colder ones occurring in winter. The within-the-year temporal pattern of the precipitation is always very smoothed. The highest values relate to November in Gospić, September in Bjelovar and June in Gradište, for the period 1981–2022.

The Pearson correlation coefficients for the sequential monthly time-series of the meteorological stations and the respective nearest grid center points were obtained in order to complement the information of the boxplots, being 0.98, 0.86 and 0.87 for precipitation in Gospić, Bjelovar and Gradište, respectively, and 0.99 for the mean temperatures in any of the stations.

2.3. Drought Index

The Standardized Precipitation Evapotranspiration Index (SPEI), developed by [6] was adopted as the climate drought index to assess the drought in mainland Croatia. The SPEI has become one of the most popular drought indices for drought monitoring and characterization [13,15,30] and is also recognized by the World Meteorological Organization (WMO) as an eligible index for those purposes [29]. It was derived from the original Standardized Precipitation Index (SPI) [5] by accounting for evapotranspiration, and it uses the water balance (WB= P – PET) time series, where P accounts for precipitation and PET the potential evapotranspiration, instead of precipitation alone, as in the SPI.

As the original formulation of the SPEI suggested, the PET method here considered was the Thornthwaite model, because of the simple data requirements for its computation which is based uniquely on monthly mean temperatures and latitude of such meteorological observations [30]. The SPEI has several advantages, such as (1) great flexibility, as it can be applied at different time-scales, allowing the representation of the response of crops, natural vegetation and hydrological systems to drought conditions; (2) less complexity, comparatively to other indices, as it requires relatively simple and well set calculations; (3) great suitability to spatial representation, allowing comparison between areas within the same region, as it is a normalized index; and, (4) in comparison with SPI, the SPEI takes into account air temperature as the most prominent variable of the water budget at the watershed level and of climate change.

The model uses a probability distribution function to fit either the monthly WB series or the aggregated WB over n months, for SPEI time-scales higher than 1 month, then those probabilities are standardized in z units (mean=0, standard deviation=1), by a equiprobability transformation and the SPEI time series are obtained. According to [6] the log-logistic distribution here adopted has provided better results than other distributions for obtaining SPEI series and has been widely used in many different drought research frameworks and countries [15,33,34,35].

The log-logistic distribution for a given variable x is given by:

$$\mathit{F}\left(\mathit{x}\right)={\left[1+{\left(\frac{\mathit{\alpha }}{\mathit{x}-\mathit{\gamma }}\right)}^{\mathit{\beta }}\right]}^{-1}$$

(1)

where α, β, and γ represent the scale, shape and location parameters that are estimated from the x sample (in this case WB).

When computing the SPEI series the drought categories adopted for analysis follow the ones proposed by [34], Table 1.

The time-scale of SPEI could be separated into three groups: short (2–4 months), medium (6–10 months) and long (10–20 months) [35]. Normally short to medium SPEI time-scales could be used for analysing the response of crops and natural vegetation to drought, and with longer the index could capture water resources availability and hence used to evaluate hydrological variations. In order to characterize the historical drought conditions in Croatia a medium (6 month) and longer (12 months) time-scales were chosen which allowed to analyse the subannual and annual variability of drought conditions (SPEI6 and SPEI12).

2.4. Drought Regional Patterns

To characterize the spatial and temporal drought patterns, principal component analysis (PCA) and nonhierarchical clustering (K-means) were applied to the SPEI time series, as considered by many other authors for drought regionalization purposes [9,10,12,37].

The PCA allows the classification of data in such a way that the spatial patterns identified for climate data could be expressed in order to highlight their similarities and differences [37]. The original intercorrelated field of SPEI could be reduced to small groups of new SPEI grid center points that are linearly uncorrelated and that explain most of the total variance. The main advantages of PCA [38] are: (1) PCAs are not affected by the lack of independence in the original variables; (2) normality of the data is preferred but not essential; and (3) only an excessive number of zeros could cause problems, which within the nature of SPEI is not a concern.

The number of PCs identified by the scree plot method to retain [39] should explain at least 75% of the accumulated variance of the SPEI field. With the main PCs selected, it is possible to produce a more stable spatial patterns by rotating the PCs with the Varimax procedure, providing a clearer division between components by the redistribution of the explained variance, preserving their orthogonality and producing more physically explainable patterns [41,42]. The new elements obtained are referred to as rotated principal components (RPCs).

The results of PCs can be based on the eigenvalues, the correlations between PCs and the original variables (factor loadings), or the observation coordinates in the PC which are linear combinations of all the grid point series (factor scores). In order to identify spatially disjunctive areas within the gridded SPEI field, the mapping of the correlations between the RPCs retained and the original SPEI series were also used to validate the classification obtained, confirming the selection of the optimal number of PCs. Accordingly, whenever a group of grid center points had the highest correlation with an RPC, a new region was delimited.

In order to validate the regional classification obtained by the PCA, a cluster analysis using the k-means method [9,10,14] was also applied to see if the number of regions obtained with K-means is equivalent to those obtained by the PCA. The K-means method, similarly to the PCA, has the ability to divide the dataset into homogeneous and distinct groups with members with similar characteristics [42].

The optimal number of clusters was evaluated by the Euclidean distances between the created clusters, which yields the lowest possible number with the greatest possible homogeneity, taking into account that the analysis requires that the number of clusters should be established beforehand. To overcome this aspect, which is considered one of the major unresolved issues in the cluster analysis, the test was repeated by considering several number of clusters, with the Euclidean distances being analysed to ensure maximum heterogeneity between groups [38].

2.5. Drought Yearly Occurrence Rate and Trend Analysis

The analysis of changes in the temporal occurrence rate of droughts is a very important aspect when characterizing their persistence and frequency. To evaluate the frequency of the occurrence of the periods under drought conditions over time a kernel occurrence rate estimator (KORE) was applied to a historical series of drought events with the aim of estimating how the mean yearly number of drought periods, λ, changes over time, that is, to characterize λ (t). Several authors have also used it for similar purposes, such as [44,45,46,46].

The kernel technique is a nonparametric model developed by [47] for smoothing point process data, in the case of the present application, the times of occurrence of the periods under drought conditions. It can be formulated as:

$$\widehat{\mathit{\lambda }}(\mathit{t})={\mathit{h}}^{-1}\sum _{\mathit{i}=1}^{\mathit{m}}\mathit{K}\left(\frac{\mathit{t}-{\mathit{T}}_{\mathit{i}}}{\mathit{h}}\right)$$

(2)

where $\widehat{\mathit{\lambda }}\left(\mathit{t}\right)$ is the estimated occurrence rate in each instant t, m is the number of events, K is the kernel function, h the bandwidth, and (t-T_i) the difference between the instant t and the i^th occurrence. In the applications carried out, a Gaussian kernel was used, according to [43].

In hydrological time series, trend analysis is an important and popular tool for better understanding the effects of climate variation and anthropogenic influences. In the present paper the temporal variability of droughts regarding the occurrence of trends was addressed based on the Modified Mann-Kendall (MMK) [22] test coupled with Sen’s slope estimator test [48] applied to the SPEI time series. The basic Mann-Kendall test [49] is a rank-based non-parametric trend test that requires data to be independent and randomly ordered; however, the SPEI inherently shows by its nature, positive autocorrelation in the time series, which becomes stronger as the time-scale of the SPEI increases.

The Modified Mann–Kendall test (MMK), on the contrary, can in fact address the issue of serial correlation through the use of the variance correction approach as first introduced and explained by [22]. According to [50] the ability to eliminate the influence of autocorrelation on test significance is the main advantage of using the MMK test over the basic MK test.

3. Results

3.1. Spatial Drought Patterns

The E-OBS grid dataset with 73 years (from 1950 to 2022) of monthly precipitation and mean monthly temperatures obtained from daily data, whose coordinates are shown in Figure 1a, allowed the calculation of the series of SPEI6 and SPEI12, with lengths of 871 and 865 months, respectively. For each of the time-scales, the SPEI values are comparable as they represent normalized values.

The PCA analysis was then applied to the 103 grid time series of SPEI at the time-scales of 6 and 12 months in order to regionalize the SPEI field into regions with different drought temporal behaviours. The amount of variance explained and the scree plot based on the Eigenvalues are shown for the first 10 principal components (D values on xx axis), on Figure 3, with the threshold adopted for selection of the optimum number of PCs.

From Figure 3 and SPEI6 the D1 explained a large percentage of the total variance (78.1%) while D2 and D3 explained 8.5% and 4.5%, respectively. According to the proposed criteria, the threshold adopted retained three PCs, with a total of 91.1% of the variance explained, that were then selected for the Varimax orthogonal rotation, and the D4, explaining only 2.30% was left out. Three new rotated principal components (RPCs) were obtained with 28.8 %, 27.6% and 34.7% of total variance explained, for RPC1, RPC2 and RPC3, respectively. For SPEI12 the D1 explained 75.6% and D2 and D3 explained 9.5% and 4.4% of total variance, respectively. The three rotate PCs explained a total variance of 89.5%, which resulted in new RPCs with 29.0%, 36.4% and 24.1%, respectively for RPC1, RPC2 and RPC3.

The mapping of the factor loadings obtained for the RPCs is shown on Figure 4. The interpolation method used to map the correlations, between the three chosen RPCs and each of the 103 SPEI series, was the Inverse Distance Weighting (IDW) available on QGIS version 3.16 Hannover, which computes the Inverse Distance to a Power gridding combined with the nearest neighbour method. For SPEI6, the contour line for a significant correlation of r = 0.6 allowed the identification of three different regions spatially disjunctive and without overlapping defined as: the Central North region (RPC D1), Eastern region (RPC D2), and Southern region along the Adriatic Coast (RPC D3), as represented in Figure 4a–c. Identical results were obtained for SPEI12, validating the initial adoption of three PCs to rotate (Figure 4d to f).

Regionalization of SPEI6 and SPEI12 presented in Figure 4 could be a consequence of the climatological characteristics of the continental part of the country – decrease of precipitation from West to the East and light increase of mean annual air temperature in the same direction. However, the southern coastal region has a typical Mediterranean climate, mainly characterized by the lowest precipitation occurring during the warm period of the year, high summer temperatures and a much higher mean annual air temperature. With the objective of validating the classification results of the PCA, regarding the three regions obtained, a non-hierarchical cluster analysis was applied through the k-means method and the optimal number of clusters for the SPEI field was also obtained.

To evaluate the K-means classification, the Euclidean distances between clusters were examined to ensure that heterogeneity between clusters prevailed, and so a clustering of 2, 3, and 4 group regions was tested. The spatial representation of the resulting cluster test is shown in Figure 5. According to the Euclidean distances for SPEI6 (Table 2) the classification with 2 groups (Figure 5a) is considered inadequate (0.546), since when grouped into 3 clusters (Figure 5b) the distance between the centroids of the new cluster 2 and cluster 3 is still above 0.5 (0.525). The 4 cluster group solution (Figure 5c) is also inadequate since the distances between clusters 3 and 4 and between clusters 1 and 3 are below 0.5 (0.417 and 0.391, respectively), giving the 3 cluster group the best and optimum solution. For SPEI12, the Euclidean distances in Table 2 (Figure 5d–f) provides quite a similar result which gives the optimal solution for the 3 groups clustering, with almost the same areas represented as for the SPEI6 (one grid center point of difference between the 3 cluster solution when comparing SPEI6 and SPEI12).

The spatial classification that resulted from the K-means cluster analysis is equivalent to the one obtained with the PCA, suggesting that both methods can be used to identify spatial areas with different temporal patterns among them. It is therefore reasonable to state that mainland Croatia could be divided by the PCA into three different regions (D1, D2 and D3 – Figure 4) with 34, 25 and 44 grid center point series for SPEI6, respectively, and with 36, 22 and 45 for SPEI12. For both SPEI time-scales, approximately the same areas were identified as D1 (Central North region), D2 (Eastern region), and D3 (Southern region), giving consistency to the drought regionalization in the country for the selected SPEI time-scales.

For some of the subsequent analysis, the SPEI6 and SPEI12 regional series, which are representative of the resulting homogeneous regions were considered, corresponding to the respective factor scores.

3.2. Temporal Evolution of Drought Areas

The areas affected by drought are an important indicator of the severity of drought events and can also indicate the frequency of various categories of drought. For that purpose, the percentages of area affected by moderate, severe, and extreme drought, according to the SPEI categories of Table 1, were computed for the SPEI6 and SPEI12 series in the three Croatian regions (Figure 6 and Figure 7). The area of each grid center point was accumulated in order to calculate the total area inside each region D1, D2 and D3 attributed to a specified drought category with values of SPEI6 and SPEI12 within those limits. To further understand the linear and nonlinear temporal evolution of the total areas under drought, linear trend lines along with a moving average considering a running length of 120 months, i.e., 10 years expressing the percent change of drought area were fitted.

In general terms, it can be seen that for both SPEI6 (Figure 6) and SPEI12 (Figure 7) fields in each of the regions (D1 to D3), the evolution of the percentage of drought affected areas shows similar trends and that from the year 2000 onwards, the frequency of occurrences of extreme droughts, with an evidently affected cumulative area, increased comparatively to the rest of the time period in both SPEI time-scales. Indeed, droughts are spread equally among each Croatian region and the analysed time window, so that the period with the smallest areas affected by drought being the years 1975 to 1983, with mainly moderate droughts. For SPEI6 and for the entire period - 1950 to 2022 - the number of months when drought affected more than 50% of the total area in each region was 180, 192, and 171, respectively, for the Central North region (D1), Eastern region (D2), and Southern region (D3), and 189, 176, and 183 for the same respective regions and period in SPEI12. The documented extraordinary drought in Croatia that took place in 2011/2012 [51] is clearly seen in each of the time-scales, especially in SPEI12, where regions D1 and D2 and almost all the area of D3 experienced an extreme drought situation. Another major drought event occurred in 1950/1951 and 2005/2006 for SPEI6 with 100% of the country experiencing an extreme drought. In years 1951/1952 and 2005/2006, the SPEI12 droughts affected almost all the regions D1, D2 and D3.

The change in the percentage of total drought affected area at both time-scales presents similar linear trends, showing an increasing trend with time in regions D1 and D2, while for region D3, it shows a flattening trend for SPEI6 and a slightly decreasing trend SPEI12 (Figure 6 and Figure 7, respectively).

3.3. Changes in Yearly Drought Occurrence Rate

In order to analyse the changes in the yearly drought occurrence rate, the KORE estimator was applied in each of the three regions, D1 to D3, to the regionalized SPEI series represented by the factor scores previously obtained by the PCA. The results, shown in Figure 8 and Figure 9 for SPEI6 and SPEI12, respectively, refer to the occurrences of moderate, severe and extreme droughts, which means SPEI values between the drought thresholds (DT) of -0.84 and -1.28, between -1.28 and -1.65, and lower than -1.65, respectively. The dates of the periods under drought conditions are identified by vertical black ticks on the x-axes. A 95% bootstrap confidence band was assigned to each λ(t) curve.

The bandwidth h that appears in equation (2) was obtained using Silverman’s ‘rule of thumb’ [52] which resulted in values ranging from ca. 30 to 110 months.

Comparison of the results of Figure 8 and Figure 9 for the same drought threshold (DT) shows that the general temporal behaviour of the moderate, severe and extreme drought occurrence rates is relatively independent of the SPEI time-scale (6 or 12 months). For SPEI12 in region D1 (Central North region) there is a strong general increase in the occurrence of moderate droughts starting in the 1990’s, followed by a decrease initiated around 2010 till the end of the analysed period, with similar behaviour for severe and extreme droughts. For SPEI6, the same fluctuation pattern in the rate of occurrences of drought is seen but much less pronounced, being that for severe drought category (DT=-1.28) after the highest rate of annual droughts occurs in 2010 is achieved, there is only a slightly decrease towards the end of the analysed period.

For region D2 (Eastern region) the same general increase pattern in the drought occurrence rate among all the drought thresholds and for both SPEI time-scales is seen in both Figure 7 and Figure 8. Starting around the beginning of the 1980’s for SPEI6 and the 1970’s for SPEI12, a persistent increase in the yearly drought occurrence rate is clear, especially for moderate droughts (DT=-0.84) where the increase rate is more pronounced. Region D3 (Southern region) is the one with a different pattern, with an identified general decrease in the drought frequency towards the present, for both SPEI6 and SPEI12, though much less evident, a general slightly decrease over time can be seen for extreme drought (DT=-1.65) whereas for moderate and severe droughts, an increase in the drought occurrences started around 1975 and reached a peak in 1990 to immediately start decreasing again till the end of the period analysed.

The obtained results show that the areas most prone to drought occurrence in Croatia, where an upward trend in the number of droughts per year can be seen since 1950, are the regions D1 and D2 in both SPEI time-scales. The Southern coastal region (D3) with decreasing rates of drought occurrences has a different pattern during the entire observed period, showing some progressively drought weakening (Figure 8 and Figure 9).

3.4. Temporal Trend Analysis

The temporal variability of droughts from 1950–2022 was analysed by exploring significant trends, either positive or negative, in the SPEI6 and SPEI12 time series fields with the modified Mann-Kendall test (MMK) [22]. Figure 10 shows the drought trends in the different homogenous regions using Sen’s slope and the results of the MMK test. The color scale represents the magnitude of the slope, and the black dots represents the statistical significant trend at the 0.05 significance level (p < 0.05).

From the results, it is clear the overall upward trend (blueish) in D3, suggesting progressively less drought susceptibility along the study period, strictly following the relief of the Southern region (D3) which is somehow consistent with the downward trend in the number of drought occurrences of Figure 8 and Figure 9 in both time-scales, especially on SPEI12. Only a few grid center points have statistical significant upward trends: 3 for SPEI6 and 9 for SPEI12, respectively in Figures 10a and b. The region of Istria and some coastal parts of Zadar and Sibenik-Knin are the exceptions to the upward trends of D3 with some reddish areas experiencing downward trends for both SPEI timescales, however only one of them denoting statistical significance (Zadar for SPEI12).

In both Central North region (D1) and Eastern region (D2) an opposite general downward trend (reddish) showing more susceptibility to droughts along the study period, was obtained, meaning an intensification of the drought phenomena. This is more pronounced for SPEI12 due to the higher trend slope (a more reddish area) and a slightly higher number of grid center points with statistical significant trends. In fact, almost all the trends identified in regions D1 and D2 have a statistical significance, with marked grid center points covering most of the regions (black dots). These results are somehow in agreement with those of Figure 8 and Figure 9 in terms of the frequency of different categories of drought occurrences that show an upward trend for both SPEI timescales, except for Figure 9 in region D1, which is not so clear.

4. Discussion

According to Seventh National Communication of the Republic of Croatia under the UNFCCC (2018) and Strategy on Adaptation to Climate Change in Croatia (2020) [54,55], several priorities have been highlighted. One of them is strengthening the professional, scientific and management capacities to deal with risks of climate change impacts on water resources. Within the proposed measures and activities are the development of monitoring systems, data acquisition and research as fundamental for estimation of vulnerabilities in agricultural and environmental sectors. Within this framework the obtained regionalization, by dividing mainland Croatia into three distinct homogenous regions, will certainly give new prospects and help future activities, especially relevant for the most susceptible areas of Croatia.

Recent investigations on extreme hydrological events were mostly focused on observed changes in the seasonal and annual air temperature and precipitation amounts in Croatia [26], or in studying the spatiotemporal characteristics of mean and maximum dry spells in Croatia [55], but they did not deal with drought as a complex hydro-climatological event inter-related with many morphoclimatic characteristics of a certain area.

Research involving drought phenomena in Croatia, its severity in terms of duration and intensity was mostly conducted on regional or local scales [20,21,22], or focused on specific drought events such as a meteorological analysis of an extreme drought in 2011/2012 which seriously affected the territory of Croatia [51]. In this research, the spatial and temporal characteristics of droughts obtained by using the SPEI6 and SPEI12, calculated upon the two main climate drivers (precipitation and temperature), seems to be aligned with the geomorphological characteristics of the regions identified. Region D3, which covers the mountain system of the Dinarides and the Adriatic basin, has the highest average elevation in the country and the highest average annual precipitation, while in regions D1 and D2, representing the Pannonian basin, the average annual precipitation progressively decreases towards the eastern areas. The Eastern region (D2) cover all the Lowland areas, with elevations up to 200 m and represent 50% of the total area of the country.

Although these results should be taken with caution due to the complexity of the climate system, they further highlight the need for the development of future drought research and in-depth studies on drought variability and the respective main drivers in the regions identified. Similar results were found in some literature examples around the world such as the work of [8,11,12,50] were the identification of homogeneous regions were consistent with the morphological and climate characteristics of those territories, which seems to play an important role in shaping drought areal patterns.

From the results obtained for SPEI6, time-scale typically associated with agricultural drought occurrence [6], the linear trend and the 120 month moving average in the drought areal evolution of regions D1 and D2 (Figure 6 and Figure 7) are illustrative of the upward trend on the areas affected by drought, which is particularly relevant since the Pannonian region (covered by D1 and D2) is the most important and the largest agricultural region of Croatia with highly developed intensive arable farming and high yields of most of the crops [56]. If the drought affected areas in these regions occurred within the months of May to October droughts tend to be more critical since they reflect the hottest period of the year in Croatia (Bjelovar and Gradište from Figure 2a) where the atmospheric evaporative demand is higher. In most of the region D3, the trend in areal evolution is almost negligible and even slightly negative for the SPEI12, effective time-scale at characterizing the country’s historical drought at an annual time-scale.

In general terms, for both regions D1 and D2 the number of drought occurrences through KORE is also increasing, showing some oscillatory periods along time, independently of the drought category and SPEI time-scale. For region D3 and both SPEI6 and SPEI12 some evident downward trend especially on moderate and severe droughts and a slightly decrease on extreme droughts is seen. The MMK trends applied to the SPEI grid center point series seems to have some relation with the drought areal evolution and occurrences rate analysis, since the trend signals are predominantly significant in regions D1 and D2, showing more locations with significant negative trends when comparing to just a few mainly non-significant positive trends in region D3. In region D3, for both SPEI time-scales the results seems to be in accordance with the results of [26] for the period 1960-2020, since for the mountainous system of the Dinarides, the authors identified a significant increase in precipitation due to very wet days and by using the daily intensity index, as well as for the Imotski station (Southern Adriatic coast) were an non-significant upward trend in the annual precipitation was detected [57]. However, the temperature effect on the aggravation of drought conditions seems to be less complex since it seems clear by several authors that statistically significant increases in air temperature especially in the mean annual temperatures, have been occurring consistently all over the country in the last 50 years [28,59].

5. Conclusions

To better deal with climate change and its impacts, it is crucial, as a first step, to understand climate variability at a sub-regional scale, so a generalized assessment of historical drought behaviour emerges as one of the main priorities to characterize the existence of water deficits at a subannual and annual scale throughout any country.

The present study accomplishes part of this objective by giving a quite complete insight into the drought phenomena in Croatia with such spatial resolution. Based on the E-OBS dataset, the SPEI was calculated on 6 and 12 month time-scales and for the period from 1950 to 2022. The PCA and K-means methods were used for drought regionalization and to explore the spatial patterns of drought distribution, describe and plot drought characteristics such as drought occurrence frequency and drought severity by the calculation of the percentage of drought affected areas. The drought trends were also analysed using the MMK method coupled with the Sen’s slope. The main findings of the study are as follows:

(1)

Based on PCA and K-means validation Croatia was divided into three homogeneous regions: Central North region (D1), Eastern region (D2) and Southern region (D3).

(2)

Central North region (D1) and Eastern region (D2) showed an upward trend in the percentage of areas affected by drought in the whole study period for both SPEI6 and SPEI12, but in Southern region (D3) a negligible trend was obtained for SPEI6 and a downward trend, meaning fewer areas progressively affected by drought were obtained for SPEI12. Both D1 and D2 areas have large non-irrigated agricultural land and grassland, resulting in high ecological vulnerability.

(3)

Region D1 (Central North region) experiences an increase in the drought occurrence rate from 1950 until around 2010 and some decrease in the last 10 years period, especially pronounced in SPE12. The Eastern region (D2) experienced a generalized continuous increase in the drought occurrence rate from 1950 to 2022 in all drought categories and SPEI time-scales. In the Southern region (D3), a decrease in the drought occurrence rate was obtained with one interruption peak. According to the nature of the SPEI calculation procedure, an increase in the number of drought occurrences over the years means progressively more periods of time with negative water balances which induces necessarily to bigger challenges on water management practices.

(4)

A generalized change towards more susceptibility to drought conditions in most of the areas of D1 and D2 was obtained with the MMK test with strong statistical significance in both SPEI6 and SPEI12. Given the Sen´s slope values obtained from the trend analysis applied to the SPEI series, more intense drought events are expected in those areas. In Southern region (D3), the trend into less susceptibility to drought conditions spatially follow the mountainous areas of the Dinarides but with less statistical significance. The region of Istria and some coastal parts of Zadar and Sibenik-Knin are the exceptions to the general pattern found in D3 since some localized areas will become a bit more susceptible to drought, which is seen in both SPEI time-scales.

(5)

In general terms, the west (Mediterranean climate) is becoming less susceptible to drought while the east (continental climate) is tending to become more prone to have an intensification of drought events, showing a greater increase in the areas affected by drought over the years and an increasing rate of occurrence of the number of annual droughts. Although the Mediterranean region is usually at the center of drought research, it is in the mainly agricultural mainland of Croatia that drought conditions seems to have worsened.

Due to Croatia's complex morphoclimatic characteristics, there are still some limitations when it comes to using the results obtained, such as delving into the underlying causes and implications of the regional spatio-temporal drought patterns found in the homogeneous regions, or how these patterns relate to climate change and other environmental factors, such as soil moisture or vegetation, and what their roles are in shaping drought susceptibility areas. Clearer evidence from this study suggests that, in addition to climate drivers, topographical features greatly influence the spatial patterns here identified. However, the spatio-temporal variability of droughts identified, given the limitations of using only SPEI6 and SPEI12, suggests that other time-scales should be considered in future studies in order to obtain more information on the drought behaviour in each homogenous region.

Following the results here presented, future studies should also investigate the main teleconnection patterns that could explain drought variability in the country, while, since temperature plays an important role on drought conditions, the effect of considering more physically realistic PET methods in drought characterization, such as the FAO-56 Penman-Monteith [59], should also be considered.

To reduce and mitigate drought impacts, the development and implementation of risk-based drought management plans is encouraged in order to support policy makers and water resources managers in developing the best strategies. The authors hope that the results thus achieved could provide valuable support information by highlighting more susceptible regions to drought than others, were the severity, the yearly number of occurrences and the affected areas has been increasing over the years.