ABSTRACT
The escalating impacts of global climate change significantly affect regional hydrological systems, particularly in northern areas such as Estonia. This study investigates the hydro l ogical sensitivity of Estonian catchments to climatic variability, focusing on the interplay between surface water and groundwater. Using data from 42 river catchments, it employs various statistical methods in hydrology, emphasizing the autocorrelation function, cross-cor relation function, baseflow index, and flow duration curve. The analysis spans the years 2012-2022, integrating hydrological, spatial, and water quality parameters. The research identifies four distinct hydrological behavior clusters: plateau, sandstone upland, carbonate upland, and lowland. Key findings include diverse catchment sensitivities to groundwater recharge, the role of baseflow in streamflow stabilization, the memory effect in catchment responses, and insights from the flow duration curve on flow variability and extremes. The LightGBM model, predicting focus parameters, highlights the critical influence of air tempera ture and snowpack on streamflow characteristics. This study underscores the diverse hydro logical sensitivities of Estonian catchments to hydroclimatic changes, emphasizing the impor-tance of considering catchment-specific characteristics in water resource management and policy-making. Con -tributing to the broader understanding of hydrological processes, it provides valuable insights for future research and environmental planning in the face of climate variability and change.
Keywords:
hydrological sensitivity, surface water/groundwater interaction, autocorrelation function (ACF), cross-correlation function (CCF), baseflow index (BFI), flow duration curve (FDC), climate change
Globaalne kliimamuutus mõjutab märkimisväärselt piirkondlikke hüdroloogilisi süsteeme, eriti põhjapoolsetes parasvöötme riikides nagu Eesti. Uurimus käsitleb Eesti valgalade hüdroloogilist tundlikkust kliimamuutuste suhtes, keskendudes pinna- ja põhjavee koosmõju ilmingutele. Töös kasutatakse 42 jõe hüdromeetriajaama äravoolu ja valgala andmeid ning rakendatakse erinevaid hüdroloogilis-statistilisi meetodeid. Analüüs hõlmab 2012.-2022. veeaasta hüdroloogilisi, ruumilisi ja veekvaliteedi andmeid. Suuremat tähelepanu pööratakse ära-voolu autokorrelatsioonifunktsiooni (ACF), ristkorrelatsioonifunktsiooni (CCF), baasvooluindeksi (BFI) ja voolu kestuse kõvera (FDC) tulemustele ning nende seotusele teiste valgala parameetritega. Valgalade hüdroloogilise iseloomu järgi saab need jagada nelja klastrisse: platoo, liivakivi kõrgustik, karbonaatkivimi kõrgustik ja madalik. Töö tulemustest selgub, kuidas erinevad valgalade klastrites äravoolu tundlikkus toite- sündmuste suhtes, baasvoolu dünaamika, autokorrelatsioonifunktsiooni kirjeldatav "mäluefekt" ja näiteks ka põuapäevade arv. LightGBM gradienti suurendava masinõppe mudeli abil tuuakse esile õhutemperatuuri ja lumi katte suur mõju hüdroloogilisele tundlikkusele. Uurimus näitab Eesti valgalade erinevat hüdroloogilist tundlikkust hüdrokliimaa-tiliste muutuste suhtes, rõhutades valgala spetsiifiliste omaduste arvestamise tähtsust veevarude majandamisel ja poliitika kujundamisel. Töö panustab laiemasse arusaama hüdroloogilistest prot sessidest ning annab väärtus - likku teavet tulevasteks teadusuuringuteks ja keskkonnaplaneerimiseks kliima muutuse tingimustes.
1. Introduction
Global climate change is anticipated to accelerate, with significant changes in pre -cipitation patterns and large spatial heterogeneity in temperature increases projected across different regions (IPCC 2022). Over the last century, changes in temperature and precipitation patterns have disrupted regional hydrological systems, leading to alterations in the intensity, frequency, and duration of peak and low flows, with increasing occurrences of severe streamflow droughts (Teutschbein et al. 2022, 2015). Climate change significantly affects northern hydrology, particularly the accumulation and melt of seasonal snow cover, which is crucial for catchment runoff and ground -water recharge (Rodhe 1998; Earman and Dettinger 2011; Jenicek et al. 2016). Observed and projected changes include reduced snow depth and snow water equivalent, altered snowmelt timing, and a shift from snowfall to rain in winter (Berghuijs et al. 2014; Meriö et al. 2019; Grogan et al. 2020; Ranasinghe et al. 2021). These changes have already decreased peak stream flows and shifted them earlier in the year in northern Europe, resulting in lower summer-autumn flows and increased spring soil moisture deficits (Blöschl et al. 2019; Douville et al. 2021; Jaagus et al. 2017; Ruosteenoja et al. 2018; Viru and Jaagus 2020). Winter precipitation is projected to increase by up to 10% by 2100, although summer changes remain uncertain (increase <5%) (Coppola et al. 2021; Jaagus and Mändla 2014; Ruosteenoja and Jylhä 2021). Rising temperatures (2.5-3.8 °C in winter and 2.0-3.4 °C in summer) are expected to increase evapotranspiration, reduce summer runoff and recharge, and increase water surplus in colder seasons (Barnett et al. 2005; Okkonen et al. 2011). Con sequently, both short- and long-term changes in stream discharge and ground water storage are anticipated (Okkonen and Kløve 2011; Donnelly et al. 2017; Costantini et al. 2023).
Groundwater processes, which tend to be slower and less uniform than surface water events, lead to prolonged summer recessions and decreased late-summer and autumn discharges (Vorobevskii et al. 2022). Climate change is expected to sig -nificantly alter groundwater levels, aquifer storage, and base -flow, mainly due to changes in groundwater recharge timing and intensity (Taylor et al. 2013; Smerdon 2017; Wu et al. 2020; Nygren et al. 2020). In northern Europe and North America, studies have demonstrated groundwater systems to be particularly sensitive to climate change. For instance, groundwater streams in northwestern USA are more sensitive to climate change than runoff-dominated streams (Jefferson et al. 2008). In northern Finland, warmer climates are pre -dicted to increase groundwater recharge and shift peak water levels earlier, affecting groundwater-surface water inter-actions (Okkonen and Kløve 2011). Research on the Baltic States' groundwater anticipates a general rise in levels with min imal seasonal variation (Babre et al. 2023), while in Estonia, a 20-40% increase in groundwater recharge is pro -jected under warmer climates (Vallner 1998).
There was a necessity to offer a general evaluation of the hydrological sensitivity of Estonian surface water and as -sociated groundwater bodies, including their interactions, in the context of patchy and scarce dedicated monitoring data. This challenge led to the development of the approach used in this study, which draws partly on the findings of Koit (2022) and Koit et al. (2022) regarding the auto- and cross-correlation functions of Estonian lowland rivers. The study's approach is based on the assumption that the characteristics of the interplay between surface water and groundwater can be identified by analyzing the catchment's streamflow hydro -graphs (Killian et al. 2019; Teutschbein et al. 2015). We set out to examine how different gauged catchments in Estonia responded to hydroclimatic variability on a water-yearly scale during the period of 2012-2022 by quantifying various aspects of streamflow, such as its memory (Schuler et al. 2022; Sutanto and van Lanen 2022), the proportion of ground -water runoff (baseflow), and the sensitivity of streamflow to groundwater recharge events.
The autocorrelation function (ACF) characterizes the degree of linearity in a time series, indicating the extent to which a data point is affected by previous data points, often referred to as the memory effect (Bailly-Comte et al. 2008; Mangin 1984; Schuler et al. 2022). The concept of memory effect, as reflected through the ACF, is influenced by the effective storage of particular surface water and the associated groundwater catchment, with higher storage leading to greater memory or inertia in discharge. Catchments featuring greater memory can be considered less susceptible to meteorological drought compared to those with lower memory. As ground -water storage plays a crucial role in streamflow memory (Cochand et al. 2019), the contribution of groundwater to total streamflow should be estimated. The ratio of baseflow to total streamflow, defined as the baseflow index (BFI), represents the slow or delayed contribution to streamflow that is influenced by catchment hydrogeology. Complementary to the ACF, the BFI helps validate the presence of periodicity or persistence observed in the autocorrelation function, providing additional evidence of groundwater dynamics influencing streamflow patterns.
Furthermore, the reaction and sensitivity of streamflow dynamics to recharge events could be quantified to assess the interconnection of surface water and groundwater in a catch -ment. We applied the cross-correlation function (CCF) to quantify the linear dependency of discharge on the recharge signal. To overcome the bottlenecks associated with using precipitation as an input to the CCF in a climate where some precipitation falls as snow and is therefore tied up in the snow -pack for some time, as pointed out by Koit et al. (2022) and Pärn et al. (2024), we experimentally used a modeled ground -water recharge signal as the input in this study. To support the findings from the ACF, CCF, and BFI, we also used the flow duration curve (FDC), another robust hydrological method, for reflecting the variability within the stream flow regime. FDCs enable the delineation of flow permanence and the identification of hydrological extremes. In the context of cli -mate change, the FDC is instrumental in revealing the alter -ations in flow regimes and aiding in the understanding of how catchment characteristics modulate the response to cli -matic variability (Vogel and Fennessey 1994) . To evaluate the performance of the focus parameters - ACF, CCF, BFI, and FDC -, we applied a variety of other hydro logical and spatial analysis-derived statistics.
2. Materials and methods
2.1. Studied catchments
In this study, the hydrological sensitivity of 42 gauged river catchments in Estonia, ranging from 52 to 1813 km2 in size, was assessed (Fig. 1; Table S1). Estonia is located in the bor -eal biogeographical region between high latitudes 57°30' N and 59°49' N, and lies in the transition zone of the Baltic Sea maritime and continental climate regions (i.e., Cfb and Dfb Köppen-Geiger climate regions; Kottek et al. 2006; Beck et al. 2018). Based on the ERA5-Land dataset (Muñoz Sabater 2019), the mean annual precipitation during 2012-2022 was 747 mm, with about 18% falling as snow. On average, around 537 mm of precipitation was removed by total evaporation, leaving 210 mm for the generation of surface water and groundwater runoff. The most important hydrological events of the water year are the spring snowmelt floods (Järvet 1998). Furthermore, seasonal hydrological sensitivity is in -fluenced by snow deposition in the winter and evapotran -spiration during the summer season (Koit et al. 2022).
Estonia's geological and hydrogeological context is de -fined by its location in the NW part of the East European Platform. In Estonia, the low-lying monoclinal platform (mean elevation 50 m asl; max 317 m asl in SE uplands) consists of Neoproterozoic and Paleozoic sedimentary rocks, which host multiple aquifer systems belonging to the Baltic Artesian Basin. Silurian and Ordovician carbonate rocks predom i -nantly outcrop on the N-NW half and Devonian sandstones in the S-SE half of the territory. The sedimentary rocks are overlain by Quaternary deposits, formed mainly during the Late Weichselian glaciation. The territory of Estonia was degla -ciated between 15 to 13 ka BP (Kalm 2006). The Quaternary cover, 5-10 m thick on average (up to 100-200 m thick in the SE uplands), consists of different types of sediments by gen -esis, among which the most widespread are glacial tills, gla -ciolacustrine sands, silts and clays, and glaciofluvial sands and gravels (Raukas and Kajak 1997). Throughout the Holocene, cool and humid climate has been favorable to the devel -opment of vast peatlands, which cover up to 20% of the ter -ritory, usually overlying Quaternary sediments with poor per -meability.
2.2. Data
Data on streamflow, land cover, geology, hydrogeology, water quality, and hydroclimate for 42 Estonian catchments (Table S1) were collected from various public databases. The full data -set contained 97 parameters (Tables S2 and S6), from which further selections were made according to the purpose of the specific analysis. Daily streamflow and water temperature time series for the 2012-2022 water years were provided by the Estonian Environment Agency. Water quality data were ac quired from the Environmental Monitoring Information System database. Spatial data of catchment topography, geo -logy, hydrogeology, and land use/cover were extracted from various maps made publicly available at Geoportal by the Estonian Land Board. Monthly climatic data (2 m air tem -pera ture, total precipitation (P), total evaporation (E), snow depth (SDE), and snow water equivalent (SWE)) from the ERA5-Land dataset by Muñoz Sabater (2019) were also used. The time-variable data were averaged by water years (1 October - 30 September). All analyzed parameters and their abbreviations are listed and defined in Table S6.
2.3. Methods
2.3.1. Time series analysis
The analysis of daily streamflow time series recorded at the 42 gauging stations forms the basis of the study. In the fol -lowing, we describe the methodology for calculating the focus parameters and other relevant metrics. Detailed descriptions of all parameters discussed in this work are provided in the supplementary content. Time series analysis was performed using Microsoft Excel software (Microsoft Corporation USA), unless noted otherwise.
To characterize streamflow memory, ACF lags (k) were calculated using daily streamflow data. The ACF calculations were carried out using R, v. 4.2.3 (R Core Team, Austria), and RStudio, v. 2023.03.0+386 (Posit Team, USA), applying a threshold value of rk < 0.2, widely used in earlier studies (Koit et al. 2022; Mangin 1984; Schuler et al. 2022). The formula for ACF according to Larocque et al. (1998) is as follows: where rk is the autocorrelation at lag k, k is the time delay from 0 to maximum lag time, Ck is the covariance at lag k, C0 is the variance at lag k, n is the length of the studied time series, xt is the series value at time t, and x- is the mean value of the series.
The BFI was calculated using the approach of Ladson et al. (2013), applying the Lyne and Hollick (1979) recursive digital filter:
where qf (i) is the quickflow response at the i-th time step, q(i) is the total streamflow at the i-th time step, qb(i) is the base -flow at the i-th time step, and a (0.91) is the filter para meter that alters the shape of separation. The filter is applied mul -tiple times on the dataset. In the case of daily data, three passes (forward, backward, and forward) are commonly used. The BFI is defined as the ratio of baseflow volume to total streamflow volume over a specified period (water year).
The sensitivity of streamflow to modeled groundwater recharge was analyzed by calculating the CCF. The CCF quan tifies the linear dependency between two time series, xt and yt (Larocque et al. 1998). The interpretation of the CCF allows to assess the transfer of pressure impulse of ground -water recharge to the surface water catchment (Bailly-Comte et al. 2008; Mayaud et al. 2014; Worthington 2019). Accord -ing to Larocque et al. (1998), the formula for cross-correlation is as follows:
where Cxy(k) is the cross-correlogram, and sx and sy are the standard deviations of the time series xt and yt , respectively.
The groundwater recharge estimates were modeled using the Precipitation-Runoff Modeling System (PRMS-IV), devel -oped by Markstrom et al. (2015), as implemented by Hunt (2021) and Pärn et al. (2024). The PRMS was constructed for the Selja River catchment located in northern Estonia (Varangu gauging station in Fig. 1), and the modeled recharge estimates were used for all the assessed river catchments. A detailed description of the model's structure, parameters, and cali -bration results are described in Hunt (2021).
CCF correlation coefficients (r) and lag times (k) between the modeled groundwater recharge and streamflow were calculated for four different time periods: the entire dataset period (2012-2022; CCF mean r/k), each water year (CCF wy r/k), autumn floods (CCF autumn r/k), and spring floods (CCF spring r/k). The CCF calculations were performed using Python 3.9.13 (Rossum and Drake 2009) and the following libraries: Pandas (McKinney 2010), NumPy (Harris et al. 2020), scikit-learn (Pedregosa et al. 2011), Matplotlib (Hunter 2007), and SciPy (Virtanen et al. 2020).
The KarstID R Shiny application, developed by Cinkus et al. (2023), was used to carry out further discharge time series analyses on the whole period of 2012-2022 to calculate various parameters (regulation time, k max, a mean, IR, and SVC) for characterizing the hydrological functioning of the studied catchments. The regulation time, as defined by Larocque et al. (1998), was obtained from the spectral density function and quantifies the duration of the input signal's influence, providing insight into the impulse response length within the hydrological system. The indicator k max was ex -tracted from the analyzed recession curves and char acterizes the capacity of dynamic storage in the catch ment. The re -cession coefficient a mean characterizes the drainage of the catchment storage. The indicator IR allows to estimate a system's capacity to filter and attenuate the pre cipitation signal. The spring variability coefficient (SVC) is the ratio between Q90 and Q10 (Q10/Q90) that serves as another indicator of flow variability.
We calcu l a t e d t h e s e v e n -day mean annual minimum discharge (Q MAM7) as a more robust low-flow statistic in com parison to the one-day minimum. We used it in com -bination with the mean discharge (Q mean) to calculate a ratio (Q MAM7/Q mean), serving as an indicator of low-flow sensi tivity, similar to Stoelzle et al. (2020).
The definition of the number of drought days by Pärn and Mander (2012) was used as the count of days during a water year when the flow rate falls below 10% of the average flow rate during the summer half-year (April-September). The anal ysis was carried out in Python 3.9.13 (Rossum and Drake 2009), using the Pandas library (McKinney 2010).
To evaluate the sensitivity of streamflow to changes in precipitation and total evaporation, the median precipitation (eP-Q) and total evaporation (eE-Q) elasticity of runoff were calculated after Sankarasubramanian et al. (2001). The equa -tions for the precipitation elasticity of runoff (7) and total evaporation elasticity of runoff (8) are as follows:
where eP-Q is the precipitation elasticity of runoff (unitless), Q (mm) is runoff, P (mm) is precipitation, and E (mm) is evaporation. In the case of the latter, the total evaporation values from the ERA5-Land dataset were used. Water-yearly median elasticity values were calculated using monthly precipitation, evaporation, and streamflow amounts.
To assess the efficiency of the studied catchments, the effective catchment index (ECI) was calculated as described by Liu et al. (2020). The ECI is defined as follows:
where Q is the long-term mean discharge, P is the sum of precipitation, and AET is the actual evapotranspiration. In the case of the latter, the total evaporation values from the ERA5-Land dataset were used. Positive and negative ECI values indicate a streamflow excess or deficit, respectively, in relation to the climatic input. The average net quantity of ground -water imports and exports in the catchment is shown by the water balance deviation.
To further characterize the dominant trends in the stream -flow regime of the studied catchments, we calculated the linear slopes of the <Q10 (FDC <Q10) and Q50-Q100 (FDC Q50-Q100) segments of the FDC, using standardized daily discharge time series of 2012-2022. The FDC analysis was carried out in Python 3.9.13 (Rossum and Drake 2009), using the Pandas (McKinney 2010) and NumPy (Harris et al. 2020) libraries.
2.3.2. Catchment characteristics
All analyses with spatial data were performed using ArcGIS Pro 3.1.0 software (ESRI 2023). The mean and standard deviation (SD) of catchment elevation, and the stream gra dient were extracted from the 25-meter digital surface elevation model of Estonia provided by the Estonian Land Board. The percentage of land cover in the catchments was extracted from the Estonian Basic Map 1:10000 (Estonian Land Board).
The percentage coverage of the dominant Quaternary sedimentary cover types (alvar, glaciofluvial, glaciola cus -trine, peat, glacial, marine) was extracted from the 1:400000 Geological Map of Estonia (Estonian Land Board). Average Quaternary cover thickness values for all the catch ments were obtained from the Quaternary cover thickness layer of the hydrogeological model of the Baltic Artesian Basin by Virbulis et al. (2013).
The aquifer specific yield score was calculated based on the specific yield values reported in the 1:400000 hydro geo logical map of Estonia (Estonian Land Board). First, the areal coverage percentage of each specific yield zone in the catch ment was calculated. Each specific yield zone was then as -signed a weight and multiplied by its percentage coverage in a particular catchment (the weights and yield ranges are shown in Table S7). The resulting weight fractions were sum -med up as a unitless parameter. A higher value indicates better transmissivity of the catchment's aquifers, and vice versa. The groundwater protection score was calculated in the same way, based on the 1:400000 hydrogeological map of Estonia (Estonian Land Board). A higher value indicates a higher degree of protection of catchment's aquifers and vice versa.
By subtracting the groundwater head elevation, obtained from the 1:400000 hydrogeological map of Estonia (Estonian Land Board), from the surface elevation of the 25-meter digital surface elevation model of Estonia (Estonian Land Board), the mean and SD depth to the groundwater level of the first bedrock aquifer from the ground surface (GWL depth mean/SD) was calculated.
2.3.3. Multivariate analysis
The multivariate relationships of the mean values (2012-2022 water year mean) of 70 hydrological, spatial, and water quality parameters (listed in Table 1) in the 42 catchments were evaluated through factor analysis (n = 70 × 42). For further analysis, the catchments were clustered using agglomerative hierarchical clustering (AHC) based on multivariate similarities. The AHC employed Euclidean distance as the measure of dissimilarity between data points. The agglomer -ation technique of choice was Ward's method, which priori -tizes minimal increases in total within-cluster variance during cluster merging. To pinpoint the optimal number of clusters, we relied on the adapted Calinski and Harabasz index, which assesses clustering quality based on the ratio of between-cluster to within-cluster variance. Varimax rotation was ap -plied to improve the interpretability by maximizing the vari -ance of squared factors loadings by column. For a given factor, high loadings become higher and low loadings be -come lower. The analysis was carried out using the XLSTAT Forecast 2023.1.3.1407 software (Lumivero 2023).
2.3.4. LightGBM modeling
To understand the influence of selected hydroclimatic factors on our focus parameters (ACF wy k, BFI, CCF mean k, CCF mean r, CCF wy k, CCF wy r, CCF autumn k, CCF autumn r, CCF spring k, CCF spring r, FDC <Q10, and FDC Q50-Q100), the LightGBM (light gradient-boosting machine) was employed. LightGBM is a distributed, high-performance gradient boosting framework based on decision tree algo -rithms, designed specifically to be efficient and scalable (Ke et al. 2017). Prior to modeling, the dataset was standardized to have zero mean and unit variance. This ensures that all fea tures have equal scales, which can enhance the stability and interpretability of the model.
The LightGBM model was configured with a specific set of hyperparameters, including numleaves = 31, learningrate = 0.05, and nestimators = 1000. These were selected based on preliminary testing. The model was trained using the squared error (l2) as the objective function. A random subset con -taining 20% of the data was set aside for validation. The performance of the model was assessed on this validation set after training. Mean squared error (MSE) and R-squared (R²) were computed to evaluate the model's accuracy and pre -dictive capability.
To interpret the influence of each hydroclimatic factor, feature importance was assessed by evaluating the gain and split count metrics. Furthermore, to provide a more detailed breakdown of how each feature impacts each prediction, SHAP (SHapley Additive exPlanations) values were com puted. SHAP values are grounded in game theory and offer insights into how each feature contributes, either positively or negatively, to the predictive outcome (Lundberg and Lee 2017).
The entire analysis was conducted in Python 3.9.13 (Rossum and Drake 2009), leveraging the capabilities of LightGBM (Ke et al. 2017), scikit-learn (Pedregosa et al. 2011), and SHAP (Lundberg and Lee 2017) packages.
3. Results and discussion
3.1. General multivariate characterization of studied catchments
Two Varimax-rotated factor components, D1 and D2 (ex -plain ing 26.4% and 21% variability, respectively), were ex -tracted from a selection of 70 streamflow and catchment-derived parameters listed in Table 1. Notably, factor D1 exhibits strong loadings for parameters related to streamflow dy namics, such as Q CV, Q MAM7/Q mean, SVC, CCF mean r, 8EQ, drought days, SEC CV, k max, BFI, a mean, etc. The catch ments positively correlating with factor D1 are at lower elevations, flatter, and underlain by low- to medium-yielding carbonate rock aquifers with shallow groundwater levels. These corresponding catchments also have a higher propor -tion of peatlands, forests, peat soils, and artificial drainage, which likely contribute to higher mean values of chemical oxygen demand (CODMn) due to enhanced leaching from organic soils. Such catchments are characterized by the most sensitive hydrology (high values for a mean, CCF mean r, Q CV, SVC, eE-Q, drought days, SEC CV, FDC <Q10, etc.), and modest storage capacity (low k max, BFI). The catch -ments corresponding to D1 can be relatively efficient (high ECI), and in this respect, the catchments of the West Estonian Archipelago stand out in particular (they are more likely to gain additional groundwater contribution from neighboring surface catchments). The catchments with a strong negative correlation to D1 are generally situated at higher elevations, underlain by high-yielding sandstone aquifers, and feature a higher probability of glaciofluvial sediments and a larger proportion of BFI.
The factor D2 exhibits significant positive representations from variables such as ACF wy k, regulation time, FDC <Q10, CCF spring k, and baseflow (BF) runoff, and the occurrence of high-yielding carbonate rock and sandstone aquifers, which is indicative of greater regulation, storage capacity, and retention of water in the catchment. These char -acteristics contribute to greater streamflow memory. Due to the presence of more fertile soils (e.g., Luvisols, Mollic Cambrisols) in areas underlain by high-yielding carbonate aquifers, the catchments positively correlating with D2 usually feature a higher percentage of arable land, along with higher and more stable loads of dissolved solids (including total nitrogen, Ntot mean). On the other hand, the catchments with negative correlation have a thicker Quaternary sedi -mentary cover and are underlain by sandstone aquifers with a greater depth to the groundwater level (i.e., a thicker vadose zone). Quaternary cover thickness and deeper vadose zone also contribute to higher groundwater protection values.
The above-described 70 parameters (Table 1) were further analyzed using the AHC to classify the studied catchments. The AHC classified the 42 catchments into four clusters (Fig. 2; Table S1), which were named as follows: plateau (1), sand -stone upland (2), carbonate upland (3), and lowland (4).
The plateau (1) cluster comprises the catchments with relatively flat topography located on the Harju and Viru pla -teaus, and the neighboring plains (Fig. 2). In terms of average elevations, these watersheds fall between highlands and lowlands (Fig. 3a). The common denominator of this cluster is the articulation of the landscape (cuesta-like scarps, hil l -ocks, alvars, and mire basins, as described in Fig. 3), formed during the gradual retreat of the earlier stages of the Baltic Sea (Fig. 3e, g). The cluster, dominated by medium-yielding carbonate aquifers (Fig. 3h), exhibits a combination of hydrol -ogical characteristics resulting from the interaction between mire basins with poor infiltration capacity and elevated karstified bedrock hillocks (Fig. 3d). Due to the relatively thin Quaternary cover (Fig. 3c), groundwater is generally rather weakly protected (Fig. 3f). The important role of wetlands is also indicated by CODMn mean values that are slightly higher than in the upland clusters (as seen from the significant contribution from the D1 factor in Tables 1 and S1).
The sandstone upland (2) cluster groups together the more elevated (Fig. 3a) SE Estonian catchments (Fig. 2) with pro -nounced gradients (Fig. 3b), underlain by sandstone aquifers (Fig. 3h), and a relatively thick Quaternary cover (Fig. 3c). These catchments exhibit a somewhat more damped hydrol -ogical response due to the thick vadose zone and dominant intergranular porosity (Fig. 3d, h). Thus, these aquifers have relatively well protected groundwater (Fig. 3f).
The carbonate upland (3) cluster comprises the catch -ments of the Pandivere Upland (Fig. 2) in NE Estonia. These catchments are characterized by a significant contribution of high-yielding karst aquifers (Fig. 3d, h). Unlike the sandstone upland (2) cluster, the thickness of the Quaternary cover is significantly thinner, usually 2-10 m (Fig. 3c). Due to the higher mean altitude of the catchments (Fig. 3a), there are fewer glaciolacustrine sediments deposited, and the shares of poorly drained depressions and peatlands are lower com -pared to the plateau (1) and lowland (4) catchments. This leads to more efficient infiltration but also makes the aquifers more vulnerable, as manifested in groundwater quality (Fig. 3f). As the upland has preferable conditions for agriculture, there are problems with high nitrogen levels, as seen from the significant contribution from the D2 factor (Tables 1 and S1).
The lowland (4) cluster represents the flat and low-lying catchments (Fig. 3a, b) in W Estonian lowlands and the West Estonian archipelago (Fig. 2). These areas have only rel -atively recently emerged above the Baltic Sea level due to post-glacial isostatic land uplift. In this cluster, carbonate aquifers are dominant, especially in the higher areas of the islands, where they can be fairly high-yielding (Fig. 3d, h). Albeit with a similar thickness as in clusters 1 and 3 (Fig. 3c), the Quaternary sedimentary cover here is more clayey than in the higher areas, as evidenced by higher groundwater pro -tection scores (Fig. 3f), and marine sediments also occur more frequently (Fig. 3g). The land cover in the lowland (4) cluster is dominated by peatlands and forests on artificially drained peat soils (Fig. 3e, g). Therefore, the highest CODMn mean values occur in this cluster, as seen from the significant contribution from the D1 factor (Tables 1 and S1).
Nõges et al. (2022) also clustered 16 Estonian surface watersheds, using parameters such as catchment area, land use, land cover, population density, etc., resulting in three clusters. The authors then characterized the clusters with an em -phasis on the manifestations of anthropogenic impact through relevant nutrient loadings. The clustering was significantly different from ours: rivers in the same cluster in some cases did not have a clear resemblance, which would have resulted from the similar hydrogeology or hydrological behavior of the catchments. For example, some carbonate upland (3) and sandstone upland (2) rivers fell into the same cluster.
3.2. The main statistical relationships of the focus parameters
This chapter covers the general findings based on the 2012- 2022 mean values of the focus parameters (ACF, BFI, CCF, and FDC). The emphasis is set on the streamflow regime, its relationship with groundwater contributions to streamflow, and the sensitivity of streamflow to groundwater recharge signals. Spearman's correlation analysis was performed to evaluate the relationships of ACF, BFI, CCF, and FDC metrics with other streamflow time series and catchment-derived parameters, based on the mean values (n = 42) for the 2012-2022 water year period. The correlation analysis results are shown in Table S3.
3.2.1. Autocorrelation function
Catchments exhibit a memory effect in streamflow, closely linked to how they regulate discharge, as evidenced by the strong correlation between ACF wy k and regulation time (day) (p = 0.94, p < 0.05; Table S3) (Larocque et al. 1998). This effect, indicative of catchment memory, is amplified in systems with substantial effective storage and notable groundwater contributions, as demonstrated by positive correlations with BF runoff (p = 0.7, p < 0.05) and flatter slopes of FDC <Q10 (p = 0.83, p < 0.05). Interestingly, longer ACF lags align with extended spring CCF lags (p = 0.83, p < 0.05), yet are inversely related to CCF r (? = -0.68; ? = -0.63, p < 0.05 for both), pointing to intricate seasonal dynamics. Stable catchments exhibit a more pronounced streamflow memory, emphasized by a negative correlation with ACF lag CV (? = -0.67, p < 0.05). In addition, relationships with specific geological features are seen. A significant positive relationship was found between ACF lag and very high-yielding carbonate aquifers (? = 0.52, p < 0.05), while negative relationships were observed with thick Quaternary cover (? = -0.56, p < 0.05) and the presence of sandstone aquifers (? = -0.51, p < 0.05).
3.2.2. Baseflow index
Baseflow mitigates streamflow variability, as shown by the strong negative correlation of the BFI with measures such as Q CV (p = -0.97, p < 0.05) and SVC (p = -0.75, p < 0.05; Table S3). Higher BFI values, indicating increased ground -water contributions, buffer against extreme low flows, with minimum streamflow values approaching the mean (Q MAM7/ Q mean: p = 0.9, p < 0.05; Q min/Q mean: p = 0.88, p < 0.05). This stability, often a function of substantial catchment storage (positive correlation with k max, p = 0.83, p < 0.05), manifests in damped responses to external stimuli. The dampening effect is further reflected in correlations with CCF metrics (p = 0.64; p = 0.42; p = 0.42; p < 0.05 for all) and fewer drought days (p = -0.82, p < 0.05).
3.2.3. Cross-correlation function
There was a general inverse relationship between ACF, BFI, and BF runoff with respect to all four (mean, water year, autumn, and spring) CCF r metrics (Table S3). Thus, in catch -ments where the memory effect and baseflow contribution are greater (i.e., streamflow regime is more stable), there exists lower sensitivity to groundwater recharge events.
Similarly, FDC <Q10 and regulation time also showed negative correlations with all four CCF coefficient metrics (Table S3). These negative correlations suggest that catch -ments with damped high-flow responses have less pro -nounced overall hydrological responses.
The general relationship between CCF coefficient metrics and lags is inverse, indicating that as the strength of the as -sociation between the recharge and streamflow signals in -creases, the CCF lag decreases. ACF lags also tend to vary more as the CCF lag time decreases (Table S3). The water-yearly CCF wy r shows a significant negative correlation with spring floods, CCF spring r (p = -0.84, p < 0.05), followed by autumn floods, CCF autumn r (p = -0.68, p < 0.05). This supports the claim that spring floods, generally caused by snowmelt, are the most important hydrological events of the water year, as also showed by Koit (2022) and Koit et al. (2022).
3.2.4. Flow duration curve
Catchments with greater memory and longer regulation times exhibited flatter slopes in high-flow conditions (FDC <Q10) and relatively steeper slopes in median- to low-flow con ditions (FDC Q50-Q100). Such connections were evident in the upland catchments (clusters 2 and 3), suggesting a pro -longed release of water from relatively homogeneous reser voirs, such as productive aquifer systems, which sustain groundwater release over prolonged periods. In catchments with less permeable surfaces and limited storage capacity, steep FDC <Q10 slopes were observed, indicating rapid changes and greater variability in streamflow. However, in these catchments, the FDC "flattens out" in low-flow con -ditions (FDC Q50-Q100) because, after rapid runoff during high flows, the yield drops significantly, as the catchments lack sufficient storage to sustain baseflow. Seasonal CCF metrics (CCF spring r: ? = -0.71, p < 0.05; CCF autumn r: ? = -0.78, p < 0.05) showed strong negative correlations with FDC <Q10 and positive correlations with FDC Q50-Q100. CCF time lags, such as CCF wy k (water-yearly lag), showed a similar pattern, with a positive correlation with FDC <Q10 and varied correlations with FDC Q50-Q100. This suggests that catchments with longer CCF lags tend to have less vari -ability in streamflow.
3.3. Differences between the clusters
The catchments in the plateau (1) cluster exhibited pro -nounced memory effects in hydrological processes, with a streamflow ACF lag averaging around 41.3 days (Fig. 4a). The mean BFI value (0.65) suggests a relatively balanced contribution from both surface runoff and baseflow, while the contribution of the latter is still greater and more stable than in lowland (4) catchments (Fig. 4b). This consistency is reflected in the CCF coefficients, where the plateau (1) cluster shows moderate correlation strengths (Fig. 4d, f, h, j), par -ticularly in autumn and spring, and moderate lag times (k) across all seasons (Fig. 4c, e, g, i). This is because the catchments in this cluster are located at intermediate ele -vations, featuring a balanced land cover: a mix of perched peatland basins and thin-surfaced karstified carbonate hil l -ocks (Fig. 3e, g). The mean CCF lag is 2.7 days, second only to the lowland (4) cluster (Fig. 4c, d), which shows the relatively high hydrological sensitivity of the plateau catch -ments. Additionally, the plateau (1) cluster ranks second only to the lowland (4) cluster, based on the number of drought days (Fig. 5c).
The sandstone upland (2) catchments featured a relatively short mean memory effect, with an average streamflow ACF lag of 34.4 days (Fig. 4a). However, they featured a high BFI (Fig. 4b) of around 0.73, and high k max values (Fig. 6a), implying significant groundwater contribu tions from relatively permeable bedrock formations with great storage capacity. The latter somewhat contradict the unexpectedly low ACF lag values of the sandstone upland (2) cluster and the general theory behind the streamflow memory effect. This anomaly results from the peculiar hydrological regime of the sandstone upland (2) catchments (Fig. 6). During the spring floods, the biggest peaks of the water year, the response in these catchments can be as rapid as in the lowland (4) catchments. This is reflected by similar means in FDC <Q10 values (Fig. 4k) and the shortest spring CCF lags (Fig. 4i), which in the uplands is most likely caused by the greater gradients (Fig. 3b) and steeper slopes, as indicated by the great devi -ation in catchment elevation (Fig. 3a). In some cases, the relative spring quickflow flashiness is further amplified by the poor permeability of the thick Quaternary cover layer.
However, the catchments in the sandstone upland (2) exhibit by far the highest and most stable Q MAM7 values and the highest Q MAM7/Q mean ratio values (Figs 5d and 6b), even surpassing the water-rich carbonate upland (3) catch ments. This reflects the significant buffering effect of the sandstone aquifers, resulting in a stable and sustained baseflow drainage regime (Fig. 5d). At the same time, the sandstone upland (2) catchments are relatively ineffective (low ECI) and least affected by evapotranspiration during the summer seasons (Fig. 5b, e, f).
The carbonate upland (3) cluster stands out with the longest streamflow ACF lags, averaging around 48.5 days (Fig. 4a). The highest BFI value (Fig. 4b) among the clusters suggests catchments with significant groundwater storage capacity. Longer regulation times (see Table S2) signify an extended period for the system to reach equilibrium after a recharge event. This emphasizes the buffered dynamics of the catchments in the cluster. This cluster also shows the longest CCF lags (Fig. 4c, e, g, i) and the lowest CCF coefficients (Fig. 4d, f, h, j) across all seasons, signifying a delayed system response. These catchments receive an important ground -water contribution from high- to very high-yielding carbonate aquifers, covered by a relatively thin layer of Quaternary sediments (Fig. 3h, d). Infiltration is diffuse in these relatively flat but elevated and well-drained watersheds, which in turn contributes to the modest CCF response but also ensures en -hanced groundwater recharge capacity in the catch ments, reaching over 0.44 mm/d, as noted by Jaagus et al. (1998). Due to the stable baseflow regime, there were prac ti cally no drought days in either of the upland clusters (Fig. 5c).
Catchments in the lowland (4) cluster exhibit short stream -flow A C F lags averaging 34 5 days (F ig 4a) T he dom i nance of surface runoff is indicated by the lowest BFI, around 0.55 (Fig. 4b). The importance of surface runoff and the lower intensity of groundwater recharge in the catchments of the lowland (4) cluster, e.g., in Vihterpalu and Laadi, was also shown by Jaagus et al. (1998) in their model ing results. The lowland (4) cluster also showed the highest median CCF coefficients (r) across all seasons (Fig. 4d, f, h, j), demon -strating a strong correlation between the streamflow and recharge signal. The lowest average CCF lag (except in spring, when sandstone upland (2) catchments can be more responsive) suggests rapid catchment responses (Fig. 4c, e, g, i) due to the more impermeable geology and/or land cover domi nated by peatland basins. The surface runoff dominance is also indicated by the high efficiency of the catchments in the lowland (4) cluster (Fig. 5e, f) and the significantly strong correlation between the mean CCF coefficient and precipi tation elasticity to runoff (? = 0.88, p < 0.05; Table S4). The hydrological sensitivity of the low land (4) cluster is summed up by the highest mean number of drought days (Fig. 5c).
3.4. Determining hydroclimatic drivers of focus parameter variability
To evaluate the variability of the focus parameters between water years and their dependence on selected hydroclimatic parameters listed in Table 2, the LightGBM model was used.
The standardized (rescaled 0-1) SHAP values are shown in Table 3, and the observed versus predicted values of the focus parameters are shown in Fig. 7.
General hydroclimatic parameters, such as the sum of precipitation, evaporation, and associated hydrological ratios, showed low impacts on predictions for ACF lag and CCF metrics (Table 3). This suggests that while these parameters are the main drivers of hydrological processes, their direct short-term impact on streamflow characteristics may be less pronounced. The model placed considerable emphasis on air temperature and water temperature, as the latter is directly dependent on the former (? = 0.45, p < 0.05; Table S5), for capturing hydrological processes related to the focus para -meters (Tables 3 and S5). Air temperature, particularly its mini mum, maximum, and variation values, emerged as a pivotal factor across multiple focus parameters (Table 3).
Significant deviations in air temperature from mean values, trending towards warmer conditions, are associated with ACF lags, which is also reflected, e.g. in the overall importance of water temperature variation (Water temp CV in Table 3). Mean air temperature is closely related to the various aspects of snowpack: higher mean air tem pera -tures correspond to decreases in SDE and SWE (? = -0.67; ? = -0.71, p < 0.05 for both; Table S5). Higher air temperature (min, max, mean, and deviation) was positively related to longer ACF lags (? = 0.42, p < 0.05; Table S5), as clearly seen in the 2020 water year (Fig. 8). However, a water year with significant snowpack and abrupt spring snowmelt, such as 2013, produced shorter ACF lags (Fig. 8). The importance of snow-related parameters in Table 3, such as the mean SDE, SWE, and the previous water year's SWE, highlights the great influence of snowpack dynamics, compared to annual precipitation amounts, on predicting the streamflow char -acteristics quantified by the ACF, CCF, FDC, and BFI metrics. A weak but significant negative correlation (? = -0.11, p < 0.05; Table S5) hints that higher precipitation amounts may be associated with shorter ACF lags, possibly due to increased runoff and more frequent and quicker hydrological responses. At the same time, precipitation conditions of the previous water year, which affect antecedent moisture conditions, appear to be even more important than current precipitation amounts (? = -0.15; ? = -0.29, p < 0.05 for both; Table S5).
The results showed that the CCF metrics of the water year (CCF wy) are most closely related to the CCF properties (strength and lag) of the spring half-year (CCF spring), fol -lowed by autumn floods (CCF autumn) (Tables 3 and S5). This suggest that significant snowmelt events, primarily in spring but sometimes also in the autumn season (early snowmelt), are decisive in shaping the hydrological regime.
In the example of the wet but snowless 2020 water year (Fig. 8), the scattered CCF and ACF metrics indicate greater variability in hydrological responses. This is especially evi -dent in the CCF spring responses (Fig. 8), reflecting the broader dispersion of responses due to the lack of a strong and uniform snowmelt signal. As winter precipitation is projected to increase (Coppola et al. 2021; Jaagus and Mändla 2014; Ruosteenoja and Jylhä 2021), and rising air tempera -tures (2.5-3.8 °C in winter and 2.0-3.4 °C in summer) are expected to reduce snowpack while increasing water surplus in colder seasons (Barnett et al. 2005; Okkonen et al. 2011), hydrological variability is expected to increase in the coming decades.
Preceding moisture conditions (previous water year SWE and precipitation) can also have a significant effect (Table 3). Snowpack and snowmelt dynamics can affect a wide variety of parameters in the following water year. Koit et al. (2022) showed that autumn CCFs best highlighted and characterized the sensitivity and responsiveness of aquifers interacting with surface catchments, as supported by the relatively higher CCF autumn r values (Fig. 8). This is because in autumn, aquifers and associated catchments recover from the summer low flow/drought period without the masking effect of snowmelt (see 2013 and 2018 vs. 2020 in Fig. 8). While seasonal CCFs provide valuable insights into hydroclimatic variability, the mean CCFs (CCF mean), calculated for the entire period of 2012-2022, were less sensitive to calculation anomalies, capturing the overall hydrological sensitivity of a catchment with sufficient detail. This makes them more suitable for catchment classification.
As seen from the SHAP values for FDC <Q10 (Table 3), high-flow conditions are predominantly influenced by snow -pack conditions and air temperature. The more negative FDC <Q10 values during the 2013 water year (Fig. 8), which correspond to steeper limbs of flood peaks, can be attributed to abrupt spring snowmelt and ice-damming (Estonian Environment Agency 2014). The sandstone upland (2) catch -ments stand out in particular, exhibiting the most negative values.
During the snow-rich 2013 water year, the mean BFI value approached 0.5, indicating a greater contribution of surface runoff (Fig. 8). The decrease in BFI values was par -ticularly significant in the lowland (4) cluster, which testifies to the greater role of surface runoff during snowy years (in the case of abrupt snowmelt). Conversely, in the dry 2018 water year, baseflow became more prominent (Fig. 8), as seen from elevated BFI values. Low to mid-flow con ditions (FDC Q50-Q100 in Table 3) were significantly in fluenced by the hydro logical balance and extreme con ditions, such as drought. In dry years (2013 and 2018), FDC Q50-Q100 slopes were more likely to be flat (fewer negative values), especially in the lowland (4) cluster (Fig. 8). Flatter FDC Q50-Q100 slopes occur during longer periods when water bodies rely mainly on baseflow, or even dry up completely. This was also reflected in BFI values (Table 3), which were influenced by extreme conditions (drought days) and water temperature variability.
The hydrological sensitivity of Estonian catchments is most clearly manifested in the influence exerted by winter-spring snow dynamics and droughts, as demonstrated by the influence of air temperature and precipitation deviations, along with the number of drought days, on focus parameters. This aligns with previous studies in Estonia, which have analyzed the water balance of Estonian catchments and long-term interdependence between atmospheric conditions and river runoff (Roosaare et al. 1998; Jaagus et al. 2017; Kotta et al. 2018). These studies found that historical changes in Estonian river runoff generally coincided with variations in air temperature and snow cover duration, caused by changes in atmospheric circulation during the winter season. In addition, it has been shown that long-term (multi-decadal) regime shifts in river runoff corresponded to the alternation of wet and dry periods (Jaagus et al. 2017). As snow cover duration in the region is projected to decrease and air tem -peratures increase, climate change will thus lead to significant alterations in the hydrological regime in the region, including a decrease in spring floods and higher discharge in autumn and winter due to winter rains (Meier et al. 2022). However, as shown by this study, the direct effects of temperature and precipitation changes on specific catchment runoff remain unclear, as they will act in combination and depend strongly on the local factors (Stahl et al. 2010).
4. Conclusions
We a ssessed the i ntricate i n t e rplay betw e en h y dro c l i matic factors and hydrological sensitivity across different catchment types in Estonia. By employing advanced statistical and ma -chine learning techniques, such as LightGBM modeling sup -ported by SHAP value interpretations, the study revealed the nuanced influences of air temperature, snowpack dynamics, and preceding hydrological conditions on stream flow char -acteristics.
Catchments with significant groundwater contributions and effective storage, particularly in the upland clusters, showed slightly different behavior in terms of memory effects and hydrological responses depending on geological con -ditions, while other parameters rather supported their high buffering capacity and drought resistance. Conversely, low -land cluster catchments, where soils with poor infiltration capacity, peatlands, and artificial drainage are widespread, displayed higher sensitivity and quicker responses to hydro -climatic variability.
The study also underscores the importance of air tem -perature and snowpack conditions in modulating streamflow, particularly in shaping high-flow conditions and baseflow contributions. It suggests that air temperature, through its influence on water temperature and snowpack dynamics, has a strong bearing on streamflow responses.
Additionally, the results demonstrate that the short- to medium-scale hydrological sensitivity of catchments is not solely dictated by immediate precipitation and evaporation dynamics but is also significantly influenced by antecedent moisture conditions from the previous water year. This under -standing has profound implications for water resource man -age ment, especially in the face of climate variability and change.
This study contributes to the broader understanding of hydrological processes by revealing the varying degrees of catchment sensitivity to hydroclimatic changes, serving as a valuable resource for future research and policy-making in hydrological and environmental planning. The nuanced char -acterizations of catchment responses derived from this re -search can inform targeted strategies to mitigate the impacts of extreme hydrological events, thereby enhancing the re -silience of water resource systems against the backdrop of a changing climate. These findings will also help in designing measures to preserve and achieve the good status of ground -water and surface water bodies, as set out in the EU Water Framework Directive (European Commission 2000), under evolving climatic conditions. Some detailed considerations:
o ACF lag is a robust and informative parameter that in -tegrates hydroclimatic, hydrologic, and hydrogeologic influences, while being relatively easy to calculate. How -ever, other supporting parameters are needed to fully understand the hydrological behavior of a watershed.
o Although the length of the ACF lag can be generally interpreted as a measure of hydrological inertia, there may be deviations from the general understanding in certain cases, as was evident in the sandstone upland (2) catch -ments. In these cases, other parameters must also be used to support the interpretation of the memory effect and to highlight the characteristics of the hydrograph, which may remain hidden when relying solely on ACF data.
o Contrary to intuition, the longest ACF lags were observed during dry and droughty water years, not necessarily in water years with lots of precipitation, storm events, or snowmelt that typically generate the most groundwater recharge. The long lags would be caused by the uninter -rupted length of streamflow recession during extended low flow periods. In contrast, wet years with repeated flood events alternate with shorter recession periods, producing shorter overall ACF lags.
o CCF is excellent in showing streamflow response/ sensitivity to recharge events; however, their automatic calculation requires precaution. Calculation results be -come more uncertain as the observation period shortens. The timing of individual flood events and their position -ing in the calculation "window" has a decisive impact on the calculation of seasonal CCF.
o The FDC segments, especially the extremes such as FDC <Q10 and the median-high flows (FDC Q50-Q100), provide a stark contrast in catchment responses to high-and low-flow conditions. These parameters indicate the variability and stability of streamflow, with high-flow conditions influenced by factors such as snowmelt and air temperature, while low- to mid-flow conditions reflect the hydrological balance and the catchment's ability to sustain baseflow during dry periods.
o It is constructive to note that the sensitivity of these parameters to hydroclimatic factors can significantly af fect water resource management strategies. For instance, a catchment with a high BFI and prolonged ACF lags suggests a system with significant groundwater influence, likely to be resilient to short-term climatic variations but potentially vulnerable to long-term climate shifts that could alter groundwater recharge patterns. In contrast, catchments with short CCF lags and steep FDC slopes during high flows may require careful monitoring for flood risks, while those with flatter FDC curves during low flows will need strategies to combat drought impacts.
Data availability statement
Data not already included in the paper and its supplementary materials will be made available upon request.
Acknowledgments
This study was financed by the project LIFE21-IPC-EE-LIFE-SIP AdaptEst/101069566 "Implementation of national climate change adaptation activities in Estonia". We also thank the reviewers and editors of the manuscript for their insightful feedback that greatly improved it. The publication costs of this article were partially covered by the Estonian Academy of Sciences.
Supplementary online data
Supplementary online data to this article can be found at
https://doi.org/10.3176/earth.2025.01S and include tables detailing the hydrol ogical characteristics and statistical analyses of the assessed river catchments:
Tab le S1. Hydrological characteristics, flow metrics, and factor scores (D1 & D2) for clustering catchments Table S2. Descriptive statistics of parameters for the total sample and by cluster
Table S3. Spearman correlation matrix of mean parameter values (2012-2022)
Table S4. Spearman correlation matrix of mean parameter values by cluster (2012-2022)
Table S5. Spearman correlation matrix of water year means for focus and climatic parameters (2012-2022) Table S6. Definitions of assessed parameters
Table S7. Input parameters for specific yield score calculations
These tables support the classification, correlation, and interpretation of the hydrological characteristics of the studied river catchments.
References
Babre, A., Popovs, K., Kalvans, A., Jemeljanova, M. and Dëlina, A. 2023. Forecasting the groundwater levels in the Baltic through standardized index analysis. Weather and Climate Extremes, 45, 100728. https://doi.org/10.1016/j.wace.2024.100728
Bailly-Comte, V, Jourde, H, Roesch, A., Pistre, S. and Batiot-Guilhe, C. 2008. Time series analyses for karst/river interactions assess -ment: case of the Coulazou river (southern France). Journal of Hydrology, 349(1-2), 98-114. https://doi.org/10.1016/jjhydrol. 2007.10.028
Barnett, T. P, Adam, J. C. and Lettenmaier, D. P 2005. Potential impacts of a warming climate on water availability in snow-dominated regions. Nature, 438, 303-309.
Beck, H. E., Zimmermann, N. E., McVicar, T. R., Vergopolan, N., Berg, A. and Wood, E. F. 2018. Present and future Köppen-Geiger climate classification maps at 1 -km resolution. Scientific Data, 5, 180214. https://doi.org/10.1038/sdata.2018.214
Berghuijs, W. R., Woods, R. A. and Hrachowitz, M. 2014. A pre -cipitation shift from snow towards rain leads to a decrease in streamflow. Nature Climate Change, 4, 583-586. https://doi.org/ 10.1038/nclimate2246
Blöschl, G., Hall, J., Viglione, A., Perdigão, R. A. P, Parajka, J., Merz, B. et al. 2019. Changing climate both increases and de -creases European river floods. Nature, 573, 108-111. https://doi.org/ 10.1038/s41586-019-1495-6
Cinkus, G., Mazzilli, N. and Jourde, H. 2023. KarstID: an R Shiny application for the analysis of karst spring discharge time series and the classification of karst system hydrological functioning. Environmental Earth Sciences, 82, 136. https://doi.org/10.1007/ s12665-023-10830-5
Cochand, M., Christe, P, Ornstein, P and Hunkeler, D. 2019. Groundwater storage in high alpine catchments and its con -tribution to streamflow. Water Resources Research, 55(4), 2613-2630. https://doi.org/10.1029/2018WR022989
Coppola, E., Nogherotto, R., Ciarlo', J. M., Giorgi, F, van Meijgaard, E., Kadygrov, N. et al. 2021. Assessment of the European climate projections as simulated by the large EURO-CORDEX regional and global climate model ensemble. Journal of Geophysical Research: Atmospheres, 126(4). https://doi.org/10.1029/2019JD03 2356
Costantini, M., Colin, J. and Decharme, B. 2023. Projected climate-driven changes of water table depth in the world's major groundwater basins. Earth's Future, 11(3). https://doi.org/10.10 29/2022EF003068
Donnelly, C, Greuell, W., Andersson, J., Gerten, D., Pisacane, G., Roudier, P and Ludwig, F. 2017. Impacts of climate change on European hydrology at 1.5, 2 and 3 degrees mean global warming above preindustrial level. Climatic Change, 143, 13-26. https:// doi.org/10.1007/s10584-017-1971-7
Douville, H, Raghavan, K., Renwick, J., Allan, R. P, Arias, P A., Barlow, M. et al. 2021. Water cycle changes. In Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change (Masson-Delmotte, V, Zhai, P, Pirani, A., Connors, S. L., Péan, C, Berger, S. et al., eds). Cambridge University Press, Cambridge, New York, 1055-1210. https://doi.org/doi:10. 1017/9 781009157896.010
Earman, S. and Dettinger, M. 2011. Potential impacts of climate change on groundwater resources - a global review. Journal of W a t e r & C l i m a t e C h a n g e , 2(4), 213-229. https://doi.org/10.21 66/wcc.2011.034
ESRI. 2023. ArcGIS Pro. https://www.esri.com/en-us/arcgis/products/ arcgis-pro/overview (accessed 2024-10-02).
Estonian Environment Agency. Environmental Monitoring Information System. https://kese.envir.ee/kese/welcome.action (accessed 2024-10-02).
Estonian Environment Agency. https://keskkonnaagentuur.ee/ (accessed 2024-10-02).
Estonian Environment Agency. 2014. Hydrological Yearbook 2013. Estonian Environment Agency, Tallinn.
Estonian Land Board. Estonian Land Board Geoportal. https:// geoportaal.maaamet.ee (accessed 2024-10-02).
European Commission. 2000. Directive 2000/60/EC of the European Parliament and of the Council of 23 October 2000 establishing a framework for community action in the field of water policy. Official Journal of the European Communities, L 327, 1-72.
Grogan, D. S., Burakowski, E. A. and Contosta, A. R. 2020. Snowmelt control on spring hydrology declines as the vernal window lengthens. Environmental Research Letters, 15, 114040. https://doi.org/10.1088/1748-9326/abbd00
Harris, C. R., Millman, K. J., van der Walt, S. J., Gommers, R., Vir tanen, P. , Co ur na peau, D. et a l. 2 020. Arra y progr ammi ng w ith NumPy. Nature, 585, 357-362. https://doi.org/10.1038/s41586-020-2649-2
Hunt, M. 2021. Modeling of the water balance in the Selja River basin with the PRMS hydrological model. Master's thesis. University of Tartu, Estonia.
Hunter, J. D. 2007. Matplotlib: a 2D graphics environment. Computing in Science & Engineering, 9(3), 90-95. https://doi.org/ 10.1109/MCSE.2007.55
IPCC (Intergovernmental Panel on Climate Change). 2022. Climate Change 2022: Impacts, Adaptation, and Vulnerability. Con -tribution of Working Group II to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge University Press, Cambridge, New York.
Jaagus, J. and Mändla, K. 2014. Climate change scenarios for Estonia based on climate models from the IPCC Fourth Assessment Report. Estonian Journal of Earth Sciences, 63(3),
Jaagus, J., Järvet, A., Roosaare, J., Tamm, T. and Vallner, L. 1998. Integrated assessment of climate change impact on water re -sources in Estonia. In Country Case Study on Climate Change Impacts and Adaptation Assessments in the Republic of Estonia (Tarand, A. and Kallaste, T., eds). SEI-T & UNEP, Tallinn, 59-76.
Jaagus, J., Sepp, M., Tamm, T., Järvet, A. and Mõisja, K. 2017. Trends and regime shifts in climatic conditions and river runoff in Estonia during 1951-2015. Earth System Dynamics, 8(4), 963- 976. https://doi.org/10.5194/esd-8-963-2017
Järvet, A. 1998. Long-term changes in time series of water balance elements. In Country Case Study on Climate Change Impacts and Adaptation Assessments in the Republic of Estonia (Tarand, A. and Kallaste, T., eds). SEI-T & UNEP, Tallinn, 69-71.
Jefferson, A., Nolin, A., Lewis, S. and Tague, C. 2008. Hydro -geologic controls on streamflow sensitivity to climate variation. Hydrological Processes, 22(22), 4371-4385. https://doi.org/10. 1002/hyp.7041
Jenicek, M., Seibert, J., Zappa, M., Staudinger, M. and Jonas, T. 2016. Importance of maximum snow accumulation for summer low flows in humid catchments. Hydrology and Earth System Sciences, 20(2), 859-874. https://doi.org/10.5194/hess-20-859-2016
Kalm, V. 2006. Pleistocene chronostratigraphy in Estonia, south -eastern sector of the Scandinavian glaciation. Quaternary Science Reviews, 25(9-10), 960-975. https://doi.org/10.1016/j.quascirev. 2005.08.005
Ke, G., Meng, Q., Finley, T., Wang, T., Chen, W., Ma, W. et al. 2017. LightGBM: a highly efficient gradient boosting decision tree. In 31st International Conference on Neural Information Processing Systems, Long Beach, California, USA, 4-9 December 2017 (Luxburg, U. von, Guyon, I., Bengio, S., Wallach, H. and Fergus, R., eds). Curran Associates Inc., New York, 3149-3157.
Killian, C. D., Asquith, W. H., Barlow, J. R. B., Bent, G. C., Kress, W. H., Barlow, P. M. and Schmitz, D. W. 2019. Char -acterizing groundwater and surface-water interaction using hydro -graph-separation techniques and groundwater-level data through -out the Mississippi Delta, USA. Hydrogeology Journal, 27, 2167- 2179. http://dx.doi.org/10.1007/s10040-019-01981-6
Koit, O. 2022. Surface water and groundwater interaction in shallow karst aquifers of Lower Estonia. PhD thesis. Tallinn University, Estonia. https://doi.org/10.13140/RG.2.2.13634.86726
Koit, O., Mayaud, C., Kogovšek, B., Vainu, M., Terasmaa, J. and Marandi, A. 2022. Surface water and groundwater hydraulics of lowland karst aquifers of Estonia. Journal of Hydrology, 610, 127908. https://doi.org/10.1016/j.jhydrol.2022.127908
Kotta, J., Herkül, K., Jaagus, J., Kaasik, A., Raudsepp, U., Alari, V. et al. 2018. Linking atmospheric, terrestrial and aquatic e n -viron ments: regime shifts in the Estonian climate over the past 50 years. PLoS ONE, 13(12). https://doi.org/10.1371/journal. pone.0209568
Kottek, M., Grieser, J., Beck, C., Rudolf, B. and Rubel, F. 2006. World map of the Köppen-Geiger climate classification updated. Meteorologische Zeitschrift, 15(3), 259-263. http://dx.doi.org/10. 1127/0941-2948/2006/0130
Ladson, A. R., Brown, R., Neal, B. and Nathan, R. 2013. A standard approach to baseflow separation using the Lyne and Hollick filter. Australasian Journal of Water Resources, 17(1), 25-34. https:// doi.org/10.7158/13241583.2013.11465417
Larocque, M., Mangin, A., Razack, M. and Banton, O. 1998. Con -tribution of correlation and spectral analyses to the regional study of a large karst aquifer (Charente, France). Journal of Hydrology, 205(3-4), 217-231. https://doi.org/10.1016/S0022-1694(97)00155-8
Liu, Y., Wagener, T., Beck, H. E. and Hartmann, A. 2020. What is the hydrologically effective area of a catchment? Environmental Research Letters, 15, 104024. https://doi.org/10.1088/1748-9326/aba7e5
Lumivero. 2023. XLSTAT statistical and data analysis solution. https://www.xlstat.com/en (accessed 2024-10-02).
Lundberg, S. M. and Lee, S.-I. 2017. A unified approach to inter -preting model predictions. In 31st Inter national Conference on Neural Information Processing Systems, Long Beach, California, USA, 4-9 December 2017 (Luxburg, U. von, Guyon, I., Bengio, S., Wallach, H. and Fergus, R., eds). Curran Associates Inc., New York, 4765-4774.
Lyne, V. D. and Hollick, M. 1979. Stochastic time-variable rainfall-runoff modeling. In Hydrology and Water Resources Symposium, Perth, Australia, 10-12 September 1979. Institution of Engineers National Conference, 89-92.
Mangin, A. 1984. The use of autocorrelation and spectral analyses to obtain a better understanding of hydrological systems. Journal of Hydrology, 67(1-4), 25-43.
Markstrom, S. L., Regan, R. S., Hay, L. E., Viger, R. J., Webb, R. M. T., Payn, R. A. and LaFontaine, J. H. 2015. PRMS-IV, the precipitation-runoff modeling system, version 4. In Te c h n i q u e s a n d M e t h o d s 6 - B 7 . US Geological Survey, Reston. https://doi.org/10.3133/tm6B7
Mayaud, C., Wagner, T., Benischke, R. and Birk, S. 2014. Single event time series analysis in a binary karst catchment evaluated using a groundwater model (Lurbach system, Austria). Journal of Hydrology, 511, 628-639. https://doi.org/10.1016/j.jhydrol.2014.02.024
McKinney, W. 2010. Data structures for statistical computing in Python. In 9th Python in Science Conference, Austin, USA, 28 June - 3 July 2010 (van der Walt, S. and Millman, J., eds). scipy.org, 51-56.
Meier, H. E. M., Kniebusch, M., Dieterich, C., Gröger, M., Zorita, E., Elmgren, R. et al. 2022. Climate change in the Baltic Sea region: a summary. Earth System Dynamics, 13(1), 457-593. https://doi.org/ 10.5194/esd-13-457-2022
Meriö, L.-J., Ala-aho, P., Linjama, J., Hjort, J., Kløve, B. and Marttila, H. 2019. Snow to precipitation ratio controls catchment storage and summer flows in boreal headwater catchments. Wa t e r Resources Research, 55(5), 4096-4109. https://doi.org/10.1029/ 2018WR023031
Muñoz Sabater, J. 2019. ERA5-Land monthly averaged data from 1950 to present. Copernicus Climate Change Service (C3S) Climate Data Store (CDS). https://doi.org/10.24381/cds.68d2bb30
Nõges, T., Vilbaste, S., McCarthy, M. J., Tamm, M. and Nõges, P. 2022. Long-term data reflect nitrogen pollution in Estonian rivers. Hydrology Research, 53(12), 1468-1479. https://doi.org/10.2166/ nh.2022.057
Nygren, M., Giese, M., Kløve, B., Haaf, E., Rossi, P. M. and Barthel, R. 2020. Changes in seasonality of groundwater level fluctuations in a temperate-cold climate transition zone. Journal of Hydrology X, 8, 100062. https://doi.org/10.1016/j.hydroa.2020. 100062
Okkonen, J. and Kløve, B. 2011. A sequential modelling approach to assess groundwater-surface water resources in a snow domi -nated region of Finland. Journal of Hydrology, 411(1-2), 91-107. https://doi.org/10.1016/j.jhydrol.2011.09.038
Okkonen, J., Jyrkama, M. and Kløve, B. 2011. A conceptual ap -proach for assessing the impact of climate change on groundwater and related surface waters in cold regions (Finland). Hydro -geology Journal, 18, 429-439. http://dx.doi.org/10.1007/s10040-009-0529-9
Pärn, J. and Mander, Ü. 2012. Increased organic carbon concentra -tions in Estonian rivers in the period 1992-2007 as affected by deepening droughts. Biogeochemistry, 108, 351-358. https://doi.org/ 10.1007/s10533-011-9604-0
Pärn, J., Walraevens, K., Hunt, M., Koit, O., van Camp, M., Ivask, J. et al. 2024. Unveiling the hydrological response of NO3-rich springs to seasonal snowmelt in a karstic carbonate upland. Journal of Hydrology, 641, 131724. https://doi.org/10.1016/j. jhydrol.2024.131724
Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V, Thirion, B., Grisel, O. et al. 2011. Scikit-learn: machine learning in Python. Journal of Machine Learning Research, 12(85), 2825-2830. Protasjeva, M. S. and Eipre, T. (eds). 1972.PecypcunoeepxHocmHbix eod CCCP (Surface Water Resources of USSR). Gidrometeoizdat, Leningrad. Ranasinghe, R., Ruane, A. C, Vautard, R., Arnell, N., Coppola, E., Cruz, F. A. et al. 2021. Climate change information for regional impact and for risk assessment. In Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change (Masson-Delmotte, V, Zhai, P, Pirani, A., Connors, S. L., Péan, C, Berger, S. et al, eds). Cambridge University Press, Cambridge, New York, 1767-1926. https://doi.org/10.1017/978 1009157896.014 Raukas, A. and Kajak, K. 1997. Quaternary cover. In Geology and Mineral Resources of Estonia (Raukas, A. and Teedumäe, A., eds). Estonian Academy Publishers, Tallinn, 125-136.
Rodhe, A. 1998. Snowmelt-dominated systems. In Isotope Tracers in Catchment Hydrology (Kendall, C. and McDonnell, J. J., eds). Elsevier, Amsterdam, 391-434.
Roosaare, J., Jaagus, J. and Järvet, A. 1998. Modelling the influence of climate change on river runoff. In Country Case Study on Climate Change Impacts and Adaptation Assessments in the Republic of Estonia (Tarand, A. and Kallaste, T., eds). SEI-T & UNEP, Tallinn, 75-83.
Va n R o s s u m , G . a n d D r a k e , F . L . 2 0 0 9 . Python 3 Reference Manual. CreateSpace, Scotts Valley.
Ruosteenoja, K. and Jylhä, K. 2021. Projected climate change in Finland during the 21st century calculated from CMIP6 model simulations. Geophysica, 56(1-2), 39-69.
Ruosteenoja, K., Markkanen, T., Venäläinen, A., Räisänen, P. and Peltola, H. 2018. Seasonal soil moisture and drought occurrence in Europe in CMIP5 projections for the 21st century. Climate Dynamics, 50 (3-4), 1177-1192. https://doi.org/10.1007/s00382-017-3671-4
Sankarasubramanian, A., Vogel, R. M. and Limbrunner, J. F. 2001. Climate elasticity of streamflow in the United States. Wa t e r Resources Research, 37(6), 1771-1781. https://doi.org/10.1029/ 2000WR900330
Schuler, P., Campanyà, J., Moe, H., Doherty, D., Williams, N. H. and McCormack, T. 2022. Mapping the groundwater memory across Ireland: a step towards a groundwater drought sus -ceptibility assessment. Journal of Hydrology, 612, 128277. https://doi.org/10.1016/j.jhydrol.2022.128277
Smerdon, B. D. 2017. A synopsis of climate change effects on ground water recharge. Journal of Hydrology, 555, 125-128. https://doi.org/10.1016/j.jhydrol.2017.09.047
Stahl, K., Hisdal, H., Hannaford, J., Tallaksen, L. M., van Lanen, H. A. J., Sauquet, E. et al. 2010. Stream flow trends in Europe: evidence from a dataset of near-natural catchments. Hydrology and Earth System Sciences, 14(12), 2367-2382. https://doi.org/10.5194/ hess-14-2367-2010
Statistics Estonia. 2023. Population Demographics 2023. https://www. stat.ee (accessed 2024-10-02).
Stoelzle, M., Schuetz, T., Weiler, M., Stahl, K. and Tallaksen, L. M. 2020. Beyond binary baseflow separation: a delayed-flow index for multiple streamflow contributions. Hydrology and Earth System Sciences, 24(2), 849-867. https://doi.org/10.5194/hess-24-849-2020
Sutanto, S. J. and van Lanen, H. A. J. 2022. Catchment memory explains hydrological drought forecast performance. Scientific Reports, 12, 2689. https://doi.org/10.1038/s41598-022-06553-5
Tay l o r, R . G ., S c anlon , B., D öl l , P. , R o dell , M . , B e e k , R . van, Wa da, Y. et al. 2013. Ground water and climate change. Nature Climate Change, 3, 322-329. http://dx.doi.org/10.1038/nclimate1744
Teutschbein, C., Grabs, T., Karlsen, R. H., Laudon, H. and Bishop, K. 2015. Hydrological response to changing climate conditions: spatial streamflow variability in the boreal region. Wa t e r R e s o u r c e s R e s e a r c h , 51(12), 9425-9446. https://doi.org/ 10.1002/2015WR017337
Teutschbein, C., Quesada Montano, B., Todorovic, A. and Grabs, T. 2022. Streamflow droughts in Sweden: spatiotemporal patterns emerging from six decades of observations. Journal of Hydrology: Regional Studies, 42, 101171. https://doi.org/10.1016/j.ejrh.2022. 101171
Vallner, L. 1998. Assessment of climate change impact on ground -water. In Country Case Study on Climate Change Impacts and Adaptation Assessments in the Republic of Estonia (Tarand, A. and Kallaste, T., eds). SEI-T & UNEP, Tallinn, 83-89.
Vir b ulis, J. , Be thers, U., Sak s, T. , Se nnik o vs, J. and Timuhins, A. 2013. Hydrogeological model of the Baltic Artesian Basin. Hydro -geology Journal, 21, 845-862. https://doi.org/10.1007/s10040-013-0970-7
Vir tanen , P. , Gomme rs, R., Olipha nt, T. E., Haber lan d, M., Reddy T., Cournapeu, D. et al. 2020. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature Methods, 17, 261-272.
Vi ru , B. a n d Jaagus, J. 2020. Spatio-temporal variability and seasonal dynamics of snow cover regime in Estonia. Theoretical and Applied Climatology, 139(1-2), 759-771. https://doi.org/10. 1007/s00704-019-03013-5
Vogel, R. M. and Fennessey, N. M. 1994. Flow duration curves. I: new interpretation and confidence intervals. Journal of Water Resources Planning and Management, 120(4), 485-504.
Vorobevskii, I ., Kronenberg, R. and Bernhofer, C. 2022. Linking different drought types in a small catchment from a statistical perspective - case study of the Wernersbach catchment, Germany. Journal of Hydrology X, 15, 100122. https://doi.org/10.1016/j. hydroa.2022.100122
Worthington, S. R. H. 2019. How preferential flow delivers pre-event groundwater rapidly to streams. Hydrological Processes, 33(17), 2373-2380. https://doi.org/10.1002/hyp.13520
You have requested "on-the-fly" machine translation of selected content from our databases. This functionality is provided solely for your convenience and is in no way intended to replace human translation. Show full disclaimer
Neither ProQuest nor its licensors make any representations or warranties with respect to the translations. The translations are automatically generated "AS IS" and "AS AVAILABLE" and are not retained in our systems. PROQUEST AND ITS LICENSORS SPECIFICALLY DISCLAIM ANY AND ALL EXPRESS OR IMPLIED WARRANTIES, INCLUDING WITHOUT LIMITATION, ANY WARRANTIES FOR AVAILABILITY, ACCURACY, TIMELINESS, COMPLETENESS, NON-INFRINGMENT, MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. Your use of the translations is subject to all use restrictions contained in your Electronic Products License Agreement and by using the translation functionality you agree to forgo any and all claims against ProQuest or its licensors for your use of the translation functionality and any output derived there from. Hide full disclaimer
© 2025. This work is published under https://creativecommons.org/licenses/by/4.0/ (the “License”). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.
Abstract
The escalating impacts of global climate change significantly affect regional hydrological systems, particularly in northern areas such as Estonia. This study investigates the hydro l ogical sensitivity of Estonian catchments to climatic variability, focusing on the interplay between surface water and groundwater. Using data from 42 river catchments, it employs various statistical methods in hydrology, emphasizing the autocorrelation function, cross-cor relation function, baseflow index, and flow duration curve. The analysis spans the years 2012-2022, integrating hydrological, spatial, and water quality parameters. The research identifies four distinct hydrological behavior clusters: plateau, sandstone upland, carbonate upland, and lowland. Key findings include diverse catchment sensitivities to groundwater recharge, the role of baseflow in streamflow stabilization, the memory effect in catchment responses, and insights from the flow duration curve on flow variability and extremes. The LightGBM model, predicting focus parameters, highlights the critical influence of air tempera ture and snowpack on streamflow characteristics. This study underscores the diverse hydro logical sensitivities of Estonian catchments to hydroclimatic changes, emphasizing the impor-tance of considering catchment-specific characteristics in water resource management and policy-making. Con -tributing to the broader understanding of hydrological processes, it provides valuable insights for future research and environmental planning in the face of climate variability and change.