Modelling the Impacts of Climate Change on Shallow Groundwater Conditions in Hungary

: The purpose of the present study was to develop a methodology for the evaluation of direct climate impacts on shallow groundwater resources and its country-scale application in Hungary. A modular methodology was applied. It comprised the deﬁnition of climate zones and recharge zones, recharge calculation by hydrological models, and the numerical modelling of the groundwater table. Projections of regional climate models for three different time intervals were applied for the simulation of predictive scenarios. The investigated regional climate model projections predict rising annual average temperature and generally dropping annual rainfall rates throughout the following decades. Based on predictive modelling, recharge rates and groundwater levels are expected to drop in elevated geographic areas such as the Alpokalja, the Eastern parts of the Transdanubian Mountains, the Mecsek, and Northern Mountain Ranges. Less signiﬁcant groundwater level drops are predicted in foothill areas, and across the Western part of the Tisz á nt ú l, the Duna-Tisza Interﬂuve, and the Szigetköz areas. Slightly increasing recharge and groundwater levels are predicted in the Transdanubian Hills and the Western part of the Transdanubian Mountains. Simulation results represent groundwater conditions at the country scale. However, the applied methodology is suitable for simulating climate change impacts at various scales. points calibration preferred value Tikhonov regularization two-step calibration calibration


Introduction
Groundwater is the world's largest distributed store of fresh water, and as such, it plays a central part in sustaining ecosystems and human population. Groundwater comprises almost one-third of the Planet's freshwater reserves, and 96% of liquid-phase freshwater [1].
The climate projections of the Intergovernmental Panel on Climate Change (IPCC) indicate a significant temperature rise and alterations in the amount and frequency of precipitation throughout the rest of the 21st century [2][3][4]. Changing climate variables influence the hydrological cycle through impacting on [5] Precipitation has a direct impact on groundwater recharge and an indirect impact on human groundwater withdrawals. Recharge is very sensitive to precipitation, and small changes in rainfall can result in more significant changes as far as recharge is concerned [6]. Higher temperature increases evapotranspiration rates, which may in turn Water 2021, 13, x FOR PEER REVIEW 3 of 25 Hungary has a continental climate, with hot summer seasons. The winters are cold and snowy. The average value of annual temperature is 9.7 °C. The extreme values of temperature range between −29 and 42 °C. The average temperature in the summer is 27-35 °C and in the winter it is 0 to −15 °C. The average value of annual rainfall is approximately 600 mm. The annual temperature distribution of Hungary averaged between 1961 and 1990 is indicated in Figure 2. In central Europe, both the maximum and minimum temperatures warmed over the 20th century. The most significant rise was manifested at minimum temperature. There was more significant summer warming in the first half of the century than in the second half [13].
The analysis of climate data indicates that heat waves became more frequent, longer, and more intense in every season in the entire Carpathian Region. Cold waves show the  The distribution of annual rainfall of Hungary averaged between 1961 and 1990 is indicated in Figure 3. According to observed data, there was a significant decrease in the amount of annual precipitation during the 20th century. The most significant drop in precipitation occurred during the spring seasons, when precipitation dropped by 25% since the beginning of the 20th century. The summer precipitation remained stagnant during the past 100 years. The autumn and winter precipitations showed a decrease by approximately 12-14%. In parallel with decreasing trends, the precipitation events show a more intensive pattern, which is likely to increase run-off rates and the risk of floods. Furthermore, it decreases recharge rates and groundwater resources. The results of the ALADIN and RegCM regional climate modelling studies [15] indicate that annual average temperature in Hungary is expected to rise by 3-5 °C by the end of the century. Highest warming rates are expected in the summer season. The warming will be more intense in Eastern Hungary than in the western part of the country. The predicted annual temperature change between the reference 1961-1990 and 2071-2100 periods is indicated in Figure 4. The results of the ALADIN and RegCM regional climate modelling studies [15] indicate that annual average temperature in Hungary is expected to rise by 3-5 • C by the end of the century. Highest warming rates are expected in the summer season. The warming will be more intense in Eastern Hungary than in the western part of the country. The predicted annual temperature change between the reference 1961-1990 and 2071-2100 periods is indicated in Figure 4.
The above climate projections indicate the decrease of annual rainfall amount until the end of the century. At the same time, the temporal distribution of rainfall is predicted to become more balanced with the wettest summer season, which become dryer, while fall seasons become wetter during the following decades. The predicted annual rainfall change between the reference periods of 1961-1990 and 2071-2100 is indicated in Figure 5.  The above climate projections indicate the decrease of annual rainfall amount until the end of the century. At the same time, the temporal distribution of rainfall is predicted to become more balanced with the wettest summer season, which become dryer, while fall seasons become wetter during the following decades. The predicted annual rainfall change between the reference periods of 1961-1990 and 2071-2100 is indicated in Figure  5.

