Introduction

Globally, groundwater plays a pivotal role in meeting the water needs of nearly 2 billion people and contributes over 40% to irrigation water consumption (Kundzewicz and Döll 2009; Famiglietti 2014; Rodell et al. 2018). Despite its significance, groundwater depletion has become a recognized global phenomenon (Jasechko and Perrone, 2021). The impact of climate change on groundwater resources, influenced by both natural and anthropogenic processes, and the reciprocal influence of groundwater on the climate system have been subject to examination (Taylor et al. 2013, Kløve et al. 2014). As climate change is anticipated to amplify the frequency of droughts, groundwater emerges as a dependable backup water source to fulfill diverse sectoral demands during periods of water scarcity (Langridge and Daniels 2017).

Drought, a complex climate-related occurrence with multifaceted economic, social, and environmental repercussions, has been the subject of past investigations. Spinoni et al. (2015) identified the twenty-two most significant drought events in Europe between 1950 and 2012, while subsequent years, such as 2003, 2015, 2018, and 2019, witnessed major droughts in Central Europe (Boergens et al. 2020). The summer of 2022 in Europe, evaluated across four geographical regions, experienced below-average soil moisture, precipitation, and relative humidity (C3S 2022a).

The climate projection-based drought trend for Central Europe shows the smallest changes, especially in terms of drought severity, under Regional Concentration Pathway (RCP) 4.5, while the increase in frequency is greater under RCP8.5 (Spinoni et al. 2018). In the case of Central and Eastern Europe, which is characterized by a moderate and uniform water supply, all scenarios show negligible changes. The expected monthly distribution of available water for the wet seasons (months 8–10) shows that the wet season is becoming wetter. A higher increase in evaporation is observed, especially in the RCP 8.5 scenario. These changes are significant in the wet season but not in the dry season. Higher emission scenarios show a more pronounced increase in evaporation (Konapala et al. 2020).

However, this cannot be dealt with as simply as a "wet gets wetter" scenario. Catchy phrases like “wet gets wetter, dry gets drier” or “global desiccation” should be avoided, and it is more useful to refer directly to the processes and variables that underlie these frameworks (Zaitchik et al. 2023). Future changes in groundwater storage (GWS) are determined not only by projected changes in precipitation, but also by other hydrological processes (such as evaporation and snowmelt) (Wu et al. 2020). In Central–Eastern Europe, the Carpathian chain plays an important role in the climate of the region, with variations in each parameter (Cheval et al. 2014). In addition, the Carpathian Basin is a significant regional water supply system (Taylor et al. 2013).

Global and regional scale studies (Konapala et al. 2020, Spinoni et al. 2018; Pokhrel et al. 2021; Vicente-Serrano et al. 2023) shows roughly concordant small or moderate changing climate estimations for the future of Central–Eastern Europe. Local scale studies for the Carpathian Basin (Mezősi et al. 2014; Fehér and Rakonczai 2019; Garamhegyi et al. 2018; 2020; Kovács and Jakab 2021) highlight the sensitivity of landscapes based on climate change induced shallow groundwater (SGW) fluctuation. Based on the Gravity Recovery and Climate Experiment (GRACE) and the GRACE Follow‐On (GRACE-FO) satellite data in Central Europe, a mass deficit of terrestrial water storage (TWS) was observed (Boergens et al. 2020). With using of GRACE and GRACE-FO data, the TWS or GWS changes due to natural and human drivers can be assessed (Thomas and Famiglietti 2019, Boergens et al. 2020; Wu et al. 2020, Zhang et al. 2022). If observed time series of the physical parameters required for modeling are available, then application of hydrological models, like Soil and Water Assessment Tool (SWAT), can help with evaluation (Zhang et al. 2022). These methods are capable of showing the observed changes but insufficient for analysis of the future changes. Given the non-stationary nature of hydrological statistics, relying solely on statistical parameters derived from historical data is insufficient (Milley et al. 2008, Szöllősi-Nagy 2022). There are numerous examples in the international literature about how to investigate the directly influence of climate change on groundwater resources. Two main approaches outlined, first one, where for the estimation of climate change effects on SGW is based on machine learning, like long short-term memory (LSTM), convolutional neural networks (CNNs), nonlinear autoregressive networks with exogenous input (NARX), and artificial neural networks (ANNs) (Jeong and Park 2019; Guzman et al. 2019; Fang et al. 2020; Wunsch et al. 2021, 2022, Sfîcă et al. 2022), and second one, where SGW response is calculated based on physical parameters with a hydrological model, like SWAT, MODFLOW, and HYDRUS-1D (Hosseinizadeh et al. 2019, Wang et al. 2021, Fallahi et al. 2023, Riedel et al. 2023, Sheikha-BagemGhaleh et al. 2023), also, as a third approach, instead of the SGW responses drought indices are calculated (El Asri et al. 2019, Afzal and Ragab 2020); however, the use of the downscaled general circulation models (GCMs) and regional climate models (RCMs) as input data is the same in all approaches.

The former study has proven that in the twenty-first century arable lands on the Great Hungarian Plain, which is a part of the Carpathian Basin, will face increasing environmental hazards owing to the unfavorable trends of the climate change on the basis of regional climate model simulations, and the most important change is the intensive increase of the drought hazard, which makes this hazard probably the most serious hazard of the region (Mezősi et al. 2014). Contradictions between global and local scale forecasts must be resolved. It can be assumed that, for the decision makers and the stakeholders, the exact quantified changes of SGW are more comprehensible than drought indices. The present study aims to forecast the changes in SGW levels (GWLs)—this is vital for the forthcoming decades—and the foundation of risk management strategy. The study examines the outcome of a "no intervention" water management scenario for GWL by 2050 and 2100. A hydrological process-based approach was chosen for the tests. The model's climatic inputs are scenarios 4.5 and RCP 8.5 from the bias-corrected RCM dataset. Only direct meteorological inputs are considered, while highly uncertain anthropogenic factors such as groundwater extraction or irrigation are excluded. For the analysis, the use of a robust and fast hydrological model is favorable, especially with a simply but good description of the physical processes. This paper follows "Occam’s razor" philosophy.

