Introduction

Climate change is projected to impact groundwater resources globally (IPCC 2014, 2021); many countries and regions may face declining groundwater levels such as the US, India, and Germany (Amanambu et al. 2020; Wunsch et al. 2022). Compared with porous aquifers, karst aquifers are potentially more sensitive to climate change, as their fast flow paths respond quickly to reduced recharge due to decreasing precipitation and increasing temperature (Hartmann et al. 2017). Being a vital freshwater resource, karst aquifers and springs play a vital role in fulfilling the drinking water needs of ~9.2% of the global population (Stevanović 2019). Understanding their changes is thus crucial for sustainably managing the usage while meeting the increasing demands of climate change and population growth (Mahler et al. 2021).

Many studies have attempted to project groundwater changes in the future climate by using groundwater models (e.g. MODFLOW, MIKE SHE) which take future meteorological variables either derived from global climate models (GCMs) and regional climate models (RCMs; Toews and Allen 2009; Kumar 2012; Doummar et al. 2018; Amanambu et al. 2020) or by applying a change factor to the historic climate (Chen et al. 2018; Dubois et al. 2020). However, the exact magnitudes of model projections are often highly uncertain due to (1) the influential drivers of groundwater, which are considered static in the model but are subjected to change, such as agricultural irrigation, land use and land cover changes, natural climate variations, and anthropogenic climate change (Green et al. 2011; Hartmann et al. 2014; Goderniaux et al. 2015); (2) the delayed response of groundwater to climate change (Panagopoulos and Lambrakis 2006; Hosseini et al. 2017); (3) the highly variable precipitation under climate change, and the increasing temperature that impacts not only evapotranspiration but also the timing of snowmelt (Chen et al. 2018; Lorenzi et al. 2022).

Quantifying the historic impact of climate change on groundwater contributes to understanding the causes of the change and increases confidence in the model projections. Analysing long-term hydrographs at natural catchments potentially allows for the assessment of this impact (Peterson and Western 2014; Ravbar et al. 2018). However, very few studies have quantified the historic impact due to the heterogeneous subsurface properties and flow processes (e.g. exchange between groundwater and surface water; Hartmann et al. 2014; Bonotto et al. 2022), multiple influential drivers of groundwater (e.g. pumping, land use changes, human activities; Shapoori et al. 2015; Kovačič et al. 2020), and lack of long-term groundwater observations (Taylor et al. 2013; Olarinoye et al. 2020). These issues are particularly outstanding for karst aquifers. Due to limited lengths of records and complicated hydrogeological characteristics, a few studies have investigated the water availability in the karst regions within a time frame of 20 years (Xanke and Liesch 2022; Lorenzi et al. 2022). Only a handful of studies have assessed the trends over a half-century (Kimmeier and Bouzelboudjen 2001; Fiorillo et al. 2015, 2021; Leone et al. 2021), and some of them struggled to attribute the trend of changes to climate change (Bonacci 2007). The long-term historic impact of climate change and variations on karst aquifers and springs is thus still largely unknown.

Climate change impacts the global water cycle and its impact is spatially variable because of the varying precipitation patterns and changing land surface-atmosphere feedbacks (IPCC 2021). Downscaling the global impact to the local scale potentially allows water managers and policymakers to understand the impact on their managing sites and regulate groundwater usage and licence the pumping limits accordingly. Besides that, new methodologies can be best developed at the local scale through analysis in great detail, then generalised and transferred to other sites. Complementing research at the global scale with detailed research at the local scale is thus necessary.

This study aims to establish and test an integrated approach to quantify the impact of climate variability and change on karst spring discharge from both historic and future perspectives. The proposed workflow of this approach is to (1) quantify the hydrodynamic characteristics of the karst system and the historic trends of the climate variables and discharge by using statistical approaches; (2) build up a lumped parameter model for the karst site and evaluate it with the knowledge learned from the historic analysis; and (3) use the model to simulate future discharge changes under climate change. Step 1 is widely lacking in the studies which attempted to project future water resource changes (Toews and Allen 2009; Kumar 2012; Amanambu et al. 2020), but it is crucial for understanding the hydrological behaviour of the system and the projected future changes. The Blautopf spring in southern Germany is studied, as it is one of the largest and most important karst springs in central Europe (Ufrecht et al. 2016). It drains a large catchment and belongs to Germany’s largest karst aquifer system (Swabian-Franconian Alb), which is essential to the drinking water supply of millions of people (Fahrmeier et al. 2022). Quantifying the response of the Blautopf system contributes to understanding the changes of karst aquifers and springs in the temperate climates where around 19% of the land surface consists of carbonate rocks (Goldscheider et al. 2020). The research questions below are answered in this article:

  1. 1.

    What are the historic climatic changes at the study site, in terms of temperature, evapotranspiration, precipitation, intensity and duration of snowmelt?

  2. 2.

    Did the observed climatic changes impact the karst discharge? If yes, what was the impact on the mean, low, and high flow, as well as the intensity and duration of hydrological extremes?

  3. 3.

    Assuming the IPCC climate change scenarios, how will the climate and discharge (minimum, mean, maximum) change at this type of climatic and hydrogeological setting (temperate, snow-influenced karst system) by the end of the twenty-first century?

Study site

