Suitability of a Coupled Hydrologic and Hydraulic Model to Simulate Surface Water and Groundwater Hydrology in a Typical North-Eastern Germany Lowland Catchment

: Lowland river basins are characterised by complex hydrologic and hydraulic interactions between the di ﬀ erent subsystems (aerated zone, groundwater, surface water), which may require physically-based dynamically-coupled surface water and groundwater hydrological models to reliably describe these processes. Exemplarily, for a typical north-eastern Germany lowland catchment (Tollense river with about 400 km 2 ), an integrated hydrological model, MIKE SHE, coupled with a hydrodynamic model, MIKE 11, was developed and assessed. Hydrological and hydraulic processes were simulated from 2010 to 2018, covering strongly varying meteorological conditions. To achieve a highly reliable model, the calibration was performed in parallel for groundwater levels and river ﬂows at the available monitoring sites in the deﬁned catchment. Based on sensitivity analysis, saturated hydraulic conductivity, leakage coe ﬃ cients, Manning’s roughness, and boundary conditions (BCs) were used as main calibration parameters. Despite the extreme soil heterogeneity of the glacial terrain, the model performance was quite reasonable in the di ﬀ erent sub-catchments with an error of less than 2% for water balance estimation. The resulted water balance showed a strong dependency on land use intensity and meteorological conditions. During relatively dry hydrological years, actual evapotranspiration (ETa) becomes the main water loss component, with an average of 60%–65% of total precipitation and decreases to 55%–60% during comparatively wet hydrological years during the simulation period. Base ﬂow via subsurface and drainage ﬂow accounts for an approximate average of 30%–35% during wet years and rises up to 35%–45% of the total water budget during the dry hydrological years. This means, groundwater is in lowland river systems the decisive compensator of varying meteorological conditions. The coupled hydrologic and hydraulic model is valuable for detailed water balance estimation and seasonal dynamics of groundwater levels and surface water discharges, and, due to its physical foundation, can be extrapolated to analyse meteorological and land use scenarios. Future work will focus on coupling with nutrient transport and river water quality models.


Introduction
Hydrological processes in moderate climate lowland catchments can be complicated due to complex surface water (SW) and groundwater (GW) interactions, and due to this complexity, precise hydrological water balance information is a prerequisite to develop management practices for the sustainable use of 1.
Possibility to satisfyingly incorporate all relevant processes and boundary conditions and to holistically calibrate the model.

2.
Quantification of the meteorological effects on water balance components and the exchange of GW and SW during dry and wet hydrological years.

3.
Effects of weather variations on seasonal dynamics of GW levels and river discharges in lowlands.

4.
Use of physically-based distributed coupled models and their ability to simulate lowlands with moderate climatic conditions.
were compared on the basis of their data requirement, desired complexity level, and their ability to simulate relevant hydrologic and hydraulic processes. Both SWAT and MIKE SHE are able to simulate the desired hydrological processes shown in Figure 1, but in case of simulating the hydraulic structures and their operations, a coupled MIKE SHE and MIKE 11 model stands out as a suitable simulation tool, as the SWAT model does not take into account hydraulic structures and their operational and control strategies. Both SWAT and the MIKE SHE model require more input data in comparison to HSPF and SWIM. SWIM and SWAT are based on similar governing equations and input data requirements, but SWIM does not simulate ponds, lakes, wetlands, and drainage systems. Both HSPF and SWIM do not provide a satisfying solution to include pump flow and drainage. HSPF is more suitable to study the fate and transport of nitrogen in unsaturated, saturated, and surface flow zones. In order to simulate the backwater effect and control structures and their operational strategies, all these models need to be coupled with hydrodynamic models.
As the above described hydrologic, hydraulic, and infrastructural conditions apply large areas in the North and West of Europe, North Western, North America, and (with different meteorological conditions) also in Asia, the chosen model combination should be interesting for these areas too.
The research main questions of the present study include: 1. Possibility to satisfyingly incorporate all relevant processes and boundary conditions and to holistically calibrate the model. 2. Quantification of the meteorological effects on water balance components and the exchange of GW and SW during dry and wet hydrological years. 3. Effects of weather variations on seasonal dynamics of GW levels and river discharges in lowlands. 4. Use of physically-based distributed coupled models and their ability to simulate lowlands with moderate climatic conditions.

Study Area
Tollense river basin is an example of a typical lowland catchment located in north-eastern Germany ( Figure 2). The Tollense river starts as an outflow of the lake Tollense (at the city Neubrandenburg) Appl. Sci. 2020, 10, 1281 4 of 23 and flows about 68 km through a glacier terrain to its confluence with river Peene in the small town Demmin. The presented study was performed on the downstream section of the Tollense river with an approximate length of 30 km, staring from Klempenow to Demmin, with an average catchment area of about 400 km 2 . The Tollense river catchment is provided with artificial drainage and is primarily used for agricultural activities. Land use maps were developed in ArcGIS based on aerospace images provided by the Rapid Eye Science Archive platform. Supervised image classification performed in ArcGIS resulted that the study area consisted mainly of 18% forests and 70% arable land and pastures [40].

Study Area
Tollense river basin is an example of a typical lowland catchment located in north-eastern Germany ( Figure 2). The Tollense river starts as an outflow of the lake Tollense (at the city Neubrandenburg) and flows about 68 km through a glacier terrain to its confluence with river Peene in the small town Demmin. The presented study was performed on the downstream section of the Tollense river with an approximate length of 30 km, staring from Klempenow to Demmin, with an average catchment area of about 400 km 2 . The Tollense river catchment is provided with artificial drainage and is primarily used for agricultural activities. Land use maps were developed in ArcGIS based on aerospace images provided by the Rapid Eye Science Archive platform. Supervised image classification performed in ArcGIS resulted that the study area consisted mainly of 18% forests and 70% arable land and pastures [40]. Tollense River originating from Tollense lake in Neubrandenburg towards its downstream end in Demmin, and land use in Tollese river catchment from Klempenow to Demmin. © Rivers lakes and land use data was provided by "State office for the Environment, Nature Conservation and Geology Mecklenburg-Vorpommern" (LUNG-MV) [40].

