Surface Depression and Wetland Water Storage Improves Major River Basin Hydrologic Predictions

Surface water storage in small yet abundant landscape depressions—including wetlands and other small waterbodies—is largely disregarded in conventional hydrologic modeling practices. No quantitative evidence exists of how their exclusion may lead to potentially inaccurate model projections and understanding of hydrologic dynamics across the world’s major river basins. To fill this knowledge gap, we developed the first-ever major river basin-scale modeling approach integrating surface depressions and focusing on the 450,000-km2 Upper Mississippi River Basin (UMRB) in the United States. We applied a novel topography-based algorithm to estimate areas and volumes of ~455,000 surface depressions (>1 ha) across the UMRB (in addition to lakes and reservoirs) and subsequently aggregated their effects per subbasin. Compared to a “no depression” conventional model, our depression-integrated model (a) improved streamflow simulation accuracy with increasing upstream abundance of depression storage, (b) significantly altered the spatial patterns and magnitudes of water yields across 315,000 km2 (70%) of the basin area, and (c) provided realistic spatial distributions of rootzone wetness conditions corresponding to satellite-based data. Results further suggest that storage capacity (i.e., volume) alone does not fully explain depressions’ cumulative effects on landscape hydrologic responses. Local (i.e., subbasin level) climatic and geophysical drivers and downstream flowpath-regulating structures (e.g., reservoirs and dams) influence the extent to which depression storage volume in a subbasin causes hydrologic effects. With these new insights, our study supports the integration of surface depression storage and thereby catalyzes a reassessment of current hydrological modeling and management practices for basin-scale studies.


