A high-resolution simulation of groundwater and surface water over most of the continental US with the integrated hydrologic model ParFlow v3

Abstract

Interactions between surface and groundwater systems are well-established theoretically and observationally. While numerical models that solve both surface and subsurface flow equations in a single framework (matrix) are increasingly being applied, computational limitations have restricted their use to local and regional studies. Regional or watershed-scale simulations have been effective tools for understanding hydrologic processes; however, there are still many questions, such as the adaptation of water resources to anthropogenic stressors and climate variability, that can only be answered across large spatial extents at high resolution. In response to this grand challenge in hydrology, we present the results of a parallel, integrated hydrologic model simulating surface and subsurface flow at high spatial resolution (1 km) over much of continental North America ($\sim$ 6.3 M km${}^{\mathrm{2}})$. These simulations provide integrated predictions of hydrologic states and fluxes, namely, water table depth and streamflow, at very large scale and high resolution. The physics-based modeling approach used here requires limited parameterizations and relies only on more fundamental inputs such as topography, hydrogeologic properties and climate forcing. Results are compared to observations and provide mechanistic insight into hydrologic process interaction. This study demonstrates both the feasibility of continental-scale integrated models and their utility for improving our understanding of large-scale hydrologic systems; the combination of high resolution and large spatial extent facilitates analysis of scaling relationships using model outputs.

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; Anyah et al., 2008; Kollet and Maxwell, 2008; Maxwell and Kollet, 2008a; Jiang et al., 2009; Rihani et al., 2010; Maxwell et al., 2011; Williams and Maxwell, 2011; Condon et al., 2013; Taylor et al., 2013). While local observations and remote sensing can now detect changes in the hydrologic cycle from small to very large spatial scales (e.g., Rodell et al., 2009), theoretical approaches to connect and scale hydrologic states and fluxes from point measurements to the continental scales are incomplete. In this work, we present integrated modeling as one means to address this need via numerical experiments.

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 Maxwell et al., 2014) 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 stream–aquifer and land–energy interactions (Maxwell et al., 2007; Kollet and Maxwell, 2008; Rihani et al., 2010; Condon et al., 2013; Camporese et al., 2014), and building our understanding of the propagation of perturbations such as land cover and anthropogenic climate change throughout the hydrologic system (Maxwell and Kollet, 2008a; 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 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; Maxwell, 2013). As such, existing large-scale studies of the subsurface have focused on modeling groundwater independently (Fan et al., 2007; Miguez-Macho et al., 2007; Fan et al., 2013) and classifying behavior with analytical functions (Gleeson et al., 2011a). Similarly, continental-scale modeling of 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 groundwater has been shown to be important across scales (Krakauer 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 that can only be addressed with an integrated approach.

The importance of groundwater–surface-water interactions in governing scaling behavior 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 the challenges observing 3-D groundwater flow across a wide range of scales. 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 nonlinear, and heterogeneity in hydraulic properties 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 nonlinear flow equations at high spatial resolution while numerically scaling these physical processes up to a large spatial extent (i.e., continental scale).

In this study, we simulate surface and subsurface flow at high spatial resolution (1 km) over much of continental North America (6.3 M km${}^{\mathrm{2}})$, which is itself considered a grand challenge in hydrology (e.g., Wood et al., 2011; Gleeson and Cardiff, 2013). The domain is constructed entirely of available data sets including topography, soil texture, and hydrogeology. This simulation solves surface and subsurface flow simultaneously and takes full advantage of massively parallel, high-performance computing. The results presented here should be viewed as a sophisticated numerical experiment, designed to diagnose physical behavior and evaluate scaling relationships. While this is not a calibrated model that is intended to match observations perfectly, we do verify that behavior is realistic by comparing to both groundwater and surface-water observations.

The paper is organized as follows: first, a brief description of the model equations are provided including a description of the input variables and observational data sets 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.

Methods

The model was constructed using the integrated simulation platform ParFlow (PARallel Flow) (Ashby and Falgout, 1996; Jones and Woodward, 2001; Kollet and Maxwell, 2006) utilizing the terrain following grid capability (Maxwell, 2013). 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 (Richards, 1931) in three spatial dimensions given as $$S_{\mathrm{s}}{S}_{\mathrm{w}}\left(h\right){\displaystyle \frac{\partial h}{\partial t}}+\mathit{\varphi }{S}_{\mathrm{w}}\left(h\right){\displaystyle \frac{\partial S_{\mathrm{w}}\left(h\right)}{\partial t}}=\mathrm{\nabla }\cdot \mathit{q}+q_{\mathrm{r}},$$ where the flux term $\mathit{q}$ [LT${}^{-\mathrm{1}}$] is based on Darcy's law: $$\mathit{q}=-{\mathbf{K}}_{\mathrm{s}}\left(\mathit{x}\right)k_{\mathrm{r}}\left(h\right)\left[\mathrm{\nabla },\left(h,+,z\right),\cos ,{\mathit{\theta }}_x,+,\sin ,{\mathit{\theta }}_x\right].$$

In these expressions, $h$ is the pressure head [L]; $z$ is the elevation with the $z$ axis specified as upward [L]; K${}_{\mathrm{s}}\mathrm{(}x$) is the saturated hydraulic conductivity tensor [LT${}^{-\mathrm{1}}$]; $k_{\mathrm{r}}$ is the relative permeability [–]; $S_{\mathrm{s}}$ is the specific storage [L${}^{-\mathrm{1}}$]; $\mathit{\phi }$ is the porosity [–]; $S_{\mathrm{w}}$ is the relative saturation [–]; $q_{\mathrm{r}}$ is a general source/sink term that represents transpiration, wells, and other fluxes including the potential recharge flux, which is enforced at the ground surface [T${}^{-\mathrm{1}}$]; and $\mathit{\theta }$ [–] is the local angle of topographic slope, $S_x$ and $S_y$ in the $x$ and $y$ directions and may be written as ${\mathit{\theta }}_x={\tan }^{-\mathrm{1}}{S}_x$ and ${\mathit{\theta }}_y={\tan }^{-\mathrm{1}}{S}_y$, respectively. Note that we assume that density and viscosity are both constant, although ParFlow can simulate density and viscosity-dependent flow (Kollet et al., 2009). The van Genuchten (1980) relationships are used to describe the relative saturation and permeability functions ($S_{\mathrm{w}}\left(h\right)$ and $k_{\mathrm{r}}\left(h\right)$, 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): $$\begin{array}{cc}\mathit{k}& \cdot \left(-,{\mathbf{K}}_{\mathrm{s}},\left(\mathit{x}\right),k_{\mathrm{r}},\left(h\right),\cdot ,\mathrm{\nabla },\left(h,+,z\right)\right)=\\ & & \frac{\partial ∥h,,,\mathrm{0}∥}{\partial t}-\mathrm{\nabla }\cdot ∥h,,,\mathrm{0}∥{\mathit{v}}_{\text{sw}}+\mathit{\lambda }{q}_{\mathrm{r}}\left(\mathit{x}\right).\end{array}$$

In this equation ${\mathit{v}}_{\text{sw}}$ is the two-dimensional, depth-averaged surface-water velocity [LT${}^{-\mathrm{1}}$] given by Manning's equation; $h$ is the surface ponding depth [L] the same $h$ as is shown in Eq. (1). Note that $∥h,,,\mathrm{0}∥$ indicates the greater value of the two quantities in Eq. (3). This means that if $h<\mathrm{0}$ the left hand side of this equation represents vertical fluxes (e.g., in-/exfiltration) across the land-surface boundary and is equal to $q_{\mathrm{r}}(\mathit{x}$) and a general source/sink (e.g., rainfall, ET) rate [LT${}^{-\mathrm{1}}$] with $\mathit{\lambda }$ being a constant equal to the inverse of the vertical grid spacing [L${}^{-\mathrm{1}}$]. This term is then entirely equivalent to the source/sink term shown in Eq. (1) at the ground surface, where $\mathit{k}$ is the unit vector in the vertical, again defining positive upward coordinates. If $h>\mathrm{0}$ then the terms on the right hand side of Eq. (3) are active water that is routed according to surface topography (Kollet and Maxwell, 2006).

Maps of top soil type (applied over the top 2 m of the model) (a), elevation (m a.s.l.) (b), bottom soil type (c), potential recharge, $P$–$E$, (m yr${}^{-\mathrm{1}}$) (d), saturated hydraulic conductivity (m h${}^{-\mathrm{1}}$; applied over the bottom 100 m of the model) (e) over the model domain (f).

[Figure omitted. See PDF]

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; Maxwell, 2013). Utilizing a globally implicit solution allows for 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.$\mathrm{6}\times {10}^{\mathrm{4}}$ processors (Kollet et al., 2010; Maxwell, 2013).