Applied Methodology
We developed a modular approach to the quantitative investigation of the influence of climate conditions on the groundwater table: While the simulation of past conditions was undertaken based on measured climate parameters, prediction of future groundwater conditions was undertaken making use of the ALADIN climate model outputs produced by the Hungarian Meteorological Service [15]. A more detailed explanation of the applied methodology is provided in [20].
Any inaccuracies and uncertainties included in climate model outputs are inherited by the groundwater models and thus reflected in model uncertainty.
Within this study it was assumed that the averaged groundwater conditions of 1961-1965 represent the natural state of the shallow groundwater system, without any significant influence from mining or water extraction operations. This assumption was based on the fact that major mine dewatering operations were accelerated after this period. This time period served as a reference condition, for which model parameters were calibrated against observed groundwater levels. Each predictive simulation was undertaken making use of the initial calibrated parameter set.

Climate Data
Throughout this study, the CARPATCLIM [21] database was applied as a source of main input parameters. CARPATCLIM is a homogenized raster dataset interpolated from climate observations within the Pannonian Basin [16]. The dataset was derived from weather observations at 258 climate stations and 727 precipitation stations. In Hungary 37 climatological and 176 precipitation stations were applied [14]. The CARPAT-CLIM project area included 9 countries (Croatia, Austria, Ukraine, Romania, Serbia, the Czech Republic, Slovakia, Poland, and Hungary).
The database includes 0.1 • (about 10 km × 10 km) resolution homogenized and gridded datasets of daily meteorological variables and climate indicators measured between Water 2021, 13, 668 7 of 23 1961 and 2010. The data grid was obtained through the multiple analysis of series for homogenization software (MASH version 3.03; [22]) and the meteorological interpolation based on surface homogenized databases (MISH, version 1.03; [23]). Meteorological data (daily temperature and precipitation, global radiation, data for evapotranspiration input, average wind speed, and relative humidity) of the CARPATCLIM project, averaged for each recharge polygon, were used as input parameters in hydrological models (HELPs) for recharge calculation [24,25].
Outputs of the ALADIN regional climate model were used as boundary conditions for predictive simulations. The ALADIN climate model was developed by an international consortium led by Meteo-France. The basic model concept is described in [25][26][27][28][29]. ALADIN is widely used by Central and Eastern European meteorological services, as it is computationally efficient. Originally, it was designed for weather forecasts as a downscaling tool of the ARPEGE global model.
ALADIN is a spectral model with grid-based physical parametrisation computation. Technical details of the ALADIN model are described in [30]; regionalmodelprojections are available at [31]. Model simulations were undertaken with the ALADIN-CLIMATE model version for the Carpathian basin with a 10 km discretization. The model used the Special Report on Emissions Scenarios, SRES [32] A1B emission scenario. The SRES emission scenario A1B corresponds to the RCP6.0 scenario [33].

Climate Classification
Zonation of interpolated climate data was necessary, as one-dimensional modelling tools were applied to calculate soil water balance necessary for the assessment of groundwater conditions. Climate zones separate geographic areas with different climate characteristics. These zones were analysed separately in the course of hydrological modelling.
Out of the internationally accepted biophysical climate classification methods, the Köppen [34], the Holdridge [35] and the Thornthwaite [17] methods were applied in Hungary. The comparative analysis of these methods were made by [36] and indicated that these methods were consistent. The Thornthwaite's method was qualified to be appropriate for the mezo-scale characterization of the climatic diversity of Hungary [37]. The methodology described in [38] was applied for the calculation of Thorntwaite climate zonation.
Climate zones were determined for various time periods using average monthly values of climate parameters. The calculation was undertaken on time-averaged parameter grids on cell-by-cell basis. The resulting Thorntwaite grids were converted to polygons for further data processing and visualisation. Thorntwaite zones were applied for the delineation of recharge zones together with geology, landuse, and slope conditions. Recharge values were calculated for each time period by applying the averaged climate parameters of the corresponding recharge zone. The climate zonation for the 1961-1990 investigation period based on CARPATCLIM data is indicated in Figure 6.

Recharge Zones
Recharge zones used in this study are hydrogeological units, in which recharge conditions are assumed to show an insignificant variability. Recharge zones are also called hydrological response units according to the SWAT modelling methodology [38].
Recharge zones were delineated as a superposition of four data layers including climate zones, surface geology, landuse, and slope conditions. The surface geological map constructed by [39] was applied in the first data layer. Geological formations were reclassified into six lithological categories such as fractured, dolomite, limestone, fine porous, coarse porous, and surface waters. Landuse polygons were derived from the CORINE land cover map [40]. The large number of original landuse categories was regrouped into six main classes such as urban areas, arable land, pastures, permanent crops, forests, and water bodies. Slope categories were determined based on the 50-m resolution digital elevation model of Hungary. Two slope categories were applied such as flat areas (0-5%) and slopes (>5%). The resulting map of recharge zones is indicated in Figure 7. odology described in [38] was applied for the calculation of Thorntwaite climate zonation.
Climate zones were determined for various time periods using average monthly values of climate parameters. The calculation was undertaken on time-averaged parameter grids on cell-by-cell basis. The resulting Thorntwaite grids were converted to polygons for further data processing and visualisation. Thorntwaite zones were applied for the delineation of recharge zones together with geology, landuse, and slope conditions. Recharge values were calculated for each time period by applying the averaged climate parameters of the corresponding recharge zone. The climate zonation for the  investigation period based on CARPATCLIM data is indicated in Figure 6.

