Lateral terrestrial water fluxes in the LSM of WRF-Hydro: Benefits of a 2D groundwater representation

The interactions between the atmosphere and the land surface are characterized by complex, non-linear processes on varying time scales. The Noah-MP is a medium complexity land-surface model (LSM), which was recently selected as the new default LSM for the hydrologically enhanced Weather Research and Forecasting modelling system (WRF-Hydro). Compared to its predecessor, several parameterizations were considerably improved and new ones added, inter alia more sophisticated groundwater descriptions, which aim to replace the traditional free-drainage lower boundary condition. This study investigates the benefits that can be obtained from a two-dimensional groundwater representation within the WRF-Hydro modelling system by performing two offline simulations for the upper Danube river basin. In comparison to the free-drainage reference simulation, the lateral routing of groundwater and the two-way interaction with the water table greatly enhances small scale variability in simulated fields of soil moisture content and evapotranspiration (ET). The representation of upward fluxes from the aquifer helps to maintain higher soil moisture contents and thus ET during prolonged dry periods. These differences are rather small though (<2%) and explained by the fact that the study region is considered to be limited by radiative energy and not water availability. The most striking difference however is the performance gap in simulating streamflow. WRF-Hydro with 2d groundwater scheme clearly outperforms the reference simulation in terms of performance metrics. A comparison with hourly streamflow observations for the water year of 2016 yields average Kling-Gupta efficiencies of 0.79 versus 0.57 for the reference. Given that both model configurations were not calibrated beforehand, we conclude that the two-dimensional groundwater option is especially beneficial for applications in poorly or even ungauged catchments. Furthermore, the inclusion of a so far missing compartment of the water cycle in the WRF-Hydro modelling system allows for a more holistic representation of interactions between atmosphere land surface and subsurface, which will be advanta-geous in feedback studies with the fully coupled WRF-Hydro.


| INTRODUCTION
Besides representations for atmosphere, oceans, and cryosphere, land surface schemes are an integral part in models describing the earth system on a global or regional scale. More precisely, they are providing the lower boundary to an atmospheric model by representing a variety of hydrological, biogeophysical, and biogeochemical processes and thus information on fluxes of heat, moisture, and momentum to the atmospheric model. Taking this background into account, many of the land surface models (LSMs) were originally developed with the coupled application in focus. Consequently, these models were aligned to the spatial resolutions commonly used by atmospheric models and are limited to a one-dimensional approach in describing vertical fluxes across the compartmental interfaces. Another characteristic-contrarily to classical hydrological models-is the fact that process descriptions with relevance for the fluxes at the upper boundary are typically more detailed than those of the lower boundary, that is, the bottom of the simulated soil column in a depth of 2-10 m (Wood et al., 2011).