Introduction
With the evolution of powerful computational resources and globally available geospatial data sets (e.g., weather and topography; Bauer et al., 2015;Gasper et al., 2014;Yang et al., 2018), current-generation hydrologic models are increasingly being used across large spatial extents. Consequently, an abundance of process-based models has been developed to predict current and potential future states of water availability, water hazards, and human-water interactions across the world's major river basins (e.g., Du et al., 2018;Maxwell et al., 2015;Pokhrel et al., 2018;Prudhomme et al., 2014;Salas et al., 2018;Schewe et al., 2014;Veldkamp et al., 2017;Wood et al., 2011).
Though state-of-the-science basin-scale models encapsulate hydrologic dynamics by effectively curve-fitting observed data (Clark et al., 2017;Fatichi et al., 2016), it is commonly accepted that current major river basin-scale (hereafter referred to as basin scale) models lack the specificity and fidelity required to adequately represent system patterns and emergent processes. Today's models may therefore be providing the right answers for the wrong reasons (Kirchner, 2006;Ruddell et al., 2019). This problem persists due to limited process representations within hydrologic models along with input data uncertainty and parameter equifinality (Beven & Cloke, 2012;Kauffeldt et al., 2013). Accordingly, researchers are tackling these challenges through improved soil moisture accounting, surface-subsurface interactions, vegetation dynamics, and river routing mechanisms, coupled with spatially distributed assimilation/calibration schemes using remotely sensed data (e.g., Camporese et al., 2015;Endrizzi et al., 2014;Kollet & Maxwell, 2008;Rajib, Merwade, & Yu, 2018;Shen & Phanikumar, 2010). However, there has been surprisingly little emphasis on the hydrologic influence of surface depressions and how including their storage capacities could further enhance models' physical realism and prediction accuracy across the world's major river basins.
Surface depressions on the landscape vary greatly in size, ranging from relatively small, unmanaged water storage systems to large, regularly managed waterbodies such as lakes and reservoirs. Small surface depressions include wetlands embedded within uplands or those along river corridors, ponds, and other similar small waterbodies (Biggs et al., 2017). Water storage is a primary function of these small aquatic systems, yet lakes and reservoirs are frequently the primary focus when quantifying water availability because of their perceived large water storage capacities. Small surface depressions (hereafter referred to as surface depressions) are therefore traditionally not integrated into estimates and models of basinscale hydrologic dynamics (Dang et al., 2020;Gao et al., 2012;Liu et al., 2018;Zhou et al., 2016). Surface depressions are largely overlooked in basin-scale hydrologic model applications, despite being abundant landscape features (Downing, 2010). Yet their potential influence on watershed hydrology has recently become well recognized (e.g., Cohen et al., 2016;Creed et al., 2017;Golden et al., 2017;Lane et al., 2018;Yu & Harbor, 2019).
Surface depressions throughout watersheds perform important sink-lag-source functions affecting downstream waters (Brooks et al., 2018;Gleason et al., 2007;Lane et al., 2018;Leibowitz et al., 2008;Shook & Pomeroy, 2011;USEPA, 2015;Wilcox et al., 2011). An individual depression receives water from an upstream area and may store this water for long periods (sink), allow temporary retention and a delayed flux release (lag), or cycle through evapotranspiration, seepage, overland, and lateral outflow processes to produce a net outgoing flux that reaches other waterbodies (source) (McLaughlin et al., 2014;Rains et al., 2016;Rains, 2011). The presence of numerous surface depressions can therefore considerably affect the cumulative hydrologic response of a watershed to different drivers (changes in climate and land management), collectively impacting the magnitude, timing, and spatial pattern of water flow and distribution (Chu, 2017;Golden et al., 2016;Huang et al., 2011;Jones et al., 2019;Vanderhoof et al., 2016).
With the improved process-level understanding of how surface depressions affect watershed hydrology and increased opportunities for quantifying their water storage capacities through recent geospatial data developments (Jan et al., 2018;Jones et al., 2019), the interdisciplinary hydro-geoscience community has recently started including surface depression storage in spatially explicit models with varying levels of spatial heterogeneity, physical interactions, and connectivity (e.g., Ameli & Creed, 2017;Blanchette et al., 2019;Evenson, Jones, et al., 2018;Fossey & Rousseau, 2016;Golden et al., 2019;Nasab & Chu, 2020). However, these model applications have been limited to small or meso-watershed scales, with drainage areas ranging from a few hectares to several thousand square kilometers. A depression-integrated basin-scale model does not exist in the literature nor does evidence of how these depressional storage systems affect basin-scale hydrologic dynamics.
Despite ubiquitous scientific evidence and the emerging community consensus on the hydrological, biophysical, and biogeochemical importance of surface depressions, calls for future directions in water predictions (Bierkens et al., 2015;Clark et al., 2015Clark et al., , 2017Fatichi et al., 2016) have not explicitly acknowledged the potentially substantive effects of integrating surface depression storage in basin-scale models. Disregarding surface depression storage thus remains a persistent issue in hydro-geoscience (Golden et al., 2014. While the effects of depression storage on model results depend on depressions' abundance and spatial distribution over the landscape (Ameli & Creed, 2019;, the widespread practice of disregarding these natural landscape features elicits a serious concern: What if our model-based understanding of basin-scale hydrologic responses to future environmental conditions are potentially flawed? We submit that because of the omission of surface depressions, management and sustainability decisions based on the existing basin-scale hydrologic models may be imprecise at best and ineffective at worst. The objective of our study is to quantitatively assess whether the inclusion of surface depressions in a process-based model improves hydrologic simulations across the world's major river basins compared to the conventional approach of disregarding them. We conduct our analyses by developing the first-ever depression-integrated basin-scale modeling approach. We focus on the 450,000-km 2 Upper Mississippi River Basin (UMRB) in the United States as a test case and ask the following questions: (a) How do depressionintegrated model simulations respond to an increased upstream abundance of surface depressions? (b) Do surface depressions with small storage capacities substantially alter subbasin hydrologic responses? (c) Do changes imposed by integrating surface depressions into a basin-scale watershed model suggest enhanced physical realism in simulated water balances? We address these questions by (a) focusing on streamflow accuracy at select sites across the URMB, (b) assessing variation in subbasin water yields, and (c) evaluating the spatial consistency of subbasin rootzone wetness conditions with satellite-based estimates. Answers to our study's questions could pave the way for a new surface depressionintegrated hydrologic modeling paradigm and provide important insights for those across the globe studying, managing, and modeling surface water resources.

Methodology
We developed two contrasting model configurations using the Soil and Water Assessment Tool (SWAT; Neitsch et al., 2011) to determine how surface depressions alter hydrologic dynamics across the UMRB:

1.
The conventional model: Surface depressions were not explicitly included in the model, an approach traditionally followed in hydrologic modeling practices.

2.
The depression-integrated model: The same setup as in the conventional model, except surface depressions and their associated water storage capacities, were explicitly included.
In the following subsections, we outline the rationale for selecting UMRB as our study area, methods used to set up the UMRB model, descriptions of model configurations (with and without surface depressions), and the calibration-verification procedure.

Study Area: The Upper Mississippi River Basin
The 450,000-km 2 UMRB ( Figure 1) is predominantly an agricultural landscape, responsible for more than 50% of the total corn and soybean production in the United States (Schnitkey, 2013). It drains 15% of the Mississippi River Basin, which covers over 40% of the contiguous land area of the United States. The UMRB also contributes nearly half of the average annual nitrogenous pollution to the Gulf of Mexico, exacerbating seasonal hypoxic conditions (Goolsby et al., 2000;Turner et al., 2008). Like other major basins across the world (e.g., the Mekong in Southeast Asia; Pokhrel et al., 2018), potential changes in UMRB's hydrologic dynamics from natural and anthropogenic influences could limit fresh water availability and agricultural production, provoke flood and drought hazards, and intensify water quality problems (Quinn et al., 2019;Rajib & Merwade, 2017;Secchi et al., 2011;Wing et al., 2018). Therefore, outcomes from using this important heterogeneous region to understand the influence of surface depression storage may have immediate management implications-both to the UMRB and to other major basins worldwide. Furthermore, as a globally significant basin, the UMRB has ample stream gage stations with long-term data availability which makes model setup, evaluation, and hypothesis testing efficient.

The Hydrologic Model
SWAT is a process-based, semi-distributed, continuous-time hydrologic model capable of simulating landscape water balances, water quality, crop yield, and best management practices related to environmental changes Gassman et al., 2007;Neitsch et al., 2011). SWAT divides a watershed into multiple subbasins, which are further discretized into Hydrologic Response Units (HRUs) with identical land use, slope, and soil characteristics (Neitsch et al., 2011). Increased research interest in the hydrological, biophysical, and biogeochemical functions of surface depressions has led to a concomitant trend of depression-integrated SWAT applications in small-or meso-scale watersheds (e.g., Almendinger et al., 2014;Evenson et al., 2016;Evenson et al., 2018,b;Liu et al., 2008;Nasab et al., 2017;Wang et al., 2008). Because SWAT is also frequently used for modeling processes within major river basins worldwide (e.g., Abbaspour et al., 2015;Du et al., 2018;Rajib et al., 2020;Schuol et al., 2008), it is an ideal tool to explore how surface depressions influence basin-scale hydrologic predictions.
Geospatial input data, primarily topography, land use, soil texture, and weather forcing required to set up the URMB model, are summarized in Table 1. In the conventional model configuration, we divided the UMRB into 279 subbasins (with an average subbasin drainage area of 1,613 km 2 ). We initially discretized this model into 6,500 HRUs during our preliminary modeling experiment. However, to reduce computational intensity, we tested whether a discretized model using the 279 subbasins alone, that is, without 6,500 HRUs, would be sufficient for our research needs. When we compared the average annual water yield outputs, the 6,500 HRU-scale experimental model did not meaningfully differ from the subbasin-scale model. Therefore, following Du et al. (2018), Faramarzi et al. (2017), and Schuol et al. (2008), we elected to minimize computational overhead by not further discretizing the subbasins into smaller spatial units (e.g., the HRUs).
Subsurface drainage in low-slope, poorly drained agricultural areas is a key factor affecting UMRB's water balance. Therefore, artificial subsurface tiles drain ~20% of the total basin area (Srinivasan et al., 2010). To explicitly quantify these effects, we delineated the tiled subbasins using a geodatabase which was consistent with our land use and soil inputs (Nakagaki & Wieczorek, 2016; Table 1; Figure S1). We adopted the tile-drain parameter values that are most commonly used in the Midwestern United States (Hutchinson & Christiansen, 2013;Moriasi et al., 2012; see Figure S1).
We used the Penman-Monteith equation for evapotranspiration estimation , the Curve Number method for infiltration-surface runoff generation, and the variable storage method for river routing (Neitsch et al., 2011). We parameterized each of the 15 major lakes and reservoirs in the UMRB (Lehner & Doll, 2004) using their unique hydraulic design and storage-discharge information obtained from the U.S. National Inventory of Dams (USACE, 2018; Table 1). Specifically, we incorporated the maximum water surface area and storage volume, area and volume in normal operating conditions, and design release rate at the dam outlet. Although this release rate was used to constrain the model's reservoir routing scheme, we fit (via direct insertion assimilation) available daily dam outflow data at five locations (supporting information Figure S2).
The simulation length for both model configurations was 10 years (2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017). This 10year simulation allowed us to analyze the role of surface depression storage based on longterm average model outputs while also minimizing the sensitivity of our interpretations to the interannual variations in climatic conditions. The significance of difference between the two output time-series (conventional and depression-integrated simulations) was quantified with a paired t test (KSU Libraries, 2017).

