Improved representation of groundwater at a regional scale – coupling of mesocale Hydrologic Model ( mHM ) with OpeneGeoSys ( OGS )

Most of the current large scale hydrological models do not contain a physically-based groundwater flow component. The main difficulties in large-scale groundwater modeling include the efficient representation of unsaturated zone flow, the characterization of dynamic groundwater-surface water interaction and the numerical stability while preserving complex physical processes and high resolution. To address these problems, we propose a highly-scalable coupled hydrologic and groundwater model (mHM#OGS) based on the integration of two open-source modeling codes: the mesoscale hydrologic 5 Model (mHM) and the finite element simulator OpenGeoSys (OGS). mHM#OGS is coupled using a boundary condition-based coupling scheme that dynamically links the surface and subsurface parts. Nested time stepping allows smaller time steps for typically faster surface runoff routing in mHM and larger time steps for slower subsurface flow in OGS. mHM#OGS features the coupling interface which can transfer the groundwater recharge and river baseflow rate between mHM and OpenGeoSys. Verification of the coupled model was conducted using the time-series of observed streamflow and groundwater levels. More10 over, we force the transient model using groundwater recharge in two scenarios: (1) spatially variable recharge based on the mHM simulations, and (2) spatially homogeneous groundwater recharge. The modeling result in first scenario has a slightly higher correlation with groundwater head time-series, which further validates the plausibility of spatial groundwater recharge distribution calculated by mHM in the mesocale. The statistical analysis of model predictions shows a promising prediction ability of the model. The offline coupling method implemented here can reproduce reasonable groundwater head time series 15 while keep a desired level of detail in the subsurface model structure with little surplus in computational cost. Our exemplary calculations show that the coupled model mHM#OGS can be a valuable tool to assess the effects of variability in land surface heterogeneity, meteorological, topographical forces and geological zonation on the groundwater flow dynamics. 1 Geosci. Model Dev. Discuss., https://doi.org/10.5194/gmd-2017-231 Manuscript under review for journal Geosci. Model Dev. Discussion started: 22 September 2017 c © Author(s) 2017. CC BY 4.0 License.


Introduction
The significance of depiction of terrestrial water cycle as an integrated system has been continuously recognized.Historically, hydrologic models and groundwater models are isolated and developed in parallel, with either near-surface water flow or subsurface water flow being considered.Furthermore, these models use simple bucket-type expressions together with several vertical water storage layers to describe near-surface water flow (Refsgaard and Storm, 1995;Wood et al., 1997;Koren et al., 2004;Samaniego et al., 2010).Due to the limitation in computational capability, all traditional hydrologic models simplify water flow processes by ignoring lateral water flow, thus inevitably fall short of explicit characterization of the subsurface groundwater head dynamics (Beven et al., 1984;Liang et al., 1994;Clark et al., 2015).
The implicit groundwater representations in traditional hydrologic models show inadequacy in many aspects.Water table depth has a strong influence on near-surface water processes such as evapotranspiration (Chen and Hu, 2004;Yeh and Eltahir, 2005;Koirala et al., 2014).Moreover, water table fluctuation has been discovered as a factor affecting runoff generation, thus affects the prediction skill of catchment runoff (Liang et al., 2003;Chen and Hu, 2004;Koirala et al., 2014).Typical hydrologic models also show inadequacy in simulating solute transport and retention at the catchment scale.Van Meter et al. (2017) found that current nitrogen fluxes in rivers can be dominated by groundwater legacies.Moreover, stream-subsurface water interactions may be significant in modulating the human and environmental effects of nitrogen pollution (Azizian et al., 2017).To assess the response of groundwater to climate change, a physically based groundwater representation including lateral subsurface flow is urgently needed (Scibek and Allen, 2006;Green et al., 2011;Ferguson et al., 2016) .
On the other hand, numerous groundwater models have been developed in parallel.Groundwater models allow for both steady-state and transient groundwater flow in three dimensions with complex boundaries and a complex representation of sources and sinks.A variety of numerical codes are available such as MODFLOW (Harbaugh et al., 2000), FEFLOW (Diersch, 2013), OpenGeoSys (Kolditz et al., 2012).A constant challenging problem in groundwater modeling is the reasonable characterization of heterogeneous and noncontinuous geological properties (Dagan, 1986;Renard and de Marsily, 1997;Attinger, 2003;Heße et al., 2013;Zech et al., 2015).Groundwater models usually contain a physically-based representation of subsurface physics, but fall short in providing good representation of surface and shallow soil processes.For example, models for predicting groundwater storage change under either climate change (e.g., global warming) or human-induced scenarios (e.g., agricultural pumping) always use a constant or linear expression to represent spatially distributed recharge (Danskin, 1999;Selle et al., 2013).The groundwater numerical models may contain some packages or interfaces to simulate surface water or unsaturated zone processes (Harbaugh et al., 2000;Kalbacher et al., 2012;Delfs et al., 2012).Those packages always need extra data and right characterization of many topographical and geological properties.Parameterization of topographical and geological parameters is a big challenge due to the strong spatial and temporal heterogeneity and lack of data (Moore and Doherty, 2006;Arnold et al., 2009).
With the development of computational capability and increasingly importance in responding to climate change, the coupled models coupling two or more hydrological components together have been attracting more and more attentions.The coupled hydrological model highlights the interactions across the shallow soil column and the deep groundwater aquifer.There exist many reviews of the approaches used in coupling surface water-groundwater processes (Brian A. Ebel, 2010;Fleckenstein et al., 2010;Barthel and Banzhaf, 2015).In recent years, there have been some efforts towards coupling surface hydrological model with detailed groundwater model.Maxwell and Miller (2005) coupled LSM (Common Land Model) with a variably saturated groundwater model ParFlow as an integrated model, and demonstrated the need for improved groundwater representation in near-surface water schemes.Sutanudjaja et al. (2011) coupled a land surface model PCR-GLOBWB with a groundwater modeling code MODFLOW using offline coupling scheme, and tested the coupled model using a case study in Rhine-Meuse basin.De Graaf et al. (2015) extended this approach to the global scale.Another highly highly developed coupled model is GSFLOW, which is based on integration of the USGS Precipitation-Runoff Modeling System (PRMS) and the USGS Groundwater Flow Model (MODFLOW and MODFLOW-NWT).GSFLOW has been successfully applied to many case studies (Markstrom et al., 2008;Huntington and Niswonger, 2012;Hunt et al., 2013).
In this study, we document the development of a coupled surface hydrological and groundwater model (mHM#OGS) with an an overall aim of bridging the gap between catchment hydrology and groundwater hydrology at a regional scale.We choose two highly-scalable, open-source codes with high reputations and wide popularities in their corresponding fields: the mesoscale Hydrologic Model mHM (Samaniego et al., 2010;Kumar et al., 2013;Samaniego et al., 2013) and the THMC simulator OpenGeoSys (Kolditz et al., 2012;Wang et al., 2009;Kalbacher et al., 2012).The coupling is achieved by mechanistically accounting for the spatio-temporal dynamics of mHM generated groundwater recharge and baseflow as boundary conditions to the OGS model as a off-line coupled mode.A nested time-stepping approach is used to account for differences in time-scales of near-surface hydrological and groundwater processes, i.e., smaller time steps for typically faster surface runoff routing in mHM and larger time steps for slower subsurface flows in OGS.We applied the coupled model mHM#OGS in a Central German mesoscale catchment (850 km 2 ), and verify the model functioning using measurements of streamflow and groundwater heads from several wells located across the study area.The, herein, illustrated coupled (surface) hydrological and groundwater model (mHM#OGS) is our first attempt towards development of a large-scale coupled modeling system to analyze the spatio-temporal variability of groundwater flow dynamics at a regional scale.
The paper is structured as follows.In the next section, we describe the model concept, model structure, coupling schematic, study area and model setup used in this study.In section 3, we present the simulation result of mHM#OGS in a catchment in central Germany, where the subsurface properties are well characterized and long-term monitoring of river stage and groundwater level exist.We also assess the effects of different spatial patterns of groundwater recharge to groundwater dynamics.In the last section, conclusion and future work are discussed.

mesoscale Hydrologic Model (mHM)
The mesoscale Hydrologic Model (mHM, www.ufz.de/mhm) is a spatially explicit distributed hydrologic model that uses grid cells as a primary modeling unit, and accounts for the following processes: canopy interception, snow accumulation and melting, soil moisture dynamics, infiltration and surface runoff, evapotranspiration, subsurface storage and discharge generation, A unique feature of mHM is the application of Multiscale Parameter Regionalization (MPR).The MPR method accounts for subgrid variabilities of catchment physical characteristics such as topography, terrain, soil and vegetation (Samaniego et al., 2010;Kumar et al., 2013).The model is flexible for hydrological simulations at various spatial scales due to applying the MPR methodology.Within mHM, three levels are differentiated to better represent the spatial variability of state and inputs variables.
Detailed description of MPR as well as formulations governing hydrological processes could be referred to in Samaniego et al. (2010) and Kumar et al. (2013).

OpenGeoSys (OGS)
OpenGeoSys ( OGS) is an open-source project with the aim of developing robust numerical methods for the simulation of Thermo-Hydro-Mechanical-Chemical (THMC) processes in porous and fractured media.OGS is written by C++ with a focus on the finite element analysis of coupled multi-field problems.Parallel versions of OGS are available based on both MPI and OpenMP concepts (Wang et al., 2009;Kolditz et al., 2012;Wang et al., 2017).OGS has been successfully applied in different fields like water resources management, hydrology, geothermal energy, energy storage, CO 2 storage, and waste deposition (Kolditz et al., 2012;Shao et al., 2013;Gräbe et al., 2013;Wang et al., 2017).In the field of hydrology / hydrogeology, OGS has been applied to various topics such as regional groundwater flow and transport (Sun et al., 2011;Selle et al., 2013), contaminant hydrology (Beyer et al., 2006;Walther et al., 2014), reactive transport (Shao et al., 2009;He et al., 2015), and sea water intrusion (Walther et al., 2012), etc.

Structure of the coupled model
The coupled model mHM#OGS was developed to simulate coupled surface-water and groundwater (SW/GW) flow in one or more catchments by simultaneously calculating flow across the land surface and within subsurface materials.mHM#OGS simulates flow within three hydrological regions.The first region is limited by the upper bound of plant canopy and the lower bound of the soil zone bottom.The second region includes open-channel water, such as streams and lakes.The third region is the saturated groundwater aquifer.mHM is used to simulate the processes in the first and second region, while OGS is used to simulate the hydrological processes in the third region.The model development is guided by the following principles: -Solve the governing equations for surface water and groundwater flow using a sequential boundary condition switching technique.-Calculate model-wide and detailed water balances in both time and space.
-Use nested time stepping method to allow different time steps in surface and subsurface regimes -Allow different cell sizes in the spatial discretization of the grid cell resolution used for mHM and the finite-element solution used for OGS.
-Allow different model boundaries definition using standard specified-head, specified-flow, or head-dependent boundary 5 conditions.
An integrated workflow for coupled modeling is illustrated in Figure 1 c.The entire modeling workflow is separated into three independent parts.The first part is the data preparation and pre-processing part marked by blue background in Figure 1 c.This part includes several codes for data preparation and model pre-processing developed by the mHM and OGS communities 5 Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2017-231Manuscript under review for journal Geosci.Model Dev. Discussion started: 22 September 2017 c Author(s) 2017.CC BY 4.0 License.(Fischer et al., 2015;Kolditz et al., 2016).The second part is model coupling, which is composed of the respective computations of mHM and OGS and the data communication between two codes (see components marked by red background in Figure 1 c).
Technically, we use a self-developed data communication code GIS2FEM to convert data format and exchange information.
The detailed illustration of this part is in the following subsection.The third part is the water budget, which is designed to calculate an overall water balance, as well as component-based water budgets for all storages simulated by mHM#OGS (marked by light green background in Figure 1 c).

Boundary condition-based coupling
The subsurface flow equation is solved in OGS.OGS applies a standard (centered) Galerkin finite element method to discretize PDEs.Here we list the governing equations of groundwater flow in saturated zones used in this study.It can be expressed as: where S s is the specific storage coefficient [1/L], ψ p is the pressure head in porous media [L], t is time[T], q is the Darcy flux (LT -1 ), q s is a specified rate source/sink (T -1 ), q e is the exchange rate with the surface water (T -1 ), K s is the saturated hydraulic on Γ, where n is the outer norm of the boundary surface.
The surface flow over the catchment in either the hillslopes or open channels can be expressed using a kinematic-wave equation ( an approximation of the Saint-Venant equations).The kinematic-wave equation for flow in streams can be expressed as: where t is time [T], v is the averaged velocity vector [LT -1 ], ψ s is the surface water depth [L], q e is the exchange rate with subsurface water [LT -1 ].In mHM, a Muskingum-Cunge method is used to solve Eq. (3) (Te Chow, 1988).
The kinematic-wave approximation assumes the gravitational forces are balanced by frictional forces such that: where S o is the slope of the channel [-], and S f is the friction slope of the channel [-].This assumption is used because the potential areas of application of this model would hardly exhibit abruptly hydrographs with supercritical flows.
Manning's equation (Te Chow, 1988) is used to calculate the depth averaged velocity using discharge: where m is the Manning roughness coefficient [L -1/3 T].As an empirical formula, the Manning's formula has been widely used in surface water flow models.
The source terms q e and q e in Eq. ( 1) and Eq. ( 3) explicitly represent the communication of water between the surface and subsurface flow compartments (Camporese et al., 2010).The surface-to-subsurface flux q e is determined after solving the surface routing equation Eq. ( 3) for the following feed to the subsurface flow equation Eq. ( 1), meanwhile the subsurface-tosurface feed q e is determined by solving the subsurface flow equation for the following feed to the surface water equation.Note that the time step of subsurface flow is always larger than that of surface flow.Note that poor water balance might occur if more than 30-50 surface time steps exist in one subsurface time step (Camporese et al., 2010).Here we use a nested time stepping method to avoid the water balance problem (see the following paragraph).Besides, we use the original linear groundwater storage in mHM to calculate daily q e .
Nested time stepping is adopted in this study in order to calculate fast surface and slow subsurface flow simultaneously.As reported by Cunge (1969), the Muskingum-Cunge used to solve surface flow equations is unconditionally stable in case some perquisites being meet, such as proper grid size and time step size.The subsurface solver in OGS is however, implicit in time and limited by less restrictive precision constraints.We use a nested time stepping method to calculate daily surface processes and monthly subsurface processes in a sequential manner.This strategy automatically fits to any stepping size difference, and avoids water balance error in flux exchange between two modules.

Conversions between volumetric flux ([L 3 /T]), specific flux ([L/T]), and water head ([L]
) are performed by adjusting different time steps or cell sizes.Specifically, the time series of groundwater recharge obtained from mHM were fed to the model interface of mHM#OGS.After reading raster file of groundwater recharge, the interface assign the proper recharge value to the top surface elements of OGS mesh by checking the coordinates of the centroid of each top surface element.For each surface element, if its centroid is within a grid cell of the raster file, the value of this grid cell is assigned to the surface element.After all top surface elements have been processed, the elements that have been assigned with the recharge values are involved in the face integration calculation, whereby the recharge data is converted into nodal source terms (see details in Figure 2).The parameter set used in this study is shown in Table 1.  12 End of simulation -Close all input files and processes.

Study area and model setup
We use a mesoscale catchment upstream of Naegelstedt catchment located in central Germany, with a drainage area of about 845 km 2 to verify our model(see Figure 3).The Naegelstedt catchment comprises the headwaters of the Unstrut river basin.The Unstrut river basin is a sedimentary basin of Unstrut river, a left tributary of the Saale.It was selected in this study because there are many groundwater monitoring wells operated by Thuringian State office for the Environment and Geology (TLUG) and the Collaborative Research Center AquaDiva (Küsel et al., 2016).Morphologically, the terrain elevation within the catchment is in a range of 164 m and 516 m, whereby the higher regions are in the west and south as part of the forested hill chain of the Hainich (see Figure 3).The Naegelstedt catchment is one of the most intensively used agricultural regions in Germany.In terms of water supply, about 70% of the water requirement is satisfied by groundwater (Wechsung, 2005).About 17% of the land in this region is forested area, 78% is covered by crop & grassland and 4% is housing and transport area.The mean annual precipitation in this area is about 660 mm.

Meteorological forcings and morphological properties
We started the modeling by performing the daily simulation of mHM to calculate near-surface and soil zone hydrological processes.Several resolutions ranging from 200 m to 2 km are applied in mHM to account different scale of spatial heterogeneity.Muschelkalk (mm) and low Muschelkalk (mu).According to previous geological survey (Seidel, 2004), even the same sub-unit of Muschelkalk have a diverse hydraulic properties depending on their positions and depths.They are further divided into subunits with higher permeability (see mo1, mm1 and mu1 in Figure 4) and sub-units with lower permeability (see mo2, mm2 and mu2 in Figure 4).The uppermost layer with a depth of 10 m is set as a soil layer, whereby the hydraulic properties are set the same with mHM setting (see soil in Figure 4).A alluvium layer is set along the mainstream and major tributaries representing 10 granite and stream deposits.

Boundary conditions
Based on the steep topography along the watershed divides, groundwater is assumed to be naturally seperated and not able to pass across the boundaries of the watershed.No-flow boundaries were imposed at the outer perimeters surrounding the basin as well as at the lower aquitard.
The stream network were delineated using a digital elevation model (DEM)-based pre-processor, and then clipped based on 5 a threshold correlated with field observations.In general, all streams are regarded as perennial in this study, except for those in the mountainous region where flow is intermittent.Stream network is determined based on long-term average of accumulated routed streamflow(see Figure 5).We cut off the small tributaries by which long-term average of accumulated routed streamflow is below a threshold.In other words, we neglect the intermittent streamflow to the upper stream reaches (see the lower left graph in Figure 5).

Calibration
Calibration of the integrated model were conducted using two different conceptual scenarios in order to test the effect of spatial heterogeneity evaluated by Multiscale parameter regionalization (MPR) method in mHM.For the first conceptual scenario (SC1), spatial variability of physiographic characteristics are characterized by MPR method in mHM.The heterogeneous groundwater recharge distribution is determined within the MPR framework.For the second scenario (SC2), we kept the total amount of groundwater recharge of the catchment the same with SC1 in every time step, but use the homogeneous distributed groundwater recharge.SC2 does not consider the complex spatial variability of groundwater recharge caused by variations in climatic conditions, land use, topography and geological heterogeneity.
The coupled model was calibrated following a two-step procedure.In the first step, mHM was calibrated independent of OGS for the period 1970 -2005 by matching observed runoff at the outlet of the catchment.The first 5 years are used as "spin-up" period to set up initial conditions in near-surface soil zone.The calibration workflow is a consecutive workflow where the parameters which affect the potential evapotranspiration, soil moisture, runoff and shallow subsurface flow were first calibrated until convergence criteria was matched.The calibrated mHM model is also verified by measurements from a single eddy-covariance measurement station in the study area.The calibration goodness was handled by means of calculating the Nash-Sutcliffe efficient (NSE).
In the second step, OGS is run independent of mHM in steady-state using long-term averaged outer forcings.Spatially distributed but long-term average recharge estimated by mHM were fed as the steady-state boundary condition.The long-term average baseflow rate estimated by mHM simulations were also used as boundary condition at stream beds.The groundwater where IQ md 7525 and IQ dt 7525 are the inter-quantile ranges of the time-series of modeling result and observations, respectively.

Spatial-temporal dynamics of recharge and baseflow
Groundwater recharge has an arbitrary behaviour depending on the sporadic, irregular, and complex features of storm rainfall occurrences, geological structure, and morphological features.The temporal and spatial variability of groundwater recharge is estimated by mHM calculation with a period of 30 years from 1975 to 2005.conductivity.The greatest point-wise monthly groundwater recharge varies from 26 mm in early spring, to 51 mm in late spring, to 14 mm in winter.We have also evaluated the plausibility of groundwater recharge simulated by mHM with other reference datasets.On the large scale, the simulated groundwater recharge from mHM agrees quite well with the estimation from the Hydrological Atlas of Germany (please refer to Kumar et al. (2016)).
Figure 7 shows the boxplot and histogram of normalized monthly groundwater income and outcome over the whole catch-5 ment.We do not include the human effects (e.g., pumping, abstractions and irrigation).Therefore, baseflow to streams is considered as the only source of groundwater discharge.The boxplot shows the degree of spread and skewness of the distribution of the monthly groundwater recharge and groundwater discharge.It shows that the long-term mean value of monthly groundwater recharge and discharge are balanced with the value of approximately 8 mm/month.Due to the numerical error, a tiny difference of 2% between groundwater recharge and baseflow is observed in the boxplot.This slight bias is within an 10 acceptable interval.The figure also shows that the spread of groundwater recharge is wider than the baseflow, which demonstrate the buffering effect of groundwater storage.the distribution pattern of monthly baseflow is unimodal.It also indicates that the monthly groundwater recharge distribution has a higher deviation than baseflow distribution.

Model evaluation using discharge & groundwater heads
We use the scenario with distributed groundwater recharge calculated by mHM (SC1) as the default scenario, and the scenario with homogeneous groundwater recharge (SC2) as the reference scenario in this study.All calibrated parameters' values are in SC1 by default.For mHM calibration, the calibration result is good with R cor >0.9 for the monthly discharge simulation at the Naegelstedt station (see Figure 10).Other fluxes like evaportransporation measured at eddy-covariance stations inside this area, also shows quite reasonable correspondence to the modeled estimation (please refer to Heße et al. (2017)).
The steady-state groundwater model calibration result shows that the groundwater model can plausibly reproduce the finite numbers of observed groundwater heads within the catchment.Figure 8 shows the 1-to-1 plot of simulated and observed groundwater heads using different recharge scenario SC1 and SC2, respectively (locations of those wells are shown in Figure 3).It can be observed that the model is capable of reproducing spatially-distributed groundwater heads in a wide range, with low RMSE values of 6.22 m in SC1 and 10.14 m in SC2, respectively.There are certain differences between simulated heads and observed heads.This difference is caused by many possible reasons, such as the limited spatial resolution, uniform meshing, or over-simplified geological zonation.A smaller RMSE in SC1 indicates that mHM is able to capture spatial heterogeneity and produce a more realistic groundwater recharge distribution.The errors in simulated heads (lower two graphs in Figure 8) show that most of simulated head errors are within an interval of ±6 m in SC1, while most of errors are within a range of ±10 m in SC2.Nevertheless, there are still some flawed points where the prediction is biased.Adding more model complexity to improve the match between simulation and observation was avoided due to the large spatial scale, the limited spatial resolution of mesh, and the noise of groundwater head data (i.e., with a various time span from 10 years to 30 years).
In general, the model is capable of capturing the historical trend of groundwater dynamics, even though the mean value of simulation and observation values may differ slightly.Due to the limitation of spatial resolution and homogeneous K in each geological unit, this difference is acceptable.
To compare the model results between two scenarios, we drew the barcharts of R cor and |QRE| at each monitoring wells using two recharge scenarios, respectively (see Figure 12).The mean value and the median value of the basin scale R cor and QRE are also calculated and shown in Figure 12. Figure 12 (a) indicates that the correlation with observations of simulations using SC1 is higher than that using SC2, with the averaged R cor of 0.703 and 0.685, respectively.The standard deviation of R cor in SC1 is 0.109, which is 13% smaller than 0.125 in SC2.Considering the only difference between SC1 and SC2 is the spatial distribution of recharge, the heterogeneous groundwater recharge estimated using mHM can be verified as a better evaluation than the homogeneous spatially distributed recharge.The relative difference of R cor between SC1 and SC2 is moderate, which indicates spatial characterization recharge distribution might has a less important influence to groundwater dynamics in the study area.This phenomenon is highly related to the coarsest resolution describing meteorological forcings (e.g., precipitation) among three spatial levels used in mHM.The coupled model also shows its potential in predicting groundwater flood and drought.Figure 13 displays the seasonality of groundwater heads over the whole catchment by means of calculating the long-term mean groundwater heads in spring, summer, autumn and winter, respectively.It indicates that in general, the possibility of groundwater flood in spring & summer is higher than in autumn & winter.However, the spatial variability within groundwater flood & drought event is significant.The groudwater head in northern, eastern and southeastern mountainous areas tend to fluctuate more wildly than the central plain areas.This phenomenon is consistent with the fact that the mountainous areas have a larger recharge rate than the plain areas.
Considering the need of predicting groundwater flood & drought in extreme climate events, we select a wet month (August 2002) and a dry month (August 2003), and show the groundwater heads variation in these months.Figure 13 e) and f) show two scenes of groundwater head variation in wet season and dry season, respectively.The groundwater heads in wet season is higher than the mean values (see Figure 13 e ).The variation of groundwater heads in dry season, however, shows a strong spatial variability.The strong spatial variability of groundwater heads variation has also been found in Kumar et al. (2016).