Groundwater Data
To conduct a successful model calibration, it was necessary to improve the existing database of observed GW levels by continuous monitoring. For this, continuous GW data loggers were installed at 7 boreholes in the catchment with an hourly log rate, starting from November 2016 to the end of April 2018. Boreholes for data logger installation were selected on the basis of their geographic location, accessibility, functionality, and protection from theft and animals (Figure 3 (left)). Borehole functionality was determined based on refill tests performed in collaboration with StALU-MS, being in charge of the institutional GW monitoring. Some of the boreholes were identified as clogged due to lack of operational management and some were not found due to misleading coordinates. Figure 3 (right) shows the constructed GW contours in the study area based on recorded GW levels via data loggers. These were further used as an initial condition during simulation.

Groundwater Data
To conduct a successful model calibration, it was necessary to improve the existing database of observed GW levels by continuous monitoring. For this, continuous GW data loggers were installed at 7 boreholes in the catchment with an hourly log rate, starting from November 2016 to the end of April 2018. Boreholes for data logger installation were selected on the basis of their geographic location, accessibility, functionality, and protection from theft and animals (Figure 3 (left)). Borehole functionality was determined based on refill tests performed in collaboration with StALU-MS, being in charge of the institutional GW monitoring. Some of the boreholes were identified as clogged due to lack of operational management and some were not found due to misleading coordinates. Figure 3 (right) shows the constructed GW contours in the study area based on recorded GW levels via data loggers. These were further used as an initial condition during simulation. where green circles are the boreholes found and working; yellow circles are the boreholes found but not working; and red circles are the boreholes not found; (right) constructed groundwater contours in study area.

Climate Data
Accumulated daily precipitation data ranging from 2010 to 2018 was collected for 6 climate stations named as Demmin; Gross KiesowSchlagtow; Anklam; Friedland; Trollenhagen; and Gross Luckow, located within or nearby the defined Tollense river catchment. The Thiessen polygon interpolation method was used to divide the catchment area under each climate station shown in Figure 4 (up). Penman Monteith method [41] was used for potential evapotranspiration (ETp) calculation and it requires relative humidity, cloud cover, wind speed, and daily maximum and minimum wind temperatures to calculate the ETp [42]. For ETp calculation, data from the climate station named as "Demmin" was used as it was the nearest climate station where all the required data was available to calculate the ETp. Figure 4 (down) shows the average monthly ET and precipitation rates for the climate monitoring station "Demmin" from 2010 to 2017. MIKE SHE uses the Kristensen and Jensen method [43] to compute the Eta on the basis of specified ETp rates and soil moisture available in the root zone. where green circles are the boreholes found and working; yellow circles are the boreholes found but not working; and red circles are the boreholes not found; (right) constructed groundwater contours in study area.

Climate Data
Accumulated daily precipitation data ranging from 2010 to 2018 was collected for 6 climate stations named as Demmin; Gross KiesowSchlagtow; Anklam; Friedland; Trollenhagen; and Gross Luckow, located within or nearby the defined Tollense river catchment. The Thiessen polygon interpolation method was used to divide the catchment area under each climate station shown in Figure 4 (up). Penman Monteith method [41] was used for potential evapotranspiration (ETp) calculation and it requires relative humidity, cloud cover, wind speed, and daily maximum and minimum wind temperatures to calculate the ETp [42]. For ETp calculation, data from the climate station named as "Demmin" was used as it was the nearest climate station where all the required data was available to calculate the ETp. Figure 4 (down) shows the average monthly ET and precipitation rates for the climate monitoring station "Demmin" from 2010 to 2017. MIKE SHE uses the Kristensen and Jensen method [43] to compute the Eta on the basis of specified ETp rates and soil moisture available in the root zone.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 5 of 24 where green circles are the boreholes found and working; yellow circles are the boreholes found but not working; and red circles are the boreholes not found; (right) constructed groundwater contours in study area.

Climate Data
Accumulated daily precipitation data ranging from 2010 to 2018 was collected for 6 climate stations named as Demmin; Gross KiesowSchlagtow; Anklam; Friedland; Trollenhagen; and Gross Luckow, located within or nearby the defined Tollense river catchment. The Thiessen polygon interpolation method was used to divide the catchment area under each climate station shown in Figure 4 (up). Penman Monteith method [41] was used for potential evapotranspiration (ETp) calculation and it requires relative humidity, cloud cover, wind speed, and daily maximum and minimum wind temperatures to calculate the ETp [42]. For ETp calculation, data from the climate station named as "Demmin" was used as it was the nearest climate station where all the required data was available to calculate the ETp. Figure 4 (down) shows the average monthly ET and precipitation rates for the climate monitoring station "Demmin" from 2010 to 2017. MIKE SHE uses the Kristensen and Jensen method [43] to compute the Eta on the basis of specified ETp rates and soil moisture available in the root zone.

Land Use Data
Land use, shown in Figure 2, is based on Rapid Eye Science Achieve images and an image classification was performed in ArcGIS to obtain the land use maps. Land use was classified as arable land, wetlands, grassland, industrial areas, forests, small gardens, settlements, concrete surfaces, roads, water facilities, and miscellaneous, based on the image classification results. Table 1 shows the summary of the data used in this study and the related sources from where the data was obtained. Leaf area index (LAI) based on project "Communal waters collective development in urban areas" (KOGGE) [44] in monthly temporal resolution is shown in Figure 5 and average root depth (RD) on seasonal temporal resolution [45] is shown in Table 2.