ParFlow solves saturated subsurface flow (i.e., groundwater), unsaturated subsurface flow (i.e., the vadose zone), and surface flow (i.e., streamflow) in a continuum approach within a single matrix. Thus, complete nonlinear 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 (surface fluxes exceed the infiltration capacity; e.g., Horton, 1933) or excess saturation (subsurface exfiltration to the surface system; e.g., Dunne, 1983) (for further discussion see, e.g., Kirkby, 1988 and Beven, 2004). 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 the HydroGeoSphere model (Therrien et al., 2012), the Penn State Integrated Hydrologic Model (PIHM) (Kumar et al., 2009), and the CATchment HYdrology model (CATHY) (Camporese et al., 2010). This is a distinct contrast to more conceptually based models that may not simulate lateral groundwater flow or 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., functional relationships that 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 robust numerical solvers (Maxwell, 2013; Osei-Kuffuor et al., 2014) and exploits high-performance computing (Kollet et al., 2010) to achieve high-resolution, large-extent simulations.

Domain setup

In this study, the model and numerical experiment was directed at the continental US (CONUS) using the terrain following a grid framework (Maxwell, 2013) for a total thickness of 102 m over five model layers. The model was implemented with a lateral resolution of 1 km with nx $=$ 3342, ny $=$ 1888 and five vertical layers with 0.1, 0.3, 0.6, 1.0 and 100 m discretization for a total model dimension of 3342 km $\times$ 1888 km $\times$ 0.102 km and 31 548 480 total compute cells. The model domain and input data sets are shown in Fig. 1. All model inputs were re-projected to have an equal cell size of 1 km $\times \mathrm{1}$ km as shown in Fig. 1. Topographic slopes ($S_x$ and $S_y)$ were calculated from the Hydrological data and maps based on the SHuttle Elevation Derivatives at multiple Scales (HydroSHEDS) digital elevation model (Fig. 1b) and were processed using the r.watershed package in the GRASS geographic information system (GIS) platform. Surface roughness values were constant (10${}^{-\mathrm{5}}$ [h m${}^{-\mathrm{1}/\mathrm{3}}$]) outside of the channels and varied within the channel as a function of average watershed slope. Over the top 2 m of the domain, hydraulic properties from soil texture information of the soil survey geographic database (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). These soil types were mapped to their corresponding category in the property database and those values were used in the model simulation (e.g., saturated hydraulic conductivity, van Genuchten relationships). 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 (Condon and Maxwell, 2013, 2014) and to reflect changes in topography using the $e$-folding relationship empirically derived in Fan et al. (2007): $\mathit{\alpha }=e^{-\frac{50}{f}}$, where $f=\frac{a}{\left(\mathrm{1},+,b,\cdot ,\sqrt{S_x^{\mathrm{2}}+S_y^{\mathrm{2}}}\right)}$. For this analysis $a=20$, $b=125$, and the value of 50 (m) were chosen to reflect the midpoint of the deeper geologic layer in the model. Larger values of $\mathit{\alpha }$ reduce 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. Below the deeper geologic layer, the presence of impermeable bedrock was assumed. This assumption oversimplifies regions that have weathered or fractured systems that contribute to regional flow and aquifer systems deeper than 100 m. These assumptions are necessitated by a lack of data at this scale, not limitations of the model simulation. Note that this complex subsurface data set is assembled from many sources and is subject to uncertainty: heterogeneity within the defined geologic types, uncertainty about the breaks between geologic types, and parameter values assigned to these types. There are breaks across data set boundaries, commonly at state or province and international political delineations. The fidelity and resolution of the source information used to formulate this data set also changes between these boundaries yielding some interfaces in property values.