One aspect of past and present endeavours to further evolve
LSMs is the improvement of hydrological process descriptions within these model systems (e.g. Clark et al., 2015;Ning et al., 2019). The underlying motivation is mainly to achieve a better understanding of the interactions between the different compartments of the hydrological cycle, and thus to ultimately improve the results obtained from coupled applications, like weather and climate predictions. These interactions are currently far from being understood, as shown by the still growing number of studies investigating feedback mechanisms between the land-surface and the atmosphere (e.g. Arnault et al., 2016;Arnault et al., 2018;Arnault et al., 2019;Barlage et al., 2015;Davison et al., 2015;Fersch, Francke, et al., 2020;Forrester & Maxwell, 2020;Keune et al., 2016;Larsen et al., 2016;Rahman et al., 2015;Rummler et al., 2018;Senatore et al., 2015;Zhang et al., 2019).
In this regard, various efforts have been made to replace the crude, unidirectional, free-drainage boundary condition commonly found at the lower boundary of a typical LSM. More realistic descriptions were added that allowed for two-way interactions between the unsaturated zone and the water table. This led to a range of groundwater implementations with varying degrees of complexity and specific use cases in mind. These include one-dimensional, column-based approaches (e.g. Liang et al., 2003;Niu et al., 2005;Niu et al., 2007;Yeh & Eltahir, 2005) that are computationally very efficient and hence can be easily applied on a global scale. Medium complexity schemes additionally consider lateral flow between neighbouring grid-cells (e.g. Batelis et al., 2020;Miguez-Macho et al., 2007;Vergnes & Decharme, 2012) and have shown their applicability on continental and global scales in various studies using offline and fully coupled model configurations (e.g. Anyah et al., 2008;Barlage et al., 2015;Decharme et al., 2019;Fan et al., 2007;Martinez et al., 2016;Miguez-Macho & Fan, 2012). Three-dimensional variably saturated groundwater models like for example ParFlow (Ashby & Falgout, 1996) have also been integrated into LSMs (e.g. Kollet & Maxwell, 2006;Maxwell & Miller, 2005) and subsequently coupled to multiple regional climate models (e.g. Maxwell et al., 2007;Maxwell et al., 2011;Shrestha et al., 2014;Sulis et al., 2017). Even though these models are considerably more expensive from a computational point of view, first applications on a continental scale (e.g. Keune et al., 2016;Maxwell et al., 2015) proved the basic feasibility of such endeavours.
However, hydrologically enhanced LSMs are not only used in process studies or coupled applications as outlined above, but also advance into areas that were previously covered by classical hydrological models, like flood forecasting (e.g. Givati et al., 2016). Recently, the WRF-Hydro modelling framework  was selected to become the National Water Model (NWM) to provide streamflow forecasts over the entire continental United States. The coupling framework extends a traditional LSM with a range of modular hydrological physics options and a subgrid approach to account for example topographic heterogeneity. Groundwater processes, on the other hand, are currently represented in a comparatively simplistic way. This, however, shows the necessity to investigate the added value of an explicit two-dimensional groundwater parameterization within the WRF-Hydro modelling system, from which our research questions arise: How does the implementation of an explicit twodimensional groundwater parameterization in the WRF-Hydro modelling system impact simulated variables of (i) soil moisture, (ii) evapotranspiration, and does it improve the model's skill in reproducing observed (iii) river streamflow? To answer these questions, we implement the groundwater scheme described in Miguez-Macho et al. (2007) into the WRF-Hydro modelling system and compare results from a simulation with explicit groundwater to those obtained from a simulation using the free-drainage boundary condition, which is currently the default option in both WRF-Hydro and the NWM.
2 | MODEL DESCRIPTION 2.1 | WRF-Hydro modelling system WRF-Hydro version 5.1.1  is used as the hydrological model for this study. At its core, WRF-Hydro is a coupling framework that enhances a traditional one-dimensional land surface parameterization by providing additional descriptions for a range of hydrological processes, that is, representations for surface overland flow, saturated subsurface flow, channel routing, baseflow processes, and lakes or reservoirs. Since the typical grid scales, LSMs operate on are too coarse to accurately route water across the landscape, WRF-Hydro uses a sub-grid approach to account for topographic heterogeneity. The model can be used as a traditional hydrological model, which is often referred to as "one-way coupling", "offline", or "standalone" mode, where the meteorological forcing data has to be supplied as gridded datasets. Alternatively, the model can be run in a fully coupled configuration, where the LSM constantly exchanges information with a limited area atmospheric model. The latter is driven by meteorological forcing at the lateral boundaries of the model domain and thus no additional hydrometeorological inputs to the LSM are required. This configuration is often referred to as "fully" or "twoway" coupled mode.
The default LSM since WRF-Hydro version 5 is the Noah-MP, which is an augmented Noah-LSM with Multi-Parameterization options Yang et al., 2011). It is largely based on its predecessor (Noah-LSM), although changes to the code structure required a complete rewrite of the model code, which now follows a modular concept. Compared to the legacy Noah LSM, the range of user-selectable process descriptions has been considerably extended, like for example by the Ball-Berry canopy stomatal resistance scheme (Ball et al., 1987), a representation of dynamic vegetation (Dickinson et al., 1998), a multi-layer snowpack parameterization (Yang & Niu, 2003) with a more permeable frozen soil scheme (Niu & Yang, 2006), and several new runoff/groundwater options (e.g. Miguez-Macho et al., 2007;Niu et al., 2007). The soil column in the Noah-MP LSM has a fixed 2-m profile that is subdivided into four layers with thicknesses of 10, 30, 60, and 100 cm, respectively ( Figure 1). Soil related model parameters are a function of textural class and read in from a lookup table.