Land Use Data
Land use, shown in Figure 2, is based on Rapid Eye Science Achieve images and an image classification was performed in ArcGIS to obtain the land use maps. Land use was classified as arable land, wetlands, grassland, industrial areas, forests, small gardens, settlements, concrete surfaces, roads, water facilities, and miscellaneous, based on the image classification results. Table 1 shows the summary of the data used in this study and the related sources from where the data was obtained. Leaf area index (LAI) based on project "Communal waters collective development in urban areas" (KOGGE) [44] in monthly temporal resolution is shown in Figure 5 and average root depth (RD) on seasonal temporal resolution [45] is shown in Table 2. Table 1. Summary of the data used in this study and the sources from where the data is obtained.

MIKE SHE Model Setup
The MIKE SHE model uses different process-oriented components to describe the hydrological cycle by coupling the UZ and SZ, as shown in Table 3. MIKE SHE uses three different methods to account for UZ flow simulations, named as two-layer UZ method, Richards equation, and gravity flow method. The MIKE SHE model uses the Kristensen-Jensen method [43] to estimate the ET and interception processes by using climatic and vegetative data. Due to the computational time constraints, a simplified two-layer water balance method was used to estimate the UZ flow based on the formulation provided by Yan and Smith [46]. The main objective of the module is to compute the ETa and the volume of GW recharge and is mainly advantageous in areas with a shallow GW table (e.g., wetlands where the ETa is nearly equal to the reference ET rate). In catchments with higher UZ depths, the two-layer UZ method does not precisely characterise the UZ flow dynamics. The two-layer UZ method considers only the average conditions and does not account for the UZ hydraulic conductivity and soil moisture content relationship, and hence does not consider the soil ability to transport water to the plant roots. That means if sufficient water is available in the root

MIKE SHE Model Setup
The MIKE SHE model uses different process-oriented components to describe the hydrological cycle by coupling the UZ and SZ, as shown in Table 3. MIKE SHE uses three different methods to account for UZ flow simulations, named as two-layer UZ method, Richards equation, and gravity flow method. The MIKE SHE model uses the Kristensen-Jensen method [43] to estimate the ET and interception processes by using climatic and vegetative data. Due to the computational time constraints, a simplified two-layer water balance method was used to estimate the UZ flow based on the formulation provided by Yan and Smith [46]. The main objective of the module is to compute the ETa and the volume of GW recharge and is mainly advantageous in areas with a shallow GW table (e.g., wetlands where the ETa is nearly equal to the reference ET rate). In catchments with higher UZ depths, the two-layer UZ method does not precisely characterise the UZ flow dynamics. The two-layer UZ method considers only the average conditions and does not account for the UZ hydraulic conductivity and soil moisture content relationship, and hence does not consider the soil ability to transport water to the plant roots. That means if sufficient water is available in the root zone it will be available for ET. However, by calibrating the input parameters, the two-layer UZ method performs quite well in most of the conditions. The defined method takes into account the processes of canopy interception, ponding, and ET and considers the whole UZ to consist of two layers representing average conditions in the UZ, where vegetation data is defined as LAI and RD. The soil properties are defined by a constant infiltration capacity, soil moisture content at wilting point, field capacity, and saturated water content. The two-layer UZ method fulfils the main objectives to account for ETa and the SZ recharge mainly required in this study. The MIKE SHE model uses the 3D-Boussinesq equation to simulate the SZ flows [36], while river hydrodynamics (stage and discharge) are defined by using the simplified diffusion wave approximation of the Saint-Venant equation. Table 3. MIKE SHE water modelling components and related numerical solutions [3,5,47,48].

Component
Numerical Methods

Evapotranspiration (ET)
Kristensen and Jensen method [43] Two-layer water balance method [46] Unsaturated zone (UZ) Flow Due to the limited ability of MIKE SHE to simulate the hydraulic structures and their operations, MIKE 11, a 1D-hydrodynamic model, was used to incorporate the control structures such as the weirs, sluice, culverts, gates, etc., and their time varying operations. The MIKE SHE model used in this study uses the "ETRS_1989_UTM_Zone_33N_8stellen" coordinate system in metric units. The performance of the hydrological model was evaluated at different grid sizes (20,50,100,200, and 300 m) and a spatial resolution of 100 m was selected as an optimum compromise between computational time and model performance. The model area of 400 km 2 was divided into 100 × 100 m computational cells to represent the spatial variations in soil, land use, and geology. Model area was divided into computational grids in order to have the numerical solutions of the governing equations [49][50][51]. DEM with a resolution of 5 m was used to represent the topography in the study area and was provided by LUNG-MV. In case of temporal discretization, the computational time steps were chosen in such a way that OL-flow computational time step is always less or equal to the UZ time step, while the UZ time step is always less or equal to the SZ time step. During the temporal discretization, selected time steps are always critical for minimizing the water balance errors [3]. Aquifer BC used in this study are explained in Figure 6; zero flux BC was selected where lateral flows were expected to be minimum based on the GW contours and topographical characteristics, while in case of gradient BC, selected gradient values were estimated based on GW contours. The two-layer UZ method was used to predict the flow through or within the UZ. Soil data and Van Genuchten parameters [52,53] were obtained from LUNG-MV and were used to develop the detailed soil map. The soil column vertical discretization was kept constant in the entire model with a minimum soil thickness of 0.5 cm at the ground surface to comprehensively quantify the nonlinear evaporation and transpiration to the maximum vertical discretization of 1 m at the model outer depth. The SZ was characterised by three layers of geology and a three-dimensional Boussinesq equation was used to estimate the SZ flows. The depth of the aquifer was estimated based on geological data obtained from 60 boreholes available in the study area, and these interpolated GW levels were used as SZ initial conditions. Fine tuning of vertical and horizontal hydraulic conductivities of geological layers played a vital part in the successful model calibration.

