Hydrological Modeling to Assess the Efficiency of Groundwater Replenishment through Natural Reservoirs in the Hungarian Drava River Floodplain

Growing drought hazard and water demand for agriculture, ecosystem conservation, and tourism in the Hungarian Drava river floodplain call for novel approaches to maintain wetland habitats and enhance agricultural productivity. Floodplain rehabilitation should be viewed as a complex landscape ecological issue which, beyond water management goals to relieve water deficit, ensures a high level of provision for a broad range of ecosystem services. This paper explores the hydrological feasibility of alternative water management, i.e., the restoration of natural reservoirs (abandoned paleochannels) to mitigate water shortage problems. To predict the efficiency of the project, an integrated surface water (Wetspass-M) and groundwater model (MODFLOW-NWT) was developed and calibrated with an eight-year data series. Different management scenarios for two natural reservoirs were simulated with filling rates ranging from 0.5 m3 s−1 to 1.5 m3 s−1. In both instances, a natural reservoir with a feeding rate of 1 m3 s−1 was found to be the best scenario. In this case 14 days of filling are required to reach the possible maximum reservoir stage of +2 m. The first meter rise increases the saturation of soil pores and the second creates an open surface water body. Two filling periods per year, each lasting for around 180 days, are required. The simulated water balance shows that reservoir–groundwater interactions are mainly governed by the inflow into and outflow from the reservoir. Such an integrated management scheme is applicable for floodplain rehabilitation in other regions with similar hydromorphological conditions and hazards, too.


Introduction
Water availability is of primary importance for landscape functions. Over the last centuries, half of European wetlands and more than 95% of riverine floodplains were converted to agricultural and urban lands [1]. The human intervention has fundamentally altered the water availability of floodplain landscapes [2]. River floodplains provide various ecosystem services [3,4], including groundwater replenishment, water purification, flood protection, biodiversity/habitat, sediment and nutrient retention, recreation and tourism, cultural values, buffering capacity of climate change, and various floodplain products (fish, wood, reed, game) [5]. Changes in the water regimes of rivers are predicted to be damaging to these ecosystems under extensive human pressure, resources exploitation, three approaches for simulating lake-groundwater interactions (namely, fixed lake stages, high-K nodes, and LAK3 package) were compared by Hunt [59]. Water balance is simulated more realistically in the latest update LAK7 package [49]. The experience gathered from the MODFLOW lake packages is utilized in the modeling of water balance for the reservoir planned in the Drava floodplain.
This paper evaluates the feasibility of water replenishment through natural reservoirs using a hydrological model (Wetspass-M) and a Newton-Raphson formulation for MODFLOW-2005 (MODFLOW-NWT) with ModelMuse as a graphical interface. The consequences of applying a natural reservoir in augmenting surface water storage are investigated under variable management scenarios of recharging or feeding the reservoir. The assessment of such reservoir/groundwater interactions according to different scenarios of reservoir recharge is a precondition to establishing a new water governance in the region, which can, in turn, be the basis of exploiting economic opportunities. Sustainable ecotourism and related development would ensure a safe livelihood for the local population.

Study Area
The 15 km-wide Drava river floodplain is located on the national border between Hungary and Croatia. The Hungarian section of the floodplain is traditionally affected by floods prior to river regulation. Floods mostly take place during June-July due to an Atlantic influence and during October-November due to a Mediterranean influence. Drava river channelization began in 1750 and reduced the river length by 50% [60]. The floodplain has a total area of 500 km 2 . There are 13 tributary streams, 20 major side channels, hundreds of meander scars, and 18 oxbow lakes in the Hungarian section of the Drava River ( Figure 1) [61]. The long-term average discharge of the Drava River was 595 m 3 s −1 for the period 1896-2018 [62]. Maximum discharge of 0.1% probability is assumed to be 3070 m 3 s −1 (reached during the flood of 1827), while 170 m 3 s −1 and 900 m 3 s −1 are accounted to baseflow and bankfull discharge, respectively [62]. As a possible source of water replenishment, the Fekete-víz stream is a tributary with an average discharge of Qav = 2-2.5 m 3 s −1 . The main features of the floodplain are oxbows, scroll bar systems, abandoned channels with natural levees, backswamps, and in the north, elongated blown sand of east to west strike [63]. The digital elevation model (DEM) of the study area was obtained from the South-Transdanubian Water Management Directorate, Pécs, at 10 m resolution. The Drava floodplain is apparently flat at 83-133 m elevation, with an average relative relief 2 m km −2 ( Figure 1). Settlements were built on elevated terrain.
Water 2020, 12, x FOR PEER REVIEW 3 of 21 more realistically in the latest update LAK7 package [49]. The experience gathered from the MODFLOW lake packages is utilized in the modeling of water balance for the reservoir planned in the Drava floodplain. This paper evaluates the feasibility of water replenishment through natural reservoirs using a hydrological model (Wetspass-M) and a Newton-Raphson formulation for MODFLOW-2005 (MODFLOW-NWT) with ModelMuse as a graphical interface. The consequences of applying a natural reservoir in augmenting surface water storage are investigated under variable management scenarios of recharging or feeding the reservoir. The assessment of such reservoir/groundwater interactions according to different scenarios of reservoir recharge is a precondition to establishing a new water governance in the region, which can, in turn, be the basis of exploiting economic opportunities. Sustainable ecotourism and related development would ensure a safe livelihood for the local population.