Model Configurations
2.3.1. The Conventional "No-Depression" Model-The conventional model described in section 2.2 did not include surface depressions and their storage volume in hydrologic simulations, representing the traditional approach for watershed modeling, hypothesis testing, and management decisions. This configuration was our baseline to measure the degree of change in simulated streamflow, water yield, and rootzone soil moisture when surface depressions were included in the subsequent configuration.

The Depression-Integrated Model-
The depression-integrated model was the same as the conventional model, except with the inclusion of surface depressions and their water storage capacities (i.e., maximum fillable volume). Use of elevation (to compute depression depths) and areal extent has emerged as an efficient way to estimate surfacewater storage volume (Gao, 2015). Following this trend, we delineated potential surface depressions from topographic data (digital elevation model, DEM) and estimated the maximum surface area and storage volume therein using a novel, efficient algorithm Wu & Lane, 2016). Briefly, we first filled the original 30-m DEM using a depression-filling algorithm to create a "depressionless" DEM. Subsequently, we subtracted the original DEM from the filled DEM to identify ~455,000 depressions (each >1 ha), then we estimated the area and volume of every depression based on a statistical analysis of the DEM cells comprising that depression (see Wu & Lane, 2016;Wu et al., 2019 for details). In a previous study (i.e., Wu et al., 2019), we validated the depression volumes in selected subbasins where high-resolution bathymetric LiDAR DEMs are available.
To avoid DEM's vertical inaccuracy propagating into estimated depression volumes, some studies considered empirical area-to-volume conversion equations (e.g., Gleason et al., 2007;Jones et al., 2017). These empirical approaches may lack accuracy for depressions with complex morphometry (see, e.g., Wu & Lane, 2016). Further, most of the existing global databases (see Aires et al., 2017;Lehner & Doll, 2004;Pekel et al., 2016;Tootchi et al., 2019;Yamazaki et al., 2015) do not provide estimates of depression volumes (i.e., databases that are available are mostly useful for generating estimates of surface areas, though at coarse spatial resolutions). Our approach to efficiently estimate both area and volume of surface depressions from any given topography data Wu & Lane, 2016) fills these information gaps and affords a reproducible, first step at quantifying depressioneffects in the world's major river basins.
The individual depression-volume in the UMRB ranged between 5 × 10 −2 m 3 and 1 × 10 8 m 3 with the basin-average volume around 6 × 10 4 m 3 ( Figure 1). However, we eliminated depressions smaller than 1 ha to keep computational overhead reasonable. Further, to distinguish the cumulative effect of surface depressions from that of the large, regularly managed waterbodies, which are traditionally included in hydrologic modeling, both our model configurations explicitly integrated hydraulic geometry and storage-discharge data for the major lakes and reservoirs (section 2.2). As a result, we excluded those lakes and reservoirs from the DEM-based area-volume estimation. Any depression that fell within the channel width was also removed to avoid overestimating in-channel storage. Finally, we aggregated the depression area and volume for each subbasin to have a parsimonious simulation of depression hydrology across the UMRB (see below).