Materials and methods

Study area description

The area featured in this study is in the Great Hungarian Plain (GHP), Hungary, in the middle of the Carpathian Basin of East–Central Europe (19.38° E to 22.86° E and 46.18°N to 48.32° N) and has an extent of 40,473 km2 (Fig. 1). The GHP is a part of the Tisza River Basin. The Tisza River stands as Hungary's second-largest watercourse, boasting a comprehensive drainage basin encompassing an area of 157,186 km2 (ICPDR 2011). This watershed traverses five countries, namely Ukraine, Romania, Slovakia, Hungary, and Serbia (ICPDR 2011). Of this expanse, the Hungarian drainage area encompasses around 46,213 km2, a substantial portion of which lies within the GHP, a pivotal domain of Hungarian agriculture (ICPDR 2011; Pinke et al. 2020). Its significance is underscored by the extensive length of flood protection main dikes, spanning 4,157.7 km, with the Tisza system accounting for 2,942.3 km of this (GDOWM 2021). On a European scale, this dimension of flood protection is notable and roughly akin to the flood protection system along the Po River valley in Italy (ICOLD 2018). The elevation of the Hungarian part of the basin varies between 74.5 and 1014 m above Baltic Sea level (MASL) (Fig. 1). Our investigatory focus is directed at the GHP portion of the Tisza's watershed. The geographical positioning of this area is elucidated in Fig. 1.

Fig. 1
figure 1

The location of the study area in Central–Eastern Europe, Hungary. The white diamonds represent the spatial distribution of the SGW monitoring wells used in the study. The magenta diamond represents the well No. 2624 which also was used for modeling and will be presented later as an example. The HYDRODEM model of General Directorate of Water Management (GDOWM) was used as digital elevation model (DEM). Base map was constructed using the OpenStreetMap open database (OSMF 2023)

Full size image

Alongside flood protection embankments, the most significant interventions into the natural course of the river involved the construction of two dams as part of the “TIKEVIR” water management system (namely the hydropower plants at Tiszalök and Kisköre, see Fig. 2a), aimed primarily at securing irrigation water supply for the GHP (Koncsos 2011). However, the substantial flood waves of the early 2000s directed attention to the depletion of reserves within the flood protection system. The former system, predominantly relying on dikes, was augmented during the Further Development of the Vásárhelyi Plan (FDVP) through the incorporation of flood protection reservoirs (Fig. 2a) (Koncsos 2011). The impact of river regulation also triggered shifts in land use, with a significant portion of previously forested areas and aquatic habitats being converted into arable land (see Fig. 2a) (Demeter et al. 2022).

Fig. 2
figure 2

a Location of the water management systems and operational FDVP reservoirs along the Tisza River and the current land use and land cover (LULC) patterns based on MA 2019. 6 out of 7 FDVP reservoirs in operation, the 7th reservoir is under construction. Their summarized capacity will be 763 million m3 (GDOWM 2021). b The temporal distribution of available SGW monitoring data shows a high variability (Source GDOWM) as well as the linear trends of GWL change, c Base map was constructed using the OpenStreetMap open database (OSMF 2023)

Full size image

To preliminary analyze the groundwater resources of the GHP, the GDOWM’s database of monitoring wells was used. The time range of the provided database is 1961–2013.This database included observational time series alongside information on well positions, elevations. After error correction (typically arising from variations in well rim elevation and insufficient documentation), as well as data quality control, the discernible nature of changes was examined.

More than 600 groundwater monitoring stations are distributed across the GHP, in the study area 557 of them was considered. Most of these stations has been collecting data typically daily since the 1960s, while others were initiated between the 1970s and 2010s (Fig. 2b). Linear trend analysis (Altman and Krzywinski 2015) was conducted to uncover trends related to different observation periods (Fig. 2c). The trends and the length of the time series also vary, but no connection can be shown between the age of the well and the decreasing or increasing nature of the trend (Table 1, Supplementary Fig. S1). It is not possible to rely on trend analysis in terms of future GWL changes.

Table 1 Descriptive statistics of the GDOWM’s groundwater database for selected 557 wells in the study area
Full size table

Although it is well known that the GHP is a drought-prone region (Rakonczai 2013; Mezősi et al. 2014; Garamhegyi et al. 2018, 2020; Fehér and Rakonczai 2019; Pinke et al. 2020), moreover, the central part of the plain displaying the greatest susceptibility to prolonged and severe drought (Weidinger and Mészáros 1997; Horváth and Breuer 2023; Erdődiné Molnár and Kovács 2023), the drought of 2022 made this obvious to the general public, decision-makers, and stakeholders. The Combined Drought Indicator (CDI) of the European Drought Observatory (EDO) (Cammalleri et al. 2021) reveals that the first 10 days were exceptionally dry of July 2022. Nearly the entire GHP experienced either 'Alert' or 'Warning' conditions (EDO 2023). Periods of meteorological drought typically coincide with low streamflow, necessitating water authorities to prepare for mitigating drought risks and optimizing the utilization of existing water resources (Szamosvári 2022).

Thus, the study area was chosen because it is the most productive agricultural region of Hungary and has an important in food production (Pinke et al. 2020). The estimation of climate change induced responses of the GWLs is important, especially in a case of “no intervention” water management scenario. One of the most indicative markers of changes in groundwater resources is the variation in the groundwater table. Alongside robust measurement techniques, long-term, high-frequency datasets are available.

Methodology

As shown in Fig. 3, this research was carried out in four main phases. In the initial phase, emphasis was placed on the acquisition of relevant databases. Preliminary investigations were carried out on the GWL database maintained by the GDOWM, as outlined in the description of the study area. Subsequently, within the domain of meteorological databases, a spatial consistency test of observed time series was executed. This scrutiny resulted in the identification of the FORESEE HUN v1.0 database as a reliable foundation for subsequent investigations.

Fig. 3
figure 3

Flowchart of the present study

Full size image