The Blautopf spring is located at the southern margin of the Swabian-Franconian Jura Mountains, Germany’s largest karst landscape and karst aquifer system. The spring catchment covers an area of ~165 km2 delineated by numerous tracer tests and diving explorations (Armbruster et al. 2006; Ufrecht et al. 2016). The catchment drips gradually from northwest to southeast in a range of 878 to 609 m asl. (Fig. 1). The land surface is mainly covered by forests (~31%) and agriculture (~60%) (Köberle 2005; Lauber et al. 2014). Numerous dolines and dry valleys are visible on the land surface. Underneath the landscape is a spectacular and large cave system. The catchment is divided into two subzones by two major caves: the Blue Cave (“Blauhöhle” in German, ~14.6 km length) and the Hessenhau Cave (“Hessenhauhöhle” in German, ~6.0 km length; Schwarz et al. 2009). Each subcatchment is estimated to contribute approximately half of the spring discharge, which is the only known direct outlet of this catchment (Lauber et al. 2013, 2014). The vertical profile shows that the catchment consists of multiple layers of Upper Jurassic limestone and marl (Selg and Schwarz 2009). The thick saturated (50–120 m) and unsaturated zones (100–150 m) allow a considerable storage capacity of the karst system (Schwarz et al. 2009).

Fig. 1
figure 1

Location and map of the Blautopf spring and catchment. a Distribution of carbonate rocks in Germany and nearby countries (Chen et al. 2017), with the location of Blautopf indicated. b Map of the Blautopf catchment and spring and its nearby climate stations and river networks. c A view from the top of the Blautopf spring (photo: N. Goldscheider)

Full size image

Data and methods

Karst discharge data

The karst spring discharge (Q) recorded at station (ID 76273) from 1 Jan 1952 to 31 Dec 2021 was used in this study and obtained from the Baden-Württemberg State Office for Environment (LUBW). Records show that the discharge was estimated with the stage–discharge relationship from 1952 to 2004. Since 2005, the discharge was calculated as the discharge difference between the downstream and upstream rivers. The discharge calculated with both methods provides consistent results [Section S1 in the electronic supplementary material (ESM)]. The data were recorded in the hourly interval from 1952 to 2016 and the 15-min interval from 2016 onwards. After removing errors, the data show no daily gap. The subdaily data were then used to calculate the daily, monthly, seasonal, and annual mean discharge.

Climate data

The daily precipitation (P), mean air temperature (T), snow height (SH), and snow water equivalent (SWE) data of the nearby weather stations (Fig. 1) were used to estimate the catchment climate. The data were obtained from the Climate Data Center of Germany (Deutscher Wetterdienst 2021). Table 1 summarises the observed climate variables and the locations of the stations.

Table 1 The weather stations near the Blautopf catchment (locations shown in Fig. 1). Note, the elevation of Station 3402 was 721 m asl (48.4065° N, 9.4928° E) before 1989
Full size table

The observed precipitation at each station has a data missing rate of <1.5%. The temporal data gaps of each station were filled with the mean precipitation of the other stations on that day. The catchment precipitation was estimated with the Thiessen polygon method in which each station takes a weight of 50.3% (ID 5511), 37.6% (ID 537), 8.4% (ID 2814), and 3.7% (ID 3402). The observed daily mean air temperature at station 3402 with no temporal gap was used to represent the catchment temperature.

The Blautopf catchment was influenced by snow between October and April. Both SH and SWE data were recorded as accumulations of previous and new snow. The daily SH data had only one missing datum which was filled with the mean snow heights of the neighbouring days. Because of snow compaction, the SH data were only used to indicate the days when the catchment was covered by snow; the SWE data were used to represent the snow amount and were nearly in the weekly interval. To interpolate it to daily, the SWE data were first filled with 0 on the days when the SH was 0 and the remaining gaps were then filled by linear interpolation (note that small snow events, which resulted in 2–3 days of snow cover, between two zeros in the weekly SWE time series were not accounted). The degree-day factor (DDF) method (Martinec 1960; Çallı et al. 2022) was adopted to represent the snowmelt processes:

$${\mathrm{Snow}}_{i}= \left\{\begin{array}{c}{\mathrm{Snow}}_{i-1}+{P}_{i}, { T}_{i}\le {T}_{\mathrm{M}}\\ \mathrm{max}\left({\mathrm{Snow}}_{i-1}-{\mathrm{Melt}}_{i}, 0\right), {T}_{i}>{T}_{\mathrm{M}}\end{array}\right.$$
(1)
$${\mathrm{Lw}}_{i}=\left\{\begin{array}{c}0, {T}_{i}\le {T}_{\mathrm{M}} \\ {P}_{i}+ \mathrm{min}\left({\mathrm{Snow}}_{i-1}, {\mathrm{Melt}}_{i}\right), {T}_{i}>{T}_{\mathrm{M}}\end{array}\right.$$
(2)
$${\mathrm{Melt}}_{i}=\left\{\begin{array}{c}0, {T}_{i}\le {T}_{\mathrm{M}} \\ \mathrm{DDF}\times \left({T}_{i}-{T}_{\mathrm{M}}\right), {T}_{i}>{T}_{\mathrm{M}}\end{array}\right.$$
(3)

where Pi is the daily precipitation [mm] for day i, Ti is the daily mean temperature [°C], TM is the temperature threshold above which the snowmelt occurs, Snowi is the snow storage [mm day−1], Melti is the snowmelt amount [mm day−1], Lwi is the liquid water [mm day−1] that is available to infiltrate into the subsurface, and DDF is the degree-day factor (DDF) [mm °C−1 day−1]. The DDF (0.1–10) and TM (−1, 0, 1 °C) are the calibration parameters, and the best set of these was determined to be the one that gives the best estimates of snow storage compared with the observed SWE (Fig. S2 in the ESM).

The daily potential evapotranspiration (PET) of the catchment was estimated with a temperature-based method (Oudin et al. 2005):

$$\mathrm{PET}= \frac{{R}_{\mathrm{e}}}{\lambda \rho }\frac{{T}_{\mathrm{a}}+5}{100},\mathrm{ if}\,{ T}_{\mathrm{a}}+5>0$$
(4)
$$\mathrm{PET}=0, \mathrm{otherwise}$$
(5)

where PET is given in the unit of millimetre per day [mm day−1], Re is the extra-terrestrial radiation [MJ m−2 day−1] estimated with the FAO method by using catchment latitude and Julian day (Allen et al. 1998 and Section S2.2 in the ESM), λ is the latent heat flux (= 2.45 MJ kg−1), ρ is the water density [kg m−3], and Ta is the daily mean air temperature [°C].

Statistical analysis

The hydrodynamic characteristics of the karst system were analysed by using the correlations between the daily precipitation and discharge from 1952 to 2021. The autocorrelation and cross-correlation were calculated with the acf and ccf functions in the R stats package. The autocorrelation of discharge informs (1) the karstification degree of the system and the existence of fast and slow flow paths indicated by the slope of the function (Padilla and Pulido-Bosch 1995; Panagopoulos and Lambrakis 2006), and (2) the karst system memory which refers to the time that the water needs to flow through the karst aquifer and is determined as the time lag at which the function coefficient decreases to 0.2 (Mangin 1984; Benavente et al. 1985; Larocque et al. 1998). The cross-correlation function quantifies the response time of the discharge to precipitation at the time lag when the maximum coefficient is achieved (Panagopoulos and Lambrakis 2006).

The long-term trends of discharge and climate variables were quantified from monthly to annual scales by using the nonparametric Mann–Kendall test (Mann 1945; Kendall 1955) and Sen’s slope method (Sen 1968). A MATLAB function mannkendall v1.0.0 (Collaud Coen et al. 2020; Collaud Coen and Vogt 2021) was used to first remove the first-order autocorrelation in the data with a 3PW prewhitening procedure and then to quantify the trends. To identify any abrupt changes in the data, a Bayesian change-point detection and time-series decomposition algorithm called Rbeast (Zhao et al. 2019b) was applied. This algorithm uses a Bayesian model averaging scheme and the Markov Chain Monte Carlo stochastic sampling approach to approximate the complex relationships in the data (Zhao et al. 2019a). The results include the locations of the change points and their probability. In this study, the period was set to 1 in the algorithm for the annual time-series data.

Karst discharge modelling

Historical simulation

A reservoir model, KarstMod (Mazzilli et al. 2019), was used to simulate the time-series discharge under historic climate variations. This model was set up with three compartments (E, M, C; Fig. 2), whereby E represents the infiltration zone (soil and epikarst) which is recharged by precipitation and loses water into the atmosphere via evapotranspiration, and M and C represent the fissured rock matrix and conduits in the underground subsystems. The equations of the model are provided in the following (Mazzilli et al. 2019):

$$\frac{dE}{dt}=P-\mathrm{PET}-{Q}_{\mathrm{EM}}-{Q}_{\mathrm{EC}}, \mathrm{if}\, E\ge 0$$
(6)
$$\frac{dE}{dt}=P-\mathrm{PET}, \mathrm{if}\, {E}_{\mathrm{min}} < E<0$$
(7)
$$\frac{dM}{dt}={Q}_{\mathrm{EM}}-{Q}_{\mathrm{MS}}$$
(8)
$$\frac{dC}{dt}={Q}_{\mathrm{EC}}-{Q}_{\mathrm{CS}}$$
(9)
$${Q}_{\mathrm{EM}}={k}_{\mathrm{EM}}\times {E}_{t}, \mathrm{if}\, {E}_{t}>0,\mathrm{ otherwise}\, {Q}_{\mathrm{EM}}=0$$
(10)
$${Q}_{\mathrm{EC}}={k}_{\mathrm{EC}}\times {E}_{t},\mathrm{ if}\, {E}_{t}>0, \mathrm{otherwise}\, {Q}_{\mathrm{EC}}=0$$
(11)
$${Q}_{\mathrm{MS}}={k}_{\mathrm{MS}}\times {M}_{t}, \mathrm{if}\, {M}_{t}>0, \mathrm{otherwise}\, {Q}_{\mathrm{MS}}=0$$
(12)
$${Q}_{\mathrm{CS}}={k}_{\mathrm{CS}}\times {C}_{t}, \mathrm{if}\, {C}_{t}>0, \mathrm{otherwise}\, {Q}_{\mathrm{CS}}=0$$
(13)
$${Q}_{\mathrm{S}}={R}_{\mathrm{A}}\times \left({Q}_{\mathrm{MS}}+{Q}_{\mathrm{CS}}\right)$$
(14)

where Et, Mt, and Ct represent the water levels [mm] in the compartments E, M, and C at the time t; Emin [mm] represents the water holding capacity (of the soil and epikarst) over which the water flows from E into M and C; QAB represents the flow [mm day−1] from compartment A to B (e.g. E to M), and kAB represents the recession coefficient of the flow; QS represents the spring discharge [m3 s−1]; RA represents the recharge area [km2].

Fig. 2
figure 2

The model structure adopted to simulate the storage and flow dynamics of the Blautopf system

Full size image

The model took the daily transformed precipitation (with the DDF method) and PET as inputs to simulate the daily discharge from 1952 to 2021. Three phases were defined in the model: warm-up (1952–1956), calibration (1957–2007) and validation (2008–2021). The model was calibrated by using a quasi-Monte-Carlo approach with a Sobol sequence sampling of the parameter space (Sobol 1976; Mazzilli et al. 2019). Five parameters (kEM, kEC, kMS, kCS, and RA) were calibrated and 10,000 parameter sets were generated (parameter descriptions and values provided in Table S5 in the ESM). The parameter sensitivity was calculated with the Sobol index (Saltelli 2002) and the model performance was evaluated with the Nash-Sutcliffe efficiency (NSE) coefficient (Nash and Sutcliffe 1970):

$$\mathrm{NSE}=1- \frac{\sum {\left({Q}_{\mathrm{obs}}-{Q}_{\mathrm{sim}}\right)}^{2}}{\sum {\left({Q}_{\mathrm{obs}}-{Q}_{\mathrm{mean}}\right)}^{2}}$$
(15)

where Qobs is the daily observed spring discharge [m3 s−1], Qsim is the daily simulated discharge [m3 s−1], and Qmean is the average of the observed discharge [m3 s−1]. An NSE = 1 represents the best model performance and 0 means the model performance equivalent to using the observed mean discharge.

Simulation of future climate change impacts

The discharge was simulated with the calibrated model from 2006 to 2100 under three climate change scenarios in the 5th Phase of the Coupled Model Intercomparison Project (CMIP5): representative concentration pathways (RCPs) 2.6, 4.5, and 8.5 (most extreme scenario). Each RCP has an ensemble of climate projections generated with the GCMs which are coupled with the RCMs for estimating regional climate change. The GCM-RCM members that represent >80% of the climate spread of all model projections were adopted in this study (Deutscher Wetterdienst 2018). This criterion resulted in five GCM-RCM members being selected for RCP 2.6 and six members for RCPs 4.5 and 8.5 respectively (Table S2 in the ESM). The climate data of the adopted members were bias-corrected and downscaled to 5 km × 5 km by Brienen et al. (2020).

The daily precipitation and air temperature were extracted from each GCM-RCM member at the locations of the observed weather stations (ID 5511, 2814, 3402, 537 in Table 1). The climate data at each location was calculated as the mean of the surrounding 3 × 3 cells of the gridded climate outputs (Wunsch et al. 2022). The future catchment climates were estimated with the same methods as in section “Climate data”—that is, the catchment precipitation and temperature were estimated by spatially interpolating the extracted climate data using the Thiessen polygon method; the precipitation was transformed into the time-series liquid water with the DDF method; the PET was estimated with the method of Oudin et al. (2005). The daily transformed precipitation and PET of each GCM-RCM member were then taken by the model to simulate the discharge for each climate change scenario. The trends of the simulated annual minimum (2.5th percentile), maximum (97.5th percentile), and mean discharge and climates were quantified with the nonparametric Mann–Kendall test and Sen’s slope method (Mann 1945; Kendall 1955; Sen 1968) as in section “Statistical analysis”.

Results

Karst system hydrodynamics

The cross-correlation function of the daily precipitation and discharge (Fig. 3a) shows that the discharge responds to precipitation with a lag of 3 days. The quick response of the discharge to precipitation confirms the existence of fast flow paths in the karst system. The small peak coefficient (0.22) indicates that the precipitation signal has been strongly attenuated before flowing out to the spring. Following the maximum, no other peak is observed, indicating that one primary quick flow component exists in the system.

Fig. 3
figure 3

Cross-correlation and autocorrelation functions of the daily precipitation and discharge from 1952 to 2021. a Cross-correlation function. b Autocorrelation functions of the precipitation and discharge and the red dots on it indicate the lags when the gradient of the function changes

Full size image

The autocorrelation function of the daily precipitation decreases rapidly to 0 (Fig. 3b), indicating no autocorrelation in the data and that precipitation can be treated as a random variable. In contrast, the daily discharge shows strong autocorrelation and the system is estimated to have a memory of 40 days at the coefficient of 0.2. The changing gradients of the function between the time lags indicate the existence of both fast and slow flow components in the system. The coefficient decreases quickly from 1 to 0.5 in the first 11 days because of the fast flow and then reduces slowly, such as between the intervals of 11–16 days and 17–43 days because of the slow flow components. Overall, this karst system has highly karstified flow paths which allow the system to respond quickly to rainfall and snowmelt, and it also has a large storage capacity and fissured karst matrix as reflected by the system memory.

Discharge and climate trends

The annual mean air temperature and total PET at Blautopf show a significantly increasing trend between 1952 and 2021 with a slope of 0.03 and 0.95 °C mm year−1 (Fig. 4a). The annual total precipitation is highly variable and does not show a significant trend of changes (Fig. 4b). On average, around two-thirds of the annual total precipitation goes into the atmosphere via evapotranspiration and one-third of it recharges the aquifer (Table 2). The spring discharge shows a high variability in different time scales (monthly, seasonal, and annual). The annual mean discharge also does not show a significant trend of changes, with the maximum nearly four times the minimum.

Fig. 4
figure 4

The annual, monthly, and daily discharge and climates at Blautopf from 1952 to 2021. a–b The annual mean air temperature (T) and total potential evapotranspiration (PET); the annual total precipitation (P) and mean discharge (Q) with the average indicated (purple dashed line). c–e The monthly mean air temperature, total precipitation, and mean discharge. f–h The cumulative distribution functions (CDFs) of the daily mean air temperature, total precipitation, and mean discharge

Full size image
Table 2 Minimum, mean, and maximum of the annual total precipitation (P), annual total potential evapotranspiration (PET), the difference between the annual total P and PET, annual mean air temperature (T), and annual mean discharge (Q) from 1952 to 2021
Full size table

The monthly discharge and climate variables (Fig. 4c–e) show that June, July, and August are the warmest and also wettest months of the year. The mean discharge in these months is however lower than the annual mean rate (1.7 cf. 2.3 m3 s−1). August to October have the lowest discharge (1.2 m3 s−1 on average) of the year. March has the highest discharge but the lowest precipitation. It is also the month when the mean air temperature normally turns from negative to positive. Snowmelt thus plays a crucial role in generating the high discharge in March at the spring.

The daily mean air temperature at the catchment varies by nearly 50 °C between the minimum and maximum (from −23.0 to +26.7 °C; Fig. 4f). Nearly 20% of the days are below 0 °C and 40% of the days have no precipitation. The maximum daily precipitation is 81 mm and the highest 10% is more than 11 mm (Fig. 4g). The daily discharge varies largely between 0.1 and 32.6 m3 s−1, with an average of 2.3 m3 s−1 (Fig. 4h). The 10th and 90th percentiles of the daily discharge are 0.7 and 4.7 m3 s−1.

The annual, seasonal, and most of the monthly mean discharge and total precipitation show no significant trends of changes (Figs. S3–S5 in the ESM). The annual high (≥90th percentile) and low discharge (≤10th percentile) also show no significant annual trends of changes (Fig. S6 in the ESM). However, the discharge in January shows a significantly increasing trend, and that in April shows a significantly decreasing trend at the same rate of 0.02 m3 s−1 year−1 (Fig. 5). April is the last snow month of the catchment. The declining discharge in April is found because of the significantly decreasing precipitation, increasing temperature and evapotranspiration.

Fig. 5
figure 5

Trends of the monthly mean discharge (Q), total precipitation (P), total potential evapotranspiration (PET), and mean air temperature (T) between January and April from 1952 to 2021. The dashed lines represent the trends

Full size image

The impact of snowmelt on discharge

Timing

The occurrence of peak spring discharge highly relates to the peak snowmelt at the catchment. The changing timing and intensity of the peak snowmelt under increasing temperature and decreasing snow can hence impact the peak discharge. The peak snowmelt normally occurs between December and April at Blautopf (Fig. 6a). Results show that nearly two-thirds (n = 39) of the daily discharge maxima were observed during and after the peak snowmelt within 1 week. The peak snowmelt period is counted from the day when the peak snow accumulation (measured in SWE) is observed until it is entirely melted (e.g. Fig. 6b). The start day of the peak snowmelt shows no significant trend of changes (Fig. 6c), but the end day significantly shifts from March towards February (Fig. 6d). The duration of the peak snowmelt shows a significantly declining trend at a rate of −0.13 days year−1 (Fig. 6e). The peak snowmelt lasts 17 days on average in the first 20 years (1950–1970) but only 13 days in the last 20 years (2000–2020), showing a long-term average of 15 days. Overall, the peak snowmelt is shown to end earlier and last for a shorter duration. The peak discharge in the snow period is also found to onset earlier but this trend is not significant.

Fig. 6
figure 6

Timing of the peak snowmelt and the daily discharge maxima (Qmax) from 1952 to 2020. a The peak snowmelt period and the appearance of the daily discharge maxima. b An example of the peak snowmelt in 1953 (blue block). The y-axis is the daily observed snow accumulation (accu., measured in snow water equivalent). c–e The start day, the end day, and the duration of the peak snowmelt. The red solid line represents the trend. Note, a 1-year cycle is defined from July to June of the next year to cover a continuous snow period

Full size image

Intensity

The intensity of the peak snowmelt depends on the amount of the peak snow accumulation. Results show that around 3–23% of the annual total precipitation at Blautopf is snow, with a mean of 14% (Table 3). The peak snow accumulation is around half of the annual total snowfall on average. The highest peak snow accumulation at Blautopf was 181 mm recorded in 1988. This peak snow took only 14 days (close to the mean duration) to melt and generated the highest discharge rate (32.6 m3 s−1) at Blautopf on record. The peak snow accumulation shows a decreasing trend at all observed stations at a rate between −0.3 and −0.8 mm year−1 (Fig. S8 and Tables S3–S4 in the ESM). The most probable changing point was identified around 1988 and the peak snow accumulation has decreased by 35 mm on average since then (Fig. 7b). Around 85% of the peak snow accumulation was more than the daily precipitation maxima before 1988, but this rate has decreased to 75% after this year (Fig. S9 in the ESM). Despite that the daily precipitation maxima increasingly exceed the peak snow accumulation, they are mostly below the mean level (60 mm) of the latter. Overall, the peak snow accumulation at the catchment has shown a significant decreasing trend and the reducing magnitude has been especially severe since 1988.

Table 3 The annual peak snow accumulation (measured in snow water equivalent), total snowfall (measured in snow water equivalent), total precipitation, total number of days covered in snow, and the annual maximum discharge (Qmax) between 1952 and 2020. Note, a 1-year cycle is defined as spanning from July to June of the next year to cover a continuous snow period
Full size table
Fig. 7
figure 7

The daily peak discharge (Qmax) of the snow period (Oct–Apr) have shifted to a low state due to increasing air temperature (T) and decreasing peak snow accumulation (accu., measured in snow water equivalent) between 1952 and 2020. a The annual snow-period mean air temperature anomaly relative to its long-term average (2.3 °C). b The annual peak snow accumulation relative to its long-term average (60 mm). The interval of the legend is divided by the percentile (%ile) of the variables. c The relationship between the mean air temperature anomaly, the peak snow accumulation, and the daily peak discharge in the snow period

Full size image

The same changing point was identified in the annual snow-period (Oct–Apr) mean air temperature, which shows a significant increasing trend at a rate of +0.03 °C year−1 (Fig. 7a). The coldest snow period (–0.5 °C) was observed in 1962 and the warmest snow period (5.5 °C) was in 2006. Starting from 1988, 70% of the annual snow-period mean air temperature was above the long-term median and increased by 1 °C on average. These results show that the catchment has entered a warm and low-snow state since then (Fig. 7a,b).

With the decreasing peak snow accumulation and increasing air temperature, the peak discharge in the snow period tends to decrease, but the trend is statistically not significant. Nevertheless, it has shifted towards the low state since 1988 (Fig. 7c). Around 80% of the peak discharge rates are lower than the long-term mean (13.6 m3 s–1) and have decreased by 1.1 m3 s–1 on average since then. Before 1988, the peak snow accumulation was often above the long-term average (60 mm) and melted slowly (~17 days) in a cold temperature (<2.3 °C); after this year, the peak snow accumulation was mostly below the average and melted faster (~13 days) in a warmer temperature (>2.3 °C). The multilinear regression between the peak snow accumulation, the snow-period mean air temperature, and the peak discharge (p-value < 0.01) shows that a 1-mm change in the peak snow accumulation results in 0.003 m3 s–1 change in the peak discharge, and 1 °C change in the air temperature results in a –1.8 m3 s–1 change in the peak discharge. The low R2 (0.2), however, indicates that the multilinear regression model explains little variance in the variables.

Karst discharge simulation in the past and future

The results of model calibration show that all 10,000 parameter sets give an NSE ≥ 0.70 (Fig. S10 in the ESM). An optimum was found after exploring the whole parameter space, and the NSE of the optimum calibration is 0.74 and that of validation is 0.79 (Fig. 8). The NSE value indicates that 74% of the discharge changes at Blautopf have been explained by climate variations alone. The model has well captured the majority of the flow dynamics despite occasionally underestimating the extremely high discharges (Fig. 8 and Fig. S11 in the ESM). Overall, the model shows to be skillful in simulating the long-term discharge changes and the trend with only climate inputs, but caution is needed when using it to estimate the hydrologic high extremes.

Fig. 8
figure 8

The observed and simulated daily discharge under climate variations a in the warm-up (1952–1956) and calibration phases (1957–2007, NSE = 0.74), and b in the validation phase (2008–2021, NSE = 0.79) with simulation uncertainty, which is quantified as the differences between the maximum and minimum simulated discharges of all behavioural parameter sets (n = 10,000, NSE ≥ 0.70). Abbreviations: observed discharge (Qobs), simulated discharge (Qsim), precipitation (P), potential evapotranspiration (PET)

Full size image

The future projections simulated with the model show a substantially different warming degree between the climate change scenarios at Blautopf (Fig. 9): the annual mean air temperature is projected to increase by 0.5, 1.8, and 4.5 °C by the end of the twenty-first century for the RCPs 2.6, 4.5, and 8.5 respectively. The annual total precipitation is highly variable in the whole time series, showing no significant trend of changes for the RCPs 2.6 and 8.5, and a slightly increasing trend for the RCP 4.5 with a slope of 0.6 mm year−1 (p-value: 0.04). The long-term average of the annual precipitation between 2000 and 2100 is shown to be slightly higher (20–50 mm) than the observed annual mean (950 mm) between 1952 and 2021. The annual mean discharge also is shown to be highly variable under climate change, similar to the precipitation. In an extremely warming scenario of RCP 8.5, the annual mean discharge is projected to decrease at a rate of −0.003 m3 s−1 year−1 (p-value: 0.06). The confidence interval shows the discrepancies between the projections of the adopted GCM-RCM members. The uncertainty due to climate models is ~10% for the projected annual mean air temperature, ~14–18% for the annual total precipitation, and ~35% for the annual mean discharge.

Fig. 9
figure 9

The annual a mean air temperature (T), b total precipitation (P), and c mean discharge (Q) simulated under three climate change scenarios (RCPs 2.6, 4.5, and 8.5) between 2006 and 2100. The coloured lines represent the median of the simulation results of all GCM-RCM members for each RCP; the coloured blocks represent the 75% confidence interval (CI) of all simulations; the black lines represent the observations

Full size image

The annual minimum, mean, and maximum temperatures all increase with the climate change scenarios (Fig. 10a–c): the annual minimum temperature is projected to increase between 0.3 and 3.4 °C, the annual mean temperature may increase between 0.4 and 3.1 °C, and the annual maximum temperature may increase between 0.3 and 3.7 °C for the three RCPs by the end of the twenty-first century. For the changes of precipitation by 2100 (Fig. 10d–f), the annual minimum precipitation (Pmin) slightly decreases with the intensifying climate change scenarios: Pmin (RCP2.6) > Pmin (RCP4.5) > Pmin (RCP8.5). The annual total precipitation (indicated by the median) is projected to decrease by around 2 and 4% (−22 mm and −36 mm) for RCPs 2.6 and 8.5, but to increase by 3% (+26 mm) for RCP 4.5. The annual maximum precipitation may increase by 4 and 5% (+0.7 and +0.9 mm) for RCPs 4.5 and 8.5, but shows no trend for RCP 2.6.

Fig. 10
figure 10

The annual minimum (2.5th percentile), mean (or sum), and maximum (97.5th percentile) a–c air temperature (T), d–f precipitation (P), and g–i discharge (Q) at Blautopf simulated for three climate change scenarios (RCPs 2.6, 4.5, 8.5) in the beginning (B, 2006–2039), middle (M, 2040–2069), and end (E, 2070–2100) stages of the twenty-first century

Full size image

The annual minimum, mean, and maximum discharge all decrease with the intensifying climate change scenarios (Fig. 10g–i). The annual minimum discharge (indicated by the median) is projected to decrease by the end of the century by 43% (−0.21 m3 s−1) for RCP 2.6, 35% (−0.17 m3 s−1) for RCP 4.5, and 55% (−0.27 m3 s−1) for RCP 8.5, compared with that (0.49 m3 s−1) between 1952 and 2021. The annual mean discharge is projected to decrease by 4% (−0.10 m3 s−1) to 10% (−0.22 m3 s−1) for RCP 2.6 to RCP 8.5 by 2100. The annual maximum discharge is shown to change variably in each stage of the century, with a long-term average of ~6.7–7.0 m3 s−1 for the three RCPs. Compared with the long-term mean (13.6 m3 s−1) between 1952 and 2021, the annual maximum discharge may decrease by ~50% projected by all three RCPs.

Discussion

Karst system structure

The Blautopf karst system responds quickly to precipitation in 3 days and also has a memory of 40 days. The quick response of discharge reveals that part of the karst system is highly karstified and forms fast flow paths, which are also supported by direct observations in the active conduit network and in-cave tracer tests (Lauber et al. 2014). Besides the fast flow, some rainfall infiltrates into the small fissures in the rocks and is released slowly to the spring (Larocque et al. 1998; Panagopoulos and Lambrakis 2006). The rainfall signal is found to have significantly attenuated before the water discharges, as the newly recharged water has well mixed with the previous storage in the deep unsaturated zone (100–150 m), the large and deep phreatic conduits, and the underground lakes in the karst system. The hydrodynamic characteristics of the karst system summarised from correlation analysis (of discharge and precipitation) are also supported by the previous isotope study which traced the water flow from the surface to the discharge at this site (Schwarz et al. 2009).

Reliability of the groundwater model

In examining the adequacy of the groundwater model, the long-term observed daily discharge changes have been well simulated with only climate inputs in both calibration and validation phases (Fig. 8). The model proves to be skillful in simulating the discharge in both low and mean flow but also occasionally underestimates the extremely high discharges (Fig. S11 in the ESM), which is a commonly known issue in hydrological modelling (Cinkus et al. 2023). Future simulations of such events may thus contain higher uncertainty. Additionally, the hydrodynamic characteristics of the karst system have been well simulated by the model, which is close to the historic analysis of the observed discharge and climate data (Figs. S12 and S13 in the ESM). Specifically, the modelled discharge reacts to the rainfall and snowmelt events the same as the observed discharge with a lag of 2–3 days; however, the model has a system memory which is ~20 days longer than the observed site. This indicates that the model may need more time to return to the initial state after experiencing any disturbances; additionally, the calibrated model parameters (i.e. recession coefficients of the flows and recharge area) could potentially be useful for understanding the site and system properties. The calibrated catchment recharge area (164.70 km2) is in good consistency with the results (165 km2) from tracer tests and diving campaigns (Ufrecht et al. 2016), thus increasing confidence in the model.

The model however also shows a few caveats that may impact the future projections. For example, the model parameters show relatively high interactions with each other (Sobol sensitivity indexes in Table S6 in the ESM), theoretically indicating a lower identifiability of the parameters, which could be because of the simple and symmetric structure of the model (Fig. 2). The selection of this model structure is because (1) this structure is consistent with the previous literature (Schwarz et al. 2009) and expert knowledge of the system; (2) this structure gives the best model performance; (3) all calibrated parameters are sensitive and show an optimum (Fig. S10 in the ESM), and (4) the parameters that show relatively high interactions (i.e. kEM, kMS) could be potentially linked and interdependent by nature. For these reasons, this model structure is tentatively accepted but may contribute to uncertainty in future projections. Regarding the prediction uncertainty, some other factors can also play a role: (1) the parameters in the groundwater model and snow model are assumed as constant; however, changes in snow model parameters (e.g. DDF factor) may impact the partitioning of precipitation to snow and rainfall and thus influence the rate and timing of the peak discharge (Doummar et al. 2018); (2) the stationary assumptions of the system evolution (karstification) and the land use and land cover, which may change in the future and impact the hydrological processes; (3) and the model may struggle to simulate the extremely high and low flows because of lack of training for simulating such unobserved and challenging conditions.

Trends and mechanisms

The annual mean and minimum discharge of Blautopf do not show significant annual and seasonal trends of changes from 1952 to 2021, but rather they are projected to decrease by the end of the twenty-first century. In the scenarios of RCPs 2.6, 4.5 and 8.5, the annual mean discharge is projected to decrease by 4–10% and the annual minimum discharge may decrease by 35–55%. In contrast, the intensity of the annual high discharge has already been impacted by the warming climate and shifted to a low state (<13.6 m3 s−1) since the 1980s, and may decrease by 50% projected by all three scenarios in this century.

The changes of discharge at Blautopf could be attributed to its climatic and hydrogeologic settings. With the increasing temperature under climate change, the catchment evapotranspiration can be enhanced due to ~90% of the land surface being covered by forests and agriculture (Köberle 2005). As the annual total precipitation is projected to change slightly (−36 to +26 mm) at the catchment by the end of the century, the groundwater recharge may thus reduce due to increasing evapotranspiration (+60 to +147 mm) and the discharge may consequently decline. Besides that, the increasing temperature will not only increase the evapotranspiration but also impact the timing and intensity of snowmelt. Despite that Blautopf is located at a relatively low altitude (600–800 m asl) and experiences less snow (~14% of the precipitation) compared with the alpine regions, the peak snowmelt still plays a crucial role in generating the annual peak discharge at the spring. In the past decades, the peak snowmelt has been shown to last for a shorter duration (−4 days) and has become less intense (−35 mm in peak snow accumulation). As a result, the peak discharge has shifted towards a low state. Despite that the daily precipitation maxima increasingly exceed the peak snow accumulation since 1988, they are mostly lower than the long-term average (60 mm) of the peak snow accumulation. With the increasing annual minimum temperature in winter (+0.3 to 3.4 °C by 2100), the peak discharge may continue decreasing by nearly half of the rate due to the less intense peak snowmelt under global warming.

This study has focused on quantifying the long-term response of the karst discharge to climate variations. To understand the detailed changes of the karst discharge in the dry and wet periods, future research may identify all likely changing points in the time-series discharge and climate data and analyse the step changes (Bonacci 2007). Additionally, our model is constrained to simulate the extremely high discharge. The future simulations of such events may hence contain higher uncertainty. Despite this caveat, the model still proves to be useful in estimating the system behaviour and the trends of discharge changes under climate change.

Implications

This study demonstrated the importance of quantifying the response of karst discharge to long-term climate change and variations. Catchments like Blautopf located in the temperate climate and mid-latitude may have been thought of as safe under climate change, compared with glaciers and arid zones (IPCC 2021). This study shows that despite the not significant historic trends of changes in the annual mean and minimum discharge, the annual peak discharge at the spring has reduced due to weakened peak snowmelt. If climate change continues, all types of discharge at the spring are projected to experience long-term adverse impact by 2100. The baseflow plays a crucial role in sustaining spring discharge, especially in the dry period. The decreasing minimum discharge may indicate water scarcity risks at the spring and pose threats to the groundwater-dependent ecosystems. The decreasing high discharge may warn of potential water supply risks for local communities and infrastructures, such as in southern Germany where karst aquifers and springs are taken as one of the main drinking water supplies (Fahrmeier et al. 2022).

Quantifying the local impact of climate change also allows for comparing the responses of karst discharge at various climatic and geologic settings worldwide and understanding the variable impact of climate change. Due to limited lengths of karst discharge data and multiple influential factors (e.g. pumping), very few studies have quantified the long-term historic climate impact on karst discharge (Kimmeier and Bouzelboudjen 2001; Fiorillo et al. 2015, 2021). The karst spring discharge in central and southern Italy is found to have declined by 15–30% since 1987 (Fiorillo et al. 2015). Similarly, the Blautopf spring has experienced significant changes since 1988, but its annual mean discharge decreased by 4% (−0.1 m3 s−1) and the daily discharge maxima decreased by 8% (−1.1 m3 s−1). The annual mean discharge at the springs in Italy shows a significantly decreasing trend since the 1920s due to the drying climate (Leone et al. 2021), whereas this trend is not significant for Blautopf. These findings imply that southern Germany has experienced the adverse historic impact of climate change, although relatively slighter than the hotspots (e.g. Italy) in the Mediterranean.

When looking at the future projections in the twenty-first century, the annual low, mean, and high discharge at the karst springs in the Mediterranean are all projected to decrease, also the same for Blautopf, but their decreasing magnitudes are different. For example, the annual mean discharge is projected to decrease by 73% at a snow-governed semiarid site and by 36% at a site slightly impacted by snow in Lebanon (Doummar et al. 2018; Dubois et al. 2020), whereas at the karst springs in southern Greece and the West Bank in eastern Mediterranean, where they are almost not impacted by snow at all, the decrease may be 15–30% (Hartmann et al. 2012; Nerantzaki and Nikolaidis 2020). For Blautopf (southern Germany, temperate climate, snow-influenced), the annual mean discharge may decrease by 4–10% in this century. Despite the uncertainty, these findings suggest that the karst aquifers and springs in the Mediterranean are likely to experience severe depletion in the future climate scenario, while the sites more impacted by snow and located in more arid climates may experience larger reductions. These findings reinforce the need for maintaining and extending the groundwater monitoring networks for detection and management of groundwater changes under climate change.

Conclusions

In this study, a three-step integrated methodology has been established to quantify the long-term impact of climate variability and change on karst discharge from both historic and future perspectives, including (1) statistical analysis of the hydrodynamic characteristics of the karst system and the observed climates and discharge, (2) model calibration and evaluation with historic data analysis, and (3) simulations of karst discharge under three climate change scenarios. This methodology has been applied to quantify the response of karst discharge at a snow-influenced temperate catchment (Blautopf) in central Europe from 1952 to 2100. The key findings of this study are:

  • The Blautopf karst system, with a memory of 40 days, consists of both fast and slow flow paths that enable the discharge to respond quickly to precipitation in 3 days. The fast and slow flow components characterize the unique response of the karst aquifer to climate change and variations.

  • The annual and seasonal mean discharge do not show a significant historic trend of changes. The monthly mean discharge shows that the discharge in April (the last month of the snow period) has significantly reduced (−0.02 m3 s−1 year−1) due to decreasing precipitation and increasing air temperature and evapotranspiration. The annual mean discharge is projected to decrease by 4–10% (−0.10 to −0.22 m3 s−1) for the RCPs 2.6, 4.5 and 8.5 by 2100.

  • The peak snowmelt plays a key role in generating the discharge maxima of the spring, but it is shown that it will end earlier (shift from March to February), last shorter (−4 days), and become less intense (−35 mm in peak snow accumulation). As a result, the annual peak discharge has shifted towards a low state (<13.6 m3 s−1) under global warming since 1988. The peak discharge may continue decreasing by 50% (−7 m3 s−1) as projected by all scenarios.

  • The annual minimum discharge does not indicate a significant historic trend of changes; however, it is projected to decline by 35–55% (−0.17 to −0.27 m3 s−1) for RCPs 2.6, 4.5 and 8.5 due to increasing evapotranspiration.

Overall, this study shows the long-term historic and future impacts of climate change and variations for a snow-influenced temperate catchment in central Europe, and may suggest potential water scarcity risks at similar climatic and geologic settings (temperate climate, snow-influenced, large karst drainage area) worldwide. Despite that the impact is relatively milder compared with the Mediterranean region, quantifying and understanding the changes of karst aquifers and springs in this type of setting is necessary for mitigating the future impact of climate change.