Hydrological Modelling
The potential effects of climate change on groundwater conditions were represented via water budget calculations for each recharge unit (HRU). The HELP model [18] was applied to calculate daily water balances. The applicability of this model for regional scale assessments is well known from the literature [19] and the methodology has been applied successfully in Hungary. The simulated percolation rates (recharge values) were imported into the numerical groundwater flow model aimed at simulating the groundwater table.
The HELP (Hydrologic Evaluation of Landfill Performance) numerical model was developed by the United States Environmental Protection Agency. Originally, the model serves as an assessment tool for landfill design. The model is based on a water-balance approach and calculates evapotranspiration and drainage through layers of soil and artificial barriers. The model is often used for simulating the effects of various climate scenarios. The weather generator of the HELP model needs several meteorological variables, such as daily and monthly average mean temperature, daily and monthly accumulated total precipitation, monthly average horizontal wind speed, daily global radiation, and monthly relative humidity.
The daily weather data were extracted from the CARPATCLIM dataset and from ALADIN model outputs for the Thorntwaite climate polygons, and daily spatial averages were calculated for each climate polygon. Automated scripts were developed for the purpose of data extractions, calculations, and file conversion operations. coarse porous, and surface waters. Landuse polygons were derived from the CORINE land cover map [40]. The large number of original landuse categories was regrouped into six main classes such as urban areas, arable land, pastures, permanent crops, forests, and water bodies. Slope categories were determined based on the 50-m resolution digital elevation model of Hungary. Two slope categories were applied such as flat areas (0-5%) and slopes (>5%). The resulting map of recharge zones is indicated in Figure 7.

Hydrological Modelling
The potential effects of climate change on groundwater conditions were represented via water budget calculations for each recharge unit (HRU). The HELP model [18] was applied to calculate daily water balances. The applicability of this model for regional scale assessments is well known from the literature [19] and the methodology has been applied successfully in Hungary. The simulated percolation rates (recharge values) were imported into the numerical groundwater flow model aimed at simulating the groundwater table.
The HELP (Hydrologic Evaluation of Landfill Performance) numerical model was developed by the United States Environmental Protection Agency. Originally, the model serves as an assessment tool for landfill design. The model is based on a water-balance Besides meteorological input, the HELP code requires the definition of soil profiles for the calculation of one-dimensional transient water balance. Soil profiles were defined by analysing grain size distributions of soil samples collected systematically as part of the national soil mapping campaign, and organized in a soil logging database. A characteristic soil profile was assigned to each lithological category. Based on grain size distribution data, soil layers were classified according to the United States Department of Agriculture (USDA) soil classification triangle. Default hydraulic parameters defined in HELP were assigned to each soil category. As the uppermost three metres of observed soil profiles show negligible vertical variability, and the average depth of groundwater is within this range, homogeneous soil profiles were applied. The applicability of homogeneous profiles was verified and confirmed through extensive sensitivity analysis. The applied porosity values ranged between 0.3 and 0.46, saturated hydraulic conductivities ranged between 0.05 and 5 m/day, field capacity ranged between 0.05 and 0.25, and wilting point ranged between 0.03 and 0.12.
Simulated percolation rates (recharge) were verified against literature annual values and were also compared with monitoring well hydrographs of selected test sites. Default soil parameters stored in the HELP database were applied and fine-tuned through calibration against observed water level hydrographs, where it was necessary.
The effects of landcover and slope were simulated using a range of runoff curve numbers. The runoff curve number (also called a curve number or simply CN) is a lumped empirical parameter that is applied in hydrology for calculating direct runoff or infiltration from rainfall excess. Applied curve numbers were adjusted in order to obtain realistic recharge rates for each type profile. CN has a range from 30 to 100, where lower numbers indicate permeable soils and low runoff potential and large numbers indicate high runoff potential. The lower the curve number, the more permeable the soil is. Adjusted curve numbers ranged between 91 and 94.
Recharge rates were simulated using the finalised soil profiles for each recharge zone applying spatially averaged climate parameters for the corresponding climate zones. Thorntwaite climate zonation for each simulation was applied from the relevant time period.
The distribution of averaged recharge for the 1961-1965 reference period is indicated in Figure 8. Recharge calculations undertaken on ALADIN model results indicate decreasing recharge rates throughout the 21st century by up to 50 mm/y in the mountainous areas such as the Northern and Transdanubian Mountain Ranges and the Mecsek Mountains. These models do not predict recharge deficiency in the Alpokalja area, however slightly increasing recharge is predicted in parts of the Great Hungarian Plain and the Transdanubian Hills.