Coupling of MIKE 11 and MIKE SHE
MIKE 11 uses complete dynamic wave formulation of the Saint-Venant equation and is able to simulate the wide range of hydraulic structures such as weirs, culverts, gates, and bridges, etc. The river network in the defined catchment containing the main Tollense River, its tributaries, and control hydraulic structures are shown in Figure 7 (right). A fully-dynamic 1D-Saint-Venant equation was used to describe channel flow in the Tollense river and its tributaries. Two weirs, named as Osten and Tückhude, equipped with adjustable gates and variable operational strategy, were modelled in MIKE 11, as shown in Figure 7 (left). Channel cross-sections were based on ADCP surveys conducted at 29 locations along the Tollense river, in the river section starting from Demmin (downstream) to Klempenow (upstream). The river bed was considered fully connected with the SZ, as a result the exchange between the river and the GW aquifer was defined by the hydraulic conductivity of the aquifer instead of the river bed material. In case of the Tollense River, time series of discharge at daily temporal resolution was inserted as an upstream BC at Klempenow and average water levels were used as a downstream BC at Demmin. The data used in the present study in the MIKE11 model includes river network, control strategy of hydraulic structures, river cross sections, hydraulic structures and their geometry, seasonal Manning's roughness coefficients, river boundary conditions, and hydrodynamic parameters.

Coupling of MIKE 11 and MIKE SHE
MIKE 11 uses complete dynamic wave formulation of the Saint-Venant equation and is able to simulate the wide range of hydraulic structures such as weirs, culverts, gates, and bridges, etc. The river network in the defined catchment containing the main Tollense River, its tributaries, and control hydraulic structures are shown in Figure 7 (right). A fully-dynamic 1D-Saint-Venant equation was used to describe channel flow in the Tollense river and its tributaries. Two weirs, named as Osten and Tückhude, equipped with adjustable gates and variable operational strategy, were modelled in MIKE 11, as shown in Figure 7 (left). Channel cross-sections were based on ADCP surveys conducted at 29 locations along the Tollense river, in the river section starting from Demmin (downstream) to Klempenow (upstream). The river bed was considered fully connected with the SZ, as a result the exchange between the river and the GW aquifer was defined by the hydraulic conductivity of the aquifer instead of the river bed material. In case of the Tollense River, time series of discharge at daily temporal resolution was inserted as an upstream BC at Klempenow and average water levels were used as a downstream BC at Demmin. The data used in the present study in the MIKE11 model includes river network, control strategy of hydraulic structures, river cross sections, hydraulic structures and their geometry, seasonal Manning's roughness coefficients, river boundary conditions, and hydrodynamic parameters.
A dynamic coupling of MIKE SHE and MIKE 11 considers the exchange of data among the two models after every computational time step and their coupling is done by means of line segments, serving as river links in MIKE 11 to the adjacent MIKE SHE grids. River links are created for coupled reaches and their locations are defined automatically from MIKE 11 river point coordinates that define the river branches of a hydrodynamic model. During simulation, the water levels within the coupled reaches are transported from MIKE 11 H-points to the adjacent MIKE SHE river links, as shown in Figure 8. In return, MIKE SHE estimates the OL-flow to each river link as an exchange between the neighbouring grid squares and the river aquifer, and finally these terms are sent back to the corresponding MIKE 11 H-points as an outflow or inflow for the following computational time step. If the water level in a grid square gets higher than its topographic level, it is considered as flooded, and as soon as a grid square gets flooded MIKE SHE computes infiltration, seepage, OL-flow, and ET in the similar manner as for a grid square with SW ponding resulting from precipitation and surface runoff, or the water table intercepting the ground surface [54].
Appl. Sci. 2020, 10, x FOR PEER REVIEW 10 of 24 cross sections, hydraulic structures and their geometry, seasonal Manning's roughness coefficients, river boundary conditions, and hydrodynamic parameters. A dynamic coupling of MIKE SHE and MIKE 11 considers the exchange of data among the two models after every computational time step and their coupling is done by means of line segments, serving as river links in MIKE 11 to the adjacent MIKE SHE grids. River links are created for coupled reaches and their locations are defined automatically from MIKE 11 river point coordinates that define the river branches of a hydrodynamic model. During simulation, the water levels within the coupled reaches are transported from MIKE 11 H-points to the adjacent MIKE SHE river links, as shown in Figure 8. In return, MIKE SHE estimates the OL-flow to each river link as an exchange between the neighbouring grid squares and the river aquifer, and finally these terms are sent back to the corresponding MIKE 11 H-points as an outflow or inflow for the following computational time step. If the water level in a grid square gets higher than its topographic level, it is considered as flooded, and as soon as a grid square gets flooded MIKE SHE computes infiltration, seepage, OLflow, and ET in the similar manner as for a grid square with SW ponding resulting from precipitation and surface runoff, or the water table intercepting the ground surface [54].  A dynamic coupling of MIKE SHE and MIKE 11 considers the exchange of data among the two models after every computational time step and their coupling is done by means of line segments, serving as river links in MIKE 11 to the adjacent MIKE SHE grids. River links are created for coupled reaches and their locations are defined automatically from MIKE 11 river point coordinates that define the river branches of a hydrodynamic model. During simulation, the water levels within the coupled reaches are transported from MIKE 11 H-points to the adjacent MIKE SHE river links, as shown in Figure 8. In return, MIKE SHE estimates the OL-flow to each river link as an exchange between the neighbouring grid squares and the river aquifer, and finally these terms are sent back to the corresponding MIKE 11 H-points as an outflow or inflow for the following computational time step. If the water level in a grid square gets higher than its topographic level, it is considered as flooded, and as soon as a grid square gets flooded MIKE SHE computes infiltration, seepage, OLflow, and ET in the similar manner as for a grid square with SW ponding resulting from precipitation and surface runoff, or the water table intercepting the ground surface [54].