Discussion and conclusion
A coupled hydrologic model mHM#OGS is proposed and applied it in a mesoscale catchment in central Germany.A boundary condition-based off-line coupling method is applied to depict the dynamic flow exchanges between the surface and subsurface water regimes.This coupling method, together with nested time stepping, allow the surface and subsurface parts to be solved sequentially, keep computational surplus to a minimum, while avoid possible water balance problem in the flow exchange processes.The result shows a promising prediction capability in surface and subsurface water modeling via calibration and comparison to groundwater time series.The SC1 using spatially heterogeneous groundwater recharge distribution is more plausible than SC2 which uses spatially homogeneous groundwater recharge.The results of this study highlight the successful application of Multiscale Parameter Regionalization (MPR) method in characterizing spatial heterogeneous groundwater recharge.In the spatial scale of 10 3 km 2 (the scale in this study), the MPR method shows a moderate improvement in groundwater recharge representation.Note that MPR has been successfully applied at a larger scale over Europe (Thober et al., 2015;Kumar et al., 2013;Zink et al., 2016;Rakovec et al., 2016).The effectiveness of MPR to characterize groundwater dynamics at larger scales (e.g., 10 4 -10 6 km 2 ) even global scales is still unknown and needs to be explored.Moreover, MPR has been proved its capability to produce better runoff prediction in cross-validated locations (e.g., ungauged basins) (Samaniego et al., 2010).To date, the effectiveness of MPR in unguaged basins has not been verified by groundwater head dynamics.In the next step, we may use the groundwater time series to test the effectiveness of MPR in ungauged basins using the coupled model mHM#OGS .
The convincing results of this study provide a new possibility in improving classic large-scale disbributed hydrologic models, such as the current unmodified version of mHM (Samaniego et al., 2010;Kumar et al., 2013), VIC (Liang et al., 1994), PCR- coupling processes in large-scale hydrologic cycles, which is significant for a wide range of real-world applications, including land subsidence, agricultural irrigation, nutrient circulation, salt water intrusion, drought, and heavy metal transport.
We realize that there are several limitations in the current model.The first one is the fact that we use an off-line coupling scheme instead of a full coupling scheme between mHM and OGS.Consequently, we do not account for the explicit feedback from OGS to mHM, although it might be less important related to the large subsurface time step.This approach has certain advantages of less computational consumption and better numerical stability, and fits perfectly on the long-term large-scale groundwater modeling in this study.In the next step, we may try to incorporate the full coupling scheme in the next version of mHM#OGS model.Via the full coupling scheme, the dynamic interactions between overland flow and groundwater flow, and between soil moisture dynamics and groundwater dynamics are explicitly accounted.This approach is open to a broader spectrum of calibration options, such as calibration using remotely sensed soil moisture data.The second one is that we do not use the parallel computing.Although the whole simulation is conducted on the EVE linux cluster at UFZ, which is a high performance computing platform, we do not use the distributed computing to reduce computational efforts.In the future, a parallel version of mHM#OGS is needed to reduce computation time for the computationally-expensive full coupling procedure.