Groundwater Modelling
The goal of groundwater modelling was to simulate the water table distribution within the study area under various climate conditions. A two-dimensional steady state model was built for this purpose that was calibrated for transmissivity.
Modelling was performed with the USGS three-dimensional finite-difference groundwater model MODFLOW (MODFLOW-NWT). Calibration was performed with the modelindependent parameter estimation code PEST and Beopest (version 17.05). Visual MOD-FLOW Flex was used as the graphical user interface during most of the modelling except for some of the calibration tasks.
For simplicity, in mountainous areas of open karst terrain, where shallow aquifers are absent, karst water table was simulated, and was considered to be hydraulically connected to adjacent shallow groundwater bodies. Model extent included the political borders of the country. The model domain was horizontally discretized into a rectangular grid with uniform cell sizes of 1000 m × 1000 m. Only one model layer was used extending from 0 m a.s.l. to the land surface elevation. Recharge as calculated and estimated according to the description presented earlier in this paper was selected as a boundary condition for the model at the upper boundary. Major rivers of Hungary and lakes with a surface area larger than 10 ha were selected as constant head internal boundary conditions for the model (Figure 11). Stream and lake stages were selected as the topographic elevation resulting in a slight overestimation of the water table elevation especially in the lowlands. This, however, does not have any negative impacts on our results as we were estimating changes instead of true heads and water table elevations.

Groundwater Modelling
The goal of groundwater modelling was to simulate the water table distribution within the study area under various climate conditions. A two-dimensional steady state model was built for this purpose that was calibrated for transmissivity.
Modelling was performed with the USGS three-dimensional finite-difference ground-water model MODFLOW (MODFLOW-NWT). Calibration was performed with the model-independent parameter estimation code PEST and Beopest (version 17.05). Visual MODFLOW Flex was used as the graphical user interface during most of the modelling except for some of the calibration tasks.
For simplicity, in mountainous areas of open karst terrain, where shallow aquifers are absent, karst water table was simulated, and was considered to be hydraulically connected to adjacent shallow groundwater bodies. Model extent included the political borders of the country. The model domain was horizontally discretized into a rectangular grid with uniform cell sizes of 1000 m × 1000 m. Only one model layer was used extending from 0 m a.s.l. to the land surface elevation. Recharge as calculated and estimated according to the description presented earlier in this paper was selected as a boundary condition for the model at the upper boundary. Major rivers of Hungary and lakes with a surface area larger than 10 ha were selected as constant head internal boundary conditions for the model (Figure 11). Stream and lake stages were selected as the topographic elevation resulting in a slight overestimation of the water table elevation especially in the lowlands. This, however, does not have any negative impacts on our results as we were estimating changes instead of true heads and water table elevations.
Model calibration was performed using historical average water levels in 641 wells and 152 springs between 1961 and 1965. Wells and springs were selected to provide a relatively uniform coverage of the model area ( Figure 11).  Model calibration was performed using historical average water levels in 641 wells and 152 springs between 1961 and 1965. Wells and springs were selected to provide a relatively uniform coverage of the model area ( Figure 11).
Water abstractions were not included in the model simulations. Simulated water table conditions are thus hypothetical distributions with the purpose of demonstrating the direct impacts of climate conditions rather than predicting actual future groundwater conditions. As stated earlier, baseline conditions (natural-state) for the model were the average groundwater levels for the period of 1961-1965 as it can be reasonably assumed that compared with later periods, the shallow groundwater was not yet significantly affected by high industrial, agricultural, and residential water abstractions. This was confirmed by historical monitoring data series. Calibrated transmissivities obtained as a result of these assumptions were directly used for the simulation of the predictive scenarios.
Due to the more-or-less arbitrarily chosen vertical model dimensions, calibrated hydraulic transmissivities do not have a direct physical meaning in terms of shallow aquifer properties, however, they provide information on the coupled impact of aquifer conductivity and aquifer thicknesses. The calibrated distribution of hydraulic transmissivity is shown in Figure 12. The calibrated transmissivities were qualitatively compared to known transmissivities at those locations where the model thickness approaches the true aquifer thickness the most, and found that the calibrated values at these locations are realistic.
ble conditions are thus hypothetical distributions with the purpose of demonstrating the direct impacts of climate conditions rather than predicting actual future groundwater conditions.
As stated earlier, baseline conditions (natural-state) for the model were the average groundwater levels for the period of 1961-1965 as it can be reasonably assumed that compared with later periods, the shallow groundwater was not yet significantly affected by high industrial, agricultural, and residential water abstractions. This was confirmed by historical monitoring data series. Calibrated transmissivities obtained as a result of these assumptions were directly used for the simulation of the predictive scenarios.
Due to the more-or-less arbitrarily chosen vertical model dimensions, calibrated hydraulic transmissivities do not have a direct physical meaning in terms of shallow aquifer properties, however, they provide information on the coupled impact of aquifer conductivity and aquifer thicknesses. The calibrated distribution of hydraulic transmissivity is shown in Figure 12. The calibrated transmissivities were qualitatively compared to known transmissivities at those locations where the model thickness approaches the true aquifer thickness the most, and found that the calibrated values at these locations are realistic.