Study Area
The 15 km-wide Drava river floodplain is located on the national border between Hungary and Croatia. The Hungarian section of the floodplain is traditionally affected by floods prior to river regulation. Floods mostly take place during June-July due to an Atlantic influence and during October-November due to a Mediterranean influence. Drava river channelization began in 1750 and reduced the river length by 50% [60]. The floodplain has a total area of 500 km 2 . There are 13 tributary streams, 20 major side channels, hundreds of meander scars, and 18 oxbow lakes in the Hungarian section of the Drava River ( Figure 1) [61]. The long-term average discharge of the Drava River was 595 m 3 s −1 for the period 1896-2018 [62]. Maximum discharge of 0.1% probability is assumed to be 3070 m 3 s −1 (reached during the flood of 1827), while 170 m 3 s −1 and 900 m 3 s −1 are accounted to baseflow and bankfull discharge, respectively [62]. As a possible source of water replenishment, the Fekete-víz stream is a tributary with an average discharge of Qav = 2-2.5 m 3 s −1 . The main features of the floodplain are oxbows, scroll bar systems, abandoned channels with natural levees, backswamps, and in the north, elongated blown sand of east to west strike [63]. The digital elevation model (DEM) of the study area was obtained from the South-Transdanubian Water Management Directorate, Pécs, at 10 m resolution. The Drava floodplain is apparently flat at 83-133 m elevation, with an average relative relief 2 m km −2 ( Figure 1). Settlements were built on elevated terrain.  The meteorological data for the period 2010-2018 (time selected as this study period) were obtained from the South-Transdanubian Water Management Directorate, Pécs. The climate of the Drava Plain is moderately wet and warm in the west under Atlantic influence, and moderately dry and warm in the east under Mediterranean influence [64]. The long-term mean annual temperature is 10.8 • C in the east and 10.2 • C in the west. Absolute maximum and minimum temperatures are 35 • C and −17 • C, respectively. The number of cold days (mean daily temperature below −10 • C) is 10, while the number of hot days (mean temperature above 30 • C) is 21. The long-term average annual rainfall ranges from 383 mm y −1 to 1057 mm y −1 , with a mean value of 682 mm y −1 . The rainy season (summer and spring) has a share of 60% of the annual rainfall, while the remaining 40% falls in the drier season (winter and autumn) [65]. Daily evaporation from the lake surface was calculated using the Penman open-water equation [66] based on wind speed, hourly incoming solar radiation, air temperature, relative humidity, air temperature, and wind-speed. Evaporation from the lake is concentrated on the summer months and shows a large variation between 0.0 mm d −1 and 6.8 mm d −1 , with a mean value of 1.85 mm d −1 . Daily potential evapotranspiration (PET) was computed from meteorological data using the Food and Agriculture Organization of the United Nations (FAO) Penman-Monteith method [67] and varies from 0.0 mm to 46 mm, with an average of 19.8 mm. About 88% of the PET is attributed to summer and spring.

Soil Data
Altogether, 72 soil samples were taken from the floodplain for particle size analysis and the identification of textural classes (Figure 2a). Over around 50% of the floodplain soils are of sandy texture, in 12% clay, in 8% loam, while in over 30% of the area sandy loam and clay loam occur. Fluvisols formed under periodical water surplus and forest vegetation dominate the area. Soil thickness is highly variable, from a few centimeters to a meter, the top 0.5 m has more organic material. In addition to particle size distribution, soils were studied using GPR (ground penetrating radar) and satellite image analysis, in combination with the survey of landforms [68,69].

Meteorological Data
The meteorological data for the period 2010-2018 (time selected as this study period) were obtained from the South-Transdanubian Water Management Directorate, Pécs. The climate of the Drava Plain is moderately wet and warm in the west under Atlantic influence, and moderately dry and warm in the east under Mediterranean influence [64]. The long-term mean annual temperature is 10.8 °C in the east and 10.2 °C in the west. Absolute maximum and minimum temperatures are 35 °C and −17 °C, respectively. The number of cold days (mean daily temperature below −10 °C) is 10, while the number of hot days (mean temperature above 30 °C) is 21. The long-term average annual rainfall ranges from 383 mm y −1 to 1057 mm y −1 , with a mean value of 682 mm y −1 . The rainy season (summer and spring) has a share of 60% of the annual rainfall, while the remaining 40% falls in the drier season (winter and autumn) [65]. Daily evaporation from the lake surface was calculated using the Penman open-water equation [66] based on wind speed, hourly incoming solar radiation, air temperature, relative humidity, air temperature, and wind-speed. Evaporation from the lake is concentrated on the summer months and shows a large variation between 0.0 mm d −1 and 6.8 mm d −1 , with a mean value of 1.85 mm d −1 . Daily potential evapotranspiration (PET) was computed from meteorological data using the Food and Agriculture Organization of the United Nations (FAO) Penman-Monteith method [67] and varies from 0.0 mm to 46 mm, with an average of 19.8 mm. About 88% of the PET is attributed to summer and spring.

Soil Data
Altogether, 72 soil samples were taken from the floodplain for particle size analysis and the identification of textural classes (Figure 2a). Over around 50% of the floodplain soils are of sandy texture, in 12% clay, in 8% loam, while in over 30% of the area sandy loam and clay loam occur. Fluvisols formed under periodical water surplus and forest vegetation dominate the area. Soil thickness is highly variable, from a few centimeters to a meter, the top 0.5 m has more organic material. In addition to particle size distribution, soils were studied using GPR (ground penetrating radar) and satellite image analysis, in combination with the survey of landforms [68,69].    All texture classes turned out in the top 20 m of the modeling space. These sediments were deposited in point bars, channels, and oxbows [69]. During the late Pleistocene-Holocene in the top 100 m fluvial coarse-sand-heavy clay series deposited in the basin (Figure 2b). Subsurface water flow and capillary rise are influenced by the above-mentioned landforms and deposits.