3
scheme which routes runoff in upstream cells along river network using the Muskingum-Cunge algorithm.The model is driven by daily meteorological forcings (e.g., precipitation, temperature), and it utilizes observable basin physical properties or signals (e.g., soil textural, vegetation, and geological properties) to infer the spatial variability of the required parameters.mHM is an open-source project written in FORTRAN 2008.Parallel versions of mHM are available based on OpenMP concepts.

Figure 1 .
Figure 1.The concept of mHM#OGS model a) the conceptual representation of hydrological processes in a catchment; b) the schematic used to couple the mHM and OGS.The upper box depicts the canopy interception, atmospheric forcing, and the land surface processes represented by mHM.The lower box depicts the saturated zone represented by the OGS groundwater model; c) the complete workflow including several interfaces with external softwares for data import, format conversion, model calibration and water balance check.

Figure 2 .
Figure 2. Transfer of groundwater recharge from mHM grid cells to OGS nodes using the model interface GIS2FEM.

4
Transfer surface-to-subsurface exchange rates to OGS -Transfer surface-to-subsurface volumetric flow rates needed for computing saturated flow as Neumann boundary conditions in OGS. 5 Steady-state calibration -Run OGS-only steady-state simulation using boundary conditions given by mHM in step 4. The calibrated K field is fed to transient model.The steady-state groundwater head serves as an initial condition of the transient mHM#OGS modeling.6 Start transient simulation of mHM#OGS -Sequence through coupled mHM and OGS components.7 Compute stepwise near-surface flow and storage in mHM -Compute spatially-distributed daily near-surface processes.8 Transfer primary variables and volumetric flow rates to OGS -The same as step 4, except this step generate time-dependent raster files of flow rates.8 Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2017-231Manuscript under review for journal Geosci.Model Dev. Discussion started: 22 September 2017 c Author(s) 2017.CC BY 4.0 License.Solve the groundwater flow equation -Calculate groundwater heads and groundwater flow velocity field in study region.10 Compute budgets -Run water budget package to check overall water balance as well as timedependent water budgets in each storages.11 Write results -Output the simulations results.
Heße et al. (2017) have already established the mHM simulation over the study area.The meteorological and morphological settings in this paper are the same with his work.For the detailed settings of meteorological forcings and morphological properties, please refer toHeße et al. (2017).2.4.2 Aquifer properties and meshingTypical distributed hydrological models use shallow soil profiles or extended soil profile to present groundwater storage.Here, we use a spatially distributed aquifer model to explicitly present groundwater storage.We set up this spatially distributed aquifer model through geological modeling based on well log data and geophysical data from Thuringian State office for the Environment and Geology (TLUG).To convert data format, we use the workflow developed byFischer et al. (2015) to convert the complex 3D geological model into open-source VTU format file that can be used by OGS.Model elements of OGS were set to a 250 m × 250 m horizontal resolution and a 10 m vertical resolution over the whole model domain.