Plot of observed streamflow (a) and observed water table depth (b).

[Figure omitted. See PDF]

Map of simulated surface flow (m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$) over the CONUS domain with two insets zooming into the Ohio River basin. Colors represent surface flow in log scale and line widths vary slightly with flow for the first two panels.

[Figure omitted. See PDF]

All input data sets are a work in progress and should be continually improved. However, we feel it is important to continue numerical experiments with the data that are currently available, while keeping in mind the limitations associated with every model input. Shortcomings in hydrogeological data sets reflect the lack of detailed unified hydrogeological information that can be applied in high-resolution continental models. This constitutes a significant source of uncertainty, which needs to be assessed, quantified, and ultimately reduced in order to arrive at precise predictions. Still, it should be noted that the purpose of this work is to demonstrate the feasibility of integrated modeling to explicitly represent processes across many scales of spatial variability using best available data. By focusing on large-scale behaviors and relationships we limit the impact of uncertain inputs.

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) field, shown in Fig. 1d, was derived from products developed by Maurer et al. (2002). They developed a gridded precipitation field from observations and simulated evaporation and transpiration fluxes using the variable infiltration capacity (VIC) model. We calculate the average difference between the two from 1950 to 2000 and apply all positive values as potential recharge ($P$–$E$) (negative values were set to zero). The model was initialized dry and the $P$–$E$ forcing was applied continuously at the land-surface upper boundary ($q_{\mathrm{r}}$ in Eq. 1) until the balance of water (difference between total outflow and $P$–$E$) was less than 3 % of storage. For all simulations a nonlinear tolerance of 10${}^{-\mathrm{5}}$ and a linear tolerance of 10${}^{-10}$ were used to ensure proper model convergence.