Calibration Methodology
Model calibration was performed by means of manual and automated calibration using PEST. PEST [41] is a nonlinear parameter estimation code. Parameter optimisation is achieved using the Gauss-Marquardt-Levenberg method to drive the differences between model predictions and corresponding field data to a minimum in a weighted least squares sense. The implementation of this search algorithm in PEST is particularly robust; hence PEST can be used to estimate parameters for both simple and complex models including large numerical spatial models with distributed parameters.
The calibration process was restricted to the adjustment of hydraulic transmissivity until an acceptable match between calculated and measured water levels was obtained.

Calibration Methodology
Model calibration was performed by means of manual and automated calibration using PEST. PEST [41] is a nonlinear parameter estimation code. Parameter optimisation is achieved using the Gauss-Marquardt-Levenberg method to drive the differences between model predictions and corresponding field data to a minimum in a weighted least squares sense. The implementation of this search algorithm in PEST is particularly robust; hence PEST can be used to estimate parameters for both simple and complex models including large numerical spatial models with distributed parameters.
The calibration process was restricted to the adjustment of hydraulic transmissivity until an acceptable match between calculated and measured water levels was obtained. Model calibration was undertaken with the assumption that field measured time-averaged water levels represent steady state (equilibrium) conditions of the groundwater system. As a result of PEST calibration using pilot points, the resulting calibrated transmissivity field is spatially distributed. A number of 592 pilot points distributed over a regular mesh were used during the calibration process. Preferred homogeneity followed by preferred value Tikhonov regularization were used during the two-step calibration process. The total number of PEST calibration iterations resulting in the lowest measurement objective function was much larger (>25 iterations) than the number of iterations (7) after which the resulting transmissivity distribution was considered acceptable. After seven iterations the measurement objective function was lowered to an acceptable value (93248). The resulting mean value of non-zero weighted residuals was equal to −1.37.
A slightly better fit ( Figure 13) was achieved for the observation group consisting of springs than for the well observation group. While most upgradient dots in the calibration plot indicate locations with higher land surface elevation, wells are mainly located in the plain areas at lower elevations. This group can be identified as the dense cluster of points in the calibration plot within the downgradient section of Figure 13. The overall calibration seems satisfactory with an average discrepancy of slightly higher calculated heads than observed that was expected based on the selection of the boundary conditions. Calibrated water levels for 1961-1965 are presented in Figure 14.
Model calibration was undertaken with the assumption that field measured time-averaged water levels represent steady state (equilibrium) conditions of the groundwater system.
As a result of PEST calibration using pilot points, the resulting calibrated transmissivity field is spatially distributed. A number of 592 pilot points distributed over a regular mesh were used during the calibration process. Preferred homogeneity followed by preferred value Tikhonov regularization were used during the two-step calibration process. The total number of PEST calibration iterations resulting in the lowest measurement objective function was much larger (>25 iterations) than the number of iterations (7) after which the resulting transmissivity distribution was considered acceptable. After seven iterations the measurement objective function was lowered to an acceptable value (93248). The resulting mean value of non-zero weighted residuals was equal to −1.37.
A slightly better fit ( Figure 13) was achieved for the observation group consisting of springs than for the well observation group. While most upgradient dots in the calibration plot indicate locations with higher land surface elevation, wells are mainly located in the plain areas at lower elevations. This group can be identified as the dense cluster of points in the calibration plot within the downgradient section of Figure 13. The overall calibration seems satisfactory with an average discrepancy of slightly higher calculated heads than observed that was expected based on the selection of the boundary conditions. Calibrated water levels for 1961-1965 are presented in Figure 14.