Figure 3 .
Figure 3.The Naegelstedt catchment used as the test catchment for this model.The left map shows elevation and locations of monitoring wells used in this study.The lower right map shows the relative location of Naegelstedt catchment in Unstrut basin.The upper right map shows the location of Unstrut basin in Germany.

Figure 4 .
Figure 4. Three-dimensional and cross Section view of hydrogeologic zonation in the Naegelstedt catchment.The upper left figure shows the complete geological characterization and zonation including alluvium and soil zone.The upper right figure shows the geological characterization along two cross sections.The lower map shows the detailed zonation of geological sub-units beneath the soil zone and alluvium.

Figure 5 .
Figure 5. Illustration of stream network used in this study.a) Stream network based on long-term average of accumulated routed streamflow; b) Stream network whereby long-term averaged accumulated monthly streamflow rate is above 1000 mm; c) Stream network whereby longterm averaged accumulated monthly streamflow rate is above 1500 mm, which is also the default setting in this study; d) Stream network whereby long-term averaged accumulated monthly streamflow rate is above 2000 mm.

Figure 6 .
Figure 6.Spatial distributions of groundwater recharge in Naegelstedt catchment (unit: mm/month) (a) during early spring, (b) late spring, and (c) winter of year 2005.
Figure 6 shows the spatial variability of groundwater recharge in three months: the early spring (March) (Figure 6 a), late spring (May) (Figure 6 b), and winter (January) (Figure 6 c).The results indicate that the largest groundwater recharge may occur at mountainous areas.Greatest recharge occurs in the upstream bedrock areas where dominant sedimentary is Muschelkalk with a relatively low hydraulic