While this study employs state of the art modeling techniques, it is important to note that the numerical simulation of this problem required significant computational resources. Simulations were split 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 better than 60 % efficiency due to the non-symmetric preconditioner used (Maxwell, 2013). The reason for this is the special architecture of JUQUEEN with only 256MB of memory per core and relatively slow clock rate. Additionally, code 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 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 M 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.

Model results were compared to available observations of streamflow and hydraulic head (the sum of pressure head and gravitational potential). Observed streamflow values were extracted from a spatial data set of current and historical US Geological Survey (USGS) stream gauges mapped to the National Hydrography Dataset (NHD) (Stewart et al., 2006). The entire data set includes roughly 23 000 stations, of which just over half (13 567) fall within the CONUS domain. For each station, the data set 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 gauge records. For comparison, stations without a reported 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. (2007, 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 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 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>\mathrm{0}$ in the top layer of the model in Eqs. 1–3). Overland flow is promoted both by topographic convergence, and surface and subsurface flux; however, with this formulation there is no requirement that all potential streams support flow. Thus, the model captures the generation 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 nonlinear physical processes.

Map of water table depth (m) over the simulation domain with two insets zooming into the North and South Platte River basin, headwaters to the Mississippi River. Colors represent depth in log scale (from 0.01 to 100 m).

[Figure omitted. See PDF]

Scatterplots of simulated vs. observed hydraulic head (a) and surface flow (b).

[Figure omitted. See PDF]

Histograms of simulated and observed water table depth (a), hydraulic head (b), median observed flow (c) and 75th percentile observed flow (d).

[Figure omitted. See PDF]

Distributions of observed and simulated streamflow by basin as indicated.

[Figure omitted. See PDF]

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 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 – all of which are present in the observations. Many of these anthropogenic impacts have been implemented into the ParFlow modeling framework (Ferguson and Maxwell, 2011; Condon and Maxwell, 2013, 2014). Although anthropogenic impacts clearly influence water resources, a baseline simulation allows for a comparison between the altered and unaltered systems in future work.