Land Use and Groundwater Depth
Land use data for the Drava floodplain in 2018 were obtained from the Coordination of Information on the Environment (CORINE) data base (https://land.copernicus.eu/pan-european/corine-land-cover/ clc2018). Land use was composed of 13 classes, with a dominance of agriculture (45%), followed by deciduous forests (23.83%), whereas shrub areas, sealed surfaces, and water bodies took up 11%, 4%, and 1.1%, respectively (Figure 3a). For daily groundwater depth, 18 observation wells were used to depict the spatial distribution of groundwater level for the Drava floodplain ( Figure 3b) coarse-sand-heavy clay series deposited in the basin (Figure 2b). Subsurface water flow and capillary rise are influenced by the above-mentioned landforms and deposits.

Land Use and Groundwater Depth
Land use data for the Drava floodplain in 2018 were obtained from the Coordination of Information on the Environment (CORINE) data base (https://land.copernicus.eu/paneuropean/corine-land-cover/clc2018). Land use was composed of 13 classes, with a dominance of agriculture (45%), followed by deciduous forests (23.83%), whereas shrub areas, sealed surfaces, and water bodies took up 11%, 4%, and 1.1%, respectively (Figure 3a). For daily groundwater depth, 18 observation wells were used to depict the spatial distribution of groundwater level for the Drava floodplain ( Figure 3b

Groundwater Recharge and Evapotranspiration
The spatial distribution of groundwater recharge and evapotranspiration were estimated by the WetSpass-M model [71]. The WetSpass model has been applied to assess water balance components in different regions around the world [72][73][74][75][76][77]. The WetSpass model considers the spatial distribution of the soil, land use, slope, groundwater depth, and the meteorological data of every grid cell are divided into four fractions (vegetated, bare soil, open water surface). The seasonal water balance for each fraction is calculated as: where P is precipitation, S is surface runoff, ET is evapotranspiration, and R is groundwater recharge. The actual evapotranspiration is calculated as a function of potential evapotranspiration and a vegetation coefficient [78]. R is determined as a residual part of the water balance. The simulated annual groundwater recharge of the Drava floodplain ranged from 0 mm y −1 to 360 mm y −1 as the lowest and highest values, with an average of 241 mm y −1 (Figure 4a). The estimated evapotranspiration varied from 211 mm y −1 to a maximum of 1394 mm y −1 , with an average value of 310 mm y −1 (Figure 4b).

Groundwater Recharge and Evapotranspiration
The spatial distribution of groundwater recharge and evapotranspiration were estimated by the WetSpass-M model [71]. The WetSpass model has been applied to assess water balance components in different regions around the world [72][73][74][75][76][77]. The WetSpass model considers the spatial distribution of the soil, land use, slope, groundwater depth, and the meteorological data of every grid cell are divided into four fractions (vegetated, bare soil, open water surface). The seasonal water balance for each fraction is calculated as: where P is precipitation, S is surface runoff, ET is evapotranspiration, and R is groundwater recharge. The actual evapotranspiration is calculated as a function of potential evapotranspiration and a vegetation coefficient [78]. R is determined as a residual part of the water balance. The simulated annual groundwater recharge of the Drava floodplain ranged from 0 mm y −1 to 360 mm y −1 as the lowest and highest values, with an average of 241 mm y −1 (Figure 4a). The estimated evapotranspiration varied from 211 mm y −1 to a maximum of 1394 mm y −1 , with an average value of 310 mm y −1 (Figure 4b).

Numerical Groundwater Modeling
MODFLOW-NWT [79] with ModelMuse [80] as a graphical user interface was used to simulate groundwater flow. It was expressly developed to solve problems which are caused by nonlinearities of drying and rewetting of the unconfined groundwater flow, similar to that treated in this paper. This is in contrast with other versions of MODFLOW, which exclude dry cells from calculation. The LAK7 package was developed to simulate lake/groundwater interactions. The lake stage is determined based on the exchanges of volumetric water in and out of the lake and the overall water balance of the lake. The seepage between a lake and its groundwater system depends on lake stage, lakebed leakage in grid cell dimensions, the hydraulic conductivity of the lakebed material, the adjacent aquifer, and hydraulic heads in the groundwater system. The seepage is calculated by Darcy's law [54], as follows: where q is the seepage rate (L T −1 ); K is the hydraulic conductivity (L T −1 ) of materials between the lake and the aquifer; ℎ is the stage of the lake (L); ℎ is the aquifer head (L); d is the distance (L) between the points at which ℎ and ℎ are measured; and the volumetric flux Q (L 3 T −1 ) is established by the following equation: where c = KA/d is the conductance (L 2 T −1 ) and K/d is the leakage (T -1 ). Inputs to the groundwater system are groundwater recharge from precipitation, river seepage and lake seepage, lateral groundwater inflow from a part of the northern boundary, while output from the groundwater system includes groundwater seepage to rivers and lakes, lateral groundwater outflow from the southwestern boundary, and evaporation from groundwater.

Groundwater Flow Model Setup
According to evidence from geological boreholes, pumping tests and GPR surveys, the groundwater flow system consists of four layers. Since the sediments in the area range from coarse sand to clay, the four layers are highly heterogeneous in space. The model was discretized by a 100 m by 100 m grid and the model was refined for lake area with 30 m by 30 m, resulting in a domain of 361 rows and 742 columns, and four layers with a total number of 1,339,310 cells, 655,390 of which were active cells. The records of GPR, supplemented with auger holes, presented high heterogeneity in the hydraulic proprieties of the sediment depending on particle size (ranging from coarse sand to clay). The fall head method was applied to calculate the hydraulic conductivity for each borehole. According to the boreholes data, the first layer consisted of four K-zones of hydraulic conductivity (

Numerical Groundwater Modeling
MODFLOW-NWT [79] with ModelMuse [80] as a graphical user interface was used to simulate groundwater flow. It was expressly developed to solve problems which are caused by nonlinearities of drying and rewetting of the unconfined groundwater flow, similar to that treated in this paper. This is in contrast with other versions of MODFLOW, which exclude dry cells from calculation. The LAK7 package was developed to simulate lake/groundwater interactions. The lake stage is determined based on the exchanges of volumetric water in and out of the lake and the overall water balance of the lake. The seepage between a lake and its groundwater system depends on lake stage, lakebed leakage in grid cell dimensions, the hydraulic conductivity of the lakebed material, the adjacent aquifer, and hydraulic heads in the groundwater system. The seepage is calculated by Darcy's law [54], as follows: where q is the seepage rate (L T −1 ); K is the hydraulic conductivity (L T −1 ) of materials between the lake and the aquifer; h 1 is the stage of the lake (L); h a is the aquifer head (L); d is the distance (L) between the points at which h 1 and h a are measured; and the volumetric flux Q (L 3 T −1 ) is established by the following equation: where c = KA/d is the conductance (L 2 T −1 ) and K/d is the leakage (T −1 ). Inputs to the groundwater system are groundwater recharge from precipitation, river seepage and lake seepage, lateral groundwater inflow from a part of the northern boundary, while output from the groundwater system includes groundwater seepage to rivers and lakes, lateral groundwater outflow from the southwestern boundary, and evaporation from groundwater.

Groundwater Flow Model Setup
According to evidence from geological boreholes, pumping tests and GPR surveys, the groundwater flow system consists of four layers. Since the sediments in the area range from coarse sand to clay, the four layers are highly heterogeneous in space. The model was discretized by a 100 m by 100 m grid and the model was refined for lake area with 30 m by 30 m, resulting in a domain of 361 rows and 742 columns, and four layers with a total number of 1,339,310 cells, 655,390 of which were active cells. The records of GPR, supplemented with auger holes, presented high heterogeneity in the hydraulic proprieties of the sediment depending on particle size (ranging from coarse sand to clay). The fall head method was applied to calculate the hydraulic conductivity for each borehole. According to the boreholes data, the first layer consisted of four K-zones of hydraulic conductivity ( Figure 5), varying from 375 to 0.15 m d −1 , while the second layer was characterized by three zones with a range of 375 and 75 m d −1 . The third and fourth layers were described by two zones, with a range from between 375 and 75 m d −1 and 0.15 to 3 m d −1 , respectively ( Figure 5). Although all the hydraulic parameters were determined in the field, because of limited spatial representativeness and the influence of modeling scale in accordance with [81][82][83], they still had to be adjusted in the calibration process. Fortunately, the majority of the units had natural boundaries.
Water 2020, 12, x FOR PEER REVIEW 7 of 21 hydraulic parameters were determined in the field, because of limited spatial representativeness and the influence of modeling scale in accordance with [81][82][83], they still had to be adjusted in the calibration process. Fortunately, the majority of the units had natural boundaries. The southern boundary is represented by the Drava River and assigned as a river package. Based on data available from observation wells, the general direction of groundwater flow is from north to south. The Fekete-víz and Okor streams comprise the northern, eastern, and part of the western boundary and they are assigned by a river package ( Figure 6a). All streams in the floodplain are assigned by a river package. The conductance of the riverbed was initially set to 5 m 2 d −1 , and then adjusted during the calibration process. MODFLOW General-Head Boundary (GHB) was assigned to the southwestern and part of the northern boundary to simulate lateral groundwater inflow and outflow of the system based on the time series of groundwater level from the available observation wells. Initially, the conductance of GHB was set to be 20 m 2 d −1 , based on geological boreholes and adjusted during the calibration process. Cún-Szaporca Lake was assigned in the model through the LAK7 package (Figure 6a). The lake simulation was performed using a separate additional inactive layer [54] with variable thickness. The top of inactive layer had a 1.0 m thickness outside the lake and the top equaled the maximum lake stage within the lake area. The bottom was equal to the bathymetrically defined lake bottom. Lakebed leakage was assigned from field measurements for the vertical hydraulic conductivity and the thickness of the lakebed. Lakebed leakage was set to 0.05 d −1 . The southern boundary is represented by the Drava River and assigned as a river package. Based on data available from observation wells, the general direction of groundwater flow is from north to south. The Fekete-víz and Okor streams comprise the northern, eastern, and part of the western boundary and they are assigned by a river package (Figure 6a). All streams in the floodplain are assigned by a river package. The conductance of the riverbed was initially set to 5 m 2 d −1 , and then adjusted during the calibration process. MODFLOW General-Head Boundary (GHB) was assigned to the southwestern and part of the northern boundary to simulate lateral groundwater inflow and outflow of the system based on the time series of groundwater level from the available observation wells. Initially, the conductance of GHB was set to be 20 m 2 d −1 , based on geological boreholes and adjusted during the calibration process. Cún-Szaporca Lake was assigned in the model through the LAK7 package (Figure 6a). The lake simulation was performed using a separate additional inactive layer [54] with variable thickness. The top of inactive layer had a 1.0 m thickness outside the lake and the top equaled the maximum lake stage within the lake area. The bottom was equal to the bathymetrically defined lake bottom. Lakebed leakage was assigned from field measurements for the vertical hydraulic conductivity and the thickness of the lakebed. Lakebed leakage was set to 0.05 d −1 .

Calibration and Sensitivity Analysis
A steady-state model was developed to represent an initial case for the transient model by applying long-term average parameters, i.e., mean precipitation, mean groundwater levels, mean PET, and mean groundwater recharge for the period from 1 January 2010 to 20 September 2018. The steady-state model was calibrated using average groundwater levels for 16 observation wells. Calibrated results showed a good agreement between observed and simulated groundwater heads, with R 2 = 0.98. The mean error was −0.016 m, and the mean absolute error was 0.09 m. However, starting the transient model by steady-state model as an initial condition proved to be unsuccessful due to a large difference between the steady-state parameters and the measured daily values at the start-up of the transient model on 1 January 2010.
The first three hydrological years from which data were available, i.e., from 1 January 2010 until 31 December 2012, were used as a warm-up period, after which the model was run (calibrated) throughout the following six hydrological years from 1 January 2013 to 20 September 2018 with a daily time step. The calibrated results of the transient model against the time series of observation well heads and the stages of the Cún-Szaporca Lake are shown in Figure 7. The calibration process of groundwater head was carried out using 16 time series of daily groundwater levels from observation wells extending over 10 years. Calibrated results showed a good match between observed and simulated groundwater heads, with R 2 = 0.98. The mean error was −0.08 m, and the root mean square error was 0.4 m (Figure 7a). Additionally, the calibration of lake stages showed a good agreement between simulated and observed stages with a correlation coefficient of 0.90, a mean error of −0.07 m, mean absolute error of 0.15 m, and root mean square error of 0.2 m (Figure 7b). The calibrated

Calibration and Sensitivity Analysis
A steady-state model was developed to represent an initial case for the transient model by applying long-term average parameters, i.e., mean precipitation, mean groundwater levels, mean PET, and mean groundwater recharge for the period from 1 January 2010 to 20 September 2018. The steady-state model was calibrated using average groundwater levels for 16 observation wells. Calibrated results showed a good agreement between observed and simulated groundwater heads, with R 2 = 0.98. The mean error was −0.016 m, and the mean absolute error was 0.09 m. However, starting the transient model by steady-state model as an initial condition proved to be unsuccessful due to a large difference between the steady-state parameters and the measured daily values at the start-up of the transient model on 1 January 2010.
The first three hydrological years from which data were available, i.e., from 1 January 2010 until 31 December 2012, were used as a warm-up period, after which the model was run (calibrated) throughout the following six hydrological years from 1 January 2013 to 20 September 2018 with a daily time step. The calibrated results of the transient model against the time series of observation well heads and the stages of the Cún-Szaporca Lake are shown in Figure 7. The calibration process of groundwater head was carried out using 16 time series of daily groundwater levels from observation wells extending over 10 years. Calibrated results showed a good match between observed and simulated groundwater heads, with R 2 = 0.98. The mean error was −0.08 m, and the root mean square error was 0.4 m (Figure 7a). Additionally, the calibration of lake stages showed a good agreement between simulated and observed stages with a correlation coefficient of 0.90, a mean error of −0.07 m, mean absolute error of 0.15 m, and root mean square error of 0.2 m (Figure 7b). The calibrated parameters of hydraulic conductivity, river conductance, general head conductance, and groundwater recharge are shown in Table 1.
Water 2020, 12, x FOR PEER REVIEW 9 of 21 parameters of hydraulic conductivity, river conductance, general head conductance, and groundwater recharge are shown in Table 1. The sensitivity analysis was achieved using computer code for universal sensitivity analysis, calibration, and uncertainty evaluation (UCODE-2005) [84] with the help of ModelMate [85]. It was studied for the groundwater recharge (RCH), hydraulic conductivity (HK), which was divided into five zones, evapotranspiration (EVT), river conductance (RIVC), and general head boundary (GHB) parameters. The parameters of the groundwater recharge (RCH), hydraulic conductivity (HKZ1), and evapotranspiration (EVT) had the highest composite scaled sensitivity values and were therefore the most sensitive parameters (Figure 7c).  The sensitivity analysis was achieved using computer code for universal sensitivity analysis, calibration, and uncertainty evaluation (UCODE-2005) [84] with the help of ModelMate [85]. It was studied for the groundwater recharge (RCH), hydraulic conductivity (HK), which was divided into five zones, evapotranspiration (EVT), river conductance (RIVC), and general head boundary (GHB) parameters. The parameters of the groundwater recharge (RCH), hydraulic conductivity (HKZ1), and evapotranspiration (EVT) had the highest composite scaled sensitivity values and were therefore the most sensitive parameters (Figure 7c).

Water Budget of the Groundwater Flow Model
The water balance at the end stress of the calibrated transient model is shown in Table 2. The inflow rate from GHB boundary to the aquifer represents 42.25% of the total inflow to the aquifer through 627,789.00 m 3 d −1 . Water seepage from streams and the Cún-Szaporca Lake represents 4.68% of the total inflow. Groundwater discharge from the aquifer to streams is 1,275,244.38 m 3 d −1 , which amounts to 85.82% of the total outflow from groundwater system as the area characterized by shallow groundwater table, while evapotranspiration represents 10% of total outflow.

Simulated Scenarios
The natural water reservoirs are the former oxbows, paleomeanders, and deeper-lying floodplain sections, which are under excess water effect. These areas may be suitable for the storage of surplus water provided by the adjacent watercourses. The model space covered two subareas: 1, the Korcsina oxbow system (Figure 6b); and 2, the Okor-Fekete-víz backswamps (Figure 6c). The exact location of the potential reservoir was selected from the analysis of the topography and hydraulic soil properties explored by auger samples. The selected subareas were simulated as natural reservoirs with filling and discharge conditions based on +2 m vertical water depth increasing theoretically with different stream discharges of 0.5, 0.75, 1, and 1.5 m 3 s −1 . (The excess water hazard on agricultural land excludes water level rise above +2 m.) Based on measurements of stream discharge, to obtain feeding with stream discharge larger than 1 m 3 d −1 , constructing a dam is required. Evaporation from the open water surface of the reservoir is a main challenge of the natural reservoir solution. The forested environs of the reservoir, however, reduce evaporation. This model is the first attempt to investigate whether water management through establishing a natural reservoir can be a major contribution to floodplain rehabilitation. Results of the management scenarios are discussed below.

Korcsina Subarea
The Korcsina paleochannel system is a potential water reservoir in the western central part of the floodplain. With streams feeding 0.5 m 3 s −1 , the total upfilling (+2 m relative water level from 98 to 100 m) in the Korcsina reservoir would take 55 days, with 2,376,000 m 3 water demand. The discharge (seepage) period lasts 11 days from a level of 100 m to 98.75 m. By the end of the filling period, seepage from the reservoir to groundwater reaches 35,815 m 3 . Recharge from groundwater to the reservoir in the filling period from the water level stage of 98 to 99 m and over this level becomes zero (Figure 8a). For filling at a rate of 0.75 m 3 s −1 , the reservoir reaches the stage 100 m after 27 days with reservoir seepage of 38,900 m 3 , while the recharge from groundwater to the reservoir takes place in the first five days of the filling process at the level of 99 m (Figure 8b). In the filling with a discharge of 1 and 1.5 m 3 s −1 , 18 and 11 days are required to achieve +2 m relative water level rise in the reservoir (Figure 8c,d). For all scenarios with different discharge to the reservoir, the reservoir stage did not reach 98 m again because of the shallow groundwater table, which allows groundwater recharge larger than the seepage when reaching a level of 98.83 m (i.e., the inflow from groundwater to the reservoir, which is equal to the seepage from the reservoir at this stage) (Figure 8).
Water 2020, 12, x FOR PEER REVIEW 11 of 21 (Figure 8c,d). For all scenarios with different discharge to the reservoir, the reservoir stage did not reach 98 m again because of the shallow groundwater table, which allows groundwater recharge larger than the seepage when reaching a level of 98.83 m (i.e., the inflow from groundwater to the reservoir, which is equal to the seepage from the reservoir at this stage) (Figure 8).

Simulation of the Water Storage Opportunities in the Okor-Fekete-víz Area
These randomly distributed and irregular-shaped paleochannels and backswamps are located at the confluence of the Okor and Fekete-víz streams. The main textural types in the west are medium to coarse sands, in the middle part fine sands, and in the southeastern part loamy sands (Figure 2a). The bed level of the Okor natural reservoir is at 95.95 m. Feeding is from the two mentioned streams. With streams feeding at a rate of 0.5 m 3 s −1 , the reservoir stage does not reach +2 m relative water level (from 95.95 to 97.95 m), but it reaches 97.7 m within 61 days (Figure 9a

Simulation of the Water Storage Opportunities in the Okor-Fekete-víz Area
These randomly distributed and irregular-shaped paleochannels and backswamps are located at the confluence of the Okor and Fekete-víz streams. The main textural types in the west are medium to coarse sands, in the middle part fine sands, and in the southeastern part loamy sands (Figure 2a). The bed level of the Okor natural reservoir is at 95.95 m. Feeding is from the two mentioned streams. With streams feeding at a rate of 0.5 m 3 s −1 , the reservoir stage does not reach +2 m relative water level (from 95.95 to 97.95 m), but it reaches 97.7 m within 61 days (Figure 9a The entrenchment and recent incision of the Drava River resulted in a dropping of groundwater tables by 2.5 m. For both reservoirs, and with a feeding discharge of 0.5 and 0.75 m 3 s −1 , the feeding stream capacity does not allow the recharge of the reservoirs in 60 days and 25 days, respectively, based on the observed discharge of the feeding stream. However, feeding the reservoir at a rate of 1.5 m 3 s −1 , a +2 m rise to be achieved in five days requires the construction of a dam on the feeding stream. Therefore, this scenario would not be economically feasible. Thus, the best scenario for both reservoirs is a feeding rate of 1.0 m 3 s −1 , as 18 and 13 days are needed for filling the Korcsina and Okor reservoirs, respectively, and water recharge to the reservoirs could be implemented by two filling periods per year, each providing water storage for six months. The drop of the water level in the reservoir after the filling period is attributed to the pressure differences between the reservoir water and the adjacent groundwater. The reservoir bed seepage analysis was a crucial step to understand the interactions between the reservoir and groundwater system. The seepage from the reservoir to the groundwater system, up to 25,000 m 3 d −1 , was dominant during stages higher than ~99 m, whereas the recharge from groundwater to the reservoir, up to 5000 m 3 d −1 , was dominant during stages below 99.00 m. Water levels in the reservoirs are primarily effected by the surrounding groundwater levels, regulated by the hyporheic flow of the streams. The first 1 m of filling increases the saturation of soil pores, which raises soil moisture content in the root zone and the second meter creates an open water body, which is useful to create wetland habitats in the floodplain. The entrenchment and recent incision of the Drava River resulted in a dropping of groundwater tables by 2.5 m. For both reservoirs, and with a feeding discharge of 0.5 and 0.75 m 3 s −1 , the feeding stream capacity does not allow the recharge of the reservoirs in 60 days and 25 days, respectively, based on the observed discharge of the feeding stream. However, feeding the reservoir at a rate of 1.5 m 3 s −1 , a +2 m rise to be achieved in five days requires the construction of a dam on the feeding stream. Therefore, this scenario would not be economically feasible. Thus, the best scenario for both reservoirs is a feeding rate of 1.0 m 3 s −1 , as 18 and 13 days are needed for filling the Korcsina and Okor reservoirs, respectively, and water recharge to the reservoirs could be implemented by two filling periods per year, each providing water storage for six months. The drop of the water level in the reservoir after the filling period is attributed to the pressure differences between the reservoir water and the adjacent groundwater. The reservoir bed seepage analysis was a crucial step to understand the interactions between the reservoir and groundwater system. The seepage from the reservoir to the groundwater system, up to 25,000 m 3 d −1 , was dominant during stages higher than~99 m, whereas the recharge from groundwater to the reservoir, up to 5000 m 3 d −1 , was dominant during stages below 99.00 m. Water levels in the reservoirs are primarily effected by the surrounding groundwater levels, regulated by the hyporheic flow of the streams. The first 1 m of filling increases the saturation of soil pores, which raises soil moisture content in the root zone and the second meter creates an open water body, which is useful to create wetland habitats in the floodplain.

Spatio-Temporal Extent of Reservoir Impact
The assessment of the spatio-temporal impact of a lake (reservoir) on groundwater conditions is required to manage the groundwater systems adjacent to artificial lakes. The developed groundwater head of lakes Korcsina and Okor for a filling scenario with a discharge of 1 m 3 s −1 at different filling and discharge periods were plotted against distances at a cross-section which passes through the lakes (Figure 6b,c) for the selected stress periods (Figures 10 and 11). The influence of the lake on groundwater levels decreases with distance (Figures 10 and 11). However, there are differences between the trends of the two sections due to the variance in hydrogeological and hydraulic conditions for the two locations. The results clearly indicate that after 10 and 20 days of filling the Korcsina reservoir, the average groundwater head increased by 1 m and 1.6 m, respectively, with respect to the base case along the cross section. For the Okor reservoir, the average groundwater head increased by 0.7 m and 1.3 m after filling periods of 10 days and 20 days, respectively, with respect to the base case, while the average groundwater table decreased by 0.9 m and 1.0 m after discharge periods of 70 days and 170 days, respectively, compared to the base case ( Figure 11). As depicted in the gradual rise in groundwater table after successive filling (Figures 10 and 11), the natural reservoirs improve aquifer storage and hence maintain the level of ecosystem services and increase agricultural productivity in the Drava basin.
Lake-aquifer exchange fluxes in the Korcsina and Okor reservoirs are shown in Figures 12 and 13, identifying the losing and gaining zones of the planned reservoirs with the help of the cell-to-cell flux terms in the model output files. Exchange fluxes between the reservoirs and the groundwater system for different periods of filling and discharge were analyzed. The sign on the graph indicates the direction of fluxes between the reservoirs and the aquifer. A positive sign represents upward reservoir seepage gain and negative represents downward reservoir seepage loss. For the Korcsina reservoir with 20 days of filling (+2 m), the seepage from reservoir to groundwater system ranged from 0 to 245 m 3 d −1 with an average value of 55 m 3 d −1 , with higher values in the northeast (Figure 12a). The average exchange between reservoir and aquifer was 1.4 m 3 d −1 after 40 days of discharge period, the majority of which the reservoir has gained from the aquifer (Figure 12b).
Water 2020, 12, x FOR PEER REVIEW 13 of 21

Spatio-Temporal Extent of Reservoir Impact
The assessment of the spatio-temporal impact of a lake (reservoir) on groundwater conditions is required to manage the groundwater systems adjacent to artificial lakes. The developed groundwater head of lakes Korcsina and Okor for a filling scenario with a discharge of 1 m 3 s −1 at different filling and discharge periods were plotted against distances at a cross-section which passes through the lakes (Figure 6b,c) for the selected stress periods (Figures 10 and 11). The influence of the lake on groundwater levels decreases with distance ( Figures 10 and 11). However, there are differences between the trends of the two sections due to the variance in hydrogeological and hydraulic conditions for the two locations. The results clearly indicate that after 10 and 20 days of filling the Korcsina reservoir, the average groundwater head increased by 1 m and 1.6 m, respectively, with respect to the base case along the cross section. For the Okor reservoir, the average groundwater head increased by 0.7 m and 1.3 m after filling periods of 10 days and 20 days, respectively, with respect to the base case, while the average groundwater table decreased by 0.9 m and 1.0 m after discharge periods of 70 days and 170 days, respectively, compared to the base case ( Figure 11). As depicted in the gradual rise in groundwater table after successive filling (Figures 10 and 11), the natural reservoirs improve aquifer storage and hence maintain the level of ecosystem services and increase agricultural productivity in the Drava basin.    Lake-aquifer exchange fluxes in the Korcsina and Okor reservoirs are shown in Figures 12 and  13, identifying the losing and gaining zones of the planned reservoirs with the help of the cell-to-cell flux terms in the model output files. Exchange fluxes between the reservoirs and the groundwater system for different periods of filling and discharge were analyzed. The sign on the graph indicates the direction of fluxes between the reservoirs and the aquifer. A positive sign represents upward reservoir seepage gain and negative represents downward reservoir seepage loss. For the Korcsina reservoir with 20 days of filling (+2 m), the seepage from reservoir to groundwater system ranged from 0 to 245 m 3 d −1 with an average value of 55 m 3 d −1 , with higher values in the northeast ( Figure  12a). The average exchange between reservoir and aquifer was 1.4 m 3 d −1 after 40 days of discharge period, the majority of which the reservoir has gained from the aquifer (Figure 12b). Similar seepage patterns can be observed for the Okor reservoir (Figure 13a,b), with average values of seepage of 25.3 m 3 d −1 and 0.7 m 3 d −1 after 14 days of filling (+2 m) and 40 days of discharge (+1 m), respectively. The patterns of spatial seepage at various reservoir stages through filling and discharge periods are different, attributed to variations in reservoir-bed leakage and reservoir stages. The magnitude and direction of seepage between the reservoir and underlying groundwater system relies on the reservoir-bed leakage and the difference between reservoir stage and hydraulic head of an underlying aquifer. In the Drava Plain, of highly variable alluvial deposits and soil textural types, Similar seepage patterns can be observed for the Okor reservoir (Figure 13a,b), with average values of seepage of 25.3 m 3 d −1 and 0.7 m 3 d −1 after 14 days of filling (+2 m) and 40 days of discharge (+1 m), respectively. The patterns of spatial seepage at various reservoir stages through filling and discharge periods are different, attributed to variations in reservoir-bed leakage and reservoir stages. The magnitude and direction of seepage between the reservoir and underlying groundwater system relies on the reservoir-bed leakage and the difference between reservoir stage and hydraulic head of an underlying aquifer. In the Drava Plain, of highly variable alluvial deposits and soil textural types, the modeling of reservoir/groundwater interchanges is central to the planning of rehabilitation measures.
(b) Figure 12. Lake-aquifer exchange fluxes for the Korcsina reservoir: (a) at the end of filling period of 14 days (+2 m reservoir stage), (b) for 40 days of discharge period (+1 m) (positive represents upward reservoir seepage gain and negative represents downward reservoir seepage loss). Similar seepage patterns can be observed for the Okor reservoir (Figure 13a,b), with average values of seepage of 25.3 m 3 d −1 and 0.7 m 3 d −1 after 14 days of filling (+2 m) and 40 days of discharge (+1 m), respectively. The patterns of spatial seepage at various reservoir stages through filling and discharge periods are different, attributed to variations in reservoir-bed leakage and reservoir stages. The magnitude and direction of seepage between the reservoir and underlying groundwater system relies on the reservoir-bed leakage and the difference between reservoir stage and hydraulic head of an underlying aquifer. In the Drava Plain, of highly variable alluvial deposits and soil textural types, the modeling of reservoir/groundwater interchanges is central to the planning of rehabilitation measures. The findings point out that several factors influence the interactions of reservoir with the groundwater system. The regulated recharge to the reservoir and the large spatio-temporal variability of interactions indicate that sophisticated and complex modeling methods are required to depict such a complex hydrological system. Therefore, the integration of a spatially-distributed water balance model (WetSpass-M), and a groundwater flow model (MODFLOW-NWT) to represent the The findings point out that several factors influence the interactions of reservoir with the groundwater system. The regulated recharge to the reservoir and the large spatio-temporal variability of interactions indicate that sophisticated and complex modeling methods are required to depict such a complex hydrological system. Therefore, the integration of a spatially-distributed water balance model (WetSpass-M), and a groundwater flow model (MODFLOW-NWT) to represent the volumetric exchanges between reservoir and groundwater was applied for case studies in the Drava floodplain.

Conclusions
The regulation works and incision of Drava River have changed the water regime and water budget of the floodplain, and, as a consequence, the water stress has increased and threatens the sustainable development of different sectors in the Drava basin. This study assessed the feasibility of the natural reservoir in augmenting the water resources of the floodplain through surface water storage and to mitigate water shortage problems caused by drought and human interventions using MODFLOW-NWT with ModelMuse as a graphical interface. The transient groundwater flow model was calibrated with a correlation coefficient of 0.96, a mean error of −0.016 m, and a mean absolute error of 0.4 m. A digital elevation model and hydraulic soil properties revealed by auger samples were used to identify possible allocations of a natural reservoir. The two selected subareas were selected as potential natural reservoir sites. Simulated reservoirs had filling and discharge conditions based on +2 m vertical water depth with different discharges from adjacent streams at 0.5, 0.75, 1, and 1.5 m 3 s −1 rates. For both reservoirs, around 60 days and 25 days of stream feeding of 0.5 and 0.75 m 3 s −1 , respectively, are required to achieve the target water depth for the reservoir (+2 m), the optimal and economic scenario of filling the reservoirs is with a stream discharge of 1 m 3 s −1 for 14 days, based on the capacity of the feeding stream. Feeding at a rate of 1.5 m 3 s −1 would need around just six days, but this is not an economic scenario as it necessitates impoundment and dam construction. With feeding at 1 m 3 s −1 , filling would be repeated twice per year, as the reservoir lasts for 180 days. The effect of the reservoir on the groundwater table in the adjacent area depends on reservoir stage, leakage, and the distribution of the groundwater head using the reservoir results in raising the groundwater table by 0.8 m at the end of the filling period. Natural reservoirs can be regarded as an economically feasible management practice of water resources retention and can help protect ecological values, the restoration of wetlands, and enhance agriculture productivity on the Hungarian section of the Drava floodplain. Moreover, the findings of this research can be utilized in planning rehabilitation measures in lowlands with water shortage.