Interactions of artificial lakes with groundwater applying an integrated MODFLOW solution

Hydrogeology Journal, Sep 2017

Artificial lakes (reservoirs) are regulated water bodies with large stage fluctuations and different interactions with groundwater compared with natural lakes. A novel modelling study characterizing the dynamics of these interactions is presented for artificial Lake Turawa, Poland. The integrated surface-water/groundwater MODFLOW-NWT transient model, applying SFR7, UZF1 and LAK7 packages to account for variably-saturated flow and temporally variable lake area extent and volume, was calibrated throughout 5 years (1-year warm-up, 4-year simulation), applying daily lake stages, heads and discharges as control variables. The water budget results showed that, in contrast to natural lakes, the reservoir interactions with groundwater were primarily dependent on the balance between lake inflow and regulated outflow, while influences of precipitation and evapotranspiration played secondary roles. Also, the spatio-temporal lakebed-seepage pattern was different compared with natural lakes. The large and fast-changing stages had large influence on lakebed-seepage and water table depth and also influenced groundwater evapotranspiration and groundwater exfiltration, as their maxima coincided not with rainfall peaks but with highest stages. The mean lakebed-seepage ranged from ~0.6 mm day−1 during lowest stages (lake-water gain) to ~1.0 mm day−1 during highest stages (lake-water loss) with largest losses up to 4.6 mm day−1 in the peripheral zone. The lakebed-seepage of this study was generally low because of low lakebed leakance (0.0007–0.0015 day−1) and prevailing upward regional groundwater flow moderating it. This study discloses the complexity of artificial lake interactions with groundwater, while the proposed front-line modelling methodology can be applied to any reservoir, and also to natural lake interactions with groundwater.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

https://link.springer.com/content/pdf/10.1007%2Fs10040-017-1641-x.pdf

Interactions of artificial lakes with groundwater applying an integrated MODFLOW solution