Sensitivity Analyses and Calibration Procedure
Sensitivity analyses were performed in order to assess the most sensitive parameters for Tollense river catchment model calibration and were performed by changing different parameters within their allowable range. During the sensitivity analyses, only one parameter was adjusted at once, while others kept unchanged during each sensitivity test. Table 4 shows the considered variables and their definitions in the MIKE SHE and MIKE 11 model. Due to the close interaction of GW and SW, calibration of GW levels and SW discharges was conducted in parallel. Especially the calibration of GW level became challenging to the extreme soil heterogeneity, which is typical for glacial soils. Auto calibration was not applied in this study, as the intention was not to achieve an optimum fit between the observed and simulated data, as the model setup without auto calibration provides better process understanding and helps to identify the model deficiencies. Calibration is a process of model testing with known inputs and outputs and is used to estimate or adjust different parameters. The calibration process was initiated after the sensitivity analysis and following calibration techniques were applied for the different river flow calibration outcomes. In case 1-when the model was unable to simulate the peak flows, a careful review of the data was done for both the precipitation and river flow, as this happens normally due to the reason that either the rainfall station is not representative or due to the malfunctioning of the precipitation or flow gauges. In case 2-when the model continuously over predicted the flow, ET, soil water content, percolation, and GW recharge rates were adjusted. In case 3-where simulated flow followed the observed flow dynamics but lags the actual flow consistently, river bed leakage coefficient rates and Manning's roughness coefficient were adjusted. In case 4-when the model over predicted the peak flows but under predicted the average normal river flows, infiltration rate, interflow, and base flow recession parameters were adjusted. The coupled MIKE SHE and MIKE 11 model was calibrated by using the measured GW heads at seven selected GW observation wells shown in Figure 9, for the period ranging from November 2016 to April 2018. River flow calibration was based on flow measurements performed by StALU-MS and University of Rostock along the Tollense river at four different chainages in the river section from Klempenow to Demmin. Calibrated GW level and surface flow results are presented in Figures 10 and 11, respectively, for a period of nine years starting from 2010 to 2018. During the calibration process, hydraulic conductivity, BCs, roughness coefficient, and specific yield were adjusted to calibrate the coupled MIKE SHE and MIKE 11 model and their calibrated values are shown in Table 5.

Groundwater Dynamics
Simulated GW levels were calibrated against hourly observed GW data obtained from the data loggers installed in the study area for the period ranging from November 2016 to April 2018. In the present study, GW data from seven boreholes located in the Tollense catchment were used. Borehole HySzw 141/1989 is located near to the Tollense River, HyStav 7/2005 and HyBorr 106/1985 are directly at the bank of Gemkow tributary, while other boreholes are located at banks of small tributaries present in different sub catchments of the Tollense river basin shown in Figure 9.

Groundwater Dynamics
Simulated GW levels were calibrated against hourly observed GW data obtained from the data loggers installed in the study area for the period ranging from November 2016 to April 2018. In the present study, GW data from seven boreholes located in the Tollense catchment were used. Borehole HySzw 141/1989 is located near to the Tollense River, HyStav 7/2005 and HyBorr 106/1985 are directly at the bank of Gemkow tributary, while other boreholes are located at banks of small tributaries present in different sub catchments of the Tollense river basin shown in Figure 9. The statistical performance of the coupled MIKE SHE and MIKE 11 model in terms of predicting the GW dynamics is shown in Table 6, whereas Figure 10 shows the comparison of simulated and observed GW levels at the selected seven borehole locations located in different sub catchments of the Tollense river basin. A good agreement was observed between simulated and observed GW levels. but simulated GW reacts more dynamically. A probable reason is the selection of too course spatial and temporal discretization to simulate strong soil heterogeneity of the catchment. As expected, higher GW level occurs during the winter season, then dropping to lower levels during the summer season with increasing ET. In general, the coupled model underestimates the GW levels during the high recharge periods, but the performance of the model was not equally comparable at all the observation locations. Due to the rather shallow GW tables, there is a strong relationship spatially and temporally between the precipitation and resulted simulated GW levels. GW elevations react dynamically to local precipitation events. Boreholes located in the vicinity of the climate station "Demmin" showed that the average simulated GW elevations across the middle and downstream section varies up to around 0.15 m during an extra ordinary precipitation event The statistical performance of the coupled MIKE SHE and MIKE 11 model in terms of predicting the GW dynamics is shown in Table 6, whereas Figure 10 shows the comparison of simulated and observed GW levels at the selected seven borehole locations located in different sub catchments of the Tollense river basin. A good agreement was observed between simulated and observed GW levels. but simulated GW reacts more dynamically. A probable reason is the selection of too course spatial and temporal discretization to simulate strong soil heterogeneity of the catchment. As expected, higher GW level occurs during the winter season, then dropping to lower levels during the summer season with increasing ET. In general, the coupled model underestimates the GW levels during the high recharge periods, but the performance of the model was not equally comparable at all the observation locations. Due to the rather shallow GW tables, there is a strong relationship spatially and temporally between the precipitation and resulted simulated GW levels. GW elevations react dynamically to local precipitation events. Boreholes located in the vicinity of the climate station "Demmin" showed that the average simulated GW elevations across the middle and downstream section varies up to around 0.15 m during an extra ordinary precipitation event that occurred in October 2017. Coupled model performance was evaluated by using mean absolute error, root mean square error, correlation coefficient, and standard deviation residuals. Initially climate parameters were adjusted according to the basin climate conditions and after that the calibration process was focused mainly on horizontal and vertical conductivity and specific storage of the SZ parameters. As a result, the calibration process includes seasons with different GW dynamics due to variation in rainfall patterns and resulted varied GW levels. Specific yield and GW gradients turned out to be the most adjusted parameters, as boundary gradients normally determine the GW influx, while specific yield defines the dynamic response of the GW system, hence both have an impact on estimated water balance. * MAE = "mean absolute error in level m"; RMSE = "root mean square error in level m"; R = "correlation"; STDres = "standard deviation residuals". Appl. Sci. 2020, 10, x FOR PEER REVIEW 14 of 23 Figure 10. Observed (green) and calibrated (black) groundwater heads in the Tollense River catchment with rainfall (blue) and evapotranspiration (red) rates.

