How Complex Groundwater Flow Systems Respond to Climate Change Induced Recharge Reduction?

Abstract

Our recent knowledge about the role of different fluid driving forces in the response of groundwater flow systems to climate change is still limited. This current study aimed to evaluate possible spatial and temporal changes in complex, gravity- and overpressure-driven groundwater flow systems induced by climate change and look for general trends and characteristics of the flow field using 2D transient groundwater flow simulations. Results showed significant large-scale changes in the transient subsurface flow field and flow dynamics due to recharge reduction. Local gravity-driven flow systems are the most vulnerable to atmospheric processes, while overpressured regimes are expected to be independent of direct climatic variability. By the involvement of different degrees of overpressure, it was revealed that, as the degree of overpressure increases, the penetration depth of the topography-driven local flow systems decreases. The higher the overpressure, the lower the climate change-induced groundwater level decrease over time, suggesting the buffering effect of overpressure as a fluid driving force in the flow systems’ response to the changes in hydrologic parameters. The main novelty of the study is the involvement of different fluid driving forces in the evaluation with the combination of a regional scale investigation, which is unique in the context of climate change effects on groundwater systems.

1. Introduction and Objectives

Groundwater reserves, being a crucial water resource, have great importance in sustainable water management, which requires an effective prediction of future availability. Simultaneous changes in key hydrologic parameters, predicted by the reports of the Intergovernmental Panel on Climate Change (IPCC), can significantly modify groundwater replenishment [1,2]. Moreover, the relation between climate parameters and groundwater is considered more complicated than with surface water [3], partly because of the residence time of groundwater, which varies on a wide scale and is likely to delay, attenuate and disperse the effects of climate change [4,5,6,7].

Climate variability will likely have diverse effects on subsurface recharge rates and mechanisms [2,6,8,9,10], as recharge is a sensitive function of the hydrological factors, local geological conditions, topography and land use. While many studies have predicted reduced subsurface replenishment [11,12], the effects of climate change on recharge may not necessarily be negative in all aquifers [13,14,15].

As scientific interest regarding the influences of climate change on groundwater systems has significantly increased over the last decade, it has been confirmed that both changes (i.e., decrease and increase in recharge) can have widespread consequences on subsurface conditions and groundwater systems. Several investigations have been completed to clarify possible future groundwater level changes [16,17], modification of groundwater recharge and storage [18,19], groundwater quality changes [20], alteration in subsurface temperature conditions [21], and in surface water–groundwater interaction [22,23]. However, the majority of the prior studies focus mainly on the hydrological aspects of climate change and the groundwater systems, while studies which deal with possible spatial and temporal changes in groundwater flow pattern and flow dynamics, i.e., in the fragmentation and hierarchy of nested flow systems, are still under-represented. Additionally, our recent knowledge about the role of the different fluid driving forces in the response of groundwater flow systems to climate change is still limited. A particularly interesting situation occurs when the groundwater flow pattern of a given area is determined by various fluid driving forces, because their exposure to changes in hydrologic parameters may be quite different. Gravity-driven flow is induced by elevation differences in the water table [24], owing to the most conspicuous consequences of alteration in atmospheric processes, as it is recharged directly by the infiltrating precipitation. Beyond topography, water movement in the subsurface can also be driven by stress-related mechanisms (such as sedimentary loading, tectonic loading and stresses, erosional unloading) and buoyancy (arising from density gradients due to changes in salinity and temperature) [25,26], which usually have weak and highly delayed, or do not have a direct connection with the surficial hydrologic processes.

This study aimed to contribute to the understanding of the large-scale behaviour of such complex subsurface processes in response to climate change-induced recharge reductions by looking for general trends and characteristics of the flow field using 2D transient groundwater flow simulations.

The Duna–Tisza Interfluve in Hungary was assigned as a sample area, where two separated flow domains were identified and deeply studied [27]. The upper gravity-driven flow regime is recharged directly from the infiltrating precipitation, while the underlying overpressured domain is maintained by pore volume reduction due to compaction and tectonic compression of the basement [28,29]. Such flow regimes do not have meteoric water recharge, as the super-hydrostatic vertical pressure component prevents vertical leakage from above [30].

In the study, special emphasis was put on (i) how the relative role of the two fluid flow systems (i.e., gravity-driven and overpressured domain) may change, (ii) how the hierarchy and fragmentation of the ground flow field may alter by time, (iii) how the penetration depth of the upper, gravity-driven flow field may adjust to these changes and (iv) how groundwater-related shallow surface water bodies will be affected by these changes.

The novelty of the study is the regional scale of the investigation, which is subordinate in the literature compared to the lower scale, site-specific predictions. It requires a different approach, and in parallel, combined with a relatively long-term investigation (several decades), it could provide more extensive spatial and temporal information on the effects of climate change on groundwater flow patterns and dynamics, which is still underrepresented in the previous studies. Additionally, the involvement of different fluid driving forces in the evaluation is unique in the context of climate change effects on groundwater flow systems.

2. Sample Area

The Duna-Tisza Interfluve, located in Hungary, is characterized by flat-bottomed valleys at the River Duna and Tisza (80–95 m asl), and by a N-S orientated central ridge (100–130 m asl) along the divide (Figure 1).