Interactions of artificial lakes with groundwater applying an integrated MODFLOW solution A. A. El-Zehairy 0 1 2 3 M. W. Lubczynski 0 1 2 3 J. Gurwin 0 1 2 3 0 Irrigation and Hydraulic Engineering Department, Faculty of Engineering, Mansoura University , Mansoura 35516 , Egypt 1 Department of Water Resources, Faculty of Geo-Information Science and Earth Observation (ITC) University of Twente , P.O. Box 217, 7500 AE Enschede , The Netherlands 2 M. W. Lubczynski 3 Department of Applied Hydrogeology, Institute of Geological Sciences, University of Wroclaw , Wroclaw , Poland Artificial lakes (reservoirs) are regulated water bodies with large stage fluctuations and different interactions with groundwater compared with natural lakes. A novel modelling study characterizing the dynamics of these interactions is presented for artificial Lake Turawa, Poland. The integrated surface-water/groundwater MODFLOW-NWT transient model, applying SFR7, UZF1 and LAK7 packages to account for variably-saturated flow and temporally variable lake area extent and volume, was calibrated throughout 5 years (1-year warm-up, 4-year simulation), applying daily lake stages, heads and discharges as control variables. The water budget results showed that, in contrast to natural lakes, the reservoir interactions with groundwater were primarily dependent on the balance between lake inflow and regulated outflow, while influences of precipitation and evapotranspiration played secondary roles. Also, the spatio-temporal lakebed-seepage pattern was different compared with natural lakes. The large and fast-changing stages had large influence on lakebed-seepage and water table depth and also influenced groundwater evapotranspiration and groundwater exfiltration, as their maxima coincided not with rainfall peaks but with highest stages. The mean lakebed-seepage ranged from ~0.6 mm day−1 during lowest stages (lake-water gain) to ~1.0 mm day−1 during highest stages (lake-water loss) with largest losses up to 4.6 mm day−1 in the peripheral zone. The lakebed-seepage of this study was generally low because of low lakebed leakance (0.0007-0.0015 day−1) and prevailing upward regional groundwater flow moderating it. This study discloses the complexity of artificial lake interactions with groundwater, while the proposed front-line modelling methodology can be applied to any reservoir, and also to natural lake interactions with groundwater. Dynamics of artificial lakes; Groundwater/ surface-water relations; Numerical modelling; Water balance; Poland Introduction In most places on the Earth, groundwater and surface water are in a continuous dynamic interaction that can affect not only water quantity but also the quality of both, the groundwater and surface-water resources (Sophocleous 2002) . Therefore, water resources management moves nowadays towards the challenge of integrating groundwater and surface water as one management unit (Ala-aho et al. 2015) . This paradigm shift requires good understanding of the complexity of interactions between surface-water and groundwater and reliable methods to simulate these interactions. Among surface-water/groundwater interactions, probably the most distinct and complex are those between artificial lakes—referred to also as “reservoirs”, following the Chapman (1996) terminology—and groundwater. Artificial lakes are designed to store surface water, so their beds are expected to have low permeability and low seepage. However, they are affected by natural and particularly strong artificial (due to lake regulation) driving forces, which imply large and fast changing lake-water levels (lake-stages), implying large impact on groundwater dynamics. That impact may result, for example, in flooding of an area adjacent to a reservoir also affecting agricultural productivity or contamination of that area (Wildi 2010) . For all these reasons, the management of reservoirs and adjacent areas requires appropriate methodologies and models well accounting for surface-water/groundwater interactions. An example of such methodology, based on modelled simulation of interactions between an artificial lake and groundwater in the adjacent area, is proposed in this study. The simplest and most direct way of studying a lake interaction with groundwater, is by investigating exchange of seepage between that lake water and groundwater by point field measurements using methods such as seepage meters, thermic profiling, dye tracing, piezometers and potentio-manometers (Lee 1977; Ong and Zlotnik 2011; Otz et al. 2003; Winter et al. 1988) or a combination of methods (Anibas et al. 2009; Owor et al. 2011; Su et al. 2016) . However, such methods are vulnerable to heterogeneity and therefore difficult for spatial integration. Hydrochemical analysis and mass balance methods (Sacks et al. 1998; Stauffer 1985) are better in that respect but require a conservative tracer not influenced by land-use practices such as for example stable isotopes (Brindha et al. 2014; Kanduč et al. 2014; Krabbenhoft et al. 1990; Sacks et al. 2014) ; additionally, convenient semianalytical water balance models (Ghosh et al. 2015; Rudnick et al. 2015) involve quite a number of assumptions simplifying spatial heterogeneity. The most complete and reliable assessment of lake–groundwater interaction is by increasingly sophisticated physically based distributed and integrated hydrological models, particularly if based on solid characterization of hydrogeological conditions of a lake and its adjacent area and time series data, as proposed in this artificial lake– groundwater interaction study. The general scarcity of physically based distributed models that simulate interactions of lakes with groundwater is surprising, but even more surprising is that, till now, to the authors’ knowledge, there has not been any such study published on the interaction of artificial lakes with groundwater. Some of the few available examples of artificial lake studies, but with a different focus and methods applied than in this study, include the following: Liuzzo et al. (2015) , who used a lumped conceptual model (TOPDM) to characterize the effect of climate change on water resources availability in the Belice River Basin (Italy), in which artificial Lake Garcia is located; Fowe et al. (2015) , who used a genetic algorithm model to couple Boura reservoir (Burkina Faso) with the adjacent groundwater system and to optimize water allocation; and Chhuon et al. (2016) , who coupled the semi-distributed hydrologic SWAT model with the MODSIM decision support system to investigate hydrological consequences of a future upland reservoir on the Prek Te River, Cambodia. Simulations of interactions of natural lakes with groundwater have been carried out for decades, mainly by the widely used MODFLOW (McDonald and Harbaugh 1988) type of solutions, also applied in this study. The simplest MODFLOW lake implementations include representing the lake as a head (stage) dependent sink or source through the River Package (Leblanc et al. 2007) , General Head Boundary (GHB) Package (Mylopoulos et al. 2007) or Reservoir Package (Fenske et al. 1996) . The disadvantage of these approaches is that the simulated lake stages constrain lakebed seepages, thus a prior knowledge of their seepage amounts is required (Merritt and Konikow 2000) instead of calculating them internally. Another approach that requires cells of a model grid representing lake volume is called the “high-K” method (Lee 1996) , where extremely high hydraulic conductivity and a storage coefficient equal to 1 are used to simulate lake cells. The main disadvantages of the “high-K” method are that it is difficult to represent accurately the connections between streams and a lake (Merritt and Konikow 2000) and that the method is prone to numerical instability (Yihdego and Becht 2013) . A prototype of the current MODFLOW Lake Package (LAK1) solution that overcomes most of the aforementioned disadvantages was developed by Cheng and Anderson (1993) . The most important advantage of the Lake Package is that the lake stages are not defined as hydraulic boundaries but determined as part of the simulation. Besides, it allows for fluctuating lake levels and accordingly calculates volumetric water exchange between a lake and an aquifer following the hydraulic gradient and explicitly defined lakebed conductances. The LAK1 Package was also integrated with a primary version of the Stream Flow Routing (STR1) Package (Prudic 1989) regulating water inflow to the lake and outflow from the lake to streams. However, as compared to the current Lake Package, LAK1 did not take into account that some of the lakebed cells could become dry due to low lake-stages. Also, it did not consider the flow resistance within an aquifer, which might be substantial in magnitude and even larger than the lakebed resistance in the case of thin lakebeds (Merritt and Konikow 2000) . The most important addition in the next, LAK2 Package was its capability of simulating more than one lake. The largest improvement has taken place in the LAK3 Package (Merritt and Konikow 2000) reviewed by Hunt (2003) , where a lake is represented as a volume of space within the model grid which consists of inactive cells extending downward from the upper surface of the grid. Active model grid cells bordering this space and representing the adjacent aquifer, exchange water with the lake at a rate determined by the relative heads and by conductances that are based on grid cell dimensions, hydraulic conductivities of the aquifer material, and user-specified leakance distributions that represent the resistance to flow through the material of the lakebed. In the LAK3 Package, also, the simulation of solute transport was added. The LAK3 Package was validated through seepage meter measurements by Kidmose et al. (2011) and applied in a number of studies such as, for example, in Vaeret et al. (2009) . Another important improvement took place recently in the most up-to-date LAK7 Package, also used in this study. In that LAK7 Package, water budget can be more realistically simulated than before, thanks to its integration with the Stream Flow Routing (SFR7; Niswonger and Prudic 2005) and the Unsaturated Zone Flow (UZF1; Niswonger et al. 2006) packages of MODFLOW-NWT (Niswonger et al. 2011) , the latter package accounting for variably-saturated flow during lake-area expansion or shrinkage and replacing former the Recharge and Evapotranspiration packages of original version of MODFLOW (McDonald and Harbaugh 1988) . If interactions between a lake and groundwater are modest and with low temporal variability, all modelling solutions reviewed in the preceding text are likely applicable. However, if, as in artificial lakes, the lake stages are driven by complex, simultaneously operating natural and maninduced (lake outflow regulation) driving forces resulting in fast-changing lake stages, lake area extent, lake-water volume, hydraulic gradients and seepages, then a more sophisticated solution is needed such as integrated (fully coupled) hydrological models with variably-saturated flow, calibrated in a transient mode, which reduces more degrees of freedom than steady-state calibration (Lubczynski and Gurwin 2005) , as proposed in this study. To the authors’ knowledge, such a solution to investigate the dynamics and water budget of interactions between an artificial lake and the groundwater system adjacent to the lake, including the dynamics and pattern of lakebed seepage, has never been published, but is highly recommended for artificial lake management. Therefore, the main objectives of this study are: (1) to investigate the complexity of the effect of artificial lake stage changes upon lake– groundwater seepage exchange and hydraulic heads in the area adjacent to the lake and (2) to quantify such interactions in a spatio-temporal manner, separating natural and maninduced impacts. To address these objectives, artificial Turawa Lake (TL) in Poland was selected as a test case, mainly due to its typicality for artificial lakes, frequent and large stage fluctuations (up to ~5 m) and also because of its extensive, available database, including long time-series records of groundwater levels, lake stages, river discharges, lakebed measurements and pumping/slug tests. The proposed methodology in this study was applied to TL, although it can be extrapolated to any similar case of not only artificial but also natural lakes, so it can be considered as a practical guideline for lake management. The main novelties of this study are in: (1) first time use of a distributed, integrated, multi-layered hydrological model for the daily transient calibration of artificial lake interactions with groundwater and streams, applying volumetric water exchange, accounting for temporally variable lake area and variably-saturated water flow; (2) identification and separation of natural and human-induced impacts upon the hydrological system of the artificial lake; (3) first time assessment of the spatio-temporally variable lakebed seepage pattern of the artificial lake; (4) realistic water balancing of interactions between the artificial lake, groundwater and streams where: (a) lake stages are calculated by the model, not predefined; (b) the lake area extent is dependent on lake stage so temporally variable, implying temporally variable lake/land area ratio; (c) lake storage and related lake-water exchanges with streams (except stream inflows at the boundaries and lake regulated outflow) and groundwater are calculated by the model, not predefined; and (d) estimates of recharge, groundwater evapotranspiration and groundwater exfiltration are made applying the variably-saturated flow concept, which allows one to account for variable lake area extent and related spatially variable saturation of shallow subsurface. Data and methods Study area The Turawa Lake (TL), located in the south of Poland (Fig. 1), is an artificial lake with an average area of ~13.5 km2. The TL was constructed on the Mala Panew (MP) River, which is the right tributary of the Odra River (Fig. 1). The earthen dam of the TL (Fig. 1) was constructed in 1938. The TL is used for flow regulation downstream of the MP River, flood prevention, electric power generation and for touristic purposes (Simeonov et al. 2007) . To satisfy all these purposes, the TL management is achieved through outflow regulation implying temporal variability of the lake stages and lake area extent. The lake is shallow in the eastern part and its depth increases gradually to the west, towards its deepest point (~10 m depth, during the highest lake stage), located near the dam, at the western outlet of the lake into the MP River. The TL and the adjacent area have a good database developed within the “Odra-2006” regional cooperation project (Gurwin et al. 2004) and afterwards carried out by the University of Wroclaw (Gurwin 2010). The TL has been heavily contaminated by industrial sediments from the steel works factory in Ozimek and from the Upper Silesia coal mining industry, both located upstream of the lake in the SE direction (outside the area displayed in Fig. 1). This industrial sediment has accumulated at the lake bottom in the form of “sapropel”, which is a dark-coloured contaminated sediment of a low permeability. Besides this, there are some other organic and non-organic sediments. In 2004, the Wroclaw University investigated (Gurwin et al. 2004) spatial extent, thickness and chemical composition of the lake bottom sediment; the regular sampling schema is shown by dots within the lake area in Fig. 1. The ~100 km2 TL catchment area, selected as the study area, is presented in Fig. 1. From the north, the area is bounded by the MP catchment boundary, from the south by internal topographic divide and from the east and the west, by artificially defined boundaries. The study area has a slightly rolling, post-glacial topography with gentle hills. The majority of the study area is covered by Pinus silvestris L forest. The climatic data were recorded hourly by the hydrometeorological automated data acquisition system (ADAS; Figs. 1, 2, and 3) operating from 2003 until 2010. The ADAS installation included monitoring of precipitation, incoming solar radiation, relative humidity, wind speed, air temperature and soil heat flux. As the data set had some gaps, a gap-filling process was undertaken, using data from “Opole” meteo-station located 15 km south-west of the Turawa Lake and retrieved from NNDC climate database, “CLIMVIS” (NOAA 2004) . The study area is characterized by temperate climate with relatively cold winters (mean, min, max temperatures: −2, −19, +7 °C respectively) and warm summers (mean, min, max temperatures: 19, 10, 29 °C respectively). Large daily precipitation variability can be observed (Fig. 2) particularly during summers (June, July, Aug) characterized by frequent rain showers sometimes exceeding even 20 mm day−1. Winters (Dec, Jan, Feb) are characterized with snowfall, while during springs (March, April, May) and autumns (Sept, Oct, Nov), precipitation variability is moderate. The daily lake evaporation was computed using the Penman open-water equation (Penman 1948) on the base of hourly incoming solar radiation, relative humidity, wind speed, air temperature and using the wind-speedfunction constant values recommended by McMahon et al. (2013) . The daily lake evaporation ranged from Fig. 2 Daily precipitation (upper position) and river inflows (lower position) 0.0 to 8.2 mm day−1 (Fig. 3), with small values attributed to winter and large values to summer months. For the land surface (total model area minus lake area), the daily potential evapotranspiration (PET) was estimated using hourly incoming solar radiation, relative humidity, wind speed, air temperature, and soil heat flux following the McMahon et al. (2013) definition of potential evapotranspiration considered as “the rate at which evapotranspiration would occur from a large area completely and uniformly covered with growing vegetation, which has access to an unlimited supply of soil water, and without advection or heating effects”. In line with that definition, first, the reference evapotranspiration (ETo) using the “FAO Penman Monteith” method (Allen et al. 1998) was estimated, and then, by using the “single crop coefficient” method, the ETo was converted to daily crop potential evapotranspiration considered as potential evapotranspiration (PET). The estimated PET, used as a model input, ranged from 0.0 to 6.9 mm day−1 (Fig. 3) with small values attributed to winter and large values to summer months. The main water input to the Turawa Lake is from the MP River entering the lake from the south-eastern side (Figs. 1 and 2) gauged at the “Staniszcze Wielkie” (14.5 km upstream of the lake – not shown). Another two small rivers, Libawa and Rosa, and some drainage ditches dewatering the land surface adjacent to the lake, supply relatively little water to the lake (Fig. 2). Besides, the lake receives water from precipitation, direct runoff and groundwater inflow, the latter occurs mainly during low stages. The lake has only one outlet into the MP River through sluice gates of the hydroelectric power plant in the earthen dam, where the outflow measurements are carried out at the Turawa station. The lake also Fig. 3 Daily lake evaporation and land-surface potential evapotranspiration (both, upper position), and regulated outflow from the lake (lower position) discharges water by evaporation and lake seepage to groundwater, the latter mainly at high and average lake stages. The daily upstream and downstream river discharges were obtained from the Institute of Meteorology and Water Management-National Research Institute (IMGW-PIB) in Poland for the period from 1 November 2003 to 31 October 2008 (time selected as this study period). The temporal pattern of the natural river inflow is obviously quite different compared with the regulated outflow (compare Figs. 2 and 3). The most distinct in these patterns are large spring inflow peaks due to snow melting that are dumped by lake-water storage and outflow regulation. The large variability of the lake-water input and the management of the lake-water storage, imply periodic cumulation of water and its fast removal resulting in large and frequent watertable fluctuation. The large amplitude of the lake stages of ~5 m, as recorded within the period of this study, implied also large changes of the lake area extent, ranging from 12.10 to 16.12 km2. These area changes could have been quite well defined thanks to the detailed 5 m resolution digital elevation model of the whole catchment area and the lake bathymetry measured with sonar (Gurwin et al. 2004) . The hydrogeological investigations in the study area (Gurwin et al. 2004) focussed on the three-layer Quaternary system composed of two aquifers and an aquitard in between, underlain by impermeable Tertiary and/or Triassic sedimentary deposits (Figs. 4 and 5). The three-layer Quaternary sequence represents glacial outwash deposits. The upper unconfined aquifer, composed mainly of fine, medium, and coarse sand with boulders, gravel, and a little clay and loam, has variable thickness ranging from 1.5 to 35.0 m. The aquitard is composed of loam with some sandy clay, clay, and argillaceous material, and has thickness ranging from 1.0 to 15.5 m. The lower confined aquifer, composed mainly of gravel and sand, has variable thickness ranging between 0.7 and 36.6 m. The two aquifers’ lithology, thickness and hydraulic conductivities are available from 45 boreholes with pumping tests, most of them available for the deep confined aquifer, and from 16 slug tests carried out in the shallow unconfined aquifer in addition to some other tests outside the study area. The hydraulic conductivities of the upper aquifer ranged from 0.4 up to 44.8 m day−1 with a geometric mean of 6.2 m day−1, while of the lower aquifer, from 5.2 to 41.6 m day−1 with a geometric mean of 12.3 m day−1. In addition, a number of investigation boreholes were drilled to determine the top and bottom of each layer. The available time-series of piezometric heads and lake stages obtained from Wroclaw University comprised hourly automated logging data and manual weekly and quarterly instantaneous observations from regular measurement campaigns (Table 1). Conceptual model Based on borehole logs, pumping tests, slug tests and available hydrogeological and geophysical information, the Fig. 6 Model grid, boundary conditions and mean hydraulic head distribution groundwater flow system of the study area can be conceptualized as a system of three layers (Figs. 4 and 5), the upper unconfined aquifer, the middle aquitard and the lower confined aquifer underlain by impermeable basement. The three layers are assumed to be spatially heterogeneous and anisotropic. The regional groundwater flow in both analysed aquifers is from SE towards NW and follows the direction of the MP River that in general drains groundwater (Fig. 6). The groundwater flow system of the upper aquifer is under the influence of the lake-stage fluctuations, i.e. at the low stages, the lake drains groundwater but at the high stages, the lake induces seepage into groundwater system. The variability of the lake stages affects not only the lake seepage but also the leakage through the aquitard that in general is in upward direction, except in a few areas in the surrounding of the lake where aquitard thickness is reduced and mainly during low lakestages. The sources of water input to the groundwater system are: recharge from precipitation, lake seepage, stream seepage, and lateral groundwater inflow across the eastern and south-eastern boundaries of the study area. The sources of water output from the groundwater system are: groundwater evaporation, lateral groundwater outflow across the western boundary, groundwater seepage into the lake and streams, and groundwater exfiltration to the land surface. Numerical model The applied MODFLOW-NWT model (Niswonger et al. 2011) is particularly suitable for problems such as in this study, i.e. caused by drying and rewetting nonlinearities of the unconfined groundwater-flow equation, where groundwater heads are calculated for dry cells even if below the cell’s bottom, in contrast to other MODFLOW versions, where the dry cells are excluded from the calculation. It also improves the model convergence and computational efficiency when using nonlinear boundary conditions in LAK7, SFR7 and UZF1 packages. The LAK7 Package is designed to simulate lake–groundwater interactions and it is particularly suitable for lakes with significantly changing stages and spatial extent, because the lake stages are computed based on the volumetric water exchanges into and out of the lake and on the overall lake-water balance. Seepage between a lake and the adjacent aquifer system is quantified using Darcy’s Law (Merritt and Konikow 2000) , i.e. based on lake stage, hydraulic heads in the groundwater system, and conductances that are dependent on grid cell dimensions, hydraulic conductivities and user specified leakance distribution; inverse of the latter represents the resistance of seepage across the lakebed. The SFR7 Package (Niswonger and Prudic 2005) is designed to simulate volumetric river interactions with groundwater. It allows to implement river discharge measurements in the model and is fully integrated with LAK7 in the MODFLOW-NWT model, allowing stage-dependent, volumetric water exchanges between rivers and lakes. The UZF1 Package simulates vertical one-dimensional (1D) variably-saturated flow between land surface and the water table, applying the kinematic-wave approximation of Richard’s Equation (Niswonger et al. 2006) . In the UZF1 Package, the relation between unsaturated hydraulic conductivity and the unsaturated zone water content is defined by the Brooks and Corey function (Brooks and Corey 1966) . The UZF1 Package is fully integrated in MODFLOW-NWT, not only with SFR7, but also with LAK7, which allows the infiltration rates to be applied not only for the land area but also for the lake cells converted to dry cells when the lake shrinks. It also allows for computing unsaturated flow beneath the lake area if the groundwater head drops under the lake bottom. In standard MODFLOW transient model simulations, stress periods used to be assigned on the basis of temporal variability of: external sink/sources, state variables such as groundwater levels and/or discharges and input/output groundwater fluxes, namely gross recharge (Rg) and groundwater evapotranspiration (ETg). Such a solution was arbitrary, not flexible and often inaccurate because it did not account for water flux variability within each arbitrarily assigned stress period and had an arbitrary assumption of water fluxes, e.g. recharge (Hunt et al. 2008) . In the proposed methodology, there are no procedural limitations regarding temporal resolution of stress periods and time step lengths, while Rg and ETg are calculated internally in the UZF1 Package for each time step together with groundwater exfiltration (Exfgw, sometimes also reffered as surface leakage) based on external driving forces, i.e. rainfall and potential evapotranspiration and unsaturated zone parameterization. The proposed software combination was run under the ModelMuse modelling environment (Winston 2009) , selected because: (1) at the time of this study, it was the only pre- and post-processor that could run MODFLOW-NWT with LAK7, SFR7 and UZF1 packages; (2) it is easy and straightforward in use; (3) it is public domain software. For water balances of each layer or specific part of the model, the ZONEBUDGET (Harbaugh 1990) option under ModelMuse was used. Model setup Following the conceptual model, three MODFLOW-NWT model layers were used to simulate 3D flow through the system of two aquifers separated by an aquitard and one extra inactive layer at the top, to simulate the lake using LAK7 Package as explained in the following. The model topographic surface was assigned using the available digital elevation model and each layer was defined by subtraction of the interpolated layer thicknesses using all the available borehole data. The upper unconfined aquifer and the aquitard layer were assigned as convertible between unconfined and confined conditions, while the lower aquifer was assigned as confined. A quadratic 200 × 200 m grid, consistent with the PUWG_92 Polish coordinate system, was used. That grid size was assigned after several tests applying different grid sizes, as a compromise between model accuracy and retardant in model calibration, use of finer grid, implying extensively long model simulations. The three layers were parameterized using internally homogenous horizontal hydraulic conductivity zones (K-zones) b a s e d o n p u m p i n g t e s t s , s l u g t e s t s a n d a v a i l a b l e hydrogeological knowledge. The assigned horizontal hydraulic conductivity (Kh) for the upper and lower aquifers varied from 0.4 to 44.8 m day−1 and from 5.2 to 41.6 m day−1 respectively. The vertical hydraulic conductivity (Kv) was assigned through an anisotropy factor Kh/Kv equal to 10, following local and general literature information about glacial outwash deposits (Gurwin et al. 2005; Smerdon et al. 2007) . The specific yield (Sy) and the specific storage (Ss) were estimated based on borehole sample analysis and assigned as spatially uniform, for the first unconfined aquifer as Sy = 0.15, for the aquitard and for the second confined aquifer as Ss = 0.00001 m−1. Although all the hydraulic parameters were defined based on the field measurements, due to the modelling scale effect ( Guimerà et al. 1995 ; Zhang et al. 2006 ; Zhang et al. 2007 ) and limited spatial representativeness, they were still considered to be adjusted in the calibration process. The external boundary conditions, the same for the entire vertical extent of the model, were of two types (Fig. 6): (1) noflow boundaries along watershed divides at the N and SW boundaries, assigned based on the match between surface and groundwater divides; the same boundaries were also assigned in the lower layers, i.e. in the direction parallel to the regional groundwater flow matching the flow direction of the MP river; and (2) head dependent flow boundaries assigned using MODFLOW General Head Boundary (GHB) Package (McDonald and Harbaugh 1988) , in the E and SE, to simulate lateral groundwater inflow into the study area, and in the W, to simulate lateral groundwater outflow from the study area; the GHB heads were assigned on the basis of available piezometric heads, while the GHB conductances, initially assumed as 1.5 m2 day−1 based on the Kh measurements, were to be further adjusted in the model calibration. In case there was no time series of head measurement at a GHB boundary, instantaneous records at that boundary were used to create time series, forcing the same pattern (phase and amplitude) as the nearest time series piezometer but with the water table matching the instantaneous record. The diaphragm under the earthen dam was simulated using the Horizontal Flow Barrier (HFB) Package (Hsieh and Freckleton 1993) . The barrier thickness was assigned as 3 m thick, but its hydraulic conductivity was adjusted during the calibration process. The Turawa Lake (TL) was implemented in the model through the LAK7 Package. The lake was simulated by a separate additional (fourth) top layer, represented by inactive cells (Merritt and Konikow 2000) of variable thickness within the maximum lake extent, with the top equal to the maximum measured lake-stage and bottom equal to the bathymetrically defined lake bottom. Outside the maximum lake extent, 1.0 m thick inactive layer was applied. To represent the lakebed leakance, two uniform zones were assigned based on field measurements of the lakebed thicknesses and vertical hydraulic conductivities: 0.0007 day−1 for the western part covered with sapropel and 0.0015 day−1 for the remaining lakebed area in the eastern part. The main driving force of the model was the difference between the MP River inflow and the MP River outflow (regulated by the dam-weir). The rivers and drains were simulated by the SFR7 Package, all with rectangular cross sections and streambed thicknesses of 1.0 m. The measured daily river inflows at the external model boundary and the measured daily regulated lake outflows to the MP River, just downstream of the lake, were assigned as model inputs. The Manning coefficient for all streams was assumed as 0.035, while the assigned hydraulic conductivity of the streambeds was 0.1 m day−1, which was to be adjusted during the calibration process. The MP River width upstream of the lake was set equal to 20 m and downstream to 30 m. The widths of Libawa and Rosa rivers and all the drains were set equal to 10, 10 and 5 m respectively. As MODFLOW-NWT under the ModelMuse environment did not allow the lake outlet of the MP River to have a streambed elevation lower than the lakebed elevation, the MP River downstream of the lake was disconnected from the lake and the lake outflow measured at the Turawa gauging station downstream of the lake was assigned as free overfall input into the MP River. In the proposed modelling method, the variable in time “precipitation falling on the lake area” is directly added to the current lake volume while the “precipitation falling on the land part of the lake catchment” is partitioned in the UZF1 Package into interception, infiltration and surface runoff (RO) routed to selected streams or lakes, depending on the user choice, if either the SFR7 or LAK7 Package is active. The difference between precipitation and interception is called “applied infiltration rate”. If the applied infiltration rate is higher than the soil saturated vertical hydraulic conductivity (Kv), in addition to actual infiltration (Pa) through the soil, there is excess infiltration runoff (ROei) also known as Hortonian runoff, while if the water table goes higher than the land surface (groundwater exfiltration), then saturation excess runoff (ROsat), also known as Dunnian runoff, takes place (Niswonger et al. 2006) . The ROsat can include two components, rejected infiltration due to shallow groundwater levels, and groundwater exfiltration (Exfgw) to land surface. The UZF1 computed runoff, including groundwater exfiltration, was selected to be routed to the lake. The interception was estimated as a percentage of precipitation based on the available land use map and literature estimates of interception for the land use types available in the study area. For the dominant P. silvestris L forest land cover, interception equal to 27.3% of precipitation was assigned (Wang et al. 2007) . In the UZF1, the Pa is divided into gross groundwater recharge (Rg) and unsaturated zone evapotranspiration (ETuz) and the remaining is a change of unsaturated zone storage (ΔSuz). First, ETuz is removed from the unsaturated zone at the PET rate and if the evapotranspiration demand is not met, water is removed further from groundwater (as ETg), but only if the water table is above the predefined evapotranspiration extinction depth assigned equal to 2.0 m, in accordance with the root depth of the dominant tree species. The computed evapotranspiration rates depend also on the amount of water stored in the unsaturated zone above the predefined evapotranspiration extinction water content assigned as residual water content of 0.05 (m3 m−3); that value was also used as the initial soil water content, while the soil saturated water content was assumed as 0.30 (m3 m−3). The UZF1 Package was activated for the land area and below the lake, which allowed the applied infiltration rates to be estimated for the land area and for the lake cells converted into dry cells when the lake shrank. Model calibration and sensitivity analysis Initially, a steady-state model was developed to provide the initial condition for the transient model by simulating the long-term average water balance of the modelled system applying mean driving forces, i.e. mean precipitation, mean PET and mean state variables, i.e. mean groundwater levels and mean lake stage (173.8 m a.s.l.) for the period from 1 Nov. 2003 to 31 Oct. 2008. However, the initialization of the transient model using that steady-state model as a starting condition was unsuccessful because of a large difference between the steady-state river inflows/outflows and the measured daily discharge values at the start-up of the transient model on 1 November 2003 (later modified to be 1 November 2004). Moreover, the UZF1 Package works differently in steady-state as compared to transient conditions (Niswonger et al. 2006) , which resulted in different calibrated parameters—e.g. horizontal hydraulic conductivities (Kh), lakebed leakance, streambed hydraulic conductivities, etc. As a consequence, the steady-state model was abandoned and the transient model calibration was initialized with a warm-up (spin-up) period (Hassan et al. 2014) of 70 days, followed by the normal simulation period. However, that 70 days period was too short, as it resulted in systematic error, i.e. divergence of the simulated lake stages from the observed values. Therefore, the whole first hydrological year of the available data, i.e. from 1 November 2003 until 31 October 2004, was “sacrificed” as a warm-up period, after which the model was run (calibrated) in transient conditions throughout another four hydrological years, i.e. from 1 November 2004 to 31 October 2008. The time domain of the final transient model was discretized by assigning daily stress periods, each including one single-day time step. The model was run using the Newtonian (NWT) solver, applying the option of calculating groundwater heads “even if below cell bottoms” in the case of drying cells. The final solver head tolerance was adjusted to 0.005 m and the flux tolerance to 1,000 m3 day−1. The model complexity was set as “complex” and all the remaining solver criteria were accepted by default. The model was calibrated manually in forward mode because of its large complexity implying long simulation time when using optimisation codes such as PEST (Doherty and Hunt 2010) or UCODE (Hill and Tiedeman 2007), but also because forward calibration allows the user to understand better the model behaviour (Hassan et al. 2014) . The transient model calibration targets were to minimize the: (1) root mean square error (RMSE) of the differences between the simulated and measured groundwater heads; (2) RMSE of the differences between the measured and simulated lake stages; and (3) water balance discrepancies at each time step. The calibration process was done mainly by adjusting the number of initially assigned K-zones, their areas and the associated hydraulic conductivities (Kh). The same zones were used to adjust the specific yield (Sy) and specific storage (Ss) values; furthermore, some minor changes were made in the initially assigned hydraulic conductivities of both streambeds and the HFB, and, also, GHB conductance at the outflow boundary was slightly adjusted. The sensitivity analysis of the transient model involved assessment of: (1) the lakebed seepage variations due to changes in lakebed leakance, due to changes in hydraulic conductivity K (Kh and its corresponding value of Kv assigned as 0.1Kh) of the shallow aquifer, and due to changes in the lake inflow; (2) the magnitudes of the excursion of groundwater heads and lake stages from the calibrated values as a result of changes in: horizontal hydraulic conductivity (Kh) and specific yield (Sy) of the upper aquifer, specific storage (Ss) of the lower aquifer, the anisotropy factor (Kh/Kv), lakebed leakance, and river inflows to the lake. Specific tests were also performed in order to characterize the spatial extent of the impact of the lake stage fluctuation upon the groundwater heads in adjacent areas. For that purpose, two transects of fictitious piezometers (Fig. 6) were assigned in the upper aquifer. In each transect, the first fictitious piezometer was located 100 m away from the lake shoreline, and the distance to each next piezometer increased sequentially by 200 m. The northern transect was selected mainly to test the validity of the northern no-flow boundary assigned along the watershed divide, matching the groundwater divide and located pretty close to the lake (~1,000 m), while the southern transect was selected mainly to test the extent of the lake impact in the direction towards the GHB inflow boundary. Regarding the proximity of the northern noflow boundary from the lakeshore line, an additional experiment was carried out on the northern transect, to test the validity of that no-flow boundary. For that purpose, the northern boundary was experimentally moved 800 m to the north, i.e. up to 1,800 m distance from the shoreline, and replaced with a GHB. The new model was recalibrated and the effect of the boundary shift was analysed. Water balance Water balancing of artificial lakes, particularly when simulated with variably-saturated, transient, integrated models, is a complex matter because of many interacting surface, unsaturated zone and groundwater components (as presented in Fig. 7), the temporally changing lake versus land area ratio (so also changing flux area reference), and interplaying impacts of natural and human-imposed influences. The following water balance equations for each Fig. 7 Schematic diagram of the model setup in Turawa Lake study area; symbols the same as in Eqs. (1)–(10) model domain are indispensable to understand and quantify interactions between an artificial lake and groundwater. All these equations, are presented as model water input (IN), equal to water output (OUT) plus change of water storage (ΔS). The water balance of the whole catchment model domain (74.8 km2) can be expressed as follows: P þ Qin þ GHBin ¼ ET þ Qout þ GHBout þ ΔS ð1Þ where P is precipitation rate, Qin is stream inflows at the inlet of the modelled area, Qout is stream outflows at the outlet of the modelled area, GHBin is lateral groundwater inflow into the modelled area across the GHB boundaries, GHBout is lateral groundwater outflow from the modelled area across the GHB boundaries, ET is total evapotranspiration consisting of land surface evapotranspiration (ETLD) and lake evaporation (ELK), and ΔS is total change of storage. The total evapotranspiration (ET) and total change in storage (ΔS) can be expressed as follows: ET ¼ ETLD þ ELK ¼ I þ ETuz þ ETg þ ELK ΔS ¼ ΔSg þ ΔSuz þ ΔSLK ð2Þ ð3Þ where I is canopy interception, ETuz is unsaturated zone evapotranspiration, ETg is groundwater evapotranspiration, ΔSg is storage change of saturated zone, ΔSuz is storage change of unsaturated zone and ΔSLK is the lake storage change. The lake-water balance extracted from the Lake Package output file can be expressed as follows: P þ QLKin þ LLKin þ RO ¼ ELK þ QLKout þ LLKout ΔSLK ð4Þ where QLKin is stream inflow (rivers and drains) at the inlet of the lake, QLKout is stream outflow at the outlet of the lake, RO is total surface runoff into the lake computed by the UZF1 Package as per Eq. (5), LLKin is seepage of groundwater into the lake, LLKout is seepage from the lake into groundwater. In that equation, RO is expressed as follows: RO ¼ ROei þ ROsat ð5Þ ð6Þ ð7Þ ð8Þ ð9Þ where ROei is excess infiltration runoff (Hortonian runoff) and ROsat is saturation excess runoff (Dunnian runoff) including groundwater exfiltration to land surface (Exfgw). The land surface and unsaturated zone water balance are expressed as follows: P þ Exf gw ¼ Rg þ ETuz þ I þ RO ΔSuz which can be further divided into the following two equations: P þ Exf gw ¼ Pa þ I þ RO Pa ¼ Rg þ ETuz ΔSuz where Pa is the actual infiltration through the unsaturated zone, and Rg is gross recharge. The saturated zone water balance of the two aquifers in addition to the aquitard layer can be expressed as follows: Rg þ GHBin þ QGWin þ LLKout ¼ ETg þ GHBout þ Exf gw þ QGWout þ LLKin ΔSg where QGWin is stream seepage to groundwater, QGWout is groundwater seepage to streams. Following Hassan et al. (2014) , the net recharge (Rn) refers to the net amount of water that reaches the saturated zone after deducting groundwater evapotranspiration (ETg) and groundwater exfiltration (Exfgw) as follows: Rn ¼ Rg −Exf gw−ETg Results and discussion ð10Þ To the authors’ knowledge, there is no other published work applying integrated hydrological modelling to investigate dynamics of interactions between artificial lakes and groundwater; therefore, the results of this study are presented, discussed and compared to selected natural (not artificial) lake studies where possible. However, one should be aware that because of the dam construction and its man-induced reservoir regulation, artificial lakes not only have different main driving forces of system dynamics, different lake stage amplitude and fluctuation frequency but also different lakebed geometry, composition, permeability and thickness because of its sealing prior to reservoir operation and more active sediment deposition, both reducing lakebed leakance and, as a consequence, the lakebed seepage. As a result, reservoirs have different dynamics of interactions with groundwater than natural lakes and this will be shown and emphasized hereafter. Calibration and sensitivity analysis The results of the transient model calibration against time series of piezometric heads and lake-stages are presented in Fig. 8a,b respectively. The temporal variation shown in the Fig. 8, is mainly due to combined effects of lake outflow regulation and lake inflow, but partly also due to the relatively minor effect of climatic factors. The calibration of heads, carried out on six timeseries of daily piezometric data extending over 4 years, shows a good agreement between simulated and measured groundwater levels. That agreement is reflected by correlation coefficient 0.99, mean error −0.15 m, mean absolute error 0.27 m and RMSE 0.36 m. The calibration was additionally controlled by matching simulated heads with corresponding weekly and quarterly instantaneous piezometric observations (Fig. 1; Table 1). The groundwater mass balance discrepancy for all the stress periods ranged from −0.13 to 0.78%, and the cumulative discrepancy value at the end of the 4-year simulation period was 0.02% with an average discrepancy value equal to 0.00%. The recorded small discrepancies in the simulated groundwater heads are most likely due to errors in model parameterization, heterogeneities of the Quaternary sediments, but possibly also due to some measurement gaps and eventual measurement errors in stream discharges, climatic forces, etc. The calibration of the lake stages (Fig. 8b) also shows a good match between simulated and observed stages: correlation coefficient 0.97, mean error 0.30 m, mean absolute error 0.45 m and RMSE 0.52 m. The small biases, in Fig. 8b, are likely because of daily measurements instead of finer, e.g. hourly measurements of the lake stages and lake inflows, in addition to potential small errors in the measurements of the latter. The calibration process resulted in 20 K-zones for the shallow aquifer, 11 for the deeper aquifer and 3 for the aquitard. As compared to initially assigned pumping-test/slug-test-based spatial Kh-distribution, the calibrated Kh of the upper unconfined aquifer increased gently in the north-eastern area from 2.5 to 5 m day−1 and in the central-north from 44.8 to 60 m day−1, while in the central-south it decreased from 9 to 0.1 m day−1. The final calibrated Kh of the upper aquifer ranged from 0.1 to 60 m day−1. Larger Kh changes as compared to initial pumping test assumptions, were made in the confined aquifer; in the northern area Kh decreased from 15 to 0.5 m day−1, in the eastern area from 5 to 0.2 m day−1, in the western area from 11 to 2 m day−1, while in the central and southern areas it increased from 20 to 60 m day−1 and from 41 to 70 m day−1 respectively. The final calibrated Kh of the lower confined aquifer varied from 0.2 to 70 m day−1. Also, Kv of the aquitard decreased in the northeastern area and along the MP River from 0.001 to 0.0004 m day−1 and from 0.001 to 0.0001 m day−1 respectively, finally varying spatially from 0.0001 to 0.001 m day−1. The calibration process resulted in seven zones of different specific yield (Sy) values in the upper aquifer. The dominant Sy was 0.15 and varied from 0.07 in the northern and southern parts of the modelled area to 0.20 in the western part. In the second confined aquifer, the majority of the modelled area had specific storage (Ss) of 0.00001 m−1, except the area to the SW from the dam, n e a r p i e z o m e t e r 5 / P T- 6 ( F i g . 1 ) w h i c h h a d Ss = 0.0001 m−1. The calibrated Ss of the middle aquitard was spatially uniform and equal to 0.00001 m−1. The originally assigned lakebed leakances of 0.0007 day−1 in the western part (with sapropel) and of 0.0015 day−1 in the eastern part, were maintained. Also, the streambed thicknesses were kept equal to 1.0 m, while the adjusted streambed hydraulic conductivities ranged from 0.02 to 0.1 m day−1. The final GHB conductances at the external model boundaries varied from 1.25 to 1.50 m2 day−1. The calibrated barrier hydraulic conductivity of the diaphragm under the dam was 0.05 m day−1. The initially assumed values for the evapotranspiration extinction depth, extinction water content, soil residual water content, soil-saturated water content, stream Manning coefficients, and stream widths were kept unchanged throughout the model calibration. The key issue in studying lake–groundwater interactions is the lakebed seepage. Throughout the model calibration, it was noticed that two model parameters played a major role in determination of the lakebed seepage flux, i.e. lakebed leakance (λ) and shallow-aquifer hydraulic conductivity (K), so the sensitivities of these two were tested and the results are presented in Table 2. The application of multiplication factor 10 to λ and K resulted in respectively 2 and 5.6 times increases in the average LLKnet, ~3 times increases in the maximum daily average LLKout (similar for λ and K), and 6 and 2 times increases in maximum daily shoreline LLKout, in this case larger for λ than for K. The use of the multiplication factor 0.1 for λ and K resulted in respectively 3.9 and 2.9 times decreases in the average LLKnet, 2.9 and 1.6 times decreases in the maximum daily average LLKout, and 5 and 1.1 times decreases in the daily maximum shoreline LLKout. The model indicated substantial sensitivity of the lakebed seepage to changes of both parameters, λ and K, although generally larger for λ than for K, except of the 4 years average LLKnet when the parameters where increased 10 times. As expected, based on natural lake studies (e.g. Virdi et al. 2013) , the largest sensitivity was observed in the shoreline seepage zone. The large sensitivity of lakebed seepage to λ and K, emphasizes the importance of their reliable parameterization. The importance of the λ parameter has already been Fig. 9 Water balance of the: a whole model domain following Eq. (1), referenced to the total area (AT); b lake following Eq. (4), referenced to variable lake area (ALK); c land surface and unsaturated zone following Eq. (6), referenced to variable land area (ALD); d saturated zone emphasized, for example by Hogeboom et al. (2015) , but large LLKnet sensitivity to the increase of shallow aquifer K, three times larger than to the increase of λ, is surprising and relevant. Fortunately, in this TL study, the λ and K were well defined; considering λ, the lakebed thickness was directly measured through regular underwater sampling schema (Fig. 1), while the lakebed vertical hydraulic conductivity was estimated based on core samples (Gurwin et al. 2004) ; also the K of the shallow aquifer was well defined by the large number of pumping tests and slug tests. Also sensitivity tests of the piezometric heads and lake stages to changes of various parameters were carried out. These tests were as follow: (1) change of all Kh (so also Kv = 0.1 Kh) of the upper aquifer when multiplied by 0.1 referenced to the total area (AT) following Eq. (9); each hydrologic year starts from 1 November of the previous year and ends 31 October of the year as listed in Table 3; all values are in mm year−1 and 10, resulted in the 1.75 times increase of RMSE between o b s e r ve d a n d s i m u l a t e d he a d s ( t h e s a m e f o r b ot h multiplicators) and in the increase of RMSE between observed and simulated lake stages by factors of 2.6 and 5.8 respectively; (2) change of the anisotropy factor Kh/Kv from 10 to 2, i.e. to the value suggested by Ala-aho et al. (2015) for glacial outwash deposits, did not affect the model at all; only the change to the value 50 resulted in significant increase of RMSE of piezometric heads; (3) changes of all Sy values of the upper aquifer when multiplying them by 0.5, resulted in the rise of RMSE of the piezometric heads by factor 1.95, but did not affect the lake stage; the use of multiplier >1 did not affect the model solution at all; (4) changes of all Ss values of the lower confined aquifer, when multiplying by 0.1 or 10, resulted in negligible sensitivity of both piezometric heads and lake stages. The particularly distinct sensitivity of piezometric heads and lake stages to changes of K emphasized the importance of reliable K estimates. The calibrated model was sensitive not only to changes of parameters but also to variations of the lake inflow and outflow—for example, the use of the lake inflow multiplier of only 1.1 resulted in increase of LLKnet from 0.27 to 2.19 mm day−1 and in increase of RMSE between observed and simulated lake stages and groundwater heads by factors 7.0 and 2.95 respectively; the use of multiplier 0.9 resulted in decrease of LLKnet from 0.27 to 0.17 mm day−1 and in increase of RMSE between observed and simulated lake stages and groundwater heads by factors 5.1 and 1.4 respectively. The use of inflow multiplier 0.5 caused the lake to dry out. The high sensitivity of artificial lake models to inflow and regulated outflow, points out the critical importance of reliable estimates of these data types when modelling interactions of artificial lakes with groundwater. Water balances The water balance in the form of yearly means for the 4-year transient simulation, accounting for interactions between lake and groundwater according to Eqs. (1), (4), (6) and (9), is presented in Fig. 9a–d and in the corresponding Tables S1–S4 of the electronic supplementary material (ESM). While analysing water balances, it is important to realize that the relations between the land part of the model and the lake are dynamic, because with changing lake stages, the lake area changes (see ALK in Table 3), so also the proportion between ALK and the land surface area (ALD) changes, i.e. when one area increases the other decreases because the total model domain area (AT) remains constant (74.8 km2). Such variability creates difficulty in interpretation of the lake–groundwater interactions because water flux estimates differ depending on whether they are referenced to the temporally variable ALK, temporally variable ALD or temporally invariant AT. This can be observed in Table 3, which presents yearly IN and OUT components for the whole model domain and for the lake, the latter referenced either to AT or to ALK according to Eqs. (1) and (4), respectively. For example, the 4-year mean lake-water input referenced to the lake area, i.e. IN (ALK) = 17,864.2 mm year−1, is very different than if referenced to the total modelled area IN (AT) = 3938.0 mm year−1, while both refer to the same volumetric water transfer (Table 3); therefore, water balances of artificial lakes should always clearly state to which area they are referenced. The water balance of the whole model domain that follows Eq. (1), is presented in Fig. 9a and Table S1 of the ESM. The three water inputs as per Eq. (1) include: P (15.9% of IN), Qin (82.2% of IN), and GHBin (1.8% of IN), whereas the three water outputs include: ET (16.4% of OUT), Qout (83.3% of OUT), and GHBout (0.4% of OUT). The ΔS = IN – OUT was highly variable, ranging from 8.0 mm year−1 in the second year to 145.1 mm year−1 in the fourth year, but positive, which is attributed to the storage function of the dam reservoir and to the controlled lake outflow. The Qin and Qout were the major components of the water balance, which explains the large model sensitivity to their changes, while P and ET were several times lower. The lateral GHB-groundwater inflows/outflows were the lowest, although important components of the water balance. It is remarkable that the GHBin was 5–6 times larger than GHBout, mainly because of the sealing role of the hydraulic barrier of the dam, equipped with internal diaphragm, simulated by the HFB Package of MODFLOW. Such water balance as in this study, where lake inflow and regulated outflow contributions are much larger than rainfall or evapotranspiration, seems to be characteristic for artificial lakes constructed on large rivers. The water balance of the lake itself, accounting for temporally (e.g. yearly) variable lake area extent, is presented in Fig. 9b and Table S2 of the ESM. It follows Eq. (4), including four input fluxes: P (3.3% of IN), QLKin (95.3% of IN), LLKin (0.2% of IN), RO (1.2% of IN) and three output fluxes: ELK (4.9% of OUT), QLKout (94.3% of OUT), LLKout (0.8% of OUT)—all referenced to the lake area (ALK). The presented water balance emphasizes the dominant role of the QLKin and regulated QLKout. The losses of ELK are larger than P input to the lake. Also, the losses of LLKout are larger than LLKin, but those losses are dependent on the lake regulation—for example, forcing larger lake outflows (QLKout) than natural water inflows (QLKin) results in the lowering of lake stage, thus also in the reduction of LLKout. Regarding ΔSLK, among the 4 years analysed, only in 2006 was the ΔSLK negative, which means that only in that year did the lake output (OUT) exceed the natural lake input (IN), while in all other years the storage effect was the opposite, which was mainly due to the exceptionally large QLKout in 2006, nearly the same as QLKin which typically is much larger than QLKout. The water balance of the land surface and unsaturated zone, accounting for variable land surface area, handled by UZF1 Package as per Eq. (6), is presented in Fig. 9c and Table S3 of the ESM, and includes two input fluxes, i.e. P applied to the land part of the model domain as P – I (73% of P), which infiltrates through the soil according to the calibrated vertical hydraulic conductivity (Kv) and Exfgw (6.4% of P). On the output side, there are four fluxes: Rg (14.4% of P), ETuz (52.6% of P), I (27.3% of P), and RO (8.3% of P); note that the water balance shows a relatively small contribution of Rg which is mainly because of the large ETuz. The water balance of the saturated zone (two aquifers and aquitard) of the whole model domain as per Eq. (9), is presented in Fig. 9d and Table S4 of the ESM. The 4-year groundwater balance includes four input fluxes: Rg (11.8% P), GHBin (11.5% P), QGWin (1.6% of P), and LLKout (4.2% P); and five output Fig. 10 Monthly averages of groundwater heads (GH), groundwater evapotranspiration (ETg), groundwater exfiltration (Exfgw) and net lake seepage (LLKnet), all versus monthly lake stages within the 4-year fluxes: QGWout (10.3% P), ETg (8.9% of P), Exfgw (5.2% P), GHBout (2.2% of P), and LLKin (1.1% P). The main groundwater inputs to the study area are Rg and GHBin, while the main groundwater outputs are QGWout, ETg and Exfgw. The groundwater balance indicates that throughout the four analysed years, the TL groundwater system gained groundwater storage (ΔS = 9.0 mm year−1 = 1.5% of P), likely because of the presence of the reservoir. The only study that provides similarly detailed lake-waterbalance analysis, although over a natural (not artificial) and much smaller seepage lake (Lake Starr, only ~0.8 km2), i.e. a lake without surface inflow/outflow, was carried out by Virdi et al. (2013) . In that study, the 10-year means of rainfall and lake evaporation were 1,270 and 1,450 mm year−1, respectively, each of them being more than one order of magnitude larger than any other component of their water balance. As such, their lake dynamics were driven primarily by these two climatic driving forces, i.e. balance between rainfall and evaporation, in contrast simulation period and all referenced to the land area; a and r are the regression and correlation coefficients respectively to the artificial lake dynamics reported here, which was driven primarily by the balance between QLKin and man-induced regulation of the QLKout. It is remarkable that in artificial TL: (1) not only the lake itself but also the whole catchment domain dynamics were driven by QLKin and regulated QLKout, whereas natural driving forces such as rainfall or evapotranspiration, that typically constrain dynamics of natural lakes, in artificial TL had relatively minor importance; (2) lakebed seepages were generally lower than one would expect in a lake characterized by large stage fluctuation amplitude and extensive periods with high lake stages, implying LLKout > LLKin, which was because of solid, low-permeability and thick lakebed isolation, reflected by low lakebed leakance; (3) the contributions of ETg and particularly of Exfgw to the overall water balance were much larger relative to Rg than in natural lakes, which can be explained by reservoir influence upon the shallow water table, particularly distinct in areas close to the lake where the water table is frequently pulled up by sudden rises Fig. 11 Daily variability of land water fluxes for the 4-year model simulation period, from 1 Nov. 2004 to 31 Oct. 2008 for the following water balance components: a precipitation (P), actual infiltration rates (Pa), potential evapotranspiration (PET) and measured lake stages; b gross recharge (Rg), net recharge (Rn), groundwater exfiltration (Exfgw), groundwater evapotranspiration (ETg) and unsaturated zone evapotranspiration (ETuz); note that the water fluxes are referenced to the land area (ALD) of lake stages when outlet weir restricts QLKout; such behaviour seems to be specific for artificial lakes, creating different system dynamics than in natural lakes; (4) QGWout > QGWin due to the generally high water table, enhanced by frequent high reservoir stages, resulting in hydraulic gradient stimulating groundwater drainage by rivers and streams; (5) GHBout < GHBin because of the sealing function of the dam, i.e. mainly its subsurface diaphragm section. The preceding five groundwater balance observations seem to all be characteristic for artificial lake systems. Reservoir–groundwater interactions Water table and subsurface water fluxes versus lake stages The variability of lake stages and of the water table is presented in Fig. 8. Even a quick visual inspection allows one to notice that the water table is well correlated with the lake stages. The statistical assessment of the correlation between daily lake stages and groundwater levels in the piezometers presented in Fig. 8 indicated the best correlation, r = 0.9, for the closest to the lake, shallow aquifer piezometer 12/PT-2, while for the rest of the piezometers, r varied from 0.68 at the piezometer PT-34 to 0.49 at the piezometer PT-116, dep e n d i n g o n t h e d i s t a n c e f r o m t h e l a k e a n d t h e hydrogeological conditions of the subsurface. To assess general water-table dependence upon the lake stages around the perimeter of the lake, 14 fictitious piezometers were assigned as presented in Fig. 6 and their monthly estimates over 4-year model simulation, plotted versus lake stages as depicted in Fig. 10, provided pretty good correlation (r = 0.66). It is remarkable that despite ~5 m amplitude of the lake stage fluctuation, the observed monthly average water-table variability was quite low, on the order of 1 m (Fig. 10) and daily <1.5 m (Fig. 8), even in locations very close to the lake. Such amplitude disproportion between lake stages and groundwater levels is primarily attributed to the large hydraulic resistance (low leakance) of the lakebed. Considering subsurface water fluxes dependence on lake stages, it is not surprising that the largest correlation, r = 0.95, was observed for LLKnet because of its direct dependence on the hydraulic gradient at the interface between the lake and groundwater. The Exfgw is enhanced by sudden water-table rises; thus, the rapid changes of lake stages, typical for artificial lakes, explain rather high correlation (r = 0.87) of Exfgw with lake stages. The ETg is directly dependent on water table depth (WTD) influenced by lake stages, although within the limit defined by extinction depth, which is the reason that ETg correlation with the lake stages (r = 0.79) was lower than the other two. The Rg is the only groundwater flux not directly dependent on WTD, although the WTD determines the storage capacity of unsaturated zone, which if large, enhances ETuz reducing Rg; therefore, Rg correlation with lake stages was the lowest (r = 0.06; not displayed in Fig. 10). Temporal variability of subsurface fluxes Figure 11a shows the daily variability of Pa, P, PET and lake stages; the Pa calculated by UZF1 Package, depends on P, interception loss (I), the vertical hydraulic conductivity (Kv) of the soil and the saturation degree of the unsaturated zone. The peaks of Pa, as the peaks of P, occur in the summer months (June, July, Aug and sometimes in Sept) and coincide with the peaks of PET that enhances actual evapotranspiration, but not with the peaks of lake stages. Figure 11b illustrates the daily variability of the main subsurface water fluxes, unsaturated zone evapotranspiration (ETuz), gross recharge (Rg), grou ndwater exfiltration (Exfgw), groundwater evapotranspiration (ETg) and net recharge (Rn). The temperate climate of the study area explains that ETuz is ~5–10 times larger than ETg and that both follow the same pattern as the PET. The Rn (Eq. 10) follows the Rg in the months November–February when ETg and Exfgw are negligible, and diverges in April–August when it follows ETg, occasionally being interrupted by peaks of Rg. Exfgw represents also an interesting pattern, which in contrast to other studies such as for example Hassan et al. (2014) , has the maxima coinciding not with the largest P and Rg but with the highest lake stages, dependent mainly not on P but on the balance between MP river inflow (QLKin) and regulated reservoir outflow (QLKout), forcing the water table either to rise or decline. Such behaviour seems to be characteristic for groundwater systems interacting with artificial lakes (such as the TL), although more studies are needed to confirm that statement. The aquitard leakage, constrained by head difference between the overlying and underlying aquifers, is dependent on the lake stages influencing the aquifer heads, particularly those of the shallow aquifer. It is remarkable that throughout the 4 years of model simulation, the aquitard leakage had always an upward direction within the entire study area, being enhanced by the spatio-temporally variable head difference across the aquitard, locally reaching even up to 8 m. The largest aquitard leakage occurred during low lake stages along the earthen dam (up to 0.8 mm day−1), along the MP River in the upstream section (0.7 mm day−1) and below the lake (~0.5 mm day−1). In the remaining area the upward leakage was moderate, i.e. less than 0.3 mm day−1. During high lake stages, the upward leakage ranged from 0.0 to 0.3 mm day−1. The upward direction of the aquitard leakage seems to be specific for the TL hydrogeological system characterized by natural, upward-directed groundwater flow at the MP River Valley; thus, it must not be attributed to artificial lakes in general. Fig. 12 Spatial and temporal impact of the lake on groundwater heads at the northern transect of a series of fictitious piezometers (see Fig. 6) Spatio-temporal extent of reservoir impact One of the most important problems to solve in the groundwater systems adjacent to artificial lakes is the assessment of spatio-temporal extent of a lake impact on groundwater. Figures 12 and 13 show the observed and simulated daily lake stages and the simulated groundwater heads in the two transects of fictitious piezometers presented in Fig. 6. Both figures show that the impact of the lake on groundwater heads declines with distance from the lake; however, there are differences between the two transects attributed to different hydraulic and hydrogeological conditions. The northern transect presented in Fig. 12, is affected by lake-stage fluctuation in all fictitious piezometers tested. The strength of the lake impact on groundwater response, only slightly declined towards the watershed no-flow boundary and was characterized by a few days of lag between the lake-stage peaks or dips and the corresponding groundwater head responses, raising an issue of appropriateness of that noflow boundary. The trial of 800 m shift of that boundary to the north and its replacement with a head dependent boundary (GHB) allowed for testing and controlling eventual water transfers across the former no-flow boundary. That trial resulted in negligible flux transfer across the former no-flow boundary line, which allowed to maintain the groundwater divide along that line. Moreover, the heads outward from the groundwater divide represented different dynamics as compared to the heads of the study area, indicating that they were of different groundwater systems. That model boundary test allowed for keeping the original northern no-flow boundary as valid and also proposes a way of testing the assigned Fig. 13 Spatial and temporal impact of the lake on groundwater heads at the southern transect of a series of fictitious piezometers (see Fig. 6) boundary, in case there are doubts about its appropriate performance. The intriguing, distinct groundwater level rise in the 900m-from-the-shore curve (Fig. 12) for the cell located near the northern no-flow boundary, represents a typical MODFLOW-NWT treatment of dry cells, which may compute a head for a dry cell that is higher than that of the surrounding cells. To avoid the convergence failure due to flow discontinuity that may be caused by flow out of a dry cell, MODFLOW-NWT sets the conductance between the dry cell and the surrounding cells equal to zero. By doing so, the dry cells can be kept active without allowing water to flow out (Niswonger et al. 2011) . Figure 13 shows standard groundwater system response to the lake fluctuations at the southern transect (Fig. 6), i.e. head fluctuations declining with distance from the lake towards the GHB boundary where fluctuations are not present. Such a fluctuation pattern was expected, because the groundwater heads at the southern side are generally higher than the lake stages and the hydraulic conductivity is low. The investigation of the lake impact on groundwater at other locations of this study area, carried out by assigning more spatially distributed fictitious piezometers, indicated that the lake had almost no effect on groundwater levels to the west of the dam. At the southern side, the water table was affected by lake fluctuations within a distance of ~200–1,100 m from the lake shoreline. Moving towards the east of the lake, the lake-stage fluctuations only slightly affected the amplitude of groundwater level fluctuation; however, that effect extended over a large distance of ~1,400–1,500 m. At the north, the lake affected the water table over the entire area, i.e. from the Fig. 14 Daily variability of groundwater seepage into the lake (LLKin), seepage from the lake into groundwater (LLKout) and net lake seepage (LLKnet) during the 4-year model simulation period, from 1 Nov. 2004 shoreline to the assigned boundary along the watershed divide. Lakebed seepage The analysis of lakebed seepage is a primary characteristic of lake–groundwater interactions. Such analysis can be carried out on a temporal and/or spatio-temporal basis. Figure 14 presents temporal analysis of daily seepage flux across the entire TL bed bottom, taking into account variable lake area extent. It involves three seepage components, from aquifer to the lake (LLKin), from lake to the aquifer (LLKout) and the net seepage (LLKnet = LLKin – LLKout). The LLKout, ranging up to ~ −1.0 mm day−1, was dominant during stages higher than ~172.5 m a.s.l., occurring typically from February to September, whereas the LLKin, ranging up to ~0.6 mm day−1, was dominant during stages <172.5 m a.s.l., occurring typically from October to January. The mean LLKnet within the 4-year simulation period was approximately −0.3 mm day−1, implying that LLKout > LLKin; thus, the lake was generally loosing water towards the underlying aquifer, which is understandable because artificial lake stages are generally higher than original groundwater heads, creating additional hydraulic gradient. Applying similar methodology of seepage analysis as in this study, Virdi et al. (2013) carried out daily model simulation over a 10-year period on the previously mentioned natural Lake Starr system. However, in contrast to this study, their lake-stages exhibited multi-year trend changes, while their yearly (seasonal) lake amplitudes were substantially lower (<1 m) and the water-table amplitudes were pretty large, larger than in this study, and as an exception, even >2 m. Besides, it is remarkable that despite relatively slower system changes and smaller hydraulic gradients between lake and groundwater, their lakebed seepage peaks were still larger than in TL; this was most likely because of their much larger leakance (0.059–0.077 day−1) as compared to the TL leakance (from 0.0007 to 0.0015 day−1) and much smaller lake area (0.8 vs ~13.5 km2), thus a larger contribution of the active peripheral seepage zone (the shoreline zone between min and max lake extent) in the total lake area. The importance of the peripheral seepage zone in the lake– groundwater interaction was emphasized, MODFLOWto 31 Oct. 2008; note that the fluxes are referenced to the temporally variable lake area (ALK) modelled and field-validated by Kidmose et al. (2011) at the natural Lake Hampen (0.760 km2 surface area, 3.2 Mm3 volume and rainfall ~900 mm year−1) in western Denmark. In their steady-state model, shoreline seepage was highly spatially variable, ranging from LLKout = ~ −70 mm day−1 to LLKin = ~60 mm day−1 and was in agreement with seepage meter measurements. Such large seepage, even as for the peripheral zone seepage, can be explained by their large lakebed leakance of 0.26 day−1 and the large hydraulic conductivities of the upper aquifer ranging from 15 to 372 m day−1, both much larger than in the TL study. The importance of peripheral zone seepage at a number of natural, boreal lakes in Finland was also studied and confirmed by Ala-aho et al. (2015) , who used field measurements, airborne thermal imaging, stable isotopes and the integrated HydroGeoSphere model. Their simulated peripheral zone seepage was up to 3.3 × 10−7 m s−1 (28.5 mm day−1), while that estimated by stable isotopes was up to 5.7 × 10−7 m s−1 (49.2 mm day−1). They concluded that the lake peripheral zones represent enhanced groundwater/surface-water exchange controlled by subsurface features such as variability in lakebed thickness or hydraulic conductivity, which is to say by leakance. The role and importance of active peripheral seepage zones in natural lakes raises the question as to what the spatiotemporal seepage patterns look like in artificial lakes. For the Turawa artificial lake, such patterns together with head distributions are presented in Fig. 15. The three characteristic patterns correspond to the three different lake stages—lowest, average and highest. At the lowest lake stage (Fig. 15a), the lake extent was the smallest and the zero seepage line was the closest to the dam, with a small LLKout area (14% of the maximum lake area) at the west and a large LLKin area (61% of the maximum lake area) at the east. Outwards from the zero seepage line, the LLKout increased in the western direction towards the dam up to 1.6 mm day−1, while in the eastern direction, LLKin first increased and then decreased, while the largest LLKin was in the northern and southern peripheral zones. Similar seepage patterns can be observed for the average lake stage (Fig. 15b), with the following differences: (1) the lake area increased; (2) the zero seepage line moved further to the east and the new zero seepage line appeared in the east, maintaining LLKin in between; (3) LLKout at the western side 2007 (positive LLKin represents upward lake seepage-gain and negative LLKout represents downward lake seepage-loss) increased, while LLKin decreased as a result of the lake stage rise, reducing hydraulic gradient between the lake stage and the groundwater. A totally different seepage pattern was observed for the maximum lake stage (Fig. 15c)—the zero seepage line moved away so that the entire reservoir area was affected by LLKout with negligible value in the central part of the lake and the largest in the peripheral zone, with maximum LLKout in order of 4.6 mm day−1 at the dam side and along the northwest lake shoreline. To the authors’ knowledge, this is the first spatio-temporal seepage analysis of an artificial lake interacting with groundwater. As such there is no other artificial lake study to compare with lakebed seepages of this study. Moreover, it is also difficult to compare spatio-temporal lakebed seepage patterns of this study with seepage patterns of natural lakes, not only because these are totally different lake systems with different dynamics but also because there are very few published natural lake studies providing spatio-temporal seepage analysis that requires complex modelling methodologies of lake– groundwater interaction such as the fully coupled transient model calibration, accounting for variably-saturated flow. Probably Smerdon et al. (2007) were the first who provided such analysis applying integrated InHM model to a small boreal lake (39 ha) in Canada. The lake seepage pattern in that study: (1) varied gradually from the largest LLKin at the upstream recharge side and throughout the zero seepage line, to the largest LLKout at the downstream discharge side of the lake outlet; (2) had dominant LLKin occupying 75–83% of the lake area, thus similar to the TL seepage pattern at the lowest lake stage (Fig. 15a); but (3) was changing only by 8% within ~1.5 year of the simulation period, while in the TL study, changes were fast and dramatic—for example, the area of LLKout could change from 14% at the lowest stage to 100% at the highest stage (Fig. 15c) within less than half a year. Also, in contrast to the TL study, they attributed their spatiotemporal seepage variability mainly to summer lake evaporation, while in the TL case, the seepage is driven mainly by the balance between lake inflow and outflow, while the lake evaporation has secondary impact on the dynamics of the lake– groundwater interactions (Fig. 9b and Table S2 of the ESM). Also, for the previously mentioned Lake Starr, Virdi et al. (2013) presented daily seepage patterns for low and high lake stages acquired with similar modelling methods as presented in this study. That analysis resulted in concentric spatiotemporal lake seepage patterns, i.e. with LLKin in the peripheral parts and LLKout in the deepest central part of the lake. Their LLKin area changed from 6 to 59% and LLKout from 41 to 94%, thus more than in Smerdon et al. (2007) but still less than in the TL study. However, their seepage changes were even slower than in Smerdon et al. (2007) and none of these natural lake studies had a seepage pattern like the highest lake stage in the TL study (Fig. 15c), with the entire lake area covered by LLKout due to the high artificially forced stage in the reservoir. It is remarkable that in this artificial lake study, the spatiotemporal seepage patterns at various lake stages were different from each other and changed so quickly. These fast changing patterns and their spatio-temporal variability are attributed to large and fast lake stage variations, heterogeneity of the lakebed leakance and configuration of the lakebed bottom. Particularly the two seepage patterns of the low and average lake stages (presented in Fig. 15a,b) are very different as compared to the third pattern of the maximum lake stage presented in Fig. 15c. The fast changing and highly spatio-temporally-variableseepage patterns seem to be characteristic of artificial lakes. The direction and magnitude of seepage between the lake and underlying aquifer system depends on the difference between lake stage and hydraulic head of an aquifer underlying the lakebed and on the lakebed leakance. Natural lakes in contrast to artificial lakes have typically small (<0.5 m) amplitude of seasonal lake stage fluctuations (Ala-aho et al. 2015; Hunt et al. 2013; Smerdon et al. 2007; Taviani and Henriksen 2015) , while multi-year climate-driven stage changes (see Virdi et al. 2013) may exceptionally reach up to ~1 m year−1. For comparison, in the artificial TL, seasonal amplitude was ~5 m, but a stage change of ~2.5 m could take place within only 22 days. Such large and fast stage changes in artificial lakes, have primary impact on the lakebed-seepage direction and magnitude. Despite the large hydraulic gradient between artificial lake stages and groundwater levels, rather low seepages are expected as artificial lakes are usually well isolated at the bottom, and thus have low lakebed leakance. In the case of the TL, after the dam and reservoir construction, the lakebed leakance was additionally reduced by industrial sediment, which was transported from outside the reservoir and deposited at its bottom, increasing lakebed thickness and colmatating (chemical and mechanical sediment clogging) it, as confirmed by the regular grid of lakebed drilling and sediment sampling (Fig. 1). The Turawa lakebed leakance was indeed low (0.0007–0.0015 day−1) as compared to leakances of natural lake studies, which were one to two (0.059–0.077 day−1 in Virdi et al. 2013) or even two to three (0.26 day−1 in Kidmose et al. 2011) orders of magnitude larger. The multitude of interplaying factors affecting interactions of artificial lakes with groundwater, including maninduced lake outflow regulation and large spatio-temporal variability of interaction processes, implies that assessment of such complex hydrological systems requires complex and sophisticated modelling methods. An example of such method involving a distributed, fully coupled, transient model calibration with volumetric representation of water exchanges between artificial lake and groundwater, accounting for temporally variable lake extent and volume, and applying kinematic wave approximation to variably-saturated flow, was tested in this study-case of the artificial TL. The method proposed turned out to be highly suitable for complex artificial lake systems and is highly recommended for their management. Conclusions The dynamics of artificial lake interactions with groundwater is complex due to large and fast, man-induced changes of lake stages, which in contrast to natural lakes, are mainly dependent not on natural driving forces (e.g. P, ET) but on temporally variable balance between river inflow and regulated outflow. The complex dynamics of artificial lake systems requires realistic and sophisticated modelling solutions, i.e. (1) distributed, fully coupled model with variably-saturated f l o w a c c o u n t i n g f o r s p a t i o - t e m p o r a l l y v a r i a b l e , volumetric representation of water exchanges between the artificial lake and the groundwater, accounting for temporally variable lake extent; (2) transient model calibration at least on a daily basis, involving two or more different calibration state variables, and controlling the volumetric water budget. In this TL study, MODFLOW-NWT model with LAK7, SFR7 and UZF1 packages was used and daily transient model calibration was based on 5-year daily data, 1-year for warm-up and 4years for calibration, with lake stages and piezometric heads as calibrated state variables while controlling volumetric water budget. The complex dynamics of artificial lake interactions with groundwater requires appropriate and reliable data sets to carry reliable simulations. It is important to have at least daily time series of climatic driving forces (P, PET), lake inflows, lake outflows and state variables, including at least lake stages and some piezometric heads well-spatially distributed. Regarding system parameterization, the most important according to the carried out sensitivity analysis are reliable estimates of lakebed leakance and hydraulic conductivity of the underlying aquifer(s). Artificial lakes, in contrast to natural lakes, are net exporters of lake water into groundwater through lakebed seepage; that seepage and its spatio-temporal pattern depend not only on the lake–groundwater gradient and lakebed leakance but also on dam construction and lake bathymetry; in TL, the lakebed seepage was lower than in natural lakes despite larger hydraulic gradients, because of the much lower lakebed leakance; also the spatio-temporal seepage pattern was different than in natural lakes. The impact of artificial lakes on the water table in the adjacent to the lake area declined with distance from the lake and was dependent on the lake stage, lakebed leakance, groundwater head distribution and hydrogeological system parameterisation. This artificial lake study discovered also that groundwater evapotranspiration and groundwater e x f i l tr a t i o n , w h i le i n f l u e n c e d b y l a rg e l a k e s t a g e fluctuations, in contrast to natural lakes, had maxima coinciding not with the rainfall peaks but with the highest lake stages dependent mainly on the balance between lake inflow and regulated outflow. For operational artificial lake management models, a finer grid than used in this study (200 × 200 m) is recommended, although such grid would require powerful computational hardware. The modelling method proposed in this study can be adapted to any artificial lake system (but also for any natural lake system) to enhance its operation, planning and management, as it can provide accurate and adequate assessment of eventual consequences of regulation of lake stages upon the adjacent groundwater system. Acknowledgements We acknowledge Wroclaw University and the Institute of Meteorology and Water Management - National Research Institute (IMGW-PIB) in Poland for allowing us to use the data for this paper, Government of the Netherlands for the scholarship offered to the first author, support from Dr. Richard Niswonger and Mr. Richard Winston (the developers of MODFLOW-NWT, UZF1 Package and ModelMuse) for technical help in model development, and Dr. Zoltan Vekerdy for very helpful comments at an early stage of this paper writing. Open Access This article is distributed under the terms of the Creative C o m m o n s A t t r i b u t i o n 4 . 0 I n t e r n a t i o n a l L i c e n s e ( h t t p : / / creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Ala-aho P , Rossi PM , Isokangas E , Kløve B ( 2015 ) Fully integrated surface-subsurface flow modelling of groundwater-lake interaction in an esker aquifer: model verification with stable isotopes and airborne thermal imaging . J Hydrol 522 : 391 - 406 . doi: 10 .1016/j. jhydrol. 2014 . 12 .054 Allen RG , Pereira LS , Raes D , Smith M ( 1998 ) Crop evapotranspiration: guidelines for computing crop water requirements . FAO Irrigation and drainage paper 56 , FAO, Rome Anibas C , Fleckenstein JH , Volze N , Buis K , Verhoeven R , Meire P , Batelaan O ( 2009 ) Transient or steady-state? Using vertical temperature profiles to quantify groundwater-surface water exchange . Hydrol Process 23 : 2165 - 2177 . doi: 10 .1002/hyp.7289 Brindha K , Neena Vaman KV , Srinivasan K , Sathis Babu M , Elango L ( 2014 ) Identification of surface water-groundwater interaction by hydrogeochemical indicators and assessing its suitability for drinking and irrigational purposes in Chennai, southern India . Appl Water Sci 4 : 159 - 174 . doi: 10 .1007/s13201-013-0138-6 Brooks RH , Corey AT ( 1966 ) Properties of porous media affecting fluid flow . J Irrig Drain Eng 101 : 85 - 92 Chapman D ( 1996 ) Water quality assessments: a guide to the use of biota, sediments and water in environmental monitoring, 2nd edn . Spon , London Cheng X , Anderson MP ( 1993 ) Numerical simulation of ground-water interaction with lakes allowing for fluctuating lake levels . Ground Water 31 : 929 - 933 . doi: 10 .1111/j.1745- 6584 . 1993 .tb00866.x Chhuon K , Herrera E , Nadaoka K ( 2016 ) Application of integrated hydrologic and river basin management modeling for the optimal development of a multi-purpose reservoir project . Water Resour Manag 30 : 3143 - 3157 . doi: 10 .1007/s11269-016-1336-4 Doherty J , Hunt R ( 2010 ) Approaches to highly parameterized inversion: a guide to using PEST for groundwater-model calibration . US Geol Surv Sci Invest Rep 2010 - 5169 , 59 pp Fenske JP , Leake SA , Prudic DE ( 1996 ) Documentation of a computer program (RES1) to simulate leakage from reservoirs using the modular finite-difference ground-water flow model (MODFLOW) . US Geol Surv Open-File Rep 96-364 , 51 pp Fowe T , Nouiri I , Ibrahim B , Karambiri H , Paturel JE ( 2015 ) OPTIWAM: an intelligent tool for optimizing irrigation water management in coupled reservoir-groundwater systems . Water Resour Manag 29 : 3841 - 3861 . doi: 10 .1007/s11269-015-1032-9 Ghosh N et al ( 2015 ) Semi-analytical model for estimation of unsteady seepage from a large water body influenced by variable flows . Water Resour Manag 29 : 3111 - 3129 . doi: 10 .1007/s11269-015-0985-z Guimerà J , Vives L , Carrera J ( 1995 ) A discussion of scale effects on hydraulic conductivity at a granitic site (El Berrocal, Spain) . Geophys Res Lett 22 : 1449 - 1452 . doi: 10 .1029/ 95GL01493 Gurwin J ( 2010 ) Zagrożenie wód podziemnych w ocenie oddziaływania na środowisko planowanej renaturalizacji zbiornika retencyjnego Turawa [Groundwater hazard regarding environmental impact assessment of renaturalisation of the Turawa reservoir] . Biul PIG Hydrogeol 440 : 65 - 76 Gurwin J , Kryza J , Skowronek A (eds) ( 2004 ) Ocena stanu ekologicznego Jeziora Turawskiego w celu opracowania działań na rzecz jego poprawy [Environmental assessment of the Turawa Lake for its ecological improvement] . Geological Institute of the Wroclaw University, Wrocław, Poland Gurwin J , Kryza H , Kryza J , Poprawski L ( 2005 ) Rozpoznanie wód podziemnych w rejonie Jeziora Turawskiego dla potrzeb oceny stanu ekologicznego [Groundwater recognizing in the region of Turawskie Lake for needs of ecological state assessment] . Współczesne Problemy Hydrogeol XII: 241 - 253 Harbaugh AW ( 1990 ) A computer program for calculating subregional water budgets using results from the US Geological Survey Modular three-dimensional finite-difference ground-water flow model . US Geol Surv Open-File Rep 90-392 , 46 pp Hassan SMT , Lubczynski MW , Niswonger RG , Su ZB ( 2014 ) Surfacegroundwater interactions in hard rocks in Sardon catchment of western Spain: an integrated modeling approach . J Hydrol 517 : 390 - 410 . doi: 10 .1016/j.jhydrol. 2014 . 05 .026 Hill M , Tiedeman C ( 2007 ) Effective groundwater model calibration: with analysis of data, sensitivities, predictions, and uncertainty . Wiley, Hoboken, NJ Hogeboom RHJ , van Oel PR , Krol MS , Booij MJ ( 2015 ) Modelling the influence of groundwater abstractions on the water level of Lake Naivasha, Kenya under data-scarce conditions . Water Resour Manag 29 : 4447 - 4463 . doi: 10 .1007/s11269-015-1069-9 Hsieh PA , Freckleton JR ( 1993 ) Documentation of a computer program to simulate horizontal-flow barriers using the US Geological Survey's modular three-dimensional finite-difference ground-water flow model . US Geol Surv Open-File Rep 92-477 , 32 pp Hunt RJ ( 2003 ) Ground water-lake interaction modeling using the LAK3 Package for MODFLOW 2000 . Ground Water 41 : 114 - 118 . doi: 10 . 1111/j.1745- 6584 . 2003 .tb02575.x Hunt RJ , Prudic DE , Walker JF , Anderson MP ( 2008 ) Importance of unsaturated zone flow for simulating recharge in a humid climate . Ground Water 46 : 551 - 560 . doi: 10 .1111/j.1745- 6584 . 2007 . 00427 .x Hunt RJ , Walker JF , Selbig WR , Westenbroek SM , Regan RS ( 2013 ) Simulation of climate-change effects on streamflow, lake water budgets, and stream temperature using GSFLOW and SNTEMP, Trout Lake Watershed, Wisconsin . US Geol Surv Sci Invest Rep 2013 - 5159 , 118 pp. http://pubs.usgs.gov/sir/2013/5159/. Accessed July 2017 Kanduč T , Grassa F , McIntosh J , Stibilj V , Ulrich-Supovec M , Supovec I , Jamnikar S ( 2014 ) A geochemical and stable isotope investigation of groundwater/surface-water interactions in the Velenje Basin, Slovenia . Hydrogeol J 22 : 971 - 984 . doi: 10 .1007/s10040-014-1103-7 Kidmose J , Engesgaard P , Nilsson B , Laier T , Looms M ( 2011 ) Spatial distribution of seepage at a flow-through Lake: Lake Hampen, western Denmark . Vadose Zone J 10 : 110 - 124 . doi: 10 .2136/vzj2010.0017 Krabbenhoft DP , Bowser CJ , Anderson MP , Valley JW ( 1990 ) Estimating groundwater exchange with lakes: 1. the stable isotope mass balance method . Water Resour Res 26 : 2445 - 2453 . doi: 10 .1029/ WR026i010p02445 Leblanc M , Favreau G , Tweed S , Leduc C , Razack M , Mofor L ( 2007 ) Remote sensing for groundwater modelling in large semiarid areas: Lake Chad Basin, Africa . Hydrogeol J 15 : 97 - 100 . doi: 10 .1007/ s10040-006-0126-0 Lee DR ( 1977 ) A device for measuring seepage flux in lakes and estuaries . Limnol Oceanogr 22 : 140 - 147 . doi: 10 .4319/lo. 1977 . 22 .1. 0140 Lee TM ( 1996 ) Hydrogeologic controls on the groundwater interactions with an acidic lake in karst terrain , Lake Barco, Florida. Water Resour Res 32 : 831 - 844 Liuzzo L , Noto LV , Arnone E , Caracciolo D , La Loggia G ( 2015 ) Modifications in water resources availability under climate changes: a case study in a Sicilian Basin . Water Resour Manag 29 : 1117 - 1135 . doi: 10 .1007/s11269-014-0864-z Lubczynski MW , Gurwin J ( 2005 ) Integration of various data sources for transient groundwater modeling with spatio-temporally variable fluxes: Sardon study case , Spain. J Hydrol 306 : 71 - 96 . doi: 10 . 1016/j.jhydrol. 2004 . 08 .038 McDonald MG , Harbaugh AW ( 1988 ) A modular three-dimensional finite-difference ground-water flow todel . Techniques of WaterResources Investigations, Book 6 , Chapter A1 , US Geological Survey, Reston, VA, 586 pp McMahon TA , Peel MC , Lowe L , Srikanthan R , McVicar TR ( 2013 ) Estimating actual, potential, reference crop and pan evaporation using standard meteorological data: a pragmatic synthesis . Hydrol Earth Syst Sci 17 : 1331 - 1363 . doi: 10 .5194/hess-17- 1331 -2013 Merritt ML , Konikow LF ( 2000 ) Documentation of a computer program to simulate lake-aquifer interaction using the MODFLOW ground water flow model and the MOC3D solute-transport model . US Geol Surv Water Resour Invest Rep 00-4167 , 146 pp Mylopoulos N , Mylopoulos Y , Tolikas D , Veranis N ( 2007 ) Groundwater modeling and management in a complex lake-aquifer system . Water Resour Manag 21 : 469 - 494 . doi: 10 .1007/s11269-006-9025-3 Niswonger RG , Prudic DE ( 2005 ) Documentation of the StreamflowRouting (SFR2) Package to include unsaturated flow beneath streams: a modification to SFR1. US Geological Survey Techniques and Methods 6-A13 , USGS, Reston, VA, 50 pp Niswonger RG , Prudic DE , Regan RS ( 2006 ) Documentation of the Unsaturated-Zone Flow (UZF1) Package for modeling unsaturated flow between the land surface and the water table with MODFLOW2005 . US Geological Survey Techniques and Methods 6-A19, Book 6 , Chapter A19 , USGS, Reston, VA, 62 pp Niswonger RG , Panday S , Ibaraki M ( 2011 ) MODFLOW-NWT, a Newton formulation for MODFLOW-2005 . US Geological Survey Techniques and Methods 6-A37 , USGS, Reston, VA, 44 pp NOAA ( 2004 ) The Climate Visualization System . National Oceanic and Atmospheric Administration , US Dept. of Commerce. https://www7. ncdc.noaa.gov/CDO/cdo#TOP. Accessed 27 September 2013 Ong JB , Zlotnik VA ( 2011 ) Assessing lakebed hydraulic conductivity and seepage flux by potentiomanometer . Ground Water 49 : 270 - 274 . doi: 10 .1111/j.1745- 6584 . 2010 . 00717 .x Otz MH , Otz HK , Otz I , Siegel DI ( 2003 ) Surface water/groundwater interaction in the Piora aquifer, Switzerland: evidence from dye tracing tests . Hydrogeol J 11 : 228 - 239 . doi: 10 .1007/s10040-002-0237-1 Owor M , Taylor R , Mukwaya C , Tindimugaya C ( 2011 ) Groundwater/ surface-water interactions on deeply weathered surfaces of low relief: evidence from lakes Victoria and Kyoga, Uganda . Hydrogeol J 19 : 1403 - 1420 . doi: 10 .1007/s10040-011-0779-1 Penman HL ( 1948 ) Natural evaporation from open water, bare soil and grass . Proc Roy Soc London A193 : 120 - 146 Prudic DE ( 1989 ) Documentation of a computer program to simulate stream-aquifer relations using a modular, finite-difference, groundwater flow model . US Geol Surv Open-File Rep 88-729 Rudnick S , Lewandowski J , Nützmann G ( 2015 ) Investigating groundwater-lake interactions by hydraulic heads and a water balance . Groundwater 53 : 227 - 237 . doi: 10 .1111/gwat.12208 Sacks LA , Swancar A , Lee TM ( 1998 ) Estimating ground-water exchange with lakes using water-budget and chemical mass-balance approaches for ten lakes in ridge areas of Polk and Highlands counties , Florida. US Geol Surv Water Resour Invest Rep 98-4133 Sacks LA , Lee TM , Swancar A ( 2014 ) The suitability of a simplified isotope-balance approach to quantify transient groundwater-lake interactions over a decade with climatic extremes . J Hydrol 519 ( part D ): 3042 - 3053 . doi: 10 .1016/j.jhydrol. 2013 . 12 .012 Simeonov V , Wolska L , Kuczyńska A , Gurwin J , Tsakovski S , Protasowicki M , Namieśnik J ( 2007 ) Sediment-quality assessment by intelligent data analysis . TrAC Trends Anal Chem 26 : 323 - 331 . doi: 10 .1016/j.trac. 2006 . 12 .004 Smerdon BD , Mendoza CA , Devito KJ ( 2007 ) Simulations of fully coupled lake-groundwater exchange in a subhumid climate with an integrated hydrologic model . Water Resour Res 43 : W01416 . doi: 10 .1029/2006WR005137 Sophocleous M ( 2002 ) Interactions between groundwater and surface water: the state of the science . Hydrogeol J 10 : 348 - 348 . doi: 10 . 1007/s10040-002-0204-x Stauffer RE ( 1985 ) Use of solute tracers released by weathering to estimate groundwater inflow to seepage lakes . Environ Sci Technol 19 : 405 - 411 . doi: 10 .1021/es00135a003 Su X , Cui G , Du S , Yuan W , Wang H ( 2016 ) Using multiple environmental methods to estimate groundwater discharge into an arid lake (Dakebo Lake, Inner Mongolia , China). Hydrogeol J 24 ( 7 ): 1707 - 1722 . doi: 10 .1007/s10040-016-1439-2 Taviani S , Henriksen H ( 2015 ) The application of a groundwater/surfacewater model to test the vulnerability of Bracciano Lake (near Rome, Italy) to climatic and water-use stresses . Hydrogeol J 23 ( 7 ): 1481 - 1498 . doi: 10 .1007/s10040-015-1271-0 Vaeret L , Kelbe B , Haldorsen S , Taylor RH ( 2009 ) A modelling study of the effects of land management and climatic variations on groundwater inflow to Lake St Lucia, South Africa . Hydrogeol J 17 : 1949 - 1967 . doi: 10 .1007/s10040-009-0476-5 Virdi ML , Lee TM , Swancar A , Niswonger RG ( 2013 ) Simulating the effect of climate extremes on groundwater flow through a lakebed . Ground Water 51 : 203 - 218 . doi: 10 .1111/j.1745- 6584 . 2012 . 00969 .x Wang D , Wang G , Anagnostou EN ( 2007 ) Evaluation of canopy interception schemes in land surface models . J Hydrol 347 : 308 - 318 Wildi W ( 2010 ) Environmental hazards of dams and reservoirs. NEAR curriculum in natural environmental science . Terre Environ 88 : 187 - 197 Winston RB ( 2009 ) ModelMuse: a graphical user interface for MODFLOW-2005 and PHAST . US Geological Survey Techniques and Methods 6-A29 , USGS, Reston, VA, 52 pp. Available online at http://pubs.usgs.gov/tm/tm6A29. Accessed July 2017 Winter TC , LaBaugh JW , Rosenberry PO ( 1988 ) The design and use of a hydraulic potentiomanometer for direct measurement of differences in hydraulic head between groundwater and surface water . Limnol Oceanogr 33 : 1209 - 1214 Yihdego Y , Becht R ( 2013 ) Simulation of lake-aquifer interaction at Lake Naivasha, Kenya using a three-dimensional flow model with the high conductivity technique and a DEM with bathymetry . J Hydrol 503 : 111 - 122 . doi: 10 .1016/j.jhydrol. 2013 . 08 .034 Zhang Y , Gable CW , Person M ( 2006 ) Equivalent hydraulic conductivity of an experimental stratigraphy: implications for basin-scale flow simulations . Water Resour Res 42 : W05404 . doi: 10 .1029/ 2005WR004720 Zhang Y , Person M , Gable CW ( 2007 ) Representative hydraulic conductivity of hydrogeologic units: insights from an experimental stratigraphy . J Hydrol 339 : 65 - 78 . doi: 10 .1016/j.jhydrol. 2007 . 03 .007


This is a preview of a remote PDF: https://link.springer.com/content/pdf/10.1007%2Fs10040-017-1641-x.pdf

A. A. El-Zehairy, M. W. Lubczynski, J. Gurwin. Interactions of artificial lakes with groundwater applying an integrated MODFLOW solution, Hydrogeology Journal, 2017, 1-24, DOI: 10.1007/s10040-017-1641-x