Results
As stated above, shallow groundwater levels were first simulated for the period 1961-1965 that also served as the calibration period. Calculated (CARPATCLIM) climate parameters were used for the recharge estimates having been provided as input data for the simulation. Hydraulic properties (transmissivity) as calibrated during this simulation were used for the subsequent 1961-1990, 2021-2050, and 2071-2100 periods, and for which simulated (ALADIN) climate parameters were used as recharge estimates.
Simulated water table distribution for the 1961-1965 baseline period ( Figure 15) and corresponding unsaturated zone thickness (depth to groundwater) presented in Figure 16 indicate that water table depth is generally very close to the ground surface in the lowlands such as the Alföld and Kisalföld areas and exceed depths of 100 m under mountain areas; although these values often correspond to karstwater bodies rather than shallow groundwater table, which might be perched above deep groundwater levels.
With the purpose of demonstrating the direct impacts of climate change on shallow groundwater resources, water table difference maps were constructed. The first of the three study periods  was considered as the reference period and the water table change calculated between this reference period and the 2021-2050 average water levels is presented on Figure 17.
This indicates significant water level drops (up to about -6 m with even higher drops at some places. The latter may also be attributed to less confidence in the simulated water levels at higher altitudes) under highland areas such as the Alpokalja, Eastern part of Transdanubian Mountains, Northern Mountain Ranges, and the Mecsek. A less significant decrease of groundwater levels was simulated along the southern reaches of the Tisza River (<−1 m) and the eastern parts of the Great Hungarian Plain (<−2 m). Most parts of the Transdanubian, the Kisalföld, the western parts of the Transdanubian Mountains, and

Results
As stated above, shallow groundwater levels were first simulated for the period 1961-1965 that also served as the calibration period. Calculated (CARPATCLIM) climate parameters were used for the recharge estimates having been provided as input data for the simulation. Hydraulic properties (transmissivity) as calibrated during this simulation were used for the subsequent 1961-1990, 2021-2050, and 2071-2100 periods, and for which simulated (ALADIN) climate parameters were used as recharge estimates.
Simulated water table distribution for the 1961-1965 baseline period ( Figure 15) and corresponding unsaturated zone thickness (depth to groundwater) presented in Figure 16 indicate that water table depth is generally very close to the ground surface in the lowlands such as the Alföld and Kisalföld areas and exceed depths of 100 m under mountain areas; although these values often correspond to karstwater bodies rather than shallow groundwater table, which might be perched above deep groundwater levels.
With the purpose of demonstrating the direct impacts of climate change on shallow groundwater resources, water table difference maps were constructed. The first of the three study periods ) was considered as the reference period and the water table change calculated between this reference period and the 2021-2050 average water levels is presented on Figure 17. and the 2071-2100 averages indicates that the expected water level drops become slightly more significant at the same locations. Furthermore, the areas characterized by lowering water table elevations increase in size including some areas previously not being negatively affected such as the Duna-Tisza interfluve. The areas that previously were characterized as increasing water level zones become smaller in size and the expected water level increase was less significant.      This indicates significant water level drops (up to about -6 m with even higher drops at some places. The latter may also be attributed to less confidence in the simulated water levels at higher altitudes) under highland areas such as the Alpokalja, Eastern part of Transdanubian Mountains, Northern Mountain Ranges, and the Mecsek. A less significant decrease of groundwater levels was simulated along the southern reaches of the Tisza River (<−1 m) and the eastern parts of the Great Hungarian Plain (<−2 m). Most parts of the Transdanubian, the Kisalföld, the western parts of the Transdanubian Mountains, and the Duna-Tisza interfluve is characterized by increasing water levels, which exceed even 2 m at some locations.
The water table difference map ( Figure 18) between the reference period (1961-1990) and the 2071-2100 averages indicates that the expected water level drops become slightly more significant at the same locations. Furthermore, the areas characterized by lowering water table elevations increase in size including some areas previously not being negatively affected such as the Duna-Tisza interfluve. The areas that previously were characterized as increasing water level zones become smaller in size and the expected water level increase was less significant.

Discussion
Several studies have been published about groundwater and climate change predicting changes in the water balance. Each of these studies apply climate data for their predictions, although different approaches were used for deriving climate data series [42]. Loaiciga et al. [43] used global averages. Allen et al. [44], Brouyere et al. [45], Vaccoro [46] applied regional projections for assessing climate change impacts. Scibek et al. [47] Jyrkama and Sykes [48], Serrat-Capdevila et al. [49], Toews and Allen [50] applied downscaled climate data in their investigations. Van Roosmalen et al. [51,52], Rosenberg at al. [53], York et al. [54] applied regional climate model projections in their studies.
The majority of early modelling studies employed stochastically generated daily weather series, and applied climate change shifts for predictive modelling. Many studies applied a range of global circulation models (GCMs) or used an averaged projection calculated from various GCMs. There are only a few studies that investigated the application of different downscaling methods.
Vaccoro [46] applied two different climate scenarios to investigate recharge sensitivity under various land use scenarios. This author also investigated the changes in recharge in response to the climate variability based on historical data. Rosenberg et al. [53] applied a Hydrologic Unit Model, imposed three GCMs and simulated the changes in water yields and groundwater recharge. Loaiciga et al. [43] used historical climate data from periods of draught, average and high precipitation. These authors scaled the historical time series to create climate change scenarios. This paper also evaluated the sensitivity of water resources at the aquifer scale through the combination of groundwater pumping scenarios with climate change scenarios. York et al. [54] utilized a coupled model of land and atmosphere and investigated climate change impacts on a watershed in Kansas, USA. Yu et al. [55] developed a novel methodology for interactive coupling of hydrologic models with climate models. Allen et al. [44] used shifted climate data using a stochastic weather