Stream Flow Dynamics
The purpose behind the river flow calibration and validation was to compare simulated and observed hydrographs at the catchment outlet as adjustment of parameters for GW calibration has a direct impact on surface flows. The outflow hydrograph has the potential of representing the integrated effect of the hydrological response of the catchment as a whole. Flow measurement campaigns were launched under the BMBF project named as "Boat monitoring" [55], and flow was measured at different chainages along the river Tollense. River flow was calibrated at four chainages at 270, 1365, 7048, and 23,341 m from the downstream end of the Tollense river. The statistic quality criteria concerning river dynamics (mean absolute error, root mean square error, correlation coefficient, and standard deviation residuals) are summarised in Table 7. The simulated results were satisfactory and were able to predict the general characteristics of the discharge time series at specified chainages. Simulated hydrographs underestimate the river flows during the observed low flow periods, while overestimate the river flows during the observed peak flow discharges in the river Tollense. The probable reason behind this discrepancy is the artificial drainage in the catchment. Artificial drainage was constructed based on DEM lowest points and does not fully represent the real installed drainage system in the catchment. GW levels lower than artificial drainage makes the drainage ineffective and results in small GW contribution to the river. For periods with GW levels higher than constructed drainage in MIKE SHE, drainage becomes effective and may contribute more to the river Tollense in comparison to reality. Moreover, average rainfall over the entire area in each polygon assigned to its respective climate station, average hydraulic conductivity of an entire geological layer, and drainage time constants also effects the GW contribution to the SW flows, and cause differences between observed and simulated river discharges. Table 7. Statistics of the objective function for river flow calibration. The purpose behind the river flow calibration and validation was to compare simulated and observed hydrographs at the catchment outlet as adjustment of parameters for GW calibration has a direct impact on surface flows. The outflow hydrograph has the potential of representing the integrated effect of the hydrological response of the catchment as a whole. Flow measurement campaigns were launched under the BMBF project named as "Boat monitoring" [55], and flow was measured at different chainages along the river Tollense. River flow was calibrated at four chainages at 270, 1365, 7048, and 23,341 m from the downstream end of the Tollense river. The statistic quality criteria concerning river dynamics (mean absolute error, root mean square error, correlation coefficient, and standard deviation residuals) are summarised in Table 7. The simulated results were satisfactory and were able to predict the general characteristics of the discharge time series at specified chainages. Simulated hydrographs underestimate the river flows during the observed low flow periods, while overestimate the river flows during the observed peak flow discharges in the river Tollense. The probable reason behind this discrepancy is the artificial drainage in the catchment. Artificial drainage was constructed based on DEM lowest points and does not fully represent the real installed drainage system in the catchment. GW levels lower than artificial drainage makes the drainage ineffective and results in small GW contribution to the river. For periods with GW levels higher than constructed drainage in MIKE SHE, drainage becomes effective and may contribute more to the river Tollense in comparison to reality. Moreover, average rainfall over the entire area in each polygon assigned to its respective climate station, average hydraulic conductivity of an entire geological layer, and drainage time constants also effects the GW contribution to the SW flows, and cause differences between observed and simulated river discharges. Table 7. Statistics of the objective function for river flow calibration.

Water Balance
A water balance estimation was performed following the calibration process as the interest of this study was to get a better insight into the interaction of the different subsystems and hydrological processes in highly GW influenced lowland river basins. The water balance was performed for the whole calibration period and additionally for each hydrologic year, which is defined between 1st October to 30th September of the next year in Germany. The total water input into the model was via precipitation and surface and subsurface inflow, which was further divided into ET, runoff and change in storage, and GW and SW interactions. Water balance error was calculated after balancing all the major hydrological components that includes precipitation, ET, runoff, surface and subsurface inflow, and change in storage. With a water balance error of less than 2% during the calibration period, the total estimated water balance is satisfyingly good. Total water balance: Overall water balance results shown in Figure 12a clearly show that a considerable portion of precipitations is going back to the atmosphere as an ET loss that represents approximate average of 60% of the total precipitation. Areas of small lakes, ponds, and forests contribute to the major part of ET loss. GW contribution to the SW flow is around 30%-35% via drainage, while SZ storage has decreased over the past nine simulated years, and this is evident in the simulated GW levels in the lowland catchment. Tollense river catchment has a very small cross boundary SZ inflow and GW storage is mainly dependent on local meteorological conditions, such as precipitation and ET.
In the lowland catchment, the UZ is very shallow during the wet seasons. Due to this, infiltration and ET are vital processes that control the rate of recharge. Results illustrate that the SZ gains approximately 30%-40% of total precipitation, out of which a major portion is drained to the river via artificial drainage available in the study area.
Water balance in wet and dry hydrological years: The water balance results for a wet hydrological year (2010-2011) and dry hydrological year (2017-2018) are shown in Figure 12b,c. During the wet hydrological year, ET is 56% of total P, while drainage contributes 34% of P into the river Tollense with a positive GW storage. During the dry hydrological year, ET raises up to 65% of total P with a raised drainage contribution of 46% in to the river flows, resulting in a negative GW