Next we compare the results of the numerical experiment to observations. As noted previously, this is not a calibrated model. Therefore, the purpose of these comparisons is to evaluate model behavior and physical processes against observations not to generate input parameters. Figure 5 plots observed and simulated hydraulic head and streamflow for the data set shown in Fig. 2. Hydraulic head (Fig. 5a) is plotted (as opposed to water table depth) as 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 influenced by topography. Additional metrics and comparisons are explored below. Simulated streamflow (Fig. 5b) also agrees closely with observations. There is some bias, particularly for smaller flows (which we emphasize by plotting in log scale) which also exhibit more scatter than larger flows, and are likely 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. We see this when linear least squared statistics are computed where the $R^{\mathrm{2}}$ value increases to 0.8.

Distributions of observed and simulated water table depth by basin as indicated.

[Figure omitted. See PDF]

Figure 6 plots histograms of predicted and observed water table depth (a), hydraulic head (b), median (50th percentile) flow, and 75th percentile flows (c, 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 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 (Fan et al., 2013). The observed water tables may include anthropogenic impacts, namely, groundwater pumping, while the model 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 percentile 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 forcing, which might also explain the shallower water table depths.

Plots of topographic slope (a), hydraulic conductivity (b) potential recharge (c), water table depth (d), streamflow (e) and satellite image (f) for a region of the model covering the Platte River basin.

[Figure omitted. See PDF]

Figures 7 and 8 compare observed and simulated flows and water table depths for each of the major basins 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 River which demonstrate better agreement between model simulations 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 River watershed, where water table depths agree (Fig. 8e and g) and flows are overpredicted (Fig. 7e and g), indicates a more complex set of interactions than basic water balance driven by forcing.

Plots of topographic slope (a), hydraulic conductivity (b), potential recharge (c), water table depth (d), streamflow (e) and satellite image (f) for a region of the model covering the upper Missouri River basin.

[Figure omitted. See PDF]

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. We do this as a check to see if and how this numerical experiment compares to real observations. It is important to use a range of measures of success that might be different from that used in a model calibration where inadequacies in model parameters and process might be muted while tuning the model to better match observations. 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 a 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 these will develop into stream loci. 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 northern and southern Platte River systems join the Missouri River. 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<Y<900$ km specifically south of the Platte River) water tables increase and stream densities decrease. The satellite image (panel f) shows increases in vegetation that correspond to shallower water tables and increased stream density.

Plots of scaling relationships for simulated and median observed surface flow. Log-scale plots of surface flow as a function of contributing drainage area derived from the model simulation (a) and observations (b). Individual symbols are colored by aridity index (AI) with blue colors being humid and red colors being arid in panels (a) and (b).

[Figure omitted. See PDF]

Hydraulic conductivity also has a significant impact on water table depth and stream network density. In areas of greater recharge in the eastern portion of Fig. 9c, regions with larger hydraulic conductivity (panel b) show decreased stream network density and increased water table depths. This is more clearly demonstrated in Fig. 10 (a region in the upper Missouri River) where, except for the north-east corner, recharge is uniformly low. Slopes are also generally low (panel a), yet hydraulic conductivities show a substantial increase due to a change in data sets between state and country boundaries (panel b, $X>1250$ km, $Y>1400$ km). The relative increase in hydraulic conductivity decreases hydraulic gradients under steady-state conditions and generally increases water table depth, which in turn decreases stream network density. This change in hydraulic conductivity yields a decrease in the formation of stream networks resulting in an increase in water table depth. Thus, hydraulic conductivity has an important role in partitioning moisture between surface and subsurface flow, also 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 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 data sets.