Transitioning to the second phase, an aggregated vertical hydrological model was systematically developed. The primary objective of this model was to provide a robust and expeditious representation of changes in groundwater levels. Moving to the third phase, the model underwent a rigorous calibration and validation process, ensuring alignment with monitored GWLs.

In the fourth phase, attention was directed toward the anticipated future behavior of GWLs in the individual monitoring wells, contingent upon diverse climate scenarios. For the multi-scenario multi-model analysis, comprehensive trend analyses were conducted for both the RCP 4.5 and RCP 8.5 scenarios. These analyses provided estimates for changes in GWLs extending up to the years 2050 and 2100.

Databases

Meteorological databases

First, the expected meteorological changes in the examined area must be taken into consideration. To examine the changes and their consequences, it is necessary to calibrate and validate the descriptive model using historical data. In the research four meteorological database was used.

For the period between 1961 and 2010, the CARPATCLIM Database © European Commission—JRC, 2013 was utilized. The climatological grids cover the latitudinal range between 44°N and 50°N, and the longitudinal range between 17°E and 27°E. The spatial resolution is 0.1° × 0.1° (Szalai et al. 2013). The CARPATCLIM database is available free of charge at http://www.carpatclim-eu.org/pages/download/.

For the period between 1971 and 2019, the homogenized, gridded climate data series from the Hungarian Meteorological Service (HMS) were employed. The spatial resolution of this database is 0.1° (HMS 2021). In accordance with European Union directives, the resources provided by the Hungarian Government enabled the HMS to implement an open meteorological data policy from January 1, 2021. The database is available through https://odp.met.hu/ website.

The potential future weather trends can be quantified using regional climate models (RCMs); however, these model outcomes are inevitably accompanied by uncertainties. To estimate the anticipated hydrological changes, corrected climate model results are required. These models are based on historical observations, future projections, and consistent correction methodologies. Such datasets are available in the Open Database for Climate Change-Related Impact Studies in Central Europe—FORESEE database, first released in 2015. This database provides interpolated daily minimum and maximum temperatures, daily precipitation, and daytime average vapor pressure deficit (Dobor et al. 2015; Kern et al. 2019). The FORESEE is database is constantly being developed. FORESEE v4.0 and FORESEE HUN v1.0 released in 2022 (FORESEE 2023). The database is available after an approved access request at https://nimbus.elte.hu/FORESEE/index-data.html.

The FORESEE v4.0 offers interpolated daily meteorological fields at a spatial resolution of 0.1° × 0.1°. This interpolation covers the period from 1951 to 2020 based on observations, and for the subsequent period from 2021 to 2100, it includes 28 bias- and discontinuity-corrected projections. These projections are based on the results of the regional climate model (RCM) within the EURO-CORDEX experiment, utilizing 14 models and 2 scenarios (RCP4.5 and RCP8.5) (FORESEE 2023).

On the other hand, the FORESEE-HUN v1.0 provides interpolated daily meteorological data at the same spatial resolution (0.1° × 0.1°) but focuses specifically on Hungary. For the period from 1971 to 2021, it relies on the HUCLIM gridded dataset of the HMS, and for the subsequent period from 2022 to 2100, it presents 28 bias- and discontinuity-corrected projections. Like FORESEE v4.0, these projections are based on 14 models and 2 scenarios (RCP4.5 and RCP8.5) (FORESEE 2023).

FORESEE v4.0 consists of 87 NetCDF files in a total database of approximately 400 GB. FORESEE-HUN v1.0 also contains 87 NetCDF files, resulting in a total database size of approximately 40 GB. The calendar in both databases comprises 365 days. In leap years, 31 December is excluded. In this context, 29 February is included in the dataset (FORESEE 2023). In leap years the data for 31 December in leap years, the temperature and relative humidity data from 30 December was duplicated and zero precipitation was assumed.

The utilized meteorological data encompass daily minimum temperature (Tn) (°C), daily maximum temperature (Tx) (°C), daily precipitation (P) (mm), and daily relative humidity (RH) (−). Especially, FORESEE uses daytime average vapor pressure deficit (VPD) (Pa) instead of RH, the used conversion between VPD and RH (Murray 1967; Allen et al. 1998; Monteith and Unsworth 2013) is descripted in the Supplementary material. In calculations, the daily average temperature was employed, derived by averaging the daily minimum and maximum temperature values.

During the modeling process, an effort is made to reduce both time and storage requirements. Consequently, the entire database is not assigned to the model; only the necessary data are linked. Even with 557 monitoring wells, this still constitutes a substantial volume of data. Therefore, a query is conducted only for 76 randomly selected monitoring wells that effectively cover the study area. For wells without queried meteorological time series, data series are assigned using the Thiessen’s polygon method (Thiessen 1911). The applicability of this method has been highlighted in previous studies (Sfîcă et al. 2022), moreover, the methods suitability was tested in this study (Supplementary Figs. S2–S19). The data's spatial processing was performed using the QGIS 3.10 'A Coruña' software (http://www.qgis.org). Grids were established for each of the databases. The individual points' correspondence to grid cells was identified for querying the databases. The query was performed in the Matlab (The MathWorks Inc. 2023) software environment.

Database of groundwater monitoring wells and the digital elevation model

Groundwater monitoring database is available on request and free of charge for academic purposes from the GDOWM. Data request form available at https://www.ovf.hu/rolunk/adatigenyles. The provided database encompasses not only the time series of GWLs but also includes the identification numbers and designations of all 1504 wells, along with the ground elevation at each monitoring well, the position of the well and the well rim elevation. The GDOWM has provided the 2014 HYDRODEM 25 × 25 m resolution digital elevation model too. In accordance with this, the HD72/EOV (EPSG:23,700) coordinate system was used in this study.

Ecosystem base map of Hungary

Land use categories can be categorized based on the Ecosystem Base Map of Hungary. The mapping was conducted using data between 2015 and 2017 and offers a more detailed resolution compared to, for instance, Corine Land Cover (MA 2019, Tanács et al. 2019). Data are available at www.alapterkep.termeszetem.hu.

Soil database