Discussion
Several studies have been published about groundwater and climate change predicting changes in the water balance. Each of these studies apply climate data for their predictions, although different approaches were used for deriving climate data series [42]. Loaiciga et al. [43] used global averages. Allen et al. [44], Brouyere et al. [45], Vaccoro [46] applied regional projections for assessing climate change impacts. Scibek et al. [47] Jyrkama and Sykes [48], Serrat-Capdevila et al. [49], Toews and Allen [50] applied downscaled climate data in their investigations. Van Roosmalen et al. [51,52], Rosenberg at al. [53], York et al. [54] applied regional climate model projections in their studies.
The majority of early modelling studies employed stochastically generated daily weather series, and applied climate change shifts for predictive modelling. Many studies applied a range of global circulation models (GCMs) or used an averaged projection calculated from various GCMs. There are only a few studies that investigated the application of different downscaling methods.
Vaccoro [46] applied two different climate scenarios to investigate recharge sensitivity under various land use scenarios. This author also investigated the changes in recharge in response to the climate variability based on historical data. Rosenberg et al. [53] applied a Hydrologic Unit Model, imposed three GCMs and simulated the changes in water yields and groundwater recharge. Loaiciga et al. [43] used historical climate data from periods of draught, average and high precipitation. These authors scaled the historical time series to create climate change scenarios. This paper also evaluated the sensitivity of water resources at the aquifer scale through the combination of groundwater pumping scenarios with climate change scenarios. York et al. [54] utilized a coupled model of land and atmosphere and investigated climate change impacts on a watershed in Kansas, USA. Yu et al. [55] developed a novel methodology for interactive coupling of hydrologic models with climate models. Allen et al. [44] used shifted climate data using a stochastic weather generator. The differences were calculated based on three different GCMs and represented extreme climate conditions. This method was exploited in British Columbia, Canada. Brouyere et al. [45] extracted monthly precipitation and temperature from three GCMs in the Geer basin, Belgium. These authors calculated three future model scenarios by combining daily data with monthly change rates. These climate time series were used as inputs for an integrated model. Scibek et al. [47] used the estimates of aquifer recharge and runoff for the calculation of river discharges at the basin-scale. Jyrkama & Sykes [48] constructed various climate scenarios based on historical reference to model climate change impacts in southern Ontario, Canada. These authors calculated groundwater recharge using a distributive model. Serrat-Capdivela et al. [49] used downscaled results derived from an ensemble of GCMs and different IPCC scenarios. These authors determined spatially distributed recharge estimates from climate models. They applied a transient 3D integrated model to evaluate the hydrological conditions in the transboundary San Pedro Basin (US/Mexico). Van Roosmalen [51] used outputs from a regional climate model for two thirty-year periods and ran distributed hydrological model to simulate groundwater conditions in Denmark. Mileham et al. [56,57] harnessed a soil-water balance model to calculate water budget for several scenarios in Uganda. Rivald et al. [58] applied Regional Climate Model results as input for a catchment hydrological model (CATHY) for a catchment in Nova Scotia, Canada. Toews & Allen [50] investigated the sensitivity of recharge to climate projections of three different GCMs. These authors employed groundwater flow models for the simulation of groundwater levels in the Oliver Region in British Columbia. Tietjen et al. [59] exploited a model to investigate the dynamics of soil water under various climate conditions. Kidmose et al. [60] undertook a regional and local scale modelling study in Denmark while using the MIKE SHE integrated model environment. They applied extreme value analysis and used the ensemble approach. Raposo et al. [61] used the SWAT model to assess climate impact in Galicia-Costa, Spain. These authors utilized climate projections of various GCMs and evaluated two regional climate model scenarios. The study concluded that climate change has a significant influence on the temporal variability of recharge. Mollema & Antonelli [62] studied the seasonal variation of recharge to coastal aquifers. Potential evapotranspiration was estimated with the Thornthwaite method. The SEAWAT model was applied for simulating the influence of recharge on the size of freshwater lenses. Kumar [63] described a methodology for climate impact assessment including the downscaling of GCMs, the use of weather generators, simulation of recharge using rainfall-runoff models and calculation of groundwater conditions using MODFLOW.
Several additional studies have been published using approaches similar to those described above. Although the investigation techniques are improving, significant uncertainties are involved in climate predictions used as input parameters in climate impact assessment studies.
The simulation results achieved in this study were compared with climate impact studies undertaken in neighbouring countries. In a comprehensive study [64] the availability of water resources were analysed in 20 test areas across South East Europe (SEE region). Average climate conditions between 2021-2050 and 2071-2100 were analysed, based on regional climate model projections. Sensitivity of the different types of water resources to climate change was characterised based on the reactions of the hydrological system to changes of relevant meteorological parameters. Regional differences in the changes of water resources (recharge) were concluded to reflect the changes in precipitation and showed a clear geographical pattern. Most sites showed decreasing recharge and diminishing water resources. The test areas in the north and west showed significantly lower changes than those in the south and south-east of the SEE region. Decreases in recharge in the north and west are less than 10% until 2050 and do not exceed 20% until the end of the century. In the south and south-east of the SEE region the decrease is general and amounts to 10-40% in the first half of the century and 20-60% in the second half.