Water Balance
A water balance estimation was performed following the calibration process as the interest of this study was to get a better insight into the interaction of the different subsystems and hydrological processes in highly GW influenced lowland river basins. The water balance was performed for the whole calibration period and additionally for each hydrologic year, which is defined between 1st October to 30th September of the next year in Germany. The total water input into the model was via precipitation and surface and subsurface inflow, which was further divided into ET, runoff and change in storage, and GW and SW interactions. Water balance error was calculated after balancing all the major hydrological components that includes precipitation, ET, runoff, surface and subsurface inflow, and change in storage. With a water balance error of less than 2% during the calibration period, the total estimated water balance is satisfyingly good.
Total water balance: Overall water balance results shown in Figure 12a clearly show that a considerable portion of precipitations is going back to the atmosphere as an ET loss that represents approximate average of 60% of the total precipitation. Areas of small lakes, ponds, and forests contribute to the major part of ET loss. GW contribution to the SW flow is around 30%-35% via drainage, while SZ storage has decreased over the past nine simulated years, and this is evident in the simulated GW levels in the lowland catchment. Tollense river catchment has a very small cross boundary SZ inflow and GW storage is mainly dependent on local meteorological conditions, such as precipitation and ET.
In the lowland catchment, the UZ is very shallow during the wet seasons. Due to this, infiltration and ET are vital processes that control the rate of recharge. Results illustrate that the SZ gains approximately 30%-40% of total precipitation, out of which a major portion is drained to the river via artificial drainage available in the study area.
Water balance in wet and dry hydrological years: The water balance results for a wet hydrological year (2010-2011) and dry hydrological year (2017-2018) are shown in Figure 12b,c. During the wet hydrological year, ET is 56% of total P, while drainage contributes 34% of P into the river Tollense with a positive GW storage. During the dry hydrological year, ET raises up to 65% of total P with a raised drainage contribution of 46% in to the river flows, resulting in a negative GW storage. The water balance results show that GW is a main contributor to the surface flows during low rainfall or partial drought periods and balance SW discharges in the Tollense river catchment.
as precipitation and ET.
In the lowland catchment, the UZ is very shallow during the wet seasons. Due to this, infiltration and ET are vital processes that control the rate of recharge. Results illustrate that the SZ gains approximately 30%-40% of total precipitation, out of which a major portion is drained to the river via artificial drainage available in the study area.
Water balance in wet and dry hydrological years: The water balance results for a wet hydrological year (2010-2011) and dry hydrological year (2017-2018) are shown in Figure 12b,c. During the wet hydrological year, ET is 56% of total P, while drainage contributes 34% of P into the river Tollense with a positive GW storage. During the dry hydrological year, ET raises up to 65% of total P with a raised drainage contribution of 46% in to the river flows, resulting in a negative GW storage. The water balance results show that GW is a main contributor to the surface flows during low rainfall or partial drought periods and balance SW discharges in the Tollense river catchment.

Discussions
The bi-directional coupled hydrologic and hydraulic model was applied to the Tollense river catchment, which represents many common features of the European lowland catchments, such as shallow GW tables, interacting with SW, control structures, artificial drainage, and periodic inundation, etc. Due to these characteristics, modelling described in this study offers a huge potential to predict the lowlands response to anthropogenic activities and expected changing climatic conditions and to provide guidelines for management and conservation practices of vulnerable lowland catchments.

Coupled Model Performance
In general, the coupled model performance in the Tollense river catchment was satisfyingly able to describe the interacting hydrologic subsystems GW and SW, assessed by comparison with observed GW levels and SW discharges. Modelling results demonstrate a close association between P, ET, SW discharges, and GW levels. Despite the good representation of SW and GW dynamics, coupled model performance was not equally the same in all the sub catchments due to heterogeneity of soils and variable availability of field monitoring data. Grid resolution is very important to define the heterogeneity of the catchment according to the desired level of complexity. The finally selected resolution was based on a good compromise between simulation time, numeric robustness, and resulted accuracy. But even with a very fine grid it would not be feasible to represent small-scale variations of geology, which is only exactly known at the borehole locations from the drilling documents. Besides the required interpolation/generalization of geology, necessary estimates on the exact layout of the artificial drainage and drainage time constants are further uncertainty factors of the GW model. Due to private rights of the farmer's, drainage maps were not available during the study, thus artificial drainage was constructed based on lowest points of DEM in the defined lowland catchment. The calibration process showed that river flows are very responsive to the Manning's n, drainage time constant, and leakage coefficients. Water balance

Discussions
The bi-directional coupled hydrologic and hydraulic model was applied to the Tollense river catchment, which represents many common features of the European lowland catchments, such as shallow GW tables, interacting with SW, control structures, artificial drainage, and periodic inundation, etc. Due to these characteristics, modelling described in this study offers a huge potential to predict the lowlands response to anthropogenic activities and expected changing climatic conditions and to provide guidelines for management and conservation practices of vulnerable lowland catchments.