The eastward-dipping basement of the broader area is tectonically strongly disturbed and characterised by horst and graben topography. It is covered by Neogene semi- to unconsolidated marine, deltaic, lacustrine, fluviatile and eolian clastic sediments with a thickness from 600 m at the Duna River to over 4000 m at the Tisza River [31]. Their hydraulic properties were determined by qualitative estimation derived from lithological and facies conditions [27,29].

Hydrostratigraphically, the rock framework is divided into an upper and a lower sequence by an extensive aquitard (Algyő Aquitard) of an average hydraulic conductivity of K = 10⁻⁸ ms⁻¹. The aquitard is overlain by a highly permeable aquifer (Great-Plain Aquifer) including two geological sequences, which can be characterised together by an average hydraulic conductivity of K = 10⁻⁵ ms⁻¹ [27]. It mainly consists of fluvial, eolian and lacustrian formations interrupted by local low-permeable eolian and lacustrine ones. The lower series of units consist of Neogene siltstones, sandstones and clays (Szolnok Aquifer and Endrőd Aquitard with hydraulic conductivities varying between 10⁻⁹ and 10⁻⁶ ms⁻¹) underlain by the Pre-Pannonian Aquifer (with an estimated hydraulic conductivity of 10⁻⁶ ms⁻¹), and Pre-Neogene formations consist of carbonates, metamorphic formations and flysch with unknown hydraulic properties [32] (Figure 1).

Based on detailed hydraulic studies by previous authors, the principal features of subsurface flow conditions of the Duna-Tisza Interfluve are relatively well-known. Two different but hydraulically connected, superimposed groundwater flow domains are acting: an upper gravity-driven meteoric fresh water, which can be characterised by near-hydrostatic pressure conditions and an underlying overpressured domain of water with high total dissolved solid (TDS) content [27,29,33].

The groundwater flow in the upper aquifer is controlled by the topography, causing a downward flow beneath the highlands and an upward flow beneath the valleys. The underlying overpressured system is characterized by an upward vertical flow component, which is maintained by pore volume reduction due to compaction and tectonic compression of the basement [27,28,29,34,35]. The upward overpressured flow is forced and deflected by descending fresh meteoric waters of the central regions toward the river valleys, where they are discharging. The distribution patterns of freshwater and saline vegetation and the location of salt-affected soils, as well as the hydrochemical differences of fresh and saltwater lakes at relatively short distances from each other, reflect the different properties of the vertically superimposed gravity-driven and overpressured flow systems [27,36,37]. Lakes of the Interfluve are situated in really extensive shallow (<1 m deep) depressions and are not affected by recharging or discharging stream water.

Hydraulic communication between the upper gravity-driven and underlying overpressured systems occurs by diffusion and seepage across the aquitard through permeable, discrete structural and sedimentological discontinuities [27,29,36,38].

As the upper, gravity-driven groundwater flow system is recharged directly from the infiltrating precipitation, this flow regime is regionally renewable on a human and on a geological timescale as well [27,29]. In contrast, the deep flow regime, being a hydraulically closed system, is not replenished by precipitation, because vertical leakage from above is restricted by the excess (above hydrostatic) vertical pressure [29,30]. Consequently, this pressure regime is characterized by fluids which are not renewable [30]; the amount of water which can be released depends on the elastic properties of water and rock matrix, i.e., the regional storage capacity of the system [39,40].

Predicted Climate Change and Consequences on Groundwater Recharge

The Duna-Tisza Interfluve is characterised by a moderately continental climate with approx. 600 mm annual precipitation and a 10 °C average annual temperature. Based on regional climate predictions for Hungary, the annual precipitation is likely to decrease by about 20% by the end of the century [41], while the annual temperature is likely to increase by 4–5.4 °C [42] compared to the reference period of 1961–1990.