Conceptualizing Depressional Storage Capacity-
The subbasin-level aggregation or "lumping" of surface depressions essentially follows the hydrologic equivalent wetland (HEW) concept . It is assumed that the depressions within a subbasin exhibit similar hydrologic responses; if all these depressions are virtually aggregated to form a single HEW, their cumulative subbasin-level effect would be equivalent to what the HEW would alone produce. This concept, often implemented in process-based semi-distributed models (Feng et al., 2013;Fossey et al., 2015;Neitsch et al., 2011;Rahman et al., 2016;Voldseth et al., 2007), is associated with empirical power-law equations to simulate the hydrology of surface depressions via time-varying area-volume relationships (see, e.g., Neitsch et al., 2011). The power-law equations in our model included shape factors, defined as a function of existing water storage volume in the HEW, to represent the variable non-wetted depression geometry and hence the potential fillable volume in successive time steps. Because HEWs are aggregated virtual waterbodies without having specific geolocation within respective subbasins, it is infeasible to represent explicit depression-to-stream connectivity and hydrologic transport between a depression and its corresponding upland drainage area (see, e.g., SWAT-based algorithms developed by Evenson et al., 2016;Lee et al., 2018;Nasab et al., 2017). Therefore, we estimated inflow from the entire non-depressional portion of the subbasin at every simulation time step. According to this inflow, and the evapotranspiration and seepage losses that were extracted from the existing storage volume, we estimated the "net change" of water storage in a HEW (note, storage capacity = maximum fillable volume; section 2.3.2). Finally, this "net change" was either retained in the HEW or directly poured in the channel, using the maximum storage volume of the HEW as the threshold prior to it "spilling" (Neitsch et al., 2011). Despite such parsimony, HEW conceptualization produces acceptable results in diverse geophysical settings (e.g., Feng et al., 2013;Rahman et al., 2016).
The depression-integrated UMRB model required defining two additional parameters specific to surface depression hydrology (hydraulic conductivity and rate of evaporation). Suitable values for these two parameters were adopted from Golden et al. (2019) and kept spatially constant throughout the basin (see Table S1).

Calibration and Verification
We followed an identical calibration protocol for both model configurations (conventional and depression integrated). The calibration involved average monthly streamflow data from 25 U.S. Geological Survey gage stations across the basin (Figures 2a and S3). The calibration length was 10 years (2008-2017) following a 3-year initialization (2005)(2006)(2007). We applied a multisite optimization scheme using the Sequential Uncertainty Fitting algorithm-version 2 (SUFI-2), which is a semi-automated inverse modeling procedure available inside SWAT-CUP platform (Abbaspour, 2015). For the conventional model configuration, calibration parameters and their initial ranges (to start the SUFI-2 optimization process) were adopted from a previous UMRB SWAT setup (Rajib & Merwade, 2017; see Table S1). For consistency, the depression-integrated model used an identical set of calibration parameters. A weighted Kling-Gupta efficiency (KGE; Gupta et al., 2009;Knoben et al., 2019;Rajib, Evenson, et al., 2018) was the objective function to measure the association between simulated and gage streamflow data. KGE ranges from −∞ to 1. Model accuracy increases as KGE values move closer to 1. SUFI-2 identified the best parameter combination corresponding to the highest possible KGE values and hence the most optimal streamflow simulations across all calibration locations.
To enable a post hoc evaluation that is independent of calibration, we verified the most optimal streamflow for the same simulation period (2008-2017) Table 1). Briefly, we took the seasonal-average rootzone SMAP data (volumetric water content) for the 2016 crop growing season (June-August), aggregated estimates at subbasin level, and then performed a spatial normalization across all subbasins to create a dimensionless relative wetness index. This spatial aggregation-normalization approach reduced the effects of uncertainties common to any remotely sensed data (e.g., Chen et al., 2019). It also minimized the conceptual and structural inconsistencies between SMAP and a hydrologic model (e.g., spatial and vertical heterogeneity, and rootzone structure; Koster et al., 2018), therefore affording their one-to-one comparison at a target spatial scale (here, subbasins). To facilitate a focused assessment of depression-effects, we compared SMAP data with model simulations over a 18,300-km 2 area at the depression-dominated north-west region of the basin (see Figure 1).