As auxiliary databases, Digital, Optimized, Soil Related Maps and Information in Hungary (DOSOREMI) map for Hungary (Pásztor et al. 2017) (available at https://dosoremi.hu/maps/textura-osztaly-usda-0-5-cm/), and the physical soil diversity map of AGROTOPO database (Várallyay et al. 1979, 1980; MTA TAKI 1991) (available at https://www.elkh-taki.hu/hu/keptar/agrotopo) was used. Although both soil databases can be viewed online, data are available upon request.

Aggregated vertical hydrological model

Approach of the model construction

The effects influencing changes in groundwater resources were analyzed in two components:

  • Vertical impact arising from meteorological factors, namely precipitation and evaporation,

  • Horizontal flow caused by topographical effects, surface inhomogeneities, or anthropogenic interventions.

A model description can be associated with each impact. The vertical effect is more hydrological, while the horizontal flow is more hydraulic, which is well described by the Boussinesq equation of seepage (Boussinesq 1904). It is judicious to separate the calculation of these two impacts into two models. This is due to the relatively small inclines of the GHP, minimizing the role of horizontal flow in the water balance compared to the vertically induced meteorological impact. A comparison of the magnitudes of these two impacts reveals that the seasonal volume change per unit area (1 km2) on the GHP, calculated from the minimal and maximal groundwater levels (Supplementary Figs. S47, 48), is estimated to be between 4 × 105 and 5 × 105 m3 km−2 year−1. The replenishing effect from horizontal flow in the same area is estimated to be between 0.05 × 104 and 0.5 × 104 m3 km−2 year−1, depending on the variability of the soil hydraulic conductivity. Therefore, the vertical impact plays a dominant role in shaping water level dynamics.

The changes were described by applying and developing dynamic models reflecting time and space. The aggregated vertical hydrological model aims to describe changes in groundwater levels. Richards' equation (Richards 1931) is employed to characterize the flow of water in variably saturated porous media. Despite its widespread application in practical scenarios, obtaining a numerical solution for the equation poses a challenge. The reliable and accurate resolution of Richards' equation for a specific set of boundary conditions and soil hydraulic properties remains an ongoing challenge (Farthing and Ogden 2017). Research has already been done on a simple, robust description of the equation (Ireson et al. 2023).

The schematic diagram shows the main physical processes considered during this research (Fig. 4). During calculations, the term "groundwater level" refers to the distance between the soil surface and the groundwater table beneath it (meters below ground level, m b.g.l.). The goal is to avoid the depth-resolved description of vertical moisture transport, thereby significantly reducing computation time. However, the parametrization of the aggregated model closely resembles the Richards model, allowing the results of calibration to be transferred to models with a more detailed description in space.

Fig. 4
figure 4

Schematic diagram of the aggregated vertical hydrological model. The parameters are the following. α is the precipitation correction coefficient (−), Φh and Φv are the horizontal and vertical material fluxes from the regional groundwater flow (mm yr−1), Eeptr is the is the evapotranspiration (mm d−1), Eeva is the groundwater evaporation via capillary mechanisms (mm d−1), Eint is the evaporation from water stored on the leaf surface (mm d−1), Etrans is the transpiration (mm d−1), h is the distance from the ground surface (m), hr is the average depth of the root zone (m), I is the interceptional storage on leaf surface (mm d−1), L represents the distance between the surface and the groundwater table (m), P0 is the measured daily precipitation (mm d−1), Peff is the effective precipitation, the portion of daily precipitation reaching the soil surface (mm d−1), Pcorr is the runoff-corrected part of precipitation reaching the surface, S is the degree of saturation (−), Sr is the equilibrium saturation of the root zone (−), SS is the soil surface pore saturation (−)

Full size image

The meteorological impacts provide the upper boundary conditions for models describing subsurface hydrodynamics. Therefore, accurate formulation of these impacts is essential. The input data for the aggregated vertical hydrological model encompass daily mean air temperature (°C), daily precipitation (mm), and relative humidity (−). The model operates using a 14-day time step (while the time step's length can be adjusted at will, this interval has proven a commendable compromise for tracking dynamic changes and attaining simplification).

Considerations of the model development

The mathematical relationships describing the processes are included in the Supplementary material; only the considerations are summarized here.

At the core of regional evaporation description lies the water retention function, which characterizes the equilibrium saturation within the three-phase zone (unsaturated zone or vadose zone) by the combined effects of advection forces and the capillary mechanism. In this study, the Kovács retention curve, developed on a theoretical basis, was used (Kovács 1981). Despite documented comparisons with other relationships (Aubertin et al. 2011; Mbonipma et al. 2011), the Kovács function's advantage lies in its physical parameters (Maqsoud et al. 2017), yet it tends to overestimate equilibrium saturation in the lower centimetric logarithmic head (known as pF) range. Parameters \({h}_{c0}\) and \(n\) are subjected to calibration. The parameter search range, within which the optimal value of \({h}_{c0}\) was derived from the DOSOREMI database (Pásztor et al. 2017). The water retention curve is defined by porosity, average capillary rise height, and suction height (van Genuchten 1980). Saturation is related to suction potential through theoretical relationships encoded by Supplementary material Eqs. (5–8) (Kovács 1981).

Evapotranspiration calculations utilize a simplified function (Dunay et al. 1969; Varga-Haszonits et al. 2019) (see in Supplementary material Eqs. 9–13, 23). Choices include the Penman formula (Penman 1948) or the FAO Penman–Monteith equation (Allen et al. 1998). However, these formulas require additional parameters, like soil heat flux density, wind speed and psychrometric constant (Allen et al. 1998). Since the meteorological databases based on measurements, like HMS, do not contain these additional parameters, the use of a simplified function is favorable. For the calculation, leaf area index (LAI) was automatically calibrated within the specified value range based on LULC map of Ecosystem Base Map of Hungary (MA 2019, Tanács et al. 2019) using a global optimization procedure. For leaf surface storage, the phenomenon of interception, for description the effects of snow accumulating and melting, and for the effective precipitation the relationships from ARES model was utilized (Koncsos and Szabó 2003, Koncsos 2008). Several models were analyzed for calculating surface runoff. One of them was the Horton model describing infiltration, which allowed us to derive the fraction of precipitation above infiltration calculated in the time step (Horton 1933, 1945). Ultimately, the simplest approximation provided the best fit from the results. For the runoff-corrected precipitation, an α precipitation correction coefficient was introduced that corrects the precipitation reaching the surface with the runoff magnitude (see in Supplementary material Eqs. 19–20).

Hence, the study area located in the Carpathian Basin the regional groundwater flows was considered. Driven by differences in groundwater levels, regional subsurface water movement is present in the study area due to its basin-like nature. The conceptual model of the flow systems of the GHP’s confined aquifers was developed by Erdélyi (1975, 1981) based on Tóth's theory (Tóth 1962, 1963). In the upper portion of the GHP, within the depth range of approximately 400–1700 m, waters flow in open flow systems driven by gravity, receive input from precipitation, and return to the surface at distinct discharge areas. Thus, the hydrogeological system in the basin occasionally enables upward flow of confined groundwater, even to the water table (Tóth 2006, Mádlné 2006). Flow from highland areas of the GHP is characterized by recharge, while the valleys of the Tisza River and its tributaries primarily act as discharge areas. The correctness of this theory has been substantiated by isotope hydrology data (Deák 2013). According to 14C carbon isotope data, the vertical flow velocity of confined groundwater in the Pleistocene formation varies between approximately −45 and 56 mm yr−1 (Deák 2013). Positive values signify downward flow (or inflow), while negative values represent upward flow (or discharge). More recent research in the area between the Danube and the Tisza calculated a value of −42–9 mm yr−1, discharge areas are usually characterized by a higher GWL (Garamhegyi et al. 2018, 2020). The upward flow reaching the water table is hydrologically conditioned by the confined groundwater level surpassing the static water table and the presence of appropriate geological connections (Rakonczai and Fehér 2015). However, it should be noted that investigations in the South GHP did not detect the impact of regional upward flow. This could be attributed to the negligible proportion of upward flow relative to the groundwater changes predominantly driven by meteorological influences (Rakonczai 2013). The effect of the regional subsurface flow system at the lower boundary was considered with a constant perimeter-specific volume flux (vertical material flux) (\({\Phi }_{{\text{v}}}\)) assumed within a year. Its magnitude, according to previous research, is calculated as −20–25 mm/year considering the vertical flow rate from −45–56 mm/year (Deák 2013). The parameter optimum was determined from this range using a calibration algorithm. Horizontal material fluxes (\({\Phi }_{h}\)) were neglected. This introduces errors in the calculations; however, the error can be eliminated by linking the model and the Boussinesq equation. In this linked model, the aggregated vertical hydrological model provides the upper boundary condition material flux based on meteorological factors. The error elimination step will be analyzed in an article presenting another stage of research; the presentation of the horizontal material fluxes was omitted in this paper.

A modeling framework was established in which the HYDRODEM, AGROTOPO, DOSOREMI, GWL database and the Ecosystem Base Map of Hungary are interconnected on a 50 × 50 m resolution grid, and the meteorological time series was available at the randomly chosen GWL monitoring well positions. The model was built on this modeling platform.

Calibration and validation

The calibration and validation of the model were conducted based on meteorological databases available for the period between 1971 and 2010. The length of the time series available for individual SGW monitoring wells significantly influenced the potential for calibration. Proportional to the length of the GWL time series, approximately 40% of the data from the beginning of the observation period were employed for calibration, while the remaining 60%, not used in calibration, were reserved for validation. The use of a time series in a separate set for calibration and validation is a used hydrological modeling practice (Wunsch et al. 2022).

During calibration, the first 700 days were excluded from both the measurement and calculation periods of the objective function. This is sufficient, based on experience, to eliminate initial condition errors from the model. The applied calibration procedure is the self-developed BLIND algorithm (Koncsos et al. 1995), a global optimization gradient-free algorithm that generates a stochastic search sequence converging to the global optimum. The task is to align the model with the measurement results through the minimization of the so-called objective function.

The BLIND algorithm continuously narrows down the search interval, reducing the number of trial steps. The model's calibrated parameters are shown in the Supplementary material.

The efficiency of the calibration was evaluated with two indicators. One is the Nash–Sutcliffe efficiency (NSE), which is the objective of the calibration. The NSE index varies between \(\left]-\infty ;1\right]\) where the former indicates poor fit and the latter indicates perfect fit (Nash and Sutcliffe 1970). NSE is computed using Eq. (1).

