Interactive comment on “ Simulation of groundwater and surface water over the continental US using a hyperresolution , integrated hydrologic model ”

Introduction Conclusions References Tables Figures


Introduction
There is growing evidence of feedbacks between groundwater, surface water and soil moisture that moderate land-atmospheric energy exchanges, and impact weather and climate (Maxwell et al., 2007(Maxwell et al., , 2011Anyah et al., 2008;Kollet and Maxwell, 2008; Though introduced as a concept in the literature almost half a century ago (Freeze and Harlan, 1969), integrated hydrologic models that solve the surface and subsurface systems simultaneously have only been a reality for about a decade (VanderKwaak and Loague, 2001;Jones et al., 2006;Kollet and Maxwell, 2006). Since their implementation, integrated hydrologic models have been successfully applied to a wide range of watershed-scale studies (see Table 1 in  successfully capturing observed surface and subsurface behavior (Qu and Duffy, 2007;Jones et al., 2008;Sudicky et al., 2008;Camporese et al., 2010;Shi et al., 2013), diagnosing streamaquifer and land-energy interactions (Maxwell et al., 2007;Kollet and Maxwell, 2008;Rihani et al., 2010;Camporese et al., 2014), and building our 15 understanding of the propagation of perturbations such as land-cover and anthropogenic climate change throughout the hydrologic system (Maxwell and Kollet, 2008;Goderniaux et al., 2009;Sulis et al., 2012;Mikkelson et al., 2013).
Prior to this work, computational demands and data constraints have limited the application of integrated models to regional domains. Advances in parallel solution 20 techniques, numerical solvers, supercomputer hardware and additional data sources have only recently made large-scale, high-resolution simulation of the terrestrial hydrologic cycle technically feasible (Kollet et al., 2010;. As such, existing large scale studies of the subsurface have focused on modeling groundwater independently (Fan et al., 2007Miguez-Macho et al., 2007;) and classifying behavior 25 with analytical functions (Gleeson et al., 2011a). Similarly, continental scale modeling of the surface water has utilized tools with simplified groundwater systems that do not capture lateral groundwater flow and model catchments as isolated systems (Maurer et al., 2002;Döll et al., 2012;Xia et al., 2012). Despite the fact that lateral flow of Introduction groundwater has been shown to be important across scales (Nir et al., 2014). While much has been learned from previous studies, the focus on isolated components within what we know to be an interconnected hydrologic system is a limitation than can only be addressed with an integrated approach. The importance of groundwater surface water interactions in governing scaling be- 5 havior of surface and subsurface flow from headwaters to the continent has yet to be fully characterized. Indeed, one of the purposes for building an integrated model is to better understand and predict the nature of hydrologic connections across scales and throughout a wide array of physical and climate settings. Arguably, this is not possible utilizing observations, because of data scarcity and challenges observing 3-D ground-10 water flow across a wide range of scale. For example, the scaling behavior of river networks is well known (Rodriguez-Iturbe and Rinaldo, 2001), yet open questions remain about the quantity, movement, travel time, and spatial and temporal scaling of groundwater and surface water at the continental scale. Exchange processes and flow near the land surface are strongly non-linear, and heterogeneity in hydraulic properties 15 exist at all spatial scales. As such, a formal framework for connecting scales in hydrology (Wood, 2009) needs to account for changes in surface water and groundwater flow from the headwaters to the mouth of continental river basins. We propose that integrated, physics-based hydrologic models are a tool for providing this understanding, solving fundamental non-linear flow equations at high spatial resolution while numeri-20 cally scaling these physical processes up to a large spatial extent, i.e. the continent.
In this study, we simulate surface and subsurface flow at high spatial resolution (1 km) over much of continental North America (∼ 6 300 000 or 6.3 million km 2 ), which is itself considered a grand challenge in hydrology (e.g. Wood et al., 2011;Gleeson and Cardiff, 2014). This simulation solves surface and subsurface flow simultaneously and 25 takes full advantage of massively parallel, high-performance computing. The domain is constructed entirely of available datasets including topography, soil texture and hydrogeology. Results are compared to observations, and used to diagnose physical behavior and evaluate scaling relationships. The paper is organized as follows: first a brief Introduction description of the model equations are provided including a description of the input variables and observational datasets used for model comparison; next model simulations are compared to observations in a number of ways, and then used to understand hydrodynamic characteristics and to describe scaling.

5
The model was constructed using the integrated simulation platform ParFlow (Ashby and Falgout, 1996;Jones and Woodward, 2001;Kollet and Maxwell, 2006) utilizing the terrain following grid capability . ParFlow is a physically based model that solves both the surface and subsurface systems simultaneously. In the subsurface ParFlow solves the mixed form of Richards' equation for variably saturated flow 10 (Richards, 1931) in three spatial dimensions given as: where the flux term q [L T −1 ] is based on Darcy's law: (S w (h) and k r (h) respectively). These functions are highly nonlinear and characterize changes in saturation and permeability with pressure.
Overland flow is represented in ParFlow by the two-dimensional kinematic wave equation resulting from application of continuity conditions for pressure and flux (Kollet and Maxwell, 2006): In this equation v sw is the two-dimensional, depth-averaged surface water velocity [L T (3) are active water that is routed according to surface topography (Kollet and Maxwell, 2006). The nonlinear, coupled equations of surface and subsurface flow presented above are solved in a fully-implicit manner using a parallel Newton-Krylov approach (Jones and Woodward, 2001;Kollet and Maxwell, 2006;. Utilizing a globally-20 implicit solution allows interactions between the surface and subsurface flow system to be explicitly resolved. While this yields a very challenging computational problem, ParFlow is able to solve large complex systems by utilizing a multigrid preconditioner (Osei-Kuffuor et al., 2014;Ashby and Falgout, 1996) and taking advantage of highly scaled parallel efficiency out to more than 1.6 × 10 4 processors (Kollet et al., 2010;. Physically this means that ParFlow solves saturated subsurface flow (i.e. groundwater), unsaturated subsurface flow (i.e. the vadose zone) and surface flow (i.e. 7322 Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | streamflow) in a continuum approach within a single matrix. Thus, complete non-linear interactions between all system components are simulated without a priori specification of what types of flow occur in any given portion of the grid. Streams form purely based on hydrodynamic principles governed by recharge, topography, hydraulic conductivity and flow parameters, when water is ponded due to either excess infiltration 5 (surface fluxes exceed the infiltration capacity) or excess saturation (subsurface exfiltration to the surface system). Groundwater converges in topographic depressions and unsaturated zones may be shallow or deep depending upon recharge and lateral flows.
The physically based approach used by ParFlow is similar to other integrated hydrologic models such as Hydrogeosphere (Therrien et al., 2012), PIHM (Kumar et al., 2009) and CATHY (Camporese et al., 2010). This is a distinct contrast with more conceptually-based models that may not simulate lateral groundwater flow or may simplify the solution of surface and subsurface flow by defining regions of groundwater or the stream-network, prior to the simulation. In such models, groundwater surface water interactions are often captured as one-way exchanges (i.e. surface water loss to groundwater) or parameterized with simple relationships (i.e. curves to impose the relationship between stream head and baseflow). The integrated approach used by ParFlow eliminates the need for such assumptions and allows the interconnected groundwater surface water systems to evolve dynamically based only on the governing equations and the properties of the physical system. The approach used here requires 20 robust numerical solvers and exploits high-performance computing to achieve high resolution, large extent simulations.

Domain setup
In this study, the model for the Continental US (CONUS) was constructed using the terrain following grid framework (Maxwell, 2013) for a total thickness of 102 m Topographic slopes (S x and S y ) were calculated from the Hydrosheds digital elevation model ( Fig. 1b) and were processed using the r.watershed package in the GRASS 5 GIS platform. Over the top 2 m of the domain, hydraulic properties from soil texture information of SSURGO were applied and soil properties were obtained from Schaap and Leij (1998). Note that two sets of soil categories were available. The upper horizon was applied over the top 1 m (the top three model layers) and the bottom one over the next 1 m (the fourth model layer). Figure 1a and c show the top and bottom soil layers of the model. The deeper subsurface (i.e. below 2 m) was constructed from a global permeability map developed by Gleeson et al. (2011b). These values (Gleeson et al., 2011b) were adjusted to reduce variance Maxwell, 2013, 2014) and to reflect changes in topography using the e-folding relationship empirically-derived in Fan et al. (2007): For this analysis a = 20, b = 125 15 and the value of 50 [m] was chosen to reflect the midpoint of the deeper geologic layer in the model. Larger values of α reduced the hydraulic conductivity categorically, that is by decreasing the hydraulic conductivity indicator values in regions of steeper slope. Figure 1e maps the final conductivity values used for simulation. Note that this complex subsurface dataset is assembled from many sources. As such, there are breaks 20 across dataset boundaries, commonly at State or Province and International political delineations. The fidelity and resolution of the source information used to formulate this dataset also changes across these boundaries yielding some interfaces in property values. All input datasets are a work in progress and should be continually improved. Still, we feel it is important to continue modeling with the data that is currently available to 25 understand limitations and assess their impact on model predictions. Shortcomings in hydrogeological data sets reflect the lack of detailed and unified hydrogeological information that can be applied in high resolution continental models. No-flow boundary conditions were imposed on all sides of the model except the land surface, where the free-surface overland flow boundary condition was applied. For the surface flux, a Precipitation-Evapotranspiration (P -E , or potential recharge) product was derived from a combination of precipitation and model-simulated evaporation and transpiration fluxes for a product very similar to Maurer et al. (2002), shown in Fig. 1d. 5 The model was initialized dry and the P -E forcing was applied continuously at the land surface until the balance of water (difference between total outflow and P -E ) was less than 3 % of storage.
While this study employs state of the art modeling techniques, it is important to note that the numerical simulation of this problem is still non-trivial. Simulations were split 10 over 128 divisions in the x direction and 128 in the y direction and run on 16 384 compute-cores of an IBM BG/Q supercomputer (JUQUEEN) located at the Jülich Supercomputing Centre, Germany. These processor splits resulted in approximately 2000 unknowns per compute core; a relatively small number, yet ParFlow's scaling was still good due to the non-symmetric preconditioner used . Additionally, code 15 performance was improved using efficient preconditioning of the linear system (Osei-Kuffuor et al., 2014). The steady-state flow field was accomplished over several steps. Artificial dampening was applied to the overland flow equations early in the simulation during water table equilibration. Dampening was subsequently decreased and removed entirely as simulation time progressed. Large time steps (10 000 h) were used 20 initially and were decreased (to 1 h) as the stream network formed and overland flow became more pronounced with reduced dampening. The entire simulation utilized approximately 2.5 million core hours of compute time, which resulted in less than 1 week of wall-clock time (approximately 150 h) given the large core counts and batch submission process. 25 Model results were validated against available observations of streamflow and hydraulic head (the sum of pressure head and gravitational potential). Observed streamflow values were extracted from a spatial dataset of current and historical US Geological Survey (USGS) stream gages mapped to the National Hydrography Dataset (NHD)  Stewart et al., 2006). The entire dataset includes roughly 23 000 stations, of which just over half (13 567) fall within the CONUS domain. For each station, the dataset includes location, drainage area, sampling time period and flow characteristics including minimum, maximum, mean and a range of percentiles (1,5,10,20,25,50,75,80,90,95,99) compiled from the USGS gage records. For comparison, stations without a re-5 ported drainage area, stations not located on or adjacent to a river cell in ParFlow, and stations whose drainage area were not within twenty percent of the calculated ParFlow drainage area were filtered out. This resulted in 4736 stations for comparison. The 50th percentile values for these stations are shown in Fig. 2a. Note that these observations are not naturalized, i.e. no attempt is made to remove dams and diversions along these streams and rivers, however some of these effects will be minimized given the longer temporal averages. Hydraulic head observations of groundwater at more than 160 000 locations were assembled by Fan et al. (2007Fan et al. ( , 2013. Figure 2b plots the corresponding water table depth at each location calculated as the difference between elevation and hydraulic head. Note that these observations include groundwater pumping (most 15 wells are drilled for extraction rather than purely observation).

Results and discussion
Figures 3 and 4 plot simulated streamflow and water table depth, respectively, over much of continental North America, both on a log scale for flow ( Fig. 3) and water table depth (Fig. 4). Figure 3 shows a complex stream network with flow rates spanning 20 many orders of magnitude. Surface flows originate in the headwaters (or recharge zones) creating tributaries that join to form the major river systems in North America. Note, as discussed previously, that the locations for flowing streams are not enforced in ParFlow but form due to ponded water at the surface (i.e. values of h > 0 in the top layer of the model in Eqs. 1-3). Overland flow is promoted both by topographic convergence, of the complete stream network without specifying the presence and location of rivers in advance, but rather by allowing channelized flow to evolve as a result of explicitly simulated non-linear physical processes. The insets in Fig. 3 demonstrate multiscale detail ranging from the continental river systems to the first-order headwaters. In Fig. 4, water table depth also varies over 5 five orders of magnitude. Whereas aridity drives large-scale differences in water table depth (Fig. 1d), at smaller scales, lateral surface and subsurface flow processes clearly dominate recharge and subsurface heterogeneity (see insets to Fig. 4). Water tables are deeper in the more arid western regions, and shallower in the more humid eastern regions of the model. However, areas of shallow water table exist along arid river channels and water table depths greater than 10 m exist in more humid regions. Note that this is a pre-development simulation, thus, results do not include any anthropogenic water management features such as groundwater pumping, surface water reservoirs, irrigation or urbanization, although many of these anthropogenic impacts have been implemented into the ParFlow modeling framework (Ferguson and Maxwell, 15 2011; Maxwell, 2013, 2014). While anthropogenic impacts are clearly influential on water resources, a baseline simulation allows for a comparison between the altered and unaltered systems. Figure 5 plots observed and simulated hydraulic head and streamflow for the dataset shown in Fig. 2. Hydraulic head (Fig. 5a) is plotted (as opposed to water table depth) as 20 it is the motivating force for lateral flow in the simulation; it includes both the topography and pressure influences on the final solution. We see a very close agreement between observations and model simulations, though given the large range in hydraulic heads the goodness of fit may be somewhat obscured. Additional metrics and comparisons are explored below. Simulated streamflow (Fig. 5b) also agrees closely with observa- 25 tions. There is some bias, particularly for smaller flows, which also exhibit more scatter than larger flows, and may be due to the 1 km grid resolution employed here. Larger flows are more integrated measures of the system and might be less sensitive to resolution or local heterogeneity in model parameters.  Figure 6 plots histograms of predicted and observed water table depth (a), hydraulic head (b), median (50th percentile) flow and 75th percentile flows (c and d). The hydraulic head shows good agreement between simulated and observed (Fig. 6b). While hydraulic head is the motivation for lateral flow and has been used in prior comparisons (e.g. Fan et al., 2007) both observed and simulated values are highly dependent 5 on the local elevation. Figure 6a plots the water table depth below ground surface, or the difference between local elevation and groundwater. Here we see the simulated water table depths are shallower than the observed, something observed in prior simulations of large-scale water table depth . The observed water tables may include anthropogenic impacts, namely groundwater pumping, while the model 10 simulations do not and this is a likely cause for this difference. Also, because groundwater wells are usually installed for extraction purposes there is no guarantee that the groundwater observations are an unbiased sample of the system as a whole. Figure 6c plots the steady-state derived flow values compared to median observed flow values and Fig. 6d plots these same steady-state simulated flows compared to the 75th per- 15 centile of the observed transient flow at each station. While the ParFlow model provides a robust representation of runoff generation processes, the steady-state simulations average event flows. We see the model predicts greater flow than the 50th percentile observed flows (Fig. 6c) and good agreement between the model simulations and the 75th percentile observed flows (Fig. 6d). This indicates a potentially wet bias in the 20 forcing, which might also explain the shallower water table depths.

GMDD
Figures 7 and 8 compare observed and simulated flows and water table depths for each of the major basin encompassed by the model. Water tables are generally predicted to be shallower in the model than observations with the exception of the Upper and Lower Colorado, which demonstrate better agreement between model simulations 25 and observations than other basins. These histograms agree with a visual inspection of Figs. 2b and 4 which also indicate deeper observed water tables. Figure 8 indicates that simulated histograms of streamflow also predict more flow than the observations. This might indicate that the P -E forcing is too wet. However, a comparison of streamflow for the Colorado, where water table depths agree and flows are overpredicted, indicates a more complex set of interactions than basic water balance driven by forcing. To better diagnose model processes, model inputs are compared with model simulation outputs over example regions chosen to isolate the impact of topographic slope, forcing and hydraulic conductivity on subsurface-surface water hydrodynamics. Figure 9 juxtaposes slope, potential recharge, surface flow, water table depth, hydraulic conductivity and a satellite image composite also at 1 km resolution (the NASA Blue Marble image, Justice et al., 2002) and facilitates a visual diagnosis of control by the three primary model inputs. While the model was run to steady-state and ultimately all the potential recharge has to exit the domain as discharge, the distribution and partitioning between groundwater and streams depends on the slope and hydraulic conductivity. Likewise, while topographic lows create the potential for flow convergence, it is not a model requirement that this create a stream. Figure 9 demonstrates some of these relationships quite clearly over a portion of the model that transitions from semi-arid to more humid conditions as the North and South Platte River systems join the Missouri. 15 As expected changes in slope yield flow convergence, however, this figure also shows that as recharge increases from west to east (X > 1700 km, panel c) the model generally predicts shallower water tables and greater stream density (panels d and e, respectively). Conversely, in localized areas of decreased P -E (e.g. 700 km < Y < 900 km specifically south of the Platte River) water tables increase and stream densities de-20 crease. The satellite image (panel f) shows increases in vegetation that correspond to shallower water tables and increased stream density.
Hydraulic conductivity also has a significant impact on water table depth and, somewhat counter intuitively, stream network density. In areas of greater recharge in the eastern portion of Fig. 9c, regions with larger hydraulic conductivity (panel b) show de-25 creased stream network density and increased water table depths. This more clearly demonstrated in Fig. 10 (a region in the upper Missouri) where, except for the northeast corner, recharge is uniformly low. Slopes are also generally low (panel a), yet hydraulic conductivities show a substantial increase due to a change in datasets between state GMDD 7,2014 Simulation of groundwater and surface water over the continental an important role in partitioning moisture between surface and subsurface flow, even under steady-state conditions. While mass balance requires that overall flow must be conserved, larger conductivity values allow this flow to be maintained within the subsurface while lower conductivities force the surface stream network to maintain this flow. In turn, stream networks connect regions of varying hydrodynamic conditions and may 10 result in locally infiltrating conditions creating a losing-stream to recharge groundwater. This underscores the connection between input variables and model predictions, an equal importance of hydraulic conductivity to recharge in model states and the need to continually improve input datasets. Finally, the connection between stream flow and drainage area is a classical scaling  (Glaster, 2009). Figure 10a plots 20 simulated streamflow as a function of associated drainage area on log-log axes, and Fig. 10b plots the same variables for median observed streamflow from more than 4000 gaging stations. While no single functional relationship is evident from this plot, there is a striking maximum limit of flow as a function of drainage area with a continental scaling coefficient of n = 0.84. Both Fig. 10a and b are colored by aridity index (AI), the degree 25 of dryness of a given location. Color gradients that transition from blue (more humid) to red (more arid) show that humid basins fall along the maximum flow-discharge line, while arid basins have less discharge and fall below this line. For discharge observations (Fig. 10b)  the n = 0.9 line and more arid stations fall below this line. Essentially this means that in humid locations, where water is not a limiting factor, streamflow scales most strongly with topography and area. Conversely arid locations fall below this line because flow to streams is limited by groundwater storage.

5
Here we present the results of an integrated, multiphysics-based hydrologic simulation covering much of Continental North America at hyperresolution (1 km). This simulation provides a consistent theoretical framework for the analysis of groundwater and surface water interactions and scaling from the headwaters to continental scale (10 0 -10 7 km 2 ).
The framework exploits high performance computing to meet this grand challenge in 10 hydrology (Wood et al., 2011). We demonstrate that continental-scale, integrated hydrologic models are feasible and can reproduce observations and the essential features of streamflow and groundwater. Results show that scaling of surface flow is related to both drainage area and aridity. These results may be interrogated further to understand the role of topography, subsurface properties and climate on groundwater table 15 and streamflow, and used as a platform to diagnose scaling behavior, e.g. surface flow from the headwaters to the continent. These presented results are a first-step in high resolution, integrated, continentalscale simulation. We simulate an unaltered, or pre-development scenario of groundwater and surface water flows. As such, discussion focuses on the physical controls of 20 groundwater surface water interactions and scaling behavior; however there are obvious limitations to this scenario. Clearly reservoir management, groundwater pumping, irrigation, diversion and urban expansion all shape modern hydrology. Additionally, the steady-state simulation does not take into consideration temporal dynamics or complex land-surface processes, also important in determining the quantity and fluxes of water. 25 These limitations can all be addressed within the current modeling framework but require transient simulations and additional computational resources. Model performance is also limited by the quality of available input datasets. As noted throughout the discussion, existing datasets are subject to uncertainty and are clearly imperfect. As improved subsurface characterization becomes available this information can be used to better inform models. Though, the magnitudes of states and fluxes may change with improved datasets, the overall trends and responses predicted here are not likely to 5 change. While there are always improvements to be made, these simulations represent a critical first step in understanding coupled surface subsurface hydrologic processes and scaling at large spatial extent. This study highlights the utility of high performance computing in addressing the grand challenges in hydrological sciences and represents an important advancement 10 in our understanding of hydrologic scaling in continental river basins. By providing an integrated model we open up a useful avenue of research to bridge physical processes across spatial scales in a hydrodynamic, physics-based upscaling framework.

g. Lower Colorado (2882 Points)
Water Table Depth (