Focusing on the large-scale alterations of groundwater flow pattern and dynamics rather than site-specific absolute predictions of future groundwater levels, the most deterministic hydrologic parameters, i.e., precipitation and evapotranspiration (described by dynamic factor analysis of Király [43] were taken into account during the evaluation of groundwater recharge. Overland flow and through-flow were neglected, considering that rainfall rate rarely exceeds infiltration capacity under recent climatic conditions of the area; moreover, the runoff would also decrease in time due to decreasing precipitation.

In the case of the Interfluve, the Duna and Tisza rivers represent a regional hydraulic potential minimum of the area. Due to the groundwater flow characteristics of the sample area, the effects of the rivers on groundwater level can be detected mainly in their vicinity [27], and therefore, surface runoff was not taken into consideration during recharge calculation.

It is assumed that uncertainties generated by this site-specific simplification are not significant compared to the uncertainties inherent in the climate change predictions and the numerical investigation can, nevertheless, help to gain insight into the large-scale response of the different groundwater flow systems to the predicted climate change.

Evapotranspiration was calculated based on the simplified Turc formula [44], which requires only limited data and is simple to implement, and given as:

$$E=\frac{P}{\sqrt{0.9+\frac{P^2}{K^2}}}$$

(1)

where E is evapotranspiration (mm yr⁻¹), P is annual precipitation (mm yr⁻¹), K = 300 + 25T + 0.05T³ and T is the mean annual air temperature (°C).

Assuming P = 600 mm yr⁻¹ and T = 10 °C, which are characteristic of the area, the current evapotranspiration (E) is 435 mm yr⁻¹. Subtracting E from precipitation (P), the current groundwater recharge is 165 mm yr⁻¹. Therefore, under current conditions, only approx. 27% of precipitation can reach the saturated zone and recharge the groundwater system.

In order to easily follow the recharge reduction-induced changes in the subsurface flow processes, a pessimistic scenario was taken into account with a 20% decrease in annual precipitation and a 4 °C increase in mean annual temperature, which can lead to a significant, >60% decrease in recharge by the end of the century (2100) in the case of the sample area.

3. Numerical Simulation Method

In the current study, future modifications of groundwater flow system characteristics due to climate change were simulated using the transient numerical flow module within the Heatflow-Smoker code [45] along a 2D vertical plane. Heatflow-Smoker is a finite element code, which couples density-dependent groundwater flow and heat transport.

The model solves the continuity equation for flow [46] which can be expressed as:

$$\frac{\partial }{\partial x_i}\left[K_{i,j}\left(T\right)\left(\frac{\partial \Psi }{\partial x_j}+{\rho }_r\left(T\right)\times \overline{n_j}\right)\right]=S_s\frac{\partial \Psi }{\partial t}$$

(2)

where x_ij are the 3D spatial coordinates (m), K_ij(T) is the temperature-dependent hydraulic conductivity tensor (m s⁻¹), ψ is the equivalent freshwater head (m), ρ_r(T) is the temperature-dependent relative density of water (-), S_s is the specific storage (m⁻¹) and t is time (s).

Changes in flow patterns were followed along a heterogeneous, anisotropic, two-dimensional vertical section across the Interfluve (Figure 1 and Figure 2). The 100 km long and 1500 m deep model domain was subdivided into a regular network of 1000 × 15 elements in the x and z directions, with each rectangular prism element measuring 100 m × 100 m, respectively (with 1 element in the transverse dimension).

As a first step, a steady-state simulation was run based on site-specific hydraulic and hydrostratigraphic data, representing the current groundwater flow system. The full transient simulation was run over a period of 80 years, with a time step of 1 yr. Changes to the hydrologic parameters were based on referenced climate change predictions.

Since our interest was to highlight the possible influence of climate change on complex groundwater flow systems—driven by gravity and overpressure—and look for general characteristics of the transient flow field, some simplifications were included in the simulation approach, which was not expected to be significant in terms of the objectives of the study.

3.1. Model Dimensionality

It is obvious that representing a 3D real system as a two-dimensional cross-section overlooks some system dynamics. However, depending on the system characteristics, the scale of the investigation and the questions to be answered, in some cases a two-dimensional approach could be a suitable choice. For low-relief settings (i.e., with a flank slope < 10°), where the ratio of recharge (R) to hydraulic conductivity (K) of the uppermost subsurface layer is <0.15, 2D cross-sectional models may be justifiable since transverse flow (i.e., perpendicular to the primary regional topographic gradient) would be generally less than 10% of the total flow [47]. As in the case of the sample site, the R/K ratio is below the critical value; furthermore, the general characteristics of the dynamically changing flow regimes are the focus of the investigation, and a 2D model is considered a reasonable approach for providing the needed insight into the response of the system to the predicted changes of the key hydrologic parameters.

3.2. Driving Forces on Groundwater Flow

The current numerical investigation focuses on the influences and temporal changes of the water table gradient (i.e., gravity) and sediment compaction-induced overpressure was the focus of the investigation. In the case of the Duna-Tisza Interfluve sample area, possible effects of density-driven flow due to high total dissolved solids (TDS) concentrations and temperatures on the flow field were analysed by previous authors. It was revealed, that (i) due to the 45–50 °C km⁻¹ geothermal gradient of the area, heat can be taken into account as a driving force, but it has only a small modifying effect on the regional groundwater flow pattern of the area [48]; (ii) TDS of groundwater remains below 40 g L⁻¹ in the case of the Interfluve [27,38]; and (iii) the effects of a high geothermal gradient and a relatively low-dissolved solid content of the water can equalize each other [49]. Therefore, the influence of a buoyancy-driven force was not considered during the current study and an isothermal groundwater flow was simulated, assuming a uniform temperature of 10 °C and water density of 1000 kg m⁻³.

3.3. Representation of Overpressure

In the case of gravity-driven flow systems, subsurface hydraulic potential distribution is clearly determined by the groundwater table and flow is generated by its gradients [24]. Therefore, hydraulic head (h) values in a given basin are limited by the range of the groundwater table elevation (z_wt,min and z_wt,max, by noting that in the case of groundwater table, its elevation and hydraulic head is equal, z_wt = h_wt) and spatial distribution pattern of fluid-dynamic parameters representing normal, i.e., near-hydrostatic pressure conditions. Any alteration in subsurface pressure conditions is due to the gravitational flow of groundwater [39], i.e., pore pressures are adjusted to their hydrogeologic environment (equilibrium-type overpressure, [50]).

In general, abnormal fluid pressure is generated by transient processes, such as stress changes caused by tectonic compression and/or disequilibrium compaction; changes in fluid volume due to thermal expansion, diagenetic and metamorphic processes as well as hydrocarbon generation; or fluid movement due to density differences caused by buoyancy forces and osmosis [50,51,52,53,54,55]. In many sedimentary basins, stress-related mechanisms are the most likely causes of overpressure (OP), which is abnormally high pressure in the subsurface that exceeds the hydrostatic pressure at a certain depth.

Overpressures are created whenever the rate of pressure generation exceeds the rate of pressure dissipation. The rate of fluid-pressure dissipation across a porous medium is inversely proportional to the hydraulic conductivity of that medium. Therefore, low-conductivity regions are essential for the generation and maintenance of overpressures [50], as they could function effectively in units where sufficiently low-permeability formations are present, which inhibits the complete dewatering of the underlying formation and contributes to the excess pressure of the pore fluid [52].

According to the definition by Mádl-Szőnyi [56], overpressure conditions prevail if the deviation of the measured dynamic pore pressure (p_dyn) from the static pore pressure (p_st), characteristics for the given elevation, exceeds the limit of near-hydrostatic pressure conditions (∆p_hydrostatic). It can be given as:

$$p_{dyn}>p_{st}\pm \Delta p_{hydrostatic}=p_{st}\pm \left(\rho g\frac{z_{wt,max}-z_{wt,min}}{2}\right)\times {10}^{-6}$$

(3)

where ρ is the density of water (1000 kg m⁻³), g is the acceleration due to gravity (9.81 ms⁻²) and z_wt,max and z_wt,min are the maximum and minimum elevations of the groundwater table [m].

Based on Equation (3), abnormal pressure can be determined simply by the hydraulic head interval of groundwater level of the given basin, as h > h_wt,max represents overpressure, while h < h_wt,min is the characteristic for underpressure conditions.

In a recent study, an overpressured system is characterized by an upward vertical flow, which is maintained by pore volume reduction due to compaction and tectonic compression of the basement of the sample site [29]. Therefore, overpressure is represented by an increased hydraulic head along the bottom of the system, induced by groundwater inflow across the bottom of the studied Section 3.4. Input parameters, initial and boundary conditions of the model.

Considering the main aims of the study, i.e., the evaluation of the large-scale alterations of groundwater flow pattern and dynamics rather than a site-specific absolute prediction of future groundwater levels, a simplified hydrostratigraphic configuration was applied, focusing on the key elements of site-specific hydrogeological characteristics of the sample site. Therefore, a non-uniform thick aquitard and an overlying aquifer were tested. In order to provide a better resolution within the uppermost subsurface zone, where main changes are expected due to the recharge reduction-induced modification in the groundwater table; the aquifer was differentiated into three distinct sub-units (aquifer 1, 2 and 3, Figure 2b) in accordance with the hydrostratigraphic settings of the sample area. However, in the following interpretation, it is referred to uniformly as “aquifer”. Hydraulic parameters of the formations were defined based on the previous detailed hydrogeological evaluation of the Duna-Tisza Interfluve [27,29,36] (Figure 2b).

The numerical simulation was performed within the saturated zone, for which the groundwater table forms the upper surface, represented by the recent water table configuration of the sample site (Figure 2a). Across the upper flow boundary, outside of the river valleys and local lake areas with constant head, a specific recharge was defined, based on the recent hydrologic parameters for the steady-state simulation, as well as calculated by the model taking into consideration a linear decrease in recharge in the case of the transient model. Water inflow across the lateral boundaries represents groundwater movement toward the regional hydraulic potential minimum of the river valleys within the upper aquifer.

Upward overpressured flow from the underlying formations was represented by groundwater inflow across the lower flow boundary of a studied cross-section with different Darcy flux values in the case of the different simulated cases with a different degree of consequent hydraulic head increase. In the case of comparative simulations neglecting the effect of overpressure, the lower flow boundary was handled as no-flow, i.e., Darcy flux of q = 0 ms⁻¹ was applied (Figure 2b).

In the case of the transient simulation, the water table configuration obtained from the steady-state simulation was used as the initial condition for the transient flow model. Fixed heads assigned to represent the steady-state lake levels were eliminated; thus, the water table elevation change below the lakes was calculated by the model based on the predicted decrease in recharge across the top boundary over each time step. Comparisons of the simulated transient groundwater levels with local lakebed bottom elevations were used to follow the possible changes in surface water levels.

Rivers along the studied section represent the regional hydraulic potential minimum. As the Danube represents the second largest river in Europe, while the Tisza is its largest sub-catchment, their response to predicted climate change is difficult to predict and beyond the scope of this study. Effects of the rivers on groundwater levels in the case of the sample area can be detected mainly in their vicinity [27]. Moreover, in the case of the studied system with an elevated area between the main rivers, the relative hydraulic position of the rivers (i.e., minimum hydraulic potential) would remain in spite of any changes in river water level. Therefore, changes in rivers’ water level and streamflow under climate change were not considered. A sensitivity analysis of the effects of river water level reductions on the flow field showed no significant effects on the general conclusions of the study.

During the transient simulation, uniform specific storage of S_s = 0.001 m⁻¹ was applied, which represents a combination of relatively high storativity in shallow unconfined sediments (governed by porosity) and relatively low storativity in the deeper rock framework.

3.4. Simulated Cases

As an initial case, only gravity as a driving force was taken into account during the simulation, neglecting the effects of overpressured flow. With the help of these models, the effects of climate change-induced recharge reduction on a purely gravity-driven flow system can be analysed, and then its flow pattern and dynamic characteristics are easily comparable with the ones created by complex, gravity and overpressure-driven flow systems.

Different degrees of overpressure were also taken into account, in order to evaluate its relative role compared to gravity in canalising groundwater and also, its behaviour under the changing hydrological characteristics of the system. Considering the recent water level configuration of the sample area with h_wt,min = 83 m and h_wt,max = 110.4 m, abnormal pressure can be evaluated based on Equation 3, which denotes overpressure conditions at the bottom of the system (at 1500 m depth) with p_dyn > 15.6 MPa. Beyond a medium overpressure of 25 MPa at the bottom, characteristics for the sample area, effects of a lower (20 MPa) and a synthetic case with extreme high (65 MPa) overpressure were tested for comparison by applying a corresponding groundwater inflow across the bottom of q_in(bottom) = 5 × 10⁻¹¹ ms⁻¹ (lower OP), 10⁻¹⁰ ms⁻¹ (medium OP) and 5 × 10⁻¹⁰ ms⁻¹ (higher OP).

The numerical investigation aimed to look for general trends and characteristics of the transient flow field and focused on (i) how the relative role of the two fluid flow systems may change due to climate change-induced recharge reduction, (ii) how the hierarchy and fragmentation of the flow field may alter and (iii) how the penetration depth of upper, gravity-driven flow field may adjust to these changes. Due to the scale of the investigation, local predictions of absolute surface water level changes are beyond the scope of the study. However, as surface waters are an integral part of the groundwater flow systems, comparisons of the simulated transient groundwater levels with shallow (<1 m deep) local lakebed elevations can be used to follow the possible changes in lake water levels during the transient simulations and predict their future state: (i) will they likely be maintained or will they eventually be drying, and (ii) will they be affected by the possible changes of their interconnected groundwater systems?

4. Results

4.1. Pure Gravity-Driven Groundwater Flow Systems in a Changing Climate

Considering only the effect of gravity as the driving force, groundwater flows from the elevated ridge to the river valleys in accordance with the hydraulic head differences and hydraulic conductivities. The flow field is characterised by hierarchically nested flow with local (L1, L2, L3, L4), intermediate (I1) and regional (R1, R2) flow systems, determined by the water table (Figure 3). Groundwater recharge takes place along the elevated ridge by the infiltrating precipitation, while discharge occurs dominantly along the river valleys within relatively broad zones. The developed flow field is hydraulically separated by a subsurface water divide below the main recharge area of the ridge (at approx. 51 km distance, Figure 3b). Gravity-driven flow is settled into a more dynamic upper flow field in the aquifer with a maximum flow velocity of 2 × 10⁻⁶ ms⁻¹, and a much slower system in the aquitard with a maximum Darcy flux of 10⁻¹⁰ ms⁻¹. It is a major part of the infiltrating water moving in the upper aquifer arranged into local or intermediate systems (L1, L2, L3, L4, I1, Figure 3b), reaching a maximum penetration depth of 1300 m.

As results of the transient simulation showed, recharge reduction (from 165 mm yr⁻¹ to 54 mm yr⁻¹, calculated based on the predicted climate change) led to large-scale changes in the gravity-driven flow field. Water level change was the highest at the watershed below the recharge area, where it decreased by at most 3.6 m after 80 yr, leading to the extent of the hydraulic potential drop within the aquifer. However, due to the significantly lower diffusivity of the aquitard unit (D, ratio of hydraulic conductivity K to specific storage S_s), its hydraulic head distribution was not affected by the recharge reduction. Consequently, a hydraulic potential maximum (h_max = 109.8 m) developed within the aquitard below the subsurface water divide of the system (at approx. 51 km distance and 700 m depth), which is 3 m higher compared to the maximum hydraulic head of 106.8 m along the water table, i.e., transient processes induced by recharge reduction led to the development of a slight overpressure within a gravity-driven flow system.

This situation created an upward flow domain within the lower part of the aquifer in accordance with the hydraulic head distribution. With these changes, direct surface replenishment of the regional flow systems (R1, R1′, R2 and R2′, Figure 3c) from infiltrating precipitation ceased. In parallel, flattening of the water table led to a decreased fragmentation in the uppermost flow domain (from L1, L2, L3, L4 local and I1 intermediate flow systems to L1 and L2 local ones), and to a significant decrease in the penetration depth of the infiltrating water from 1300 m (characteristic for the initial flow field, Figure 3b) to 700 m (Figure 3c).

Simulated groundwater level changes suggest that the water level of lake 1, which takes place at the lower flank of the ridge, will not be significantly affected by the recharge reduction after 80 yr in the case of a pure gravity-driven system. However, its feeding flow system would be changed over time from a deeper, intermediate system (I1 with a penetration depth of 600 m, Figure 3b) to a shallower, local one (L1 with a penetration depth of 300 m, Figure 3c). Another clear result of the transient simulation was the continuous decrease in groundwater level below lake 2, which dropped 2.1 m after 80 yr, i.e., its lakebed was dried.

4.2. Effect of Decreasing Recharge on a Combined Flow System with Gravity-Driven and Overpressured Flow Domains

After the evaluation of the effects of recharge reduction on a pure gravity-driven flow system, the effect of an underlying overpressured flow regime was involved in the investigation, and the response of the different flow systems to climate change was analysed and compared to each other.

The flow field of the combined flow system consists of an upper, hierarchically nested gravity-driven flow with local and intermediate flow systems (L1, L2, L3, I1, Figure 4a) reaching a maximum penetration depth of 1000 m (L3), which is vertically superimposed with a domain of an upward overpressured flow. Latter is developed by a uniform groundwater inflow of 10⁻¹⁰ ms⁻¹ across the bottom of the system, leading to a significant hydraulic head increment with h_max = 1094.4 m (corresponding to p_dyn of 25 MPa). Consequently, flow pattern and dynamics significantly changed in the aquitard compared to the pure gravity-driven case; the flow direction reversed and in parallel, the maximum Darcy flux increased by two orders of magnitude (from 10⁻¹⁰ ms⁻¹ to 1.1 × 10⁻⁸ ms⁻¹). However, since groundwater table elevation, representing recent conditions, is the same as it was in the case of the pure gravity-driven flow system analysis, the main characteristics of the current steady-state flow models, including neglecting and considering the effect of the underlying overpressured flow, are very similar within the aquifer. Its flow field is hydraulically separated by the subsurface water divided below the recharge area (at approx. 51 km, Figure 4b). Groundwater in the aquifer flows toward the potential minimums represented by the river valleys with a maximum fluid flux of 2.1 × 10⁻⁶ ms⁻¹, owing to similar recharge and discharge characteristics compared to the pure gravity-driven case. One remarkable difference appears below lake 1, which is fed by the overpressured upward flow in this case (Figure 4b); on the contrary, it is recharged directly by the infiltrating precipitation, transported by an intermediate flow system in the pure gravity-driven case (Figure 3b).

As a result of the recharge reduction, the groundwater table dropped by 3.1 m after 80 yr, which is 0.5 m lower compared to the water level decrease in the pure gravity-driven case. In parallel, fragmentation of the upper gravity-driven flow system decreased to single local cells with a maximum penetration depth of 500 m (L1, L2, Figure 4c). Flow direction and dynamics within the aquitard were not affected by these changes. However, the boundary between the overpressured and the superimposed gravity-driven flow domain was shifted from the strata boundary (recent case, Figure 4b) to a shallower depth after 80 yr (Figure 4c). The flow direction of overpressure-driven upward flow, thus, deviated within the aquifer in accordance with the higher hydraulic conductivity of this unit and the subsurface hydraulic gradient forced it toward the river valleys, and consequently, the proportion of horizontal flow increased.

Results of the transient simulation also highlighted major changes, which could affect the future state of the lakes. Discharge of overpressured flow into lake 1 will be eliminated due to recharge reduction within the next 80 yr. In the case of the Duna-Tisza Interfluve sample area, where overpressured upward flow contains a high amount of total dissolved solids, it could have significant effects on the quality of the surface waters, and also on the distribution of saline vegetation and soils. In parallel, the groundwater table changes suggest a slight increase (0.4 m, Figure 4c) in lake water level after 80 yr. In contrast, the groundwater level drop of 1.6 m below lake 2 (from 108 m to 106.4 m, Figure 4b,c) is predicted by the model, i.e., it will eventually be dry after 80 yr.

As a next step, different degrees of overpressure were involved in the numerical investigation, in order to evaluate the relation of gravity and overpressure as driving forces in canalizing groundwater and its effects on the transient flow field under changing climatic conditions. The simulated case with “medium” overpressure of 25 MPa at the bottom of the system (induced by groundwater inflow of q_in(bottom) = 10⁻¹⁰ ms⁻¹ Darcy flux across the lower boundary, Figure 4c and Figure 5d) represents the conditions that prevailed in the Duna-Tisza Interfluve sample area. As a comparison, the lower overpressure of 20 MPa (induced by q_in(bottom) = 5 × 10⁻¹¹ ms⁻¹ influx) and a synthetic case with extreme high overpressure of 64 MPa (induced by q_in(bottom) = 5 × 10⁻¹⁰ ms⁻¹ influx) was also tested.

As the results showed, the degree of overpressure at the bottom of the system clearly determines the penetration depth of the upper, topography-driven flow system, which represents, at the same time, the boundary between gravity- and overpressure-driven flow domains. As the degree of overpressure (and consequent hydraulic head increment at the bottom) increases, the penetration depth of the topography-driven local flow systems decreases (e.g., from 600 m without overpressure to 200 m with high overpressure in the case of the L2 system, Figure 5). In parallel, the drop of groundwater level was obviously affected by the magnitude of upward overpressured flow. Recharge reduction led to a 3.6 m decline in the pure gravity case, while to a 3.1 m water level drop involving medium overpressure (represents characteristics similar to the sample area), and even to a smaller, 1.0 m decrease in the case of the highest overpressure (Figure 5), showing that the higher the overpressure, the lower the groundwater level decrease over time. As a consequence, changes in surface water levels are also affected, as proven by the opposite change in the case of lake 2 with the highest overpressure case (0.4 m water level rise, Figure 5b) compared to the pure gravity simulation (2.1 m water level drop, Figure 5e).

By modifications in the groundwater table in the different tested cases, alterations in flow dynamics can also be followed. By increasing the degree of overpressure, the average Darcy fluxes gradually increased compared to the fluxes of the pure-gravity case. Major changes in flux can be observed in the vertical direction in both the aquifer (aquifer 1, 2 and 3, Figure 2b) and in the underlying aquitard. However, the magnitude of flux changes is different in the high hydraulic conductivity aquifer and in the aquitard characterized by low hydraulic conductivity (Figure 6).

5. Discussion

The 2D transient groundwater flow simulation study carried out focused on understanding the effects of decreasing recharge due to climate change. The most conspicuous changes highlighted by the transient simulation were groundwater level decrease and modification of the flow system hierarchy along the studied section in both tested cases, i.e., neglecting and considering the effect of overpressure as a driving force. Water level change was highest at the watershed below the ridge, decreasing toward the river valleys. It is consistent with the findings of Zhao et al. [57], who examined the transient behaviour of nested flow systems and found that water table fluctuations are greatest near the flow divides and are less variable in the main discharge areas. Peters et al. [58] also found that long periods of dry climate had a greater effect on suppressing the groundwater level near the groundwater divide (i.e., the recharge area).

The results showed that as the groundwater level decreases, the hierarchically nested flow field gradually degrades over time, which leads to a less fragmented system with single local cells (Figure 3c and Figure 4c). This is in agreement with the conclusions of Liang et al. [59], who showed that, as groundwater recharge gradually decreases, the flow pattern changes from (i) a fragmented flow field with separated local cells, to (ii) a less fragmented flow field with nested flow systems (including local, intermediate and regional cells) and to (iii) a single regional flow cell. These changes have clear effects on lake water levels and their replenishment by groundwater. Moreover, due to the alteration of subsurface potential distribution—induced by groundwater level changes—penetration depth of infiltrating meteoric water was decreased, and flow intensity and velocity was reduced. Zhou et al. [60] also concluded that the response of the groundwater level to linearly varying climate change shows spatial variability, which obviously affects groundwater flow velocity as well.

Beyond the large-scale changes in flow system dimensions and the relationship between local, intermediate and regional flow systems, the degree of consequential modifications on the flow field is quite different in the case of combined flow systems (i.e., with gravity-driven and overpressured flow domains) compared to pure gravity-driven flow systems, highlighting the modifying effects of overpressure.

Water level decrease was highest (d_hmax = −3.6 m, Figure 5) in the case of the pure gravity-driven simulation, which was moderated by the involving overpressured flow. As the results of the numerical investigation highlighted, the higher the degree of overpressure at the bottom of the system, the lower the decrease in groundwater level after 80 yr. In parallel, the magnitude of upward overpressured flow clearly determines the penetration depth of near-surface local flow systems, which is the deepest in the pure gravity-driven case and gradually decreases with the increase in overpressure (from 600 m to 200 m, Figure 5). It was also revealed that as the flow intensity of a shallow gravity-driven flow system decreased, the relative role of overpressure started to increase, and could buffer the effects of recharge reduction on groundwater levels depending on the degree of overpressure. These results are novel, as the involvement of different fluid driving forces in the analysis of climate change effects on groundwater flow systems is unique.

Another question addressed to the effects of modified flow patterns on groundwater-related shallow surface water bodies (Figure 1). Simulated groundwater level changes suggest differences between lake 1 and lake 2 in terms of their response to recharge reduction. In the case of lake 1, taking place at the lower flank of the ridge, a slight increase in water level is expected, while lake 2, being closer to the watershed at a higher elevation, will likely become dry in the future (Figure 5). The diverse behaviour of the lakes can be originated from their different hydraulic position, i.e., fluctuation in water level is the greatest at the recharge area, while more stable in the main discharge area, described also by Zhao et al. [57]. The water level increase in lake 1—despite the decreasing recharge—could be the consequence of an increasing influence of overpressured flow while weakening the near-surface gravity-driven flow. This phenomenon was also described by Garamhegyi et al. [61] with geostatistical analysis, who experienced increasing trends in shallow groundwater levels of the overpressured flow system (considering the time interval 1961–2010), suggesting the expansion of the overpressured flow system due to the modifications in shallow fluid potential distribution, caused by water shortage of the same sample site (Duna-Tisza Interfluve).

Based on the results of the current numerical simulations, groundwater replenishment of lake 1 is expected to change from a higher ordered flow system to a lower one (from intermediate to a local one in the pure gravity-driven case, Figure 3; from overpressured flow to local gravity-driven flow in case of the combined flow system, Figure 4). These changes could also have a substantial effect on water chemistry and groundwater-dependent ecosystems. As the penetration depth of the near-surface, local flow systems and groundwater fluxes gradually decrease over time, it can be expected that subsurface temperature distribution will also be affected. Consequently, discharging groundwater could gradually become warmer due to warming air temperature, as described by previous authors [21,62,63], and could impact biogeochemical processes, as revealed by Retter et al. [64] and Barbieri et al. [65]. In the case of the Duna-Tisza Interfluve sample area, where overpressured upward flow contains a high amount of total dissolved solids, it could also have significant effects on the recent distribution of saline vegetation and soils.

In a recent study, local scale gravity-driven flow systems were found to be the most vulnerable to atmospheric processes, in agreement with the results of Kurylyk et al. [66]. In contrast, as overpressured flow regimes do not have meteoric water recharge due to the super-hydrostatic vertical pressure component which prevents vertical leakage from above, it is expected to be independent of direct climatic variability.

This study provided a modelling investigation to determine the effects of long-term recharge reduction on regional scale groundwater flow dynamics. Certain limitations are explained as follows that should be taken into consideration when understanding our results. As a large-scale model, it plays a key role in determining long-range flow pathways; however, this spatial scale is not appropriate to simulate detailed, local flow patterns offered by small-scale models. Therefore, local predictions of absolute surface water level changes are beyond the scope of the study. However, simulation results could be appropriate for predicting their future state, i.e., will they likely be maintained or will they eventually be drying in the future, and will they be affected by the possible changes in their interconnected groundwater systems?

Another issue is the dimensionality and hydrostratigraphic configuration. Representing a 3D real system as a 2D cross-section undoubtedly overlooks some system dynamics. However, due to the characteristics of the sample area, a 2D model is considered a reasonable approach for providing the needed insight into the response of the system to the predicted changes of the key hydrologic parameters. In parallel, a simplified hydrostratigraphic configuration was applied, focusing on the key elements of site-specific hydrogeological characteristics of the sample site.

In the case of the calculation of groundwater recharge, the most deterministic hydrologic parameters, i.e., precipitation and evapotranspiration were taken into account, and other parameters (overland flow, through-flow, surface runoff) were neglected in the case of the sample area. However, in other areas where these are also important, it should be included in the recharge calculation.

It is assumed that uncertainties generated by this simplification are not significant compared to the uncertainties inherent in the climate change predictions, and despite these limitations, the current results can effectively reveal general patterns of regional groundwater flow dynamics under long-term climate change.

6. Summary and Conclusions

The transient numerical groundwater flow simulation study focused on understanding the response of groundwater flow systems to climate change and aimed to evaluate the large-scale alterations of groundwater flow patterns and dynamics by looking for general trends and characteristics of the flow field with 2D transient simulations. The main novelty of the study is the large, regional scale of the investigation in contrast to lower scale, site-specific predictions and the involvement of different fluid driving forces in the evaluation, which is unique in the context of the effects of climate change on groundwater systems. The effects of recharge reduction on a pure gravity-driven flow system were analysed, and its flow pattern and dynamic characteristics were compared with the ones created by complex, gravity- and overpressure-driven flow systems. A different degree of overpressure was considered, in order to evaluate its relative role compared to pure gravity in canalizing groundwater and its behaviour under changing hydrological characteristics. It is found that significant large-scale changes are expected in the transient subsurface flow field and flow dynamics due to recharge reduction in all tested cases. From the analyses, the following main conclusions can be drawn:

• In the case of a pure gravity-driven groundwater flow system, climate change-induced recharge reduction led to (i) decreased fragmentation of the uppermost flow domain, i.e., the disappearance of some local flow systems due to the flattened groundwater table; (ii) significant decrease in its penetration depth; and (iii) the development of a slight overpressure within the aquitard unit. The latter can be explained by the significantly lower diffusivity of the aquitard, i.e., the much slower propagation of pore pressure changes in this unit.
• In the case of the combined flow system with upper gravity-driven and overlying overpressured flow domain with “medium” (25 MPa) overpressure, recharge reduction induced groundwater table drop decreased compared to the pure gravity-driven case. In parallel, the boundary between the overpressured and the superimposed gravity-driven flow domain was shifted from the strata boundary to a shallower depth after 80 yr.
• It was highlighted by the involvement of the different degrees of overpressure (with lower, 20 MPa and a synthetic case with extreme high, 64 MPa OP) that the penetration depth of the upper, topography-driven flow system (i.e., the boundary between the gravity-driven and overpressured flow domains) is inversely proportional to the degree of the overpressure. As the degree of overpressure increases, the penetration depth of the topography-driven local flow systems decreases.
• It can be concluded that the higher the overpressure, the lower the climate change-induced groundwater level decrease over time, suggesting the buffering effect of overpressure as the fluid driving force in the flow systems’ response to the changes in hydrologic parameters.

Results of this recent study provided insights into the effects of climate change on groundwater flow systems, which could provide a base to better mitigate and prepare for the consequences of predicted future changes in hydrologic parameters in regional-scale complex groundwater flow systems.