Conclusions
The present paper summarises a methodology developed for the calculation of groundwater table distributions from climate parameters. The goal of water table modelling was to develop a methodology that can be applied for the calculation of the water table under different climate conditions. This was done in order to facilitate climate impact assessment and the evaluation of climate sensitivity of groundwater aquifers.
We developed a modular approach to the quantitative investigation of the influence of climate conditions on groundwater table distribution: 1.
Climate zones were determined based on measured and simulated climate variables; 2.
Recharge zones were delineated based on slope, landuse and surface geology; 3.
Recharge was calculated for recharge zones making use of 1D analytical hydrological models; 4.
Groundwater table distribution was simulated using numerical groundwater models.
The interpolated daily climate datasets of the CARPATCLIM European database and the results of the ALADIN regional climate model were applied as input parameters in our modelling study. Climate zones were delineated quantitatively using the Thorntwaite methodology. Recharge zones (HRU's) were delineated based on surface geology, landuse, and slope maps. We applied the HELP hydrological model for the calculation of recharge. The water table was calculated numerically using the MODFLOW numerical groundwater modelling code. The applied methodology provided a quantitative link between climate conditions and shallow groundwater conditions, in addition, it was successfully used for water table modelling.
The results of recharge calculation indicate that recharge decreased by up to 50 mm/year throughout the second half of the 20th century in mountainous areas such as the Alpokalja area, the Eastern section of the Transdanubian Mountains, and the Northern Mountain Ranges. Decreasing recharge rates of similar degree are predicted until 2100 in these areas and to a lower degree (up to 20 mm/y) in the Duna-Tisza Interfluve and Mecsek Mountain area (Figures 1 and 10). Slightly increasing recharge rates (up to 20 mm/y) are predicted in the western part of the Transdanubian Mountains, and the Transdanubian Hills. Increasing recharge is a consequence of increasing annual rainfall in these areas, and also of the changing temporal pattern of rainfall.
Model predictions indicate decreasing water levels in highland areas such as the Alpokalja, Eastern Transdanubian Mountains, Mecsek, and Northern Mountain Ranges areas. The simulated water table drop can exceed 10 metres in certain locations, but it is generally less than 1 metre. Calculations based on climate model projections also indicate moderate water level drops in the Western part of the Tiszántúl area and the Duna-Tisza Interfluve, and the Szigetköz areas, which are generally less than 1 metre. Slightly rising groundwater levels (generally less than 3 metres) are predicted in the Transdanubian Hills, the western parts of the Transdanubian Mountains, and the Szamos Valley areas.
Climate change is predicted to have long-term influence on shallow groundwater resources. In large areas of the country, it is predicted to cause a moderate water table decline. This might impact surface waters, irrigation channels, and groundwater dependent ecosystems such as wetlands and riparian vegetation. Groundwater depletion also affects shallow groundwater uses such as dug wells. Groundwater level drops affect soil moisture conditions and can have adverse effects on natural vegetation and plant production. Furthermore, the significant decrease of annual rainfall over the eastern part of the country (up to 100 mm/y) is believed to cause additional indirect effects such as increased irrigational groundwater demand. This might cause additional pressure on groundwater resources and draws the attention to careful water resource management. In the western part of Transdanubia the projected increase of annual rainfall might facilitate a more effective use of groundwater resources.
Predictive simulations were based on climate simulations following the SRES A1B CO 2 emission scenario. Future emission rates different from those assumed by the A1B scenario might entail different recharge rates and groundwater conditions in the future.
The presented outputs were determined at the regional scale, and thus cannot be used for local-scale investigations. The presented methodology can be applied for modelling climate impacts at various scales for the assessment of climate vulnerability of groundwater resources.
Author Contributions: A.K. developed the methodology, performed analysis and modelling and prepared article. A.J. performed modelling and contributed to the article. All authors have read and agreed to the published version of the manuscript.

Funding:
The research reported in this paper has been supported by the National Research Development and Innovation Office of Hungary (Grant No. TKP2020 BME-IKA-VIZ).
Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Not applicable.

Conflicts of Interest:
The authors declare no conflict of interest.

Disclaimer:
Model results reflect the inaccuracies and uncertainties of the input databases. The methodology presented in this study makes the recalculation of model outputs possible making use of updated or completed datasets. The presented simulation results are representative at the country scale, and thus cannot be used for local-scale investigations. In order to increase the accuracy of the results, the resolution of the input datasets and applied models need to be increased.