Effect of Surface Depression Storage on Streamflow Accuracy
The inclusion of surface depression storage affected streamflow simulation accuracy with higher or unchanged KGE values compared to the conventional model at 67% of the target stream gage locations (increased in 16, unchanged in 4, and decreased in 10; N = 30; Figures   2a and 2b). Because the primary difference in the two models was the absence or presence of surface depressions, our results suggest that the lack of surface depression storage was the driving factor for the conventional model to largely overestimate streamflow in areas where depression storage volumes are substantive (e.g., verification site #05316580 in Figure 2c). The relative improvements in streamflow KGEs at the target locations were correlated with their respective upstream abundance of surface depression storage (Figure 2d)-a prominent signature of enhanced process representation in the depression-integrated model. In fact, at stream gages where KGE decreased slightly, upgradient surface depression storage was minimal ( Figure 2d). However, while model improvements may be reflected in KGE values closer to 1 for simulated streamflow, other improvements associated with making the model more physically realistic via increased inclusivity of hydrologic processes associated with surface depressions (e.g., spatially varying changes in water yields) cannot be fully captured by this objective function.
In the subsequent sections, we therefore evaluate the additional differences in internal hydrologic dynamics of the two contrasting model configurations (conventional and depression integrated) focusing on landscape water yield (Figures 3 and 4) and rootzone wetness conditions ( Figure 5).

Surface Depressions Alter Subbasin Water Yield Predictions
We found substantially different subbasin water yields between the conventional and depression-integrated configurations ( Figure 3). Water yield, the amount of runoff generated from the landscape after accounting for all hydrologic losses (e.g., evapotranspiration and infiltration), is the most widely used index of water availability (e.g., Ellison et al., 2012;Fekete et al., 2002;Milly et al., 2005;Zhou et al., 2015). As previously noted, studies addressing water yield and availability often overlook the hydrologic effects of surface depressions (e.g., the conventional model configuration; Figure 3a). We therefore hypothesized that inclusion of surface depression storage would change water yield simulations both in their spatial patterns and subbasin-level volumes. Our depressionintegrated UMRB model supported this hypothesis showing 10-40% reduction in subbasin (5% reduction in basin scale) annual water yields compared to the conventional model ( Figure 3b). A visual analysis between Figures 1 and 3 indicates that the spatial density and associated storage volumes of surface depressions ( Figure 1) correlate with the spatial variations in simulated water yields (Figure 3b). These results can also vary from year to year depending on the magnitude of precipitation ( Figure S4).
Given our methodological approach, the variation in water yield between the two models stemmed solely from surface depressions (e.g., floodplain and non-floodplain wetlands, ponds, and other water storage systems). This is a new finding commonly overlooked in basin-scale/global water availability assessments (e.g., Liu et al., 2018;Veldkamp et al., 2017;Zhou et al., 2016). While Lehner et al. (2011) similarly suggested the considerable potential influence of small surface water storage systems on downstream waters, our study is, to our knowledge, the first to quantify this phenomenon in one of the world's major river basins.
Simulated differences in water yield time-series between the two model configurations were statistically significant (p < 0.10) in ~70% of the basin area, which confirms the cumulative influence of surface depressions on basin-scale hydrologic dynamics (Figure 4a). However, additional evidence suggests that surface depression storage capacity (i.e., volume) alone does not fully explain their influence on downstream waters. Substantially different HEW volumes in two similarly sized subbasins can result in similar effects on the respective hydrologic responses at subbasin outlets. For example, in Figure 4a, two subbasins X and Y have similar drainage areas, yet the HEW volume is substantially larger in Y. The poorly drained soil and agricultural land use in subbasin Y suggests greater runoff potential compared to the well-drained soil and forested land use in subbasin X, but the smaller amount of precipitation in subbasin Y nullifies that potential. This, combined with the lower flow accumulation (milder slope) and less available non-depressional land (larger HEW volume), produces a runoff volume in subbasin Y that is too insufficient to fill the abundant depression storage therein (see, e.g., Nasab et al., 2017). As a result, even with different depression storage capacities, the model simulated reduction in water yields are similar for both subbasin X and Y because of the subbasins' respective structural characteristics (see Figure S5).
To further demonstrate that factors in addition to storage volumes contribute to depressions' effects on water yields across the UMRB, we mapped the subbasins where the depressionintegrated model resulted in statistically significantly different water yield outputs compared to the conventional no-depression model (Figure 4b). Clearly, there are many subbasins where surface depressions do not have a significant effect on water yield despite their relatively substantive storage volumes compared to the rest of the basin (see, e.g., the abundance of depressions in west central UMRB; Figure 1) and vice versa (e.g., the relative dearth of depressions in south/southwest UMRB; Figure 1). Therefore, how significantly depression storage influences landscape water yield/availability depends on a combination of climatic and geophysical drivers (Thorslund et al., 2017).

Surface Depressions Improve Soil Moisture Accounting
The depression-integrated model generally retained more water on the landscape ( Figure  4a). A simple justification for this response is that the runoff generated in the depressionintegrated model had to fill an extra "water-storage bucket" (a HEW; on 2.3.3). This retention correspondingly led to increased soil moisture content in the rootzone (Figure 5a) -a fundamental process associated with the "sink-lag-source" functions of surface depressions Rains et al., 2016;Rains, 2011). The correlation between surface depressions and subsurface wetness has also been recognized in recent studies (e.g., Tootchi et al., 2019). However, as we sought to calibrate the water balance by "curve-fitting" both model configurations with the same set of observed data, some subbasins had decreased water retention (hence, a decrease in rootzone moisture content) and more water runoff from the landscape despite the presence of depressions. We noted that this happened only in a few subbasins, particularly those with smaller HEW volumes (see, e.g., Figure 4a, lower left corner). Above a relative depression storage capacity of ~0.2, water is retained in the landscape. Consequently, retaining more water on the landscape and having a relatively wet rootzone was the characteristic, and anticipated, change in surface-subsurface hydrologic partitioning when depressions were included in the conventional model ( Figure 5a). However, is this change one in the right direction, meaning does it improve the overall water balance in the model? We analyzed this issue using satellite-based SMAP soil moisture estimates (Figure 5b).
The spatial pattern of rootzone wetness condition derived from SMAP revealed a largely biased water balance in the conventional model ( Figure 5b). Specifically, the conventional model simulated water balances such that subbasins in the northwest UMRB (as an example) were unrealistically "dry" (wetness index ~0.2) compared to the wettest (wetness index = 1) parts of the basin. This result occurred regardless of abundant depression storage in the basin. Those underrepresented dry subbasins in the conventional model appeared to be much wetter in the depression-integrated model and showed a strong spatial resemblance with SMAP estimates (wetness index ~0.4-0.5). These results confirm that including depression storage in conventional modeling practices can improve overall predictability of water balance and that these improvements happen for physically meaningful reasons.