| Free drainage runoff parameterization
The default runoff option in the WRF-Hydro model system is Noah-MP's free drainage condition (from now on referred to as FD), where simulated runoff is primarily governed by the three global parameters F I G U R E 1 A schematic illustration of a single soil column from the Noah-MP LSM with MMF groundwater parameterization (left) and the associated WRF-hydro routing extensions (right). Adapted from Rummler et al. (2018) REFKDT, REFDK, and SLOPE (e.g. Yucel, Onen, Yilmaz, & Gochis, 2015;Cuntz et al., 2016). The first two parameters REFKDT (.) and REFDK (.) are infiltration scaling parameters that are used to determine the surface runoff Q surface : where Q wat (m s À1 ) is the water input rate on the soil surface, and I max (m s À1 ) the soil infiltration rate, which is calculated as following: where D (m) is the liquid soil moisture deficit of the modelled soil column, δ t the time between two model time steps, Ksat (m s À1 ) is the saturated hydraulic conductivity, Kref (m s À1 ) the reference value for the saturated hydraulic conductivity (REFDK), which defaults to the value of the silty-clay-loam soil texture (2Â10À6 m s À1 ) and kdtref the REFKDT parameter. Finally, the SLOPE parameter is a scaling factor moderating the percolation Qperc (mm): where K 4 (mm) is the drainage from the bottom soil layer, and Slope (.) is a scaling factor. By default, the Noah-MP LSM assumes a global Slope value of 0.1. When using WRF-Hydro's conceptual groundwater bucket option, Q perc is collected in those buckets and subsequently redistributed to the channel network. Groundwater buckets are userdefinable and typically setup to cover the catchment area of individual river reaches or measurement locations.

| MMF runoff parameterization
Another, more complex, runoff option within the Noah-MP LSM is the two-dimensional MMF groundwater parameterization (Barlage et al., 2015;Fan et al., 2007;Miguez-Macho et al., 2007) with Recharge (m) being the flux to/from the 4th soil layer of the Noah-MP, Q lateral (m) being the sum of positive/negative fluxes to / from a grid cell and Q river (m) representing the groundwater flow to rivers: with R cond being the river conductivity (.), h the water table head (m), and z riverbed (m) the elevation of the assumed riverbed within a model grid cell. It is important to note that Q river is one-way only, so that fluxes from the riverbed to the groundwater are neglected.
The calculation of surface runoff differs from the one used by the default FD option (Equation 1). It is based on a TOPMODEL approach (Niu et al., 2005;Niu et al., 2011), describing the saturation excess of a fractionally saturated model grid cell, which is calculated by: with FCR 1 being the impermeable pore space fraction due to frozen soil within the 1st soil layer. The saturated surface fraction F sat is an exponential function of the depth to the water table and is calculated by: with f (m À1 ) being a runoff decay factor with a hardcoded value of 6 m À1 , z wt (m) the depth to the water table, and F max (.) a tunable global runoff parameter, which represents the maximum saturated surface fraction of a model grid cell and defaults to a value of 0.38 within the Noah-MP. Thus, F sat becomes F max when z wt is within the resolved layers of the LSM, that is, z wt ≤ 2 m. Finally, the change in depth of the water table can be computed by: with Θ sat (m 3 m À3 ) and Θ equ (m 3 m À3 ) being the saturated and equilibrium soil moisture content, respectively. In case that ΔS is greater than the available pore space of the auxiliary layer, the water table can rise to soil layers of the LSM. If these layers get saturated before ΔS is fully depleted, the excess amount is treated as spring discharge Q spring .
It is noted that values for Θ equ , R cond , and z riverbed are computed during the initialization of the MMF parameterization via an internal spin-up routine. This routine requires an additional dataset on climatological groundwater recharge, which is provided by the WRF

| Model adjustments
The current implementation of the MMF groundwater option within the Noah-MP LSM relies on the parallelization infrastructure provided by the Weather Research and Forecast (WRF) atmospheric model (Skamarock et al., 2019). To make use of the parameterization in the offline variant of WRF-Hydro, the relevant code passages had to be adapted accordingly. Furthermore, the MMF scheme was not originally designed to accommodate for WRF-Hydro's routing extensions and vice versa. Hence, we made arrangements to incorporate Q river and Q spring within WRF-Hydro's coupling/routing concept. In particular, Q river is now redistributed to individual river reaches via the associated groundwater buckets and Q spring is added to the infiltration excess to make it available for redistribution with WRF-Hydro's surface routing scheme. A schematic overview of the MMF parameterization within the Noah-MP and WRF-Hydro's routing extensions is given Figure 1.

| Study region and period/observations
We  The second quarter of the water year started with substantial rainfall amounts in January and February that were well above the long-term average. In May and June, a series of strong convective precipitation events lead to localized flooding. In terms of temperature, the winter months were considerably warmer than the long-term average. Spring, on the other hand, was close to average and summer only slightly warmer.