Finally, the connection between streamflow and drainage area is a classical scaling relationship (Rodriguez-Iturbe and Rinaldo, 2001), which usually takes the power law form $Q=k{A}^n$, where $Q$ is volumetric streamflow [L${}^{\mathrm{3}}$T${}^{-\mathrm{1}}$], A is the contributing upstream area [L${}^{\mathrm{2}}$], and $k$ [LT${}^{-\mathrm{1}}$] and n are empirical constants. While this relationship has been demonstrated for individual basins and certain flow conditions (Rodriguez-Iturbe and Rinaldo, 2001), generality has not been established (Glaster, 2009). Figure 11a plots simulated streamflow as a function of associated drainage area on log–log axes, and Fig. 11b plots the same variables for median observed streamflow from more than 4000 gauging 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=\mathrm{0}$.84. Both Fig. 11a and b are colored by aridity index (AI), the degree 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. 11b) the same behavior is observed, where more humid stations fall along the $n=\mathrm{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.

The model presented here represents a first, high-resolution integrated simulation over continental-scale river basins in North America using the best available data. However, primary input data sets are used (potential recharge, subsurface properties, and topography), which clearly require improvement. For example, higher-resolution simulations are feasible, given that the ParFlow model exhibits better than 80 % parallel efficiency for more than 8 billion compute cells. This could improve the surface and subsurface prediction, although we do not expect the form of the scaling relationships as shown in Fig. 11 to change with an increase in resolution. Higher-resolution simulations would require higher-resolution parameter fields that do not exist at this time. Similarly, model lower boundaries (i.e., the overall thickness of the subsurface) could be extended given information about deeper hydrogeologic formations and their properties. The model domain could be expanded to larger spatial extent, either over more of continental North America, coastlines, or even globally. Thus, the study strongly motivates improved, unified input and validation data sets for integrated hydrologic models at the continental scale, similar to data products available to the atmospheric sciences.

Conclusions

Here we present the results of an integrated, multiphysics-based hydrologic simulation covering much of continental North America at hyperresolution (1 km). This numerical experiment provides a consistent theoretical framework for the analysis of groundwater and surface water interactions and scaling from the headwaters to continental scale (10${}^{\mathrm{0}}$–10${}^{\mathrm{7}}$ km${}^{\mathrm{2}})$. The framework exploits high-performance computing to meet this grand challenge in hydrology (Wood et al., 2011; Gleeson and Cardiff, 2013; Bierkens et al., 2015). 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 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, continental-scale simulation. We simulate an unaltered, or pre-development scenario of groundwater and surface-water flows under steady-state conditions. As such, the discussion focuses on the physical controls of groundwater–surface-water interactions and scaling behavior; however, there are obvious limitations to this scenario and these simulations. Clearly reservoir management, groundwater pumping, irrigation, diversion, and urban expansion all shape modern hydrology. Work has been undertaken to include these features within the ParFlow framework at smaller scales (Ferguson and Maxwell, 2011, 2012; Condon and Maxwell, 2013, 2014) and an important next step is to scale the impacts out to the continent.

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. 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 data sets. As noted throughout the discussion, existing data sets are subject to uncertainty and are clearly imperfect. As improved subsurface characterization becomes available, this information can be used to better inform models and fully understand the propagation of uncertainty in these types of numerical experiments (e.g., Maxwell and Kollet, 2008b; Kollet, 2009). However, while the magnitudes of states and fluxes may change with improved data sets, the overall trends and responses predicted here are not likely to 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 continental scales resolving variances over 4 orders of spatial scales.

This study highlights the utility of high-performance computing in addressing the grand challenges in hydrological sciences and represents an important advancement 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.

Code availability

ParFlow is an open-source, modular, parallel integrated hydrologic platform freely available via the GNU LPGL license agreement. ParFlow is developed by a community led by the Colorado School of Mines and F-Z Jülich with contributors from a number of other institutions. Specific versions of ParFlow are archived with complete documentation and may be downloaded

http://inside.mines.edu/~rmaxwell/maxwell_software.shtml

or checked-out from a commercially hosted, free SVN repository; v3, r693 was the version used in this study. The input data and simulations presented here will be made available and may be obtained by contacting the lead author via email.

The Supplement related to this article is available online at doi:10.5194/gmd-8-923-2015-supplement.

Edited by: J. Neal