Coupled Model Performance
In general, the coupled model performance in the Tollense river catchment was satisfyingly able to describe the interacting hydrologic subsystems GW and SW, assessed by comparison with observed GW levels and SW discharges. Modelling results demonstrate a close association between P, ET, SW discharges, and GW levels. Despite the good representation of SW and GW dynamics, coupled model performance was not equally the same in all the sub catchments due to heterogeneity of soils and variable availability of field monitoring data. Grid resolution is very important to define the heterogeneity of the catchment according to the desired level of complexity. The finally selected resolution was based on a good compromise between simulation time, numeric robustness, and resulted accuracy. But even with a very fine grid it would not be feasible to represent small-scale variations of geology, which is only exactly known at the borehole locations from the drilling documents. Besides the required interpolation/generalization of geology, necessary estimates on the exact layout of the artificial drainage and drainage time constants are further uncertainty factors of the GW model. Due to private rights of the farmer's, drainage maps were not available during the study, thus artificial drainage was constructed based on lowest points of DEM in the defined lowland catchment. The calibration process showed that river flows are very responsive to the Manning's n, drainage time constant, and leakage coefficients. Water balance estimation during dry and wet years showed the interaction of different water balance components, where ET is a major water loss component, and during dry years it reaches up to 65% of total water budget and results in lowering of GW levels due to contribution of GW to SW discharges under minimum GW recharge rates. The coupled model satisfactorily represented the water balance with an error of less than 2% of total water budget.

GW and SW Interaction
The coupled model resulted in intense interactions among SW and GW. Long term, SW flows follow the pattern of GW levels in the defined catchment with the higher GW flows followed by higher SW discharges. Model calibration of SW discharges was difficult due to the limited monitoring data availability of SW discharges, and that highlights the significance of high resolution field monitoring data in hydrological modelling. GW is a major contributor to the balance of SW flows during the low flow periods, and GW contribution rises up to 45% of total SW flows during the observed partial drought in the catchment; and, with exclusion of river flows from lake Tollense, GW contributes mainly to the Tollense river nearly up to 95% of total SW flows. The simulated hydrograph shows relatively overestimated river flows during peak flow periods. A successful calibration of SZ boundary conditions and geological layers' vertical discretization, saturated hydraulic conductivities, drainage time constant, and leakage coefficient play a vital role to successfully quantify the SW and GW interactions. The above discussed differences between simulated and monitored GW levels and SW discharges are due to a combination of different sources of uncertainty: Structural (grid size, simplification of geology), input data (climate data), and parameter uncertainty (e.g., saturated hydraulic conductivity); see also Sections 3.2 and 3.3.

Key Problems Associated with Coupled Hydrologic and Hydraulic Modelling
The advantages of a physical distributed model go along with requiring extensive hydrological input data, high computing capabilities, and long computational times and thus increases the overall effort to successfully setup and calibrate the model. Input grid discretization is a priori in MIKE SHE and has to be selected wisely regarding available data, required accuracy, and computational effort. Like other hydrological models, such as SWAT, MIKE SHE can also simulate catchments size up to thousands of km 2 , but for physical distributed models this increases the spatially distributed input data even more than models based on hydrological response units (HRU) do. Like semi distributed models, MIKE SHE cannot represent sub-grid heterogeneity. In case of limited data availability, satellite data can be very helpful for topography and land use estimation, but additional terrestrial data acquisition (e.g., soil, climatic, aquifer data, etc.) needs to be gathered. Lack of detailed river cross sections, variation of seasonal Manning's n, drainage maps, and leakage coefficient rates impact the simulated river flows. Limited availability of GW boreholes and monitored GW levels to produce GW contour maps makes it harder to delineate the GW divide, and it is sometimes possible to have cross-boundary GW flows. Calibration and validation efforts increase enormously in physically-distributed models with limited observed data, risking the compensation of structural model errors with incorrect model parametrization. Sufficient amount of observed monitoring data will help to provide ease in model setup and during the calibration process.

Transfer of Methodology to Other Lowland Catchments
The coupled model MIKE SHE and MIKE 11 model has demonstrated its potential to simulate hydrological processes common within lowlands. Extensive data requirements are potential problems to apply coupled physically-distributed models in other lowland catchments. Some of the data used in this study is freely available in Germany, such as some of the geo-data and climate data provided by the DWD (German weather service) platform, and can be used for other lowlands in Germany. Manning's n values can be obtained from literature and can be used in other similar catchments. Soil properties were obtained from a local environmental protection agency; literature values or field investigations are required for sites with different soil properties. In Europe, a high resolution DEM model can be obtained from local environmental offices. Lack of detailed river cross sections can be compensated with cross sections based on DEM. Hydraulic structure dimensions can be roughly estimated with Google earth in case of missing information. Calibration of SW discharges and GW flows require monitoring data and there is no other authentic alternative rather than field investigations in the study area.

Key Contributions
An integrated hydrological model coupled with a hydrodynamic model was developed with intent to simulate the moderate climate lowland hydrology. In order to represent the surface and subsurface hydrology at a large scale, simplifications and assumptions were made in order to represent the UZ and SZ. Findings support the hypothesis that the hydrological processes in lowlands are dominated by SW and GW interactions. The key contribution of this study includes: • Development and calibration of a physically-based distributed coupled model to describe the lowlands hydrology, as integrated hydrological-hydraulic modelling has rarely been done in the north-eastern region of Germany; • The ability of physically-based coupled models to describe the lowland hydrology is demonstrated; • Simultaneous calibration of SW and GW with a physically-based integrated hydrological model proves the robustness and reliability of an integrated model. • SW and GW interaction during dry and wet hydrological years in lowlands can be reliably modelled, which is important with regard to the expected and already ongoing change of climatic conditions; • SW discharges and GW flows can be extrapolated at ungauged sites on a physical basis; • The method of coupled modelling, including the steps of model setup and calibration can be transferred to other sites with comparable characteristics.
Land use management practices and their results on SW and GW dynamics can be quantified. Further addition of a nutrient transport model is intended and will help to study the SW and GW quality under different land use scenarios.