Hydrological Drought Assessment in a Small Lowland Catchment in Croatia

: Hydrological drought is critical from both water management and ecological perspectives. Depending on its hydrological and physical features, the resilience level of a catchment to groundwater drought can differ from that of meteorological drought. This study presents a comparison of hydrological and meteorological drought indices based on groundwater levels from 1987 to 2018. A small catchment area in Croatia, consisting of two sub-catchments with a continental climate and minimum land-use changes during the observed period, was studied. The ﬁrst analysis was made on a comparison of standardized precipitation index (SPI) and standardized precipitation evapotranspiration index (SPEI). The results showed their very high correlation. The correlation between the standardized precipitation index (SPI) and standardized groundwater index (SGI) of different time scales (1, 3, 6, 12, 24 and 48 months) showed different values, but had the highest value in the longest time scale, 48 months, for all observation wells. Nevertheless, the behavior of the SPI and groundwater levels (GW) correlation showed results more related to physical catchment characteristics. The results showed that groundwater drought indices, such as SGI, should be applied judiciously because of their sensitivity to geographical, geomorphological, and topographical catchment characteristics, even in small catchment areas.


Introduction
Droughts can be studied from different perspectives. This paper discusses some aspects of hydrological drought. Meteorological drought is the first indicator of drought occurrence, duration, and severity, based only on precipitation and air temperature. Its prolonged duration can cause hydrological drought, which affects water resources, the environment, and the economy [1,2].
Meteorological droughts occur more frequently than hydrological droughts because of the physical and hydrogeological characteristics of catchments [3]. A hydrological drought is broadly defined as a negative anomaly in river discharge and groundwater levels. However, hydrological drought is more complex and has more aspects of drought impact than meteorological drought [1].
Many analyses of hydrological droughts in Europe and other continents have been published since the beginning of the 21st century.
One of the first comprehensive anlyses was the pan-European research based on a dataset of more than 600 daily streamflow records analyzed over four time periods between 1911 and 1995 (1962-1990, 1962-1995, 1930-1995, and 1911-1995) [4]. This was followed by many other studies on different aspects of drought in Europe [4][5][6][7][8][9][10][11][12] and particular European regions such as the Mediterranean [5,13,14] or Central Europe [15][16][17]. The scientists involved in these and many other studies shared conclusions regarding drought complexity, temporal and spatial diversity, and their effects on all components of the hydrological cycle. ment is predominantly agricultural, situated in the Pannonian Valley, in the Danube River basin. This area shares common features with the lowland areas of Hungary, Slovakia, and Germany, among others. Furthermore, regarding continental climate, dominance of agricultural area and relatively large available water resources, we explored the occurrence of hydrological drought over different time scales (1, 3, 6, 12, 24 and 48 months). The results obtained and the methodology discussed can be valuable for drought analysis of European agricultural regions.