Figure 7 .
Figure 7. Mean monthly water balance of groundwater over the Naegelstedt catchment.a) Boxplot indicates spread, skewness, and outliers of groundwater recharge and groundwater discharge.b) Histogram indicates the distribution of groundwater balance.c) Monthly time series of groundwater recharge and baseflow.
Figure 7 b shows the distribution of monthly groundwater recharge and monthly baseflow.The figure indicates that the distribution pattern of monthly groundwater recharge is skewed right, whereas

Figure 8 .
Figure 8. Illustration of steady-state calibration results.(a) Observed and simulated groundwater head, including RMSE and Rcor; (b) Difference between simulated and observed head related to the observaed head values.

Figure 10 .
Figure 10.Observed and modeled monthly streamflow at the outlet of Naegelstedt catchment.

Figure 11
Figure 11 presents observed and simulated groundwater heads for the period 1975-2005 in SC1.Five out of 19 monitoring wells with different geological and morphological types were chosen as samples to test the effectiveness of our model.Well 4728230786 is located at northern upland and near the mainstream, whereas well 4828230754 is located at the southwestern lowland.Both of those two wells show promising simulation results with the R cor of 0.81247 and 0.75279, and the QRE of -10

Figure 11 .
Figure 11.The comparison between measurement data (green dashed line) and model output of groundwater head anomaly in SC1 (blue solid line).(a) Monitoring well 4728230786 located at upland near stream.(b) Monitoring well 4628230773 located at mountainous area.(c) Monitoring well 4728230781 located at a hillslope at northern upland.(d) Monitoring well 4828230754 located at lowland.(e) Monitoring well 4728230783 located at northern mountain.

