Interactions of artificial lakes with groundwater applying an integrated MODFLOW solution

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.


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 Electronic supplementary material The online version of this article (doi:10.1007/s10040-017-1641-x) contains supplementary material, which is available to authorized users. 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 lakegroundwater 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 lakegroundwater 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.

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 km 2 . 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 km 2 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 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. 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 (ET o ) using the "FAO Penman Monteith" method (Allen et al. 1998) was estimated, and then, by using the "single crop coefficient" method, the ET o 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 lakenot 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 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 km 2 . 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 Hydrogeological cross-section of Lake Turawa and its surroundings (section 1-1′ in Fig. 1), modified after Gurwin et al. (2004) 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 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 (R g ) and groundwater evapotranspiration (ET g ). 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 R g and ET g are calculated internally in the UZF1 Package for each time step together with groundwater exfiltration (Exf gw , 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) based on pumping tests, slug tests and available hydrogeological knowledge. The assigned horizontal hydraulic conductivity (K h ) 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 (K v ) was assigned through an anisotropy factor K h /K v equal to 10, following local and general literature information about glacial outwash deposits Smerdon et al. 2007). The specific yield (S y ) and the specific storage (S s ) were estimated based on borehole sample analysis and assigned as spatially uniform, for the first unconfined aquifer as S y = 0.15, for the aquitard and for the second confined aquifer as S s = 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 m 2 day −1 based on the K h 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 (R O ) 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 (K v ), in addition to actual infiltration (P a ) through the soil, there is excess infiltration runoff (R Oei ) also known as Hortonian runoff, while if the water table goes higher than the land surface (groundwater exfiltration), then saturation excess runoff (R Osat ), also known as Dunnian runoff, takes place (Niswonger et al. 2006). The R Osat can include two components, rejected infiltration due to shallow groundwater levels, and groundwater exfiltration (Exf gw ) 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 P a is divided into gross groundwater recharge (R g ) and unsaturated zone evapotranspiration (ET uz ) and the remaining is a change of unsaturated zone storage (ΔS uz ). First, ET uz 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 ET g ), 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 (m 3 m −3 ); that value was also used as the initial soil water content, while the soil saturated water content was assumed as 0.30 (m 3 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 (K h ), 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 m 3 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 (K h ). The same zones were used to adjust the specific yield (S y ) and specific storage (S s ) 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 (K h and its corresponding value of K v assigned as 0.1K h ) 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 (K h ) and specific yield (S y ) of the upper aquifer, specific storage (S s ) of the lower aquifer, the anisotropy factor (K h /K v ), 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 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 km 2 ) can be expressed as follows: where P is precipitation rate, Q in is stream inflows at the inlet of the modelled area, Q out is stream outflows at the outlet of the modelled area, GHB in is lateral groundwater inflow into the modelled area across the GHB boundaries, GHB out is lateral groundwater outflow from the modelled area across the GHB boundaries, ET is total evapotranspiration consisting of land surface evapotranspiration (ET LD ) and lake evaporation (E LK ), and ΔS is total change of storage. The total evapotranspiration (ET) and total change in storage (ΔS) can be expressed as follows: where I is canopy interception, ET uz is unsaturated zone evapotranspiration, ET g is groundwater evapotranspiration, ΔS g is storage change of saturated zone, ΔS uz is storage change of unsaturated zone and ΔS LK is the lake storage change.
The lake-water balance extracted from the Lake Package output file can be expressed as follows: where Q LKin is stream inflow (rivers and drains) at the inlet of the lake, Q LKout is stream outflow at the outlet of the lake, R O is total surface runoff into the lake computed by the UZF1 Package as per Eq. (5), L LKin is seepage of groundwater into the lake, L LKout is seepage from the lake into groundwater. In that equation, R O is expressed as follows: where R Oei is excess infiltration runoff (Hortonian runoff) and R Osat is saturation excess runoff (Dunnian runoff) including groundwater exfiltration to land surface (Exf gw ). The land surface and unsaturated zone water balance are expressed as follows: which can be further divided into the following two equations: where P a is the actual infiltration through the unsaturated zone, and R g is gross recharge. The saturated zone water balance of the two aquifers in addition to the aquitard layer can be expressed as follows: Fig. 7 Schematic diagram of the model setup in Turawa Lake study area; symbols the same as in Eqs.
(1)- (10) where Q GWin is stream seepage to groundwater, Q GWout is groundwater seepage to streams. Following Hassan et al. (2014), the net recharge (R n ) refers to the net amount of water that reaches the saturated zone after deducting groundwater evapotranspiration (ET g ) and groundwater exfiltration (Exf gw ) as follows:

Results and discussion
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 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 K h -distribution, the calibrated K h 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 K h of the upper aquifer ranged from 0.1 to 60 m day −1 . Larger K h changes as compared to initial pumping test assumptions, were made in the confined aquifer; in the northern area K h 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 K h of the lower confined aquifer varied from 0.2 to 70 m day −1 . Also, K v 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 (S y ) values in the upper aquifer. The dominant S y 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 (S s ) 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 ) which had S s = 0.0001 m −1 . The calibrated S s 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 m 2 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 L LKnet ,~3 times increases in the maximum daily average L LKout (similar for λ and K), and 6 and 2 times increases in maximum daily shoreline L LKout , 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 L LKnet , 2.9 and 1.6 times decreases in the maximum daily average L LKout , and 5 and 1.1 times decreases in the daily maximum shoreline L LKout .
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 L LKnet 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 Table 2 The sensitivity of lakebed seepage to changes of lakebed leakance (λ) and shallow aquifer hydraulic conductivity K when testing: (1) 4-years average of daily net lake seepage (L LKnet ); (2) daily average of lake seepage into groundwater at the maximum lake stage (L LKout ); and (3) max shoreline seepage at maximum lake stage (L LKout ); all using 10 and 0.1 multiplication factors of λ and K  Hogeboom et al. (2015), but large L LKnet 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 K h (so also K v = 0.1 K h ) of the upper aquifer when multiplied by 0.1 and 10, resulted in the 1.75 times increase of RMSE between observed and simulated heads (the same for both 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 K h /K v 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 S y 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 S s 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 L LKnet 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 L LKnet 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 A LK in Table 3), so also the proportion between A LK and the land surface area (A LD ) changes, i.e. when one area increases the other decreases because the total model domain area (A T ) remains constant (74.8 km 2 ). 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 A LK , temporally variable A LD or temporally invariant A T . 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 A T or to A LK according to Eqs. (1) and (4), respectively. For example, the 4-year mean lake-water input referenced to the lake area, i.e. IN (A LK ) = 17,864.2 mm year −1 , is very different than if referenced to the total modelled area IN (A T ) = 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), Q in (82.2% of IN), and GHB in (1.8% of IN), whereas the three water outputs include: ET (16.4% of OUT), Q out (83.3% of OUT), and GHB out (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 Q in and Q out 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 GHB in was 5-6 times larger than GHB out , 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  . The presented water balance emphasizes the dominant role of the Q LKin and regulated Q LKout . The losses of E LK are larger than P input to the lake. Also, the losses of L LKout are larger than L LKin , but those losses are dependent on the lake regulation-for example, forcing larger lake outflows (Q LKout ) than natural water inflows (Q LKin ) results in the lowering of lake stage, thus also in the reduction of L LKout . Regarding ΔS LK , among the 4 years analysed, only in 2006 was the ΔS LK 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 Q LKout in 2006, nearly the same as Q LKin which typically is much larger than Q LKout .
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 (K v ) and Exf gw (6.4% of P). On the output side, there are four fluxes: R g (14.4% of P), ET uz (52.6% of P), I (27.3% of P), and R O (8.3% of P); note that the water balance shows a relatively small contribution of R g which is mainly because of the large ET uz .
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: R g (11.8% P), GHB in (11.5% P), Q GWin (1.6% of P), and L LKout (4.2% P); and five output fluxes: Q GWout (10.3% P), ET g (8.9% of P), Exf gw (5.2% P), GHB out (2.2% of P), and L LKin (1.1% P). The main groundwater inputs to the study area are R g and GHB in , while the main groundwater outputs are Q GWout , ET g and Exf gw . 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 km 2 ), 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 to the artificial lake dynamics reported here, which was driven primarily by the balance between Q LKin and man-induced regulation of the Q LKout .
It is remarkable that in artificial TL: (1) not only the lake itself but also the whole catchment domain dynamics were driven by Q LKin and regulated Q LKout , 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 L LKout > L LKin , which was because of solid, low-permeability and thick lakebed isolation, reflected by low lakebed leakance; (3) the contributions of ET g and particularly of Exf gw to the overall water balance were much larger relative to R g 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 (P a ), potential evapotranspiration (PET) and measured lake stages; b gross recharge (R g ), net recharge (R n ), groundwater exfiltration (Exf gw ), groundwater evapotranspiration (ET g ) and unsaturated zone evapotranspiration (ET uz ); note that the water fluxes are referenced to the land area (A LD ) Fig. 10 Monthly averages of groundwater heads (GH), groundwater evapotranspiration (ET g ), groundwater exfiltration (Exf gw ) and net lake seepage (L LKnet ), all versus monthly lake stages within the 4-year simulation period and all referenced to the land area; a and r are the regression and correlation coefficients respectively of lake stages when outlet weir restricts Q LKout ; such behaviour seems to be specific for artificial lakes, creating different system dynamics than in natural lakes; (4) Q GWout > Q GWin 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) GHB out < GHB in 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, depending on the distance from the lake and the 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 L LKnet because of its direct dependence on the hydraulic gradient at the interface between the lake and groundwater. The Exf gw 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 Exf gw with lake stages. The ET g 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 ET g correlation with the lake stages (r = 0.79) was lower than the other two. The R g is the only groundwater flux not directly dependent on WTD, although the WTD determines the storage capacity of unsaturated zone, which if large, enhances ET uz reducing R g ; therefore, R g 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 P a , P, PET and lake stages; the P a calculated by UZF1 Package, depends on P, interception loss (I), the vertical hydraulic conductivity (K v ) of the soil and the saturation degree of the unsaturated zone. The peaks of P a , 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 (ET u z ), gross recharge (R g ), groundwater exfiltration (Exf gw ), groundwater evapotranspiration (ET g ) and net recharge (R n ). The temperate climate of the study area explains that ET uz is~5-10 times larger than ET g and that both follow the same pattern as the PET. The R n (Eq. 10) follows the R g in the months November-February when ET g and Exf gw are negligible, and diverges in April-August when it follows ET g , occasionally being interrupted by peaks of R g . Exf gw 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 R g but with the highest lake stages, dependent mainly not on P but on the balance between MP river inflow (Q LKin ) and regulated reservoir outflow (Q LKout ), 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.

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 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. 13 Spatial and temporal impact of the lake on groundwater heads at the southern transect of a series of fictitious piezometers (see Fig. 6) 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) 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 (L LKin ), from lake to the aquifer (L LKout ) and the net seepage (L LKnet = L LKin -L LKout ). The L LKout , 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 L LKin , 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 L LKnet within the 4-year simulation period was approximately −0.3 mm day −1 , implying that L LKout > L LKin ; 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 km 2 ), 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 lakegroundwater interaction was emphasized, MODFLOW-modelled and field-validated by Kidmose et al. (2011) at the natural Lake Hampen (0.760 km 2 surface area, 3.2 Mm 3 volume and rainfall~900 mm year −1 ) in western Denmark. In their steady-state model, shoreline seepage was highly spatially variable, ranging from L LKout =~−70 mm day −1 to L LKin =~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 L LKout area (14% of the maximum lake area) at the west and a large L LKin area (61% of the maximum lake area) at the east. Outwards from the zero seepage line, the L LKout increased in the western direction towards the dam up to 1.6 mm day −1 , while in the eastern direction, L LKin first increased and then decreased, while the largest L LKin 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 L LKin in between; (3) L LKout at the western side increased, while L LKin 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 L LKout with negligible value in the central part of the lake and the largest in the peripheral zone, with maximum L LKout 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 lakegroundwater 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 L LKin at the upstream recharge side and throughout the zero seepage line, to the largest L LKout at the downstream discharge side of the lake outlet; (2) had dominant L LKin 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% withiñ 1.5 year of the simulation period, while in the TL study, changes were fast and dramatic-for example, the area of L LKout 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 lakegroundwater 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 L LKin in the peripheral parts and L LKout in the deepest central part of the lake. Their L LKin area changed from 6 to 59% and L LKout 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 L LKout 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 flow accounting for spatio-temporally variable, 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 exfiltration, while influenced by large lake stage 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.