Study Area
Drought analyses should be conducted in catchments with minimum anthropogenic impacts over a longer period to obtain reliable results [35]. In this study, a small catchment of the Karašica and Vučica rivers in northern Croatia was selected. The catchment is located in the continental part of Croatia, which belongs to the Danube River basin and its tributary, the Drava River sub-basin. It consists of hilly and lowland areas, mainly under agricultural usage. Based on hydrological and geological characteristics, the downstream part of the catchment is considered a perennial stream with a porous aquifer. It is divided into two sub-catchments of almost equal size ( Figure 1). The first small sub-catchment, the Karašica River catchment, has dominant lowland characteristics with a porous aquifer and a small average longitudinal slope (0.011%) (Figure 2a, Table 1). It mostly consists of agricultural land (over 50%), pastures, and forests ( Figure 2b) [35].
The data showed on Figure 2 were interpolated by IDW method (inverse distance weighting) by QGIS computer program. IDW interpolation gives weights to sample points, such that the influence of one point on another declines with distance from the new point being estimated. Interpolation results are shown as a two-dimensional raster layer (QGIS Documentation v:3.22. https://docs.qgis.org/3.22/en/docs/gentle_gis_introduction/spatial_analysis_interpolation.html, accessed on 1 February 2022)  Both sub-catchments were recognized as almost natural, with minimum human impact on the hydrological regime (small portion of irrigation land, no new constructed reservoirs or water abstractions, etc.) and steady land use and vegetation cover [35].
In this part of Croatia, more frequent drought episodes and warming started in 1988 [36]. Severe droughts occurred in 2000 and 2003, but one of the most extreme droughts was in 2011/2012. In the continental region, its duration was extremely long, with the highest magnitudes since the beginning of the twentieth century [37]. Particularly, this catchment area was characterized by an increasing trend in air temperatures and, at the same time, variability of the precipitation regime [38].
Climate change models predict future warming, especially for summer, in all Croatian regions [39].
The first small sub-catchment, the Karašica River catchment, has dominant lowland characteristics with a porous aquifer and a small average longitudinal slope (0.011%) ( Figure 2a, Table 1). It mostly consists of agricultural land (over 50%), pastures, and forests ( Figure 2b) [35]. The rest of the catchment is mostly hilly and dominated by forests. Geological characteristics of the study area are exceptionally important for the analysis of hydrological (groundwater) drought. Figure 3 presents the geological catchment map (http://webgis.hgi-cgs.hr/gk300, accessed on 18 March 2022). Most of the lowland area consists of alluvial deposits of different permeability.  The data showed on Figure 2 were interpolated by IDW method (inverse distance weighting) by QGIS computer program. IDW interpolation gives weights to sample points, such that the influence of one point on another declines with distance from the new point being estimated. Interpolation results are shown as a two-dimensional raster layer (QGIS Documentation v:3.22. https://docs.qgis.org/3.22/en/docs/gentle_gis_introduction/spatial_ analysis_interpolation.html, accessed on 1 February 2022) Figure 2c depicts the average annual precipitation from 1987 to 2018. The annual precipitation varied between 710 mm at 97.0 m above sea level, and 922 mm at 180.0 m above sea level, at the five meteorological sites. From 1987 to 2018, groundwater levels were recorded at six observation wells in the lowest part of the basin. Figure 2d depicts the average annual groundwater level during the observation period.
The rest of the catchment is mostly hilly and dominated by forests. Geological characteristics of the study area are exceptionally important for the analysis of hydrological (groundwater) drought. Figure 3 presents the geological catchment map (http://webgis.hgi-cgs.hr/gk300,

Input Data
Meteorological and hydrological data used in further analyses of drought indices were provided by the Croatian Meteorological and Hydrological Service, https://meteo.hr/index_en.php (accessed on 15 January 2022). Monthly precipitation data for the period 1987-2018 from five meteorological stations were used for the meteorological drought index. Hydrological drought indices were calculated based on the daily groundwater levels (GW) of six observation wells. These are shown in Figure 1. The positions of the observation wells and average annual groundwater levels for the period 1987-2018 are given in Table 2. Altitude "0" refers to the bottom of the observation well.

Input Data
Meteorological and hydrological data used in further analyses of drought indices were provided by the Croatian Meteorological and Hydrological Service, https://meteo. hr/index_en.php (accessed on 15 January 2022). Monthly precipitation data for the period 1987-2018 from five meteorological stations were used for the meteorological drought index. Hydrological drought indices were calculated based on the daily groundwater levels (GW) of six observation wells. These are shown in Figure 1. The positions of the observation wells and average annual groundwater levels for the period 1987-2018 are given in Table 2. Altitude "0" refers to the bottom of the observation well.

Standardized Precipitation Index (SPI)
The SPI, which is the approach most widely used in all parts of the world, regardless of meteorological or topographical conditions, was used to analyze meteorological drought. Each dataset was fitted to the gamma function to define the relationship between probability and precipitation. Once the relationship of probability to precipitation was established from the historical records, the probability of any observed precipitation data point was calculated and used, along with an estimate of the inverse normal, to calculate the precipitation deviation for a normally distributed probability density with a mean of zero and a standard deviation of unity [23]. Negative SPI and SPEI values indicated mild to extremely dry periods ( Table 3). The same index can be used to evaluate drought severity in various water resources (groundwater, open watercourses, and soil moisture) depending on the purpose of the drought analysis. Time scales of 1-6 months are appropriate for the analysis of drought in agriculture, 1-2 months are needed for meteorological drought, and 6-24 months are needed for hydrological drought. Table 3. Limit values for standardized precipitation index (SPI) and standardized precipitation evapotranspiration index (SPEI). Data from [23,40].

SPI, SPEI Values
Drought Category The results of the 1, 3, 6, 12, 24, and 48 month SPI calculations for the period 1981-2018 are given in the Supplementary Materials. Extreme droughts of short periods (1, 3, and 6 months) with SPI values < −2.0 occasionally occurred during the complete observation period. The 24 and 48-month SPIs occurred at the end of the 1990s over a longer period.

Standardized Precipitation Evapotranspiration Index (SPEI)
The SPEI was proposed by Vicente-Serrano [24]. It is based on the same calculation as the SPI, which uses the normalized gamma distribution of precipitation and provides a number of standard deviations with regard to an average value. The only difference is that precipitation is replaced by the difference between precipitation and potential evapotranspiration. The classification of SPI and SPEI values are the same and are shown in Table 3. Potential evapotranspiration was calculated based on the Thornthwaite equation. The SPEI results are provided in the Supplementary Materials.

Standardized Groundwater Level Index (SGI)
Groundwater flow is complex process, and it is difficult to obtain direct observational data related to groundwater resources [32]. Water accumulation in the soil, nonetheless, is an important water balance component that has a significant impact on agriculture and water use. SPI is a good technique for groundwater drought evaluation when applied to groundwater level data from a specific area under standard set conditions with variable accumulation times. The time periods required to achieve a maximum correlation exhibit high spatial variability (ranging from 3 to 36 months) [41]. The interpretation of the resulting SGI must reflect an appreciation of the hydrogeological environment of the observation borehole, because SGI can be strongly influenced by site-specific recharging [29]. This approach is similar to the standardized water level index (SWI), with certain adjustments [42]. Since ground water level is measured down from the surface, positive anomalies correspond to Hydrology 2022, 9, 79 7 of 15 drought and negative anomalies correspond to wet or normal conditions. To avoid this approach, we used water depth data measured from the bottom of each observation well: where G ij represents monthly water level of i-th well and j-th observation, G im is the mean value of the data series, and σ is its standard deviation.

Correlation Function
Correlations of SPI and SPEI, SPI and SGI, and SPI and GW, between 1 and 48 months were calculated using the Pearson correlation method [43]: where N represents the number of data points in each data series, and x i and y i are data points in the first and second data series, respectively. The Pearson correlation method has been recognized as the most appropriate, compared with the Kendall and Spearman correlation functions [44]. According to the physical processes of surface runoff and percolation into the soil, a temporal delay between the precipitation and water balancing components is necessary. After a certain period, rainfall deficiency and drought consequences are recognized as below-average storage conditions in surface water, reservoirs, and groundwater [41]. As a result, hydrological drought forecasting has become more complicated and unreliable.
Another study concluded that SGI values correlated well with SPI values that were observed 5-8 months ago [45].

Results and Discussion
Meteorological drought indices, standardized precipitation index and standardized precipitation evapotranspiration index were compared, in order to define their reliability for further analysis. We calculated their correlations over different time scales (1 to 48 months). The results are presented in Figure 4, showing a very high level of correlation.
All SPI/SPEI correlation values were very high, between 0.963 and 0.983. The highest values of Pearson correlation coefficients were obtained for 3 and 6 month timescales. The other basic statistic parameters, such as index of agreement (d), mean absolute error (MAE), root mean square error (RMSE), mean bias error (MBE), trendline equation (slope and intercept) and determination coefficient (R 2 ) also confirmed high correspondence between all SPI and SPEI values (Table 4). This result led us to the conclusion that application of SPI in further analyses gives acceptable accuracy.   The SPI calculation procedure only requires precipitation data, which are almost always available. Hydrological drought deficits are influenced by average catchment wetness (mean annual precipitation) and elevation (reflecting seasonal storage in the snowpack and glaciers), both of which are governed by a complex combination of climate and catchment characteristics, including catchment storage capacity [29].
The correlation between the SPI and standardized groundwater index (SGI) was highest at the longest time scale (48 months) for all observation wells. Correlations between SPI and SGI ranged between 0.52 and 0.82, as shown in Figure 6.
Calculated SGI was correlated to the SPI of the nearest meteorological station. The SPI/SGI correlation coefficients for most observation wells were continuously increasing over different time scales. Only P-23 exhibited a typical behavior. According to Table 2, it is a shallow piezometer located close to the fishpond (Figure 1), and its groundwater table was strongly impacted by surface water. It exceeded the maximum SPI/SGI correlation at the 48 month time scale, similar to the others, but the correlation coefficient at the 1 month scale was higher than on the 6 and 12 month scales. Very high correlation coefficients for P-18 were only at the 24 and 48 month scales, and for P-17 at the 48 month scale. For the other wells, except for P-23, the correlations were moderate at 24 and 48 months. The spatial distributions of the calculated SPI/SGI correlations over different time scales are presented in Figure 7. The highest values of the correlation function for all observation wells were at the longest time scale (48 months), with very high correlation coefficients.  The SPI calculation procedure only requires precipitation data, which are almost always available. Hydrological drought deficits are influenced by average catchment wetness (mean annual precipitation) and elevation (reflecting seasonal storage in the snowpack and glaciers), both of which are governed by a complex combination of climate and catchment characteristics, including catchment storage capacity [29].
The correlation between the SPI and standardized groundwater index (SGI) was highest at the longest time scale (48 months) for all observation wells. Correlations between SPI and SGI ranged between 0.52 and 0.82, as shown in Figure 6. Calculated SGI was correlated to the SPI of the nearest meteorological station. The SPI/SGI correlation coefficients for most observation wells were continuously increasing over different time scales. Only P-23 exhibited a typical behavior. According to Table 2, it As stated earlier, the catchments are situated in lowland Holocene river terraces, with the lowest parts at 80-120 m above sea level. They are built up from fluvial sediments of variable thickness, permeability, and composition (gravel, sand, silt, or clay) belonging to Gleysols, Fluvisols, and sporadically, other hydromorphic soils characterized by large clay content [46].
Input data for mapping SPI/SGI correlations are in given in the Supplementary Materials.
Further application of correlation analysis on SPI and groundwater levels (GW) showed a different relationship (Figure 8). The results showed that the curves of most correlation coefficients followed the same pattern (P-17, P-6 and P-18). Each curve in this group increased as it reached the correlation peak (0.6 to 0.66) in 12 months, and then slowly decreased to values between 0.25 and 0.48. Observation well P-58 had the highest correlation peak (0.72), which appeared much later, after 24 months. Observation well P-9 had the lowest SPI/GW correlation coefficient; the maximum value was 0.22 and it appeared in a short time period, after 6 months. After that, its correlation curve dropped to As stated earlier, the catchments are situated in lowland Holocene river terraces, with the lowest parts at 80-120 m above sea level. They are built up from fluvial sediments of variable thickness, permeability, and composition (gravel, sand, silt, or clay) belonging to Gleysols, Fluvisols, and sporadically, other hydromorphic soils characterized by large clay content [46].
Input data for mapping SPI/SGI correlations are in given in the Supplementary Materials. Further application of correlation analysis on SPI and groundwater levels (GW) showed a different relationship (Figure 8). The results showed that the curves of most correlation coefficients followed the same pattern (P-17, P-6 and P-18). Each curve in this group increased as it reached the correlation peak (0.6 to 0.66) in 12 months, and then slowly decreased to values between 0.25 and 0.48. Observation well P-58 had the highest correlation peak (0.72), which appeared much later, after 24 months. Observation well P-9 had the lowest SPI/GW correlation coefficient; the maximum value was 0.22 and it appeared in a short time period, after 6 months. After that, its correlation curve dropped to negative values. According to Table 2, it is a deep piezometer located in the urban area, close to the river (Figure 1), and its groundwater table was probably impacted by soil structure and surface water.
close to the river (Figure 1), and its groundwater table was probably impacted by soil structure and surface water.
Observation well P-23 exhibited similar behavior as in SPI/SGI correlation analysis. According to Table 2, it is a shallow piezometer located close to the fishpond (Figure 1), and its groundwater table was strongly impacted by surface water. It exceeded the maximum SPI/GW correlation at the 24 month time scale (0.33) but its correlation coefficient curve was completely different. In order to check their correlation uniformity, we introduced SPI/GW (groundwater level) correlation analysis [44]. SPI/GW correlation analysis was performed in order to examine the relationship between meteorological drought and groundwater level. The spatial distributions of the calculated SPI/GW correlations over different time scales are presented in Figure 9. Compared with Figure 7, the maximum correlation coefficients for most observation wells appeared over a shorter time period, 12 and 24 months. Input data for mapping SPI/GW correlations are in given in the Supplementary Materials. Observation well P-23 exhibited similar behavior as in SPI/SGI correlation analysis. According to Table 2, it is a shallow piezometer located close to the fishpond (Figure 1), and its groundwater table was strongly impacted by surface water. It exceeded the maximum SPI/GW correlation at the 24 month time scale (0.33) but its correlation coefficient curve was completely different.
In order to check their correlation uniformity, we introduced SPI/GW (groundwater level) correlation analysis [44]. SPI/GW correlation analysis was performed in order to examine the relationship between meteorological drought and groundwater level. The spatial distributions of the calculated SPI/GW correlations over different time scales are presented in Figure 9. Compared with Figure 7, the maximum correlation coefficients for most observation wells appeared over a shorter time period, 12 and 24 months. Input data for mapping SPI/GW correlations are in given in the Supplementary Materials.
Observed results showed lower and more variable correlation values. The correlation between the SPI and groundwater level (GW) was the highest at different time scales (6 to 24 months) for all observation wells. Maximum correlations between SPI and GW ranged between 0.22 and 0.72, as shown in Figure 8. The correlation coefficients showed much higher variability related to the position and depth of observation wells. Again, only shallow observation well P-23 had specific behavior because of its dependence on surface water (fishpond). Its SPI/GW correlation coefficients were, compared with others, very low. Additionally, P-9 had a relatively low correlation (even negative) with meteorological drought, due to its position in the urban area, very close to the Karašica River. It is a deep observation well and it was probably strongly impacted by deeper aquifers (Figure 1). Observed results showed lower and more variable correlation values. The correlation between the SPI and groundwater level (GW) was the highest at different time scales (6 to 24 months) for all observation wells. Maximum correlations between SPI and GW ranged between 0.22 and 0.72, as shown in Figure 8. The correlation coefficients showed much higher variability related to the position and depth of observation wells. Again, only shallow observation well P-23 had specific behavior because of its dependence on surface water (fishpond). Its SPI/GW correlation coefficients were, compared with others, very low. Additionally, P-9 had a relatively low correlation (even negative) with meteorological drought, due to its position in the urban area, very close to the Karašica River. It is a deep observation well and it was probably strongly impacted by deeper aquifers (Figure 1).

Conclusions
According to the analysis of meteorological and groundwater drought, it can be concluded that the SPI and SGI correlations are significant over a longer time scale. The study area is located in Danubian valley and is characterized by multi-layered alluvial lowland.
The time taken for rainwater to accumulate in the subsurface and percolate into deeper layers is important for detecting groundwater shortages. Hydrological drought

Conclusions
According to the analysis of meteorological and groundwater drought, it can be concluded that the SPI and SGI correlations are significant over a longer time scale. The study area is located in Danubian valley and is characterized by multi-layered alluvial lowland.
The time taken for rainwater to accumulate in the subsurface and percolate into deeper layers is important for detecting groundwater shortages. Hydrological drought can take up to 36 months to manifest, compared with meteorological conditions, according to published research [17].
In addition, in regions with a carbonate basement, the delay is much less, only 18 months [44,47]. However, the behavior of groundwater is complex, which confirms this finding. In this research, groundwater drought, expressed by SGI, is about 48 months' delayed. The lag impact on the SGI revealed variations, since groundwater level is influenced by a variety of hydraulic and hydrological phenomena. The correlation between the SPI and standardized groundwater index (SGI) was the highest at the longest time scale (48 months) for all observation wells. Correlations between SPI and SGI ranged between 0.52 and 0.82, as shown in Figure 6. The observed area was only 2000 km 2 , approximately, but SPI/SGI curves had almost the same shape, with the exception of P-23 ( Figure 6).
Correlation coefficients showed similar behavior related to position and depth of observation wells. Only shallow observation well P-23 showed specific behavior because of its dependence on surface water (fishpond). Its SPI/SGI correlation coefficients were the lowest.
The results confirmed that groundwater response to drought and its changes over time are very site specific, and observations are, therefore, hardly scalable from point to catchment scale [48].
Drought, particularly hydrological drought, is a complicated phenomenon. The size of the catchment area greatly affects the accuracy of drought classification. This is especially evident in groundwater drought analysis, where the location of observation wells, their height/depth, proximity to water bodies, geological properties, and other factors significantly impact drought classification. In addition, meteorological drought severity, compared with hydrological droughts, is much more complex. Groundwater drought has a long delay period and lasts for decades. Therefore, predicting hydrological drought based on meteorological data is highly risky and unreliable.  Table S1: Data base of spatial distribution of SPI/SGI correlation over different time scales presented in Figure 7; Table S2: Data base of spatial distribution of SPI/SGI correlation over different time scales presented in Figure 9.
Author Contributions: Conceptualization, T.B.; methodology, T.B. and L.T.; data analysis, L.T. and T.B.; writing-original draft preparation, T.B. and L.T.; writing-review and editing, L.T. and T.B. All authors have read and agreed to the published version of the manuscript.