Figure 12 .
Figure 12.Barplots of a) the Pearson correlation coefficient Rcor and b) absolute inter-quantile range error |QRE| in all monitoring wells in two scenarios.Each bar corresponds to an individual monitoring well in the following order: 0 -4830230779, 1 -4828230754, 2 -

Figure 12 b
Figure12b shows the distribution of the absolute value of inter-quantile range error |QRE| in SC1 and SC2.It can be found in Figure12that the distribution pattern of |QRE| is more complicated than R cor .We can see that the |QRE| in two wells are abnormally higher than the other wells.This indicates the accurate quantification of amplitude in particular locations

Figure 13 .
Figure 13.Seasonal variation of spatially-distributed groundwater heads by their anomalies after removing the long-term mean groundwater heads (unit: m).a) Long-term mean groundwater head distribution in spring; b) Long-term mean groundwater head distribution in summer; c) Long-term mean groundwater head distribution in autumn; d) Long-term mean groundwater head distribution in winter; e) Monthly mean groundwater head distribution in wet season (August 2002); f) Monthly mean groundwater head distribution in dry season (August 2003).
Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2017-231Manuscript under review for journal Geosci.Model Dev. Discussion started: 22 September 2017 c Author(s) 2017.CC BY 4.0 License.GLOBWB (Van Beek and Bierkens, 2009), WASMOD-M (Widén-Nilsson et al., 2007) .Those distributed hydrologic models do not include the function of calculating spatial-temporal groundwater heads, therefore are not able to reasonably represent groundwater head and storage dynamics in their groundwater regime.This may be insignificant in global scale hydrologic modeling, which always has a coarser resolution of 25-50 km.The physical representation of groundwater flow is needed in future global hydrologic model with finer spatial resolutions down to 1 km.Moreover, the inclusion of groundwater model OGS in the coupled model is particularly significant for areas with large sedimentary basins or deltas (e.g., the sedimentary basins of Mekong, Danube, Yangtze, Amazon, and Ganges-Brahmaputra Rivers).The coupled model mHM#OGS also provides a potential in predicting groundwater flood & drought in extreme climate events.Due to the prediction capability of mHM in ungauged basins, the coupled model is also capable of predicting groundwater flood & drought in ungauged basins, which is quite valuable due to the lack of comprehensive groundwater observations at regional scales.Providing previous work byHeße et al. (2017) in Travel Time Distributions (TTDs) using mHM, we can expand the range of their work to the complete hydrologic cycle beneath atmosphere, which is important due to the pollutant legacy in groundwater storage.The coupled model is also able to evaluate surface water and groundwater storage change under different meteorological forcings, which allows the decent study of hydrologic response to climate change (e.g.global warming).Besides, the versatility of OGS also offers the possibility to address Thermo-Hydro-Mechanical-Chemical (THMC) Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2017-231Manuscript under review for journal Geosci.Model Dev. Discussion started: 22 September 2017 c Author(s) 2017.CC BY 4.0 License.

Table 1 .
Table of hydrological parameters used in this study.The two models are coupled in a sequential manner by fed fluxes and variables from one model to another at every subsurface time step.Technically, the coupling interface converts time-series of variables and fluxes to Neumann boundary conditions, which can be directly read by OGS.The modified OGS source code can produce raster files containing the time-series of flow-dependent variables and volumetric flow rates with the same resolution of mHM grid cells which can be directly read by mHM.The detailed workflow of the coupling technique is shown in Table2.

Table 2 :
Description of computational sequence for mHM#OGS using a sequential coupling scheme Computation 1Initialize, assign, and read -Run mHM and OGS initialize procedures, OGS assign, read and prepare parameters and subroutines for later simulation.2Computenear-surfacehydrologicprocesses in mHM -Calculate near-surface processes such as snow melting, evapotranspiration, fast interflow, slow interflow, groundwater recharge, and surface runoff for each grid cell.3Compute the long-term mean of land-surface and soil-zone hydrologic processes -Compute the long-term mean of the entire simulation period and write them as a set of raster files.

Table 3 .
Estimates of hydraulic properties for calibrated steady-state groundwater model in Naegelstedt catchment.