$${\text{NSE}}=1-\frac{\sum_{j=1}^{N}{\left({L}_{o}^{j}-{L}_{m}^{j}\right)}^{2}}{\sum_{j=1}^{N}{\left({L}_{o}^{j}-\overline{{L }_{o}}\right)}^{2}}$$
(1)

where \({L}_{o}^{j}\) represents the measured groundwater level daily data at the end of the j aggregated time step (m), \({L}_{m}^{j}\) is the modeled groundwater level daily data at the end of the j aggregated time step (m), \(\overline{{L }_{o}}\) is the average of the measured groundwater level time series (m), and \(j=1\dots N\) is the number of aggregated time steps. The calibration objective function was the Nash–Sutcliffe efficiency (NSE) value.

While calibration was done based on the NSE value, the quality of the model was also examined using the root-mean-square-error (RMSE) value according to Eq. (2). RMSE (m) can vary in the range \(\left[0;\infty \right[\) and is better as it gets closer to 0 m.

$${\text{RMSE}}=\sqrt{\frac{\sum_{j=1}^{N}{\left({L}_{o}^{j}-{L}_{m}^{j}\right)}^{2}}{N}}$$
(2)

NSE and RMSE are widely used indicators in hydrological modeling (Hosseinizadeh et al. 2019, Topál et al. 2020, Sfîcă et al. 2022, Wunsch et al. 2022, Zhang et al. 2022, Sheikha-BagemGhaleh et al. 2023). In modeling practice, an NSE > 0.7 is generally considered a very good result. During validation, only monitoring wells meeting the criteria NSE > 0.4 and RMSE < 0.5 m are accepted as validated. Consequently, out of the examined 557 monitoring wells, 463 are deemed acceptable for further analyses.

Estimation of climate change effects on shallow groundwater levels

In this phase of the study, exclusive focus is placed on the impact of weather on groundwater levels. The analysis is grounded in the previously outlined simple yet robust aggregated vertical hydrological model, which serves for the estimation of groundwater levels. The meteorological input data are provided by FORESEE HUN v1.0. In total, 14 climate projections were examined for the RCP 4.5 and RCP 8.5 scenarios (Table 2). Linear trend analysis was applied to the groundwater levels modeled based on climate projections. The trend analysis was conducted over the time horizon extending to 2050 and 2100. Descriptive statistics based on the trends were compiled for all 463 validated monitoring wells.

Table 2 List of the model projections in FORESEE-HUN v1.0 as combinations of GCM and RCM simulations (FORESEE 2023)
Full size table

Results

Test of the meteorological databases

To investigate the future changes in shallow GWLs, the expected variations of individual meteorological factors (model input data) must be analyzed.

A spatial consistency test was performed on the observed time series of all of the meteorological databases to investigate their reliability. The data from the four databases were compared with each other on scatter plots with the use of the queried meteorological time series of selected 76 monitoring wells. Coefficient of determination (R2) between the individual databases based on the daily and the 10-day average values between 1971 and 2010 was calculated. To compare the time series of two databases, we calculated the minimum, average and maximum values of R2 and a linear trend line was fitted (Supplementary Figs. S20–S45 and Supplementary Tables S1, S2). All databases have good correlation in temperature (R2 = 0.99–1.00). By precipitation data, FORESEE v4.0 compared to CARPATCLIM shows an R2 = 0.68–0.93 with a mean of 0.77, while FORESEE HUN v1.0 compared to CARPATCLIM shows an R2 = 0.90–1.00 with a mean of 0.95. It can be derived from the fact, that FORESEE v4.0 uses the E-OBS database (C3S 2022b) while FORESEE HUN v1.0 was calibrated on HMS database (FORESEE 2023).

Both FORESEE databases differ from CARPATCLIM and HMS data by relative humidity, the mean value of R2 = 0.47–0.48 for daily data comparison and R2 = 0.68–0.70 for 10-day averaged data comparison, respectively (Supplementary Tables S1 and S2). FORESEE did not contain RH data, that were calculated from VPD. In case of FORESEE HUN v1.0 for the calculated RH time series, a correction coefficient was introduced with a value K = 1.279 and the 10-day averaged values was calculated. With the K correction coefficient, the R2 value did not improve significantly, but the correlation did compared to the plots (Supplementary Figs. S38–S39 and S40–S41).

Based on the test, the CARPATCLIM and the FORESEE HUN v1.0 databases were considered to further analysis.

Model calibration and evaluation

The aggregated vertical hydrological model was calibrated and validated using CARPATCLIM and FORESEE HUN v1.0 datasets. In case of CARPATCLIM, the time series are available from 1961, so the calibration can be started from 1961 if the examined GWL monitoring well has observed data too from this year. In case of the FORESEE HUN v1.0 database the start date of the calibration was 1/1/1971 if the GWL time series was given at the examined well. The calibration was automated. Calibration results for a few monitoring wells as examples depict the model’s similarity to measurements shown in Supplementary Figs. 47–48.

The validation phase entails demonstrating the model's capability to accurately replicate independent (non-utilized in calibration) data. In this context, the focus lies on assessing the model's predictive capacity over the selected observational period and evaluating its ability to anticipate future changes. The principle of validation involves simulating the calibrated parameter set on individual wells during a timeframe not employed in calibration. The observed and modeled time series show a good match, the model describes the changes in the groundwater level well (Fig. 5 and Supplementary Fig. S50).

Fig. 5
figure 5

Calibration and validation of well No. 2624 based on CARPATCLIM database. The time series of groundwater levels available between 1970 and 2012. Validation from 1986 onward. Lo (blue) is the observed, and Lm (orange) is the modeled groundwater level. NSE value of validation 0.73

Full size image

The assessed model performance for the 557 monitoring wells is the following. In case of CARPATCLIM, during calibration the median value of NSE = 0.75 with a 0.66–0.82 interquartile range (IQR) and a median value of RMSE = 0.26 m (IQR = 0.20–0.34 m), during validation the median value of NSE = 0.62 (IQR = 0.40–0.74) and a median value of RMSE = 0.35 m (IQR = 0.26–0.50 m). In case of FORESEE HUN v1.0, during calibration the median value of NSE = 0.71 (IQR = 0.60–0.79) and a median value of RMSE = 0.31 m (IQR = 0.24–0.40 m), during validation the median value of NSE = 0.71 (IQR = 0.57–0.80) and a median value of RMSE = 0.31 m (IQR = 0.24–0.54 m). Although the CARPATCLIM-based calibration shows better result, the FORESEE-based model gives hardly the same performance during calibration ad validation. Moreover, the validation results of the FORESEE-based model are better than the CARPATCLIM-based model. The results of the model performance assessment are available in the Supplementary Tables S4 and S5.

For the further analysis, a filter criterion was used based on NSE and RMSE values. The used filter result that 463 out of 557 monitoring wells are acceptable. In the case of validation of FORESEE database the mean value of NSE improved and the IQR narrowed, NSE = 0.72 (IQR = 0.62–0.80). The RMSE also improved, RMSE = 0.30 m (IQR = 0.24–0.38). With the filtering, the CARPATCLIM based model performs similarly good as the FORESEE based (Fig. 6). For the violin plots the Violin Plots for Matlab code was used (Bechtold 2016). In summary, through the analysis of a substantial number of GWL monitoring wells, it has been determined that the model has been validated and is suitable for analyzing future scenarios.

Fig. 6
figure 6

Filtered NSE (−) (left column) and RMSE (−) (right column) values of the model calibration (first row) and validation (second row) for 463 monitoring wells. The outline of the violin plot covers the probability density of the indicator values. The median of the data is indicated with a white circle; the IQR represents with shading. The blue violin plots refer to the CARPATCLIM (CC); the orange plots refer to the FORESEE HUN v1.0 (FSH) based calculations

Full size image

Impacts of future processes based on meteorological scenarios on groundwater level changes

As an example, the results for the monitoring well No. 2624 will be discussed (Fig. 7). The median GWL shows decreasing. The multi-model median of the linear trend to 2100 in case of RCP 4.5 scenario shows −0.49 cm yr−1 (IQR = −0.68–−0.29) and the mean is −0.54 cm yr−1; in case of RCP 8.5 scenario the median trend is −0.31 cm yr−1 (IQR = −0.65–−0.15) and the mean is −0.54 cm yr−1. The multi-model median of the linear trend to 2050 in case of RCP 4.5 scenario shows −0.80 cm yr−1 (IQR = −1.16–−0.60) and the mean is −0.89 cm yr−1; in case of RCP 8.5 scenario the median trend is −0.48 cm yr−1 (IQR = −1.22–−0.25) and the mean is −0.72 cm yr−1. Although in both scenarios the mean trend to 2100 shows a −0.54 cm yr−1 decreasing of the GWL, in an unpredictable way the trend to 2050 shows a greater decrease in the GWL, −0.89 cm yr−1 and −0.72 cm yr−1, respectively.

Most model projections show a decrease in both the 2050 and 2100 time horizons. In the case of the linear trend analysis until 2100, under the RCP 4.5 scenario, the NCC HIRHAM5 model indicates a slight increase, while under the RCP 8.5 scenario, the EC-EARTH HIRHAM5 and MPI-REMO2009 r2 models indicate an increase. For the time horizon up to 2050 under the RCP 4.5 scenario, EC-EARTH RACMO22E r12 is the mildest, but the IQR remains in the negative range. In the case of the RCP 8.5 scenario, the HadGEM2 RACMO22E model can also be characterized by an increasing trend. The most unfavorable forecast is given by HadGEM2 CCLM (this can be seen in Fig. 7 as well), except for the trend analysis until 2050 under the RCP 4.5 scenario where the EC-EARTH HIRHAM5 was the worst (Fig. 8).

Fig. 7
figure 7

Estimations of the changes in GWL until 2100 based on RCP 4.5 (first row) and RCP 8.5 (second row) scenarios. In the first column, the colored lines show the individual models. In the second column, the solid purple line indicates the median GWL, the blue shading is the IQR, and the light blue shading is the 5–95% percentile range

Full size image
Fig. 8
figure 8

Linear trend values of RCP 4.5 (−) (left column) and RCP 8.5 (−) (right column) scenarios of time horizon 2100 (first row) and time horizon 2050 (second row) for 463 monitoring wells. The outline of the violin plot covers the probability density of the trend values. The median of the data is indicated with a white circle; the IQR represents with shading. The order of the individual models is the same as in Table 2

Full size image

Based on the results of linear trend analysis conducted for the two time horizons across 28 climate scenarios, statistical calculations were performed to generate multi-model outcomes regarding trends for individual monitoring wells. Percentiles (25% and 75%), median, mode, expected value, variance, and standardized normal distribution were determined. Descriptive statistics can be found in Supplementary Tables S6–S9. The results of multi-model mean were also visualized on a map (Fig. 9). For the time horizon until 2100 under RCP 4.5 (Fig. 9a) and RCP 8.5 (Fig. 9b) scenarios, there is a mere 1 cm yr−1 difference in the average expected trend values. However, for the time horizon until 2050, a significant difference of 32 cm yr−1 was identified between the average expected values under RCP 4.5 and RCP 8.5 scenarios (Table 3). Estimates for the 2050 time horizon exhibit more pronounced decreasing trends. A growing trend is not evident at any monitoring well for the 2050 time horizon under the RCP 4.5 scenario (Fig. 9c). However, under the RCP 8.5 scenario, some wells above the Kiskörei Hydropower Plant’s Reservoir (‘Tisza-lake’) exhibit a positive trend (Fig. 9d). The spatial distribution of changes is illustrated on the map (Fig. 9).

Fig. 9
figure 9

Multi-model mean values based on linear trend analysis. a Time horizon up to 2100, under the RCP 4.5 scenario. b Time horizon up to 2100, under the RCP 8.5 scenario. c Time horizon up to 2050, under the RCP 4.5 scenario. d Time horizon up to 2050, under the RCP 8.5 scenario

Full size image
Table 3 Descriptive statistics of the multi-model multi-scenario analysis
Full size table

As a summary of the examined climate projections, it can be seen that the multi-median values of the linear trends all indicate a decreasing trend, and the upper limit of the IQR does not enter the positive range either. The multi-mean values are usually lower than the median value (Table 3).

Discussion

Evaluation of the model

An important observation was that neglecting horizontal water movement did not result in a significant deterioration of computational quality. This confirms that horizontal flow minimally influences groundwater level variations, with vertical transport mechanisms (precipitation, evaporation) being the primary drivers (Rakonczai 2013; Garamhegyi et al. 2018).

There is a good spatial correlation between the calibrated LAI values and the LULC base map (Fig. 10a). At the same time, in the case of soil physical parameters, there is a large discrepancy compared to the AGROTOPO map data. From the calculated hc0 parameter with the use of Wentworth scale (Wentworth 1922) the soil type can be determined (Fig. 10b, Supplementary Table S3). In the case of porosity (Fig. 10c), there is also a difference compared to the ranges known for each soil type. The variations in porosity may be attributed not only to particle size but also to factors such as land use and soil organic matter content (Robinson et al. 2022). In sum, it can be noted that unambiguous model parameters could not be assigned to categories of auxiliary maps (e.g., AGROTOPO). While different parameter values might be calculated based on these auxiliary maps in various locations, the high NSE values from calibration and validation indicate the correctness of the parameters.

Fig. 10
figure 10

Some calibrated parameters based on CARPATCLIM database. a shows the calibrated LAI values displayed on the LULC base map (MA 2019). b shows the soil types that can be derived from the hc0 calibrated parameter values on the AGROTOPO soil physics base map (MTA TAKI 1991). In c the values of porosity are plotted on the AGROTOPO soil physics base map (MTA TAKI 1991)

Full size image

Differences may arise from the significant heterogeneity of soils in reality and unknown, unquantifiable human impacts during modeling (e.g., legal or illegal groundwater extraction, changes in runoff, infiltration, and evaporation due to construction, changes in land cultivation, etc.). Moreover, multicollinearity can occur in models with multiple parameters, meaning that similar good results can be obtained for various parameter combinations.

Groundwater responses to climate change

The time horizon of 2050 can be considered authoritative. With the use of median values, in the case of the RCP 4.5 scenario, the annual water resource loss can be 368 million cubic meters projected on the size of the study area, while in the case of the RCP 8.5, it can be 243 million cubic meters. Surprisingly, under the RCP 8.5 scenario, more favorable groundwater conditions are expected between 2030 and 2080. The reason for this could be the projected rise in precipitation amounts over the study region within the context of this extreme climate scenario (Sfîcă et al. 2022). Projected changes based on FORESEE HUN v1.0 may explain the more favorable trends anticipated by the end of the century. This could be due to a 5–10% annual precipitation increase in the Great Hungarian Plain by 2100 for RCP 4.5 and 10–20% for RCP 8.5 scenarios. Simultaneously, temperature changes are projected to be 1.4–1.8 °C for RCP 4.5 and 2.8–3.2 °C for RCP 8.5 scenarios annually (FORESEE 2023). In other words, smaller temperature changes are associated with a lesser increase in precipitation, while larger temperature changes are linked to a greater precipitation increase. These two effects may result in the expected values of projected groundwater level changes being nearly identical by 2100.

In Hungary, a total of 133,126 hectares was irrigated in 2022, with the GHP accounting for 111,801 hectares. On average, around 163 million cubic meters of irrigation water are supplied across the GHP. These irrigated areas constitute a mere 2.5% of the agricultural land (HCSO 2023a, 2023b). The vulnerability of crop cultivation is well reflected by the fact that drought damage exceeded 49 billion forints (approximately 136.7 million euros) in 2022 (MA 2023). It can be seen that the irrigation water cannot replenish this amount. Moreover, the fact that human impacts were not considered during modeling implies that the potential effects of irrigation are already reflected in the model results. Therefore, the identified loss in water reserves should be interpreted beyond the replenishing effects of the current water management system.

The fluctuation of groundwater levels significantly impacts quantifiable ecosystem services, such as food production. Major arable crops like wheat or maize are highly sensitive not only to climate change but also to changes in groundwater levels, yet this sensitivity is frequently overlooked. Its importance has also been documented in the literature. A 100 mm change in groundwater level on the GHP induces an approximate yield variation of 0.23 t ha−1 for maize (Pinke et al. 2020). In a favorable case, the changes in GWLs can cause 0.01 t ha−1 average annual loss for maize.

The examined area is located within the Tisza River basin, which, due to its flood hazards and frequent droughts, is characterized by extremes in water management. The increasingly severe aridity of the area can be most effectively addressed by proper replenishment of groundwater. The amount of water needed to recharge the groundwater supply and the technical solutions involved can be considered crucial. On a smaller scale, local water retention may offer a solution if soil physical properties are suitable for these methods and decision-makers and stakeholders can agree on their application, like beaver ecosystems, beaver dam analogs, and nature-based water retention measures (NWRM) (Feiner and Lowry 2015, Dittenbrenner et al. 2018, Hartmann et al. 2019, Wolf and Hammill 2023). On a larger scale, harnessing excess water from river floods and managed aquifer recharge (MAR) could provide a solution for declining groundwater levels (Dillon et al. 2019; Hartmann et al. 2019). As a solution, managed aquifer recharge is possibility on the GHP. The equilibrium of the water balance can be improved by involving the Tisza’s flood water resources. In earlier studies, several areas were identified along the Tisza River that are suitable for storing and distributing excess floodwater (NWRM + MAR system). Since the water governance system must adhere to natural functioning, in these studies the focus was solely on gravity-based water replenishment (Murányi & Koncsos 2022a, 2022b, 2022c).

Conclusion

The central hypothesis embodied in the self-developed aggregated hydrological vertical model posits that by utilizing longer, aggregated (in our case, 14 day) time steps, precipitation infiltration to groundwater can be achieved by omitting detailed cross-sectional descriptions. Applying aggregated (multi-day) time steps enables the depiction of average changes within the time step. The results supported the working assumption that the effect of moisture transport at upper boundary values (precipitation) on groundwater level changes, regardless of the upper boundary value, demonstrates a constant proportion between the j-th and j-1-th steps, irrespective of the j value.

The other central hypothesis, verified by the model's outcomes, was that the average value of soil moisture influencing evapotranspiration can be well described within multi-day time steps both at the soil surface and at the root zone height, based on the equilibrium retention curve set by groundwater level, within the j-th aggregated interval.

The aggregated vertical hydrological model facilitated nationwide analysis involving over 500 monitoring wells, offering an acceptable computational speed and accurate depiction of reality. The evaluation of the model was also considered successful. This allows us to capture the range and unique evolution of all future climate projections.

Utilizing the FORESEE HUN v1.0 database, an effective analysis of future groundwater level trends was enabled. In this way, a tool was developed to examine the impact of possible interventions in the water management system, this study can be treated as a “no intervention” water management scenario. Changes in GWL conditions are expected to be less favorable until 2050 than in the case of trends that can be predicted until 2100.

Because wells with high forecast accuracy in the past were deliberately chosen (refer to "Model calibration and evaluation"), it can be inferred that climate forcings predominantly influence groundwater dynamics at these wells. Given that there is no fundamental alteration in system relations (such as newly installed groundwater pumping or significant changes in land use in the vicinity), reasonable results can be anticipated for our simulations. This expectation is based on our focus solely on the impact of changing climate, assuming that other boundary conditions remain constant.