Discussion
The inclusion of surface depressions in a properly parameterized hydrologic model will inevitably produce different-and as this study demonstrates, more accurate-predictions compared to the same model disregarding them. This should stimulate a reassessment of our conventional hydrologic modeling practices across large spatial domains. To gain further clarity on these issues, explore the implications of our findings, and therefore establish a meaningful rationale for depression-integrated hydrologic modeling, we address the following two questions:

1.
What concepts are critical for interpreting the cumulative hydrologic effects of surface depressions in major river basins?

2.
How would a depression-integrated model affect watershed management decisions?

What Concepts Are Critical for Interpreting the Cumulative Hydrologic Effects of Surface Depressions in Major River Basins?
To explain how surface depressions affect downstream waters, the degree of attenuation in high flow events has been the most widely acknowledged indicator (Blanchette et al., 2019;Huang et al., 2011;Shook & Pomeroy, 2011;Wilcox et al., 2011), although maintenance of baseflow is another important consideration (e.g., . While peak flow attenuation helps prioritize land management alternatives via "what-if" scenarios focused on the abundance and spatial distribution of surface depressions (e.g., Ameli & Creed, 2019;, we suggest that it cannot be a universal indicator for quantifying cumulative depression-effects. Our conjecture is based on the following "scale issues" (Blöschl & Sivapalan, 1995) not previously addressed in the literature.
First, the effect of surface depressions on peak flows is best understood in small-to mesoscale relatively unregulated watersheds (e.g., Blanchette et al., 2019;Nasab et al., 2017). In these systems, surface depressions typically demonstrate collinearity with peak flow conditions: As volumetric depression storage increases, peak flows exhibit greater attenuation (Figure 6a; also see Babbar-Sebens et al., 2013). It therefore may be expected that scaling these results further downstream and across much larger domains would result in peak flow attenuations that are proportional to the increased abundance of upstream surface depression storage. However, the converse may be true. The peak flow attenuation signal from surface depressions may dissipate when large gatekeeper lakes/reservoirs are present in the basin-scale river systems with regulated flow conditions (Figure 6b), along with numerous water withdrawal and point source discharge locations (Alexander et al., 2012;Graf, 2006;Zdankus et al., 2008). Our findings ( Figure 6) confirm that spatial-scale effects on depressions' peak flow attenuation capability is largely controlled by downstream flowpath characteristics (Quin et al., 2015), and as such, the local effects of an individual or a small group of surface depressions cannot be generalized over large spatial scales (Quin & Destouni, 2018;Thorslund et al., 2017).
Second, detecting peak flow attenuation through sub-daily or daily simulations almost invariably shows that the predominant downstream hydrologic effects come from large surface depression systems (e.g., ). Yet simulations across longer time scales (months, seasons, and years) may dampen the sub-daily and dailysimulated hydrologic signatures (Joo et al., 2018). As a result, the substantive effects from large depressions (e.g., areas with high HEW volumes in our UMRB model) represented by sub-daily or daily time scales may be insignificant across longer periods of assessment. For example, analyzing water yields using a 10-year average instead of sub-daily or daily time scale may be why some subbasins in the UMRB had relatively substantive depression storage yet a limited influence on water yield, and vice versa (Figure 4).
In this study, we spatially quantified the hydrologic changes between conventional and depression-integrated models at individual subbasin levels, instead of discrete locations across a basin-scale river network. Temporally, we used long-term average volume of water yield as the change indicator, instead of assessing peak flow attenuation in the rivers. This affords analyses of direct correlations between cause (depression storage volumes in subbasins) and effects (volumetric change in subbasin water yield), therefore explicitly informs where and how surface depressions alter hydrologic dynamics, and provides a basinscale yet locally relevant interpretation.

How Would a Depression-Integrated Model Affect Watershed Management Decisions?
It is highly likely that studies disregarding surface depression storage in their models predict future floods to be more hazardous than they may be. Similarly, drought impacts on hydrologic processes in basins with surface depression storage will be over-predicted. We suggest this based on qualitatively comparing our results to some of the major modelingbased studies in the UMRB, the entire Mississippi Basin, and the continental United States in general. For example, multiple studies applying SWAT in the UMRB (Jha et al., 2004;Wu et al., 2012) have predicted a 50% increase in water yield for a 20% increase in the amount of precipitation, suggesting an amplification of future flood magnitudes. But the magnitude of variations in water yields between our two models indicates overpredicting flood hazards in models disregarding surface depressions (Figure 3). These discrepancies in model predictions may create imprecise and/or ineffective response strategies for flood and drought mitigation. Our findings suggest that model integration of surface depression would not only improve simulated accuracy in projected future flow regimes but provide a better informed approach for managing surface depressions as low-cost "nature-based solutions" (Bridgewater, 2018;Thorslund et al., 2017;Zalewski et al., 2018), in addition to the conventional structural measures (Kundzewicz et al., 2019;Tullos, 2018).
A key driver of regional water resources management from the overall water availability standpoint lies with our ability to delineate spatial locations of potential hydrologic changes. A multi-model comparison by Cai et al. (2014) acknowledged how some of the advanced atmospheric-land surface coupled models may misrepresent the spatial variability in hydrologic processes. Four basin/continental-scale models, including the Noah land surface model, Noah with multi-parameterization scheme (Noah-MP), Variable Infiltration Capacity (VIC), and Community Land Model-4 (CLM4)-models driving the widely used North American Land Data Assimilation System (NLDAS; NASA, 2019)-were configured by Cai et al. (2014) without surface depression storage. Results unanimously showed that water across the landscape was not in the right place. The authors demonstrated how simulated water balances, such as those in the depression-dominated northwest portion of UMRB, retained less water on the landscape compared to other parts of the basin. However, they also showed that remotely sensed terrestrial water storage anomalies (derived from Gravity Recovery and Climate Experiment, GRACE; Landerer & Swenson, 2012) revealed distinctively higher water storage capacity in the northwest UMRB compared to the rest of the basin-results that closely resemble the spatial abundance of depression storage we demonstrate here ( Figure 1) and additional evidence supporting the higher rootzone wetness in our depression-integrated model ( Figure 5, section 3.3).
It is therefore reasonable to argue that prior model-based studies designed to assist in spatially explicit watershed management may not be offering the best information to support sustainable decisions. Specifically, studies quantifying the relative contribution of upstream regions to downstream flow variability, separating out the role of climate and land use and attributing potential changes to local/regional driving factors (e.g., Du et al., 2018;Frans et al., 2013;Tao et al., 2014;Wise et al., 2017), may erroneously delineate targeted hydrologic management locations (e.g., Figure 4b) if surface depression storage is not considered in model simulations. This may lead to serious implications for socio-hydrologic modeling involving crop yield, irrigation requirements, and sustainable alternatives in agroecosystems to support human demand (Feng et al., 2018;Srinivasan et al., 2010;Yaeger et al., 2014).
In addition to the altered and improved hydrologic partitioning a depression-integrated model would provide (Figures 3 and 5), it would also activate and likely improve the modeled biogeochemical functions of surface depressions, such as their nutrient attenuation properties from settling, plant uptake, or losses to the atmosphere via denitrification. Golden et al. (2019) recently demonstrated this via a depression-integrated model, suggesting a 7% average decrease in annual nitrate yield from a ~17,000-km 2 agricultural watershed compared to a conventional no-depression model. Using a depression-integrated modeling approach is particularly important for the UMRB because of its water quality problems resulting from intensive agricultural operations (Goolsby et al., 2000). Although hybrid models that integrate physical processes and data-driven statistical regression may simulate nutrient removal and decay in lakes and reservoirs (e.g., SPAtially Referenced Regression On Watershed (SPARROW) attributes model; Robertson & Saad, 2014), inclusion of millions of small depressional storage systems-even through a parsimonious approach (section 2.3.3)-would likely change existing water quality assessments. The overall changes in water quality projections from depression-integrated models would therefore likely impact environmental sustainability decisions (Fant et al., 2017;Li et al., 2016Li et al., , 2017Yen et al., 2016) including estimates of watershed-scale ecosystem services, applicable best management practices, and the economic consequences of both across the world's major river basins (Golden et al., 2017).

Summary and Conclusions
We developed the first-ever major river basin-scale hydrologic model integrating surface depressions and focusing on the ~450,000-km 2 Upper Mississippi River Basin (UMRB) in the United States. To do this, we first estimated the surface areas and storage capacities of ~455,000 surface depressions using a novel topography-based algorithm. We then integrated the aggregated depression areas and volumes at individual subbasins to most efficiently simulate depression hydrology across the UMRB. We conducted a 10-year simulation contrasting two model configurations (depression-integrated and conventional "no depression" models) and concluded the following:

1.
Streamflow simulation accuracy increased in 54%, remained unchanged in 13%, and decreased in 33% of the target locations as a result of including surface depression storage. The relative improvement of streamflow simulation accuracy at target stream locations correlated with the upstream abundance of depression storage, indicating enhanced process representation in the depression-integrated model.

2.
Subbasin water yields from the depression-integrated model significantly altered UMRB's hydrologic response compared to the conventional model across 70% of the basin area (~315,000 km 2 ). This significant variation was additive after accounting for major lakes and reservoirs typically considered in conventional modeling practices.

3.
A characteristic shift in the partitioning of hydrologic fluxes due to the inclusion of surface depressions produced a physically realistic water balance. Specifically, the depression-integrated model generally increased the overall wetness of the rootzone, showing strong spatial resemblance with satellite-based estimates.
While these results suggested improved hydrologic prediction by a depression-integrated model, we performed additional analyses to demonstrate that substantially different depression storage capacities (i.e., volumes) in similarly sized subbasins can impose similar effects on the respective hydrologic responses at subbasin outlets. To fully understand where and how surface depressions may be influential, it is therefore necessary to consider how climatic and geophysical drivers interact with the presence of surface depressions on the landscape.
We also highlighted conceptual issues, specifically how the spatial and temporal scale of analyses affect-and may impede-interpretations of depressions' cumulative hydrologic effects on downstream waters. For example, we analyzed how surface depression storage attenuated peak flows in the UMRB, and we concluded that unlike small watersheds, flowpath regulation in large river basins (i.e., gatekeeper reservoirs and dams) is a major factor affecting the quantification of cumulative depression-effects. Therefore, our results suggest that the localized hydrologic effects of an individual depression or small group of depressions cannot be generalized across basin scales.
The findings presented in this study stimulate the need for re-envisioning our conventional hydrologic modeling practices. Yet such re-envisioning would require addressing the big data challenges associated with depression-integrated modeling. Recent developments of globally accessible, sufficiently high-resolution waterbody maps, efficient algorithms to retrieve depression area and volume from hyper-resolution topographic data, and tools for seamless data-model integration across any spatial scale of interest will advance depressionintegrated modeling by improving the physical representation of the landscape and quantification of surface depression storage. The next challenge is to use these data in the most computationally efficient manner to propel the hydro-geoscience community closer to depression-integrated hydrologic modeling as the "new normal."

Supplementary Material
Refer to Web version on PubMed Central for supplementary material.   Effect of surface depressions on subbasin water yield: (a) 10-year average annual water yield (mm/year) simulated by the conventional "no-depression" model; (b) percent reduction in water yield due to the inclusion of surface depressions.  (a) Difference in water yields with varying surface depression water storage capacities across all subbasins. Each data point in this scatter-plot refers to a corresponding subbasin (N = 279); (b) subbasins with significant differences in water yield relative to the conventional model (significant at p < 0.1 level in 198 subbasins, out of which 111 were significant at p < 0.01). Surface depression storage capacity of a subbasin was estimated as the aggregated volume of all surface depressions existent therein. To express depression storage capacity using a uniform scale regardless of subbasin size, first the aggregated depression volume in each subbasin was converted into a volume per unit subbasin area (V/A) ratio; then this V/A ratio was normalized relative to the minimum and maximum V/A ratios across all subbasins. Subbasins X and Y are explicitly demarcated to show that substantially different depression storage capacities in two similarly sized subbasins ( Figure  S5) can have similar effects on respective water yield outputs. Rajib et al. Page 29 Water Resour Res. Author manuscript; available in PMC 2021 July 06.

Figure 5.
(a) Increased rootzone wetness in the depression-integrated model indicating a characteristic shift in water balance simulation compared to the conventional model. Each data point in this scatter-plot corresponds to a subbasin (N = 279). To address inconsistencies across subbasins in representing the rootzone and therefore allow a uniform comparison throughout the basin and potentially with satellite-based estimates, 10-year average annual simulated rootzone soil moisture content (mm) in each subbasin was normalized relative to the minimum and maximum soil moisture content across all subbasins; (b) spatial resemblance (bias) of subbasin rootzone wetness condition in the depression-integrated (conventional) model compared to the satellite-based estimates from Soil Moisture Active Passive (SMAP) mission. Values in the spatial maps are relative to the driest and wettest conditions throughout the UMRB.    (2018) Weather forcing Total daily precipitation, minimum, and maximum daily temperature from 1-km gridded Daymet product Thornton et al. (2018) Solar radiation, wind speed, and relative humidity from SWAT's built-in historical database (weather generator) b Using SMAP data for the entire simulation period (2008-2017) was not feasible because of data availability since 31 March 2015. Specifically, we used SMAP data for one crop growing season (1 June-31 August 2016) because this is the most "active" period in terms of vegetation growth and subsequent soil moisture dynamics. Further, data obtained for this period are not affected by snow cover and frozen soil (Reichle et al., 2017).