| Experiment design and model setup
To assess the added value of a more detailed groundwater description within the WRF-Hydro modelling system, two experiments with the two runoff parameterizations outlined in Section 2.1 (FD) and Furthermore, it is important to note that the model calibration was limited to WRF-Hydro's channel network.
To be more precise, only the settings for channel geometry and Manning's roughness coefficients were adjusted, based on results obtained from previous simulations (not shown). We opted for this approach to avoid issues introduced by overfitting model parameters, which would considerably limit the robustness of the results and thus statements on their transferability. We therefore neglect parameter and data uncertainty in this study and instead solely focus on the structural model uncertainty.
To assess the performance of the two runoff options, we employ a regional approach and compare simulated hourly streamflow to quality-controlled observations obtained from the state measurement This circumstance was accounted for in the selection process of the four study catchments, so that simulated runoff can be more easily compared with the observations. The hourly gridded meteorological forcing data used to drive both experiments is described in the following section.

| Meteorological forcing dataset
The regional reanalysis Cosmo Rea6 (Bollmeyer et al., 2014) was chosen to provide the meteorological forcing to both offline WRF-Hydro model runs. The regional reanalysis is derived from ECMWF's ERA-  Note: Values for long-term annual mean temperature and precipitation were derived from gridded observational datasets provided by the German meteorological service.
The evolution of basin-averaged, column integrated soil saturation fraction for the four study catchments is shown in Figure 4.
Across the board, the largest differences are found during the first months of the study period, where the MMF experiment maintains a higher soil saturation. This is thought to be related to a preceding drought that affected southern Germany during summer of 2015 and the subsequent months with precipitation amounts well below average. During this phase, the representation of groundwater-with the possibility to transport water upwards from the aquifer-is capable of maintaining higher soil moisture contents. The following months, however, show little differences, which is most likely related to the humid conditions with ample precipitation amounts and the fact that the FD experiment presented here uses the Noah-MP with the default scaling factor of 0.1 to constrain the percolation out of the bottommost soil layer (Equation 5). Additionally, both experiments employ WRF-Hydro's surface and subsurface routing modules, which allow for reinfiltration of surface runoff and therefore help to retain moisture within the model domain (e.g. Senatore et al., 2015;Yucel et al., 2015). This behaviour changes close to the end of the study period in August 2016, when the MMF experiment is again able to maintain wetter soils. An exception to this behaviour is the catchment of the river Woernitz, which shows lower values for most of the study period. This might be related to groundwater recharge and is further discussed in Section 4.3.

| Ramifications of the improved groundwater representation on WRF-Hydro's streamflow prediction skill
The timeseries of simulated and observed hourly streamflow for the four study catchments are presented in Figure 7 and associated performance metrics in Table 2. Overall, the metrics obtained for the MMF experiment clearly outperform those of the FD experiment in all F I G U R E 3 (left and middle panel) Average column integrated soil moisture content (mm) for the water year 2016 and the simulations with two-dimensional groundwater parameterization (MMF) and freedrainage (FD) condition, respectively. (right panel) As in the left, except for the relative differences (%) between the MMF and FD experiment. The four study catchments are marked by black solid outlines and the basin of the upper Danube by a black dotted outline four study catchments and presented measures, with the exception of PBIAS in the Iller catchment. Most striking, however, is the fact that the performance difference between the FD and MMF simulation is considerably larger in the northern than in the southern study catchments.
Indeed, the southern study catchments receive considerable higher precipitation amounts (Table 1) and therefore the soil columns in the FD experiment are also more likely to reach the higher saturation levels found with the MMF experiment (Figure 4a,b), which in turn might lead to more similar surface runoff amounts in response to a precipitation event. Additionally, the complex topography with steep slope angles ( Figure 1 and Table 1)  In contrast, the situation in the northern study catchments is quite different. Both catchments receive substantially less precipitation, and the topography is less pronounced (s. Table 1 and Figure 1).
As shown by the unsatisfactory NSE results with values of 0.12 and À 0.58, the FD experiment has serious difficulties with predicting the right timing and magnitude of peak flow events in the catchments of the rivers Wörnitz and Naab. In particular, a couple of events during the first three month of the study period are not captured at all

| DISCUSSION
The evaluation of simulated hourly streamflow at four different locations in the study region revealed that the simulation with groundwater representation was far more capable in capturing the dynamics,  Fersch, Senatore, et al., 2020;Kerandi et al., 2017;Lahmers et al., 2019;Yucel et al., 2015), Compared to the baseline WRF-Hydro configuration, the simulation with MMF groundwater scheme increases the computational time by less than 5%, which would not add significantly to computational cost in long-term simulations. This is achieved by allowing for a series of simplifying assumptions in terms of model physics, as well as the structure and properties of the subsurface. Considering that the complexity of individual model components should be balanced (Prentice et al., 2015), we argue that the presented groundwater implementation is nevertheless a valuable addition to the WRF hydromodelling system and usefully augments the already existing parameterizations for surface and subsurface routing. For more sophisticated application scenarios that for example demand for a variably saturated 3D representation of subsurface flow or a representation of human interaction through extraction of groundwater, one of the models presented in example Kollet et al. (2017) might be a more appropriate choice. Although at the cost of a considerable higher computational burden and additional data requirements (Rahman et al., 2019;Tijerina et al., 2021).
Results further demonstrated that a representation of lateral groundwater fluxes and two-way interactions between water table and soil column lead to a considerably higher spatial heterogeneity in terms of simulated soil moisture content and, to a somewhat lesser extent, evapotranspiration. Overall, a strong link to model topography was evident, which is in line with the findings of Forrester and Maxwell (2020) and Batelis et al. (2020). While a comparison on a point basis yielded considerable differences on the order of À30% to 50% for these quantities, the spatial averages for the study region (< 2%) and domain (< 1%) were very close. This can be attributed to the fact that the region is not water limited, but energy limited which is supported by analysis of long-term averaged climate datasets (Klingler et al., 2021).  (Robinson et al., 2008). Considerable differences were identified along ridges and valley structures, while only subtle differences were obtained for catchment averages. This underlines the need for observational datasets that at least (a) cover an extent on the order of tens of kilometres and (b) provide the information at suitable spatial scales.
In this context, a promising option might be the cosmic ray neuron sensing technique proposed by Zreda et al. (2008). The non-invasive measurement approach alleviates the representativeness issue of point measurements by having a measurement footprint on the order of hundreds of meters, which suits the resolutions typically used by LSMs. Measurement depths are superior to airborne or satellite-based remote sensing products and are on the order of tens of centimetres.
Consequently, the method may have potential to also capture soil moisture content in the root zone, which is of crucial importance for plant transpiration. Additionally, mobile measurement equipment (Desilets et al., 2010) can be used to retrieve measurements with greater spatial detail for specific areas of interest, for example, valley transects where our simulations showed considerable differences.
The measurement network presented in Fersch et al. (2020b), for example, could be a suitable starting point if a different strategy is used for the placement of sensors.
Finally, the results obtained for soil moisture bias and heterogeneity might also be relevant for further feedback studies with the fully coupled WRF-Hydro modelling system. As outlined by Wood et al. (2011), many warm-season precipitation events are thought to be potentially sensitive towards small-scale variabilities in land-surface properties, like for example topography, land-use and soil moisture. In particular, feedbacks are expected during the warm season and phases of weak synoptic forcing that allow for the development of mesoscale circulation (Taylor, 2015;Taylor et al., 2007). In terms of soil moisture, two mechanisms are conceivable: A positive feedback loop via the concept of precipitation recycling (e.g. Eltahir & Bras, 1996), which relates the contribution of local evapotranspiration to local precipitation. Alternatively, a negative feedback loop due to the location of convective initiation (Taylor, 2015;Baur, Keil, & Kraig, 2018). Differential heating of dry and moist patches is thought to create circulation cells within the planetary boundary layer, resulting in convergence over dry and divergence over wet patches. An increase of precipitation is thus obtained for the dry patches, since convergence supports an earlier transition from shallow to deep convection at those locations.
The findings of this study indicate that simulations with the MMF scheme may facilitate or strengthen positive feedback loops during prolonged dry phases, when the explicit groundwater representation can provide additional moisture to sustain higher evapotranspiration rates. In addition, the higher topography-induced heterogeneity of soil moisture simulated by the MMF scheme may enable or amplify negative feedback loops under favourable conditions.

| CONCLUSIONS
In this study, we replaced WRF-Hydro's empirical groundwater bucket model with a two-dimensional groundwater scheme and explored potential benefits of an explicit representation of a so far missing major compartment of the water cycle in the modelling system. The low additional computational cost and modest data requirements enable applications over large model domains-as typically found in fully coupled model setups-and over long time periods, such as the dynamical downscaling of climate scenarios. Results from a real-world application in the upper Danube catchment suggest that an explicit groundwater treatment can substantially improve the ability of the uncalibrated model to reproduce observed runoff volumes, which is particularly beneficial for model applications in poorly or even ungauged basins. The groundwater-enhanced model configuration exhibits a higher topographically-induced heterogeneity in soil moisture content and evapotranspiration, which lends itself to further investigation with the fully coupled modelling system to study potential feedbacks between soil moisture and precipitation. It is important to note, however, that the study region presented here is energy limited, which constrains the transferability of our results to other hydro-climatic regions and underlines the necessity of further model evaluation. Yet, the results of this study demonstrate that a better representation of groundwater is highly desirable in the current LSM of WRF-Hydro in particular and for LSMs with free drainage lower boundary condition in general.