The Earth has entered the Anthropocene era that is dominated by the impacts of humans (Sun et al., 2017). Today, we are living in an increasingly urbanized world with about one-half of the world population found in urban areas, and the urban population is projected to rise to 66% by 2050 (United Nations, 2014). Meanwhile, urban land uses increased by over 34% from 1980 to 2000 and is projected to double by 2030 globally, mostly in developing counties (Alig et al., 2004;Seto et al., 2012).
Rapid urbanization poses serious stresses to watershed ecosystem structure, function, and services such as water quality degradation (Grimm, Faeth, et al., 2008;Sun & Lockaby, 2012;Sun & Caldwell, 2015), localized climate impacts such as Urban Heat Island (UHI) and Urban Dry Island phenomena (Hao et al., 2018), and increased water demand in cities (Hao, Sun, Liu, & Qian, 2015;Sanchez et al., 2018). Watershed hydrology plays a critical role in regulating water quality, aquatic ecosystems, wildlife habitats, and human health (Sun et al., 2017;Sun & Lockaby, 2012). Forest hydrologists have long been interested in the hydrological consequences of converting forests to urban uses and forest management to provide water for urban populations in the eastern United States since the 1960s (Douglass, 1983;Lull & Sopper, 1969). There are renewed interest in quantifying hydrological impacts of urbanization and climate change and variability within the forest hydrology community (DeWalle, 2003;DeWalle et al., 2000;Martin et al., 2017). The impacts of urbanization on watershed water yield (Caldwell et al., 2012;Hao, Sun, Liu, Wan, et al., 2015), and specific hydrological processes such as stormflow, peakflow, and baseflow (Price, 2011) have been increasingly studied worldwide (Oudin et al., 2018;Sunde et al., 2018). However, our knowledge of the hydrological effects of urbanization at the watershed level is still limited and fragmented (Oudin et al., 2018), preventing us from developing national policies and science-based guidelines for mitigating the effects of urbanization on water resources. For example, state and federal regulatory agencies such as the U.S. EPA (2003) have long been using a "generic approximation" to describe how urban imperviousness affects stormflow, evapotranspiration (ET), and infiltration and guide stream restoration effort across the nation (Livingston & McCarron, 1992). However, lacking quantitative national data, such a simplified illustration of the water balance and its hydrological response to urbanization developed for a specific area (i.e., Florida) (Livingston & McCarron, 1992) may not be appropriate although it has been widely cited as a standard conceptual model in the literature (Arnold & Gibbons, 1996;Paul & Meyer, 2001).
Indeed, urbanization impacts on watershed hydrology and the underlying mechanisms are highly variable and complex (Caldwell et al., 2012;Martin et al., 2017). The majority of existing studies suggests that urbanization increases impervious surfaces, reduces soil infiltration (Price, 2011), and thus causes an increase in high flows and total flow (Kumar et al., 2018;Kundu et al., 2017b;Oudin et al., 2018). In addition, other hydrological processes such as vegetation ET also plays a significant role (Hao, Sun, Liu, Wan, et al., 2015). The magnitude and forms of disturbances in land use/land cover (LULC) are a major factor affecting annual water yield (Awotwi et al., 2015;Martin et al., 2017;Zipper et al., 2018). For example, converting grasslands to urban lands, or wetlands to cropland, or croplands to orchards reduced water yield (Awotwi et al., 2015;Bieger et al., 2015). In contrast, the loss of paddy fields caused a rather large rise in streamflow and groundwater level in a humid rapidly urbanizing watershed in southern China (Hao, Sun, Liu, Wan, et al., 2015). Surprisingly, some studies did not find any significant impacts of urbanization on water yield (Konrad et al., 2002;Kumar et al., 2018;Rose & Peters, 2001;Rouge & Cai, 2014). The observed variability of hydrologic response to urbanization has been attributed to the differences in the magnitude of urbanization (e.g., imperviousness; Weng, 2001), local climate (e.g., rainfall and temperature) (Ahmed et al., 2017), LULC characteristics (Kundu et al., 2017b), and temporal scale examined (Weng, 2001). However, to our knowledge, there has not been a comprehensive effort to evaluate the relative effects of these factors on hydrologic response to urbanization at a large scale. Therefore, there is a critical need to comprehensively quantify the potential impacts of future urbanization on water balances across a diverse climate, LULC, and urbanization features. Such information is extremely important for urban planning and land management at a broad scale (Grimm, Foster, et al., 2008) to allocate limited watershed ecosystem restoration resources effectively.
Our current scientific understanding of the hydrologic impacts of urbanization is mostly based on small scale theoretical modeling using traditional engineering principles (Livingston & McCarron, 1992) that often ignores the role of vegetation (Wang et al., 2008). Empirical monitoring or retrospective studies (Oudin et al., 2018) are challenged by the effects of concomitant climatic change and variability (Martin et al., 2017;Todd et al., 2007) that are often coupled with the urbanization processes (Kumar et al., 2018;Pumo et al., 2017;Putro et al., 2016;Zipper et al., 2018). The traditional "Paired Watershed" approach for detecting the hydrologic effects of a single factor of land cover change such as forest harvesting is generally not applicable to urbanization research (e.g., Baltimore Urban Long Term Ecological Research; Bhaskar & Welty, 2012), although quasi-paired watershed studies have been attempted (Boggs & Sun, 2011). Budyko-based empirical (Teuling et al., 2019;Wang & Hejazi, 2011;Zhou et al., 2015) and process-based mathematical models (Hao, Sun, Liu, & Qian, 2015;Li et al., 2016;Pumo et al., 2017;Zipper et al., 2018) have been used to project the hydrologic effects of natural and anthropogenic disturbances including urbanization and climate change and variability for individual watersheds.
The motivation of this study was to assess the combined effects of urbanization-associated land use/land cover change (LULCC) and the underlying spatially varied climate on water balances by employing a well-tested ecohydrological model at the 12-digit Hydrologic Unit Code (HUC) watershed scale across the conterminous United States (CONUS). The CONUS includes approximately 88,000 HUC12 watersheds and covers a large gradient of urbanization intensities and climates. A consistent set of climatic and biophysical data offers a unique opportunity to examine the watershed hydrologic sensitivity to urbanization under a complex climatic and disturbance gradient at the national scale.
We hypothesized that "hydrologic impacts of urbanization are not created equal." Specifically, our hypotheses were as follows: (1) water yield increases due to both increases in impervious surface area and loss of vegetation and ET (Hypothesis 1-H1) and (2) the magnitude of water yield change varies according to local climate characteristics, the types of previous land cover (e.g., grassland, shrubland, or barren with low biomass and forest with high biomass or wetland with high water availability), and the magnitude of LULC and impervious surface change (Hypothesis 2-H2). These hypotheses were used to guide our modeling analysis to understand key controls of hydrologic responses to urbanization at a national scale.
We used a process-based Water Supply Stress Index (WaSSI) model to project the effects of urbanization on watershed water balances for four time periods: 2000 (baseline), 2010, 2050, 2100. The WaSSI model has been well-validated and applied in the United States (Caldwell et al., 2012(Caldwell et al., , 2015;;Sun, Caldwell, et al., 2011;S. Sun, Sun, et al., 2016), Rwanda (Bagstad et al., 2018), China (Liu et al., 2013), and Australia (Liu et al., 2018). The model proved to be effective for understanding regional ecohydrological effects of forest thinning (Sun, Caldwell, et al., 2015), wildland fires (Hallema et al., 2018), drought (Sun, Sun, et al., 2015b, 2015c), air pollution and climate change (Duan et al., 2016), and water withdrawals (Caldwell et al., 2012), and also ecosystem service trade-off quantifications (Bagstad et al., 2018;Duan et al., 2016) in various physiographic settings. Model structure, algorithms, and inputs and outputs are found in Sun, Caldwell, et al. (2011) and Caldwell et al. (2012) and are described briefly below.
The WaSSI model simulates the water balance and performs streamflow routing at a monthly time step with a spatial resolution of an HUC12 watershed scale (~100 km 2 ). In contrast to the monthly water balance model developed by the U.S. Geological Survey (USGS) (McCabe & Markstrom, 2007;McCabe & Wolock, 2014;Wolock & McCabe, 1999), the WaSSI model considers land cover and was designed to account for the effects of land cover and impervious surface on ET and runoff compositions in addition to climate (Caldwell et al., 2012;Sun, Caldwell, et al., 2011). At its core, WaSSI quantifies ET as a function of potential evapotranspiration (PET), estimated by either temperature-based PET model or FAO Penman-Monteith Grass Reference ET method (ETo), leaf area index (LAI), and precipitation, and further constrained by soil moisture availability (Caldwell et al., 2012;Sun, Alstad, et al., 2011). Unfortunately, MODIS LAI data products exclude LAI values for urban core areas (Zhao et al., 2005). Therefore, we estimated LAI for urban areas by overlaying land use grid layers and MODIS LAI layer in this study. When LAI data were not available for certain land use 30 ×30 m cells, the LAI means of surrounding cells were adopted. The soil hydrology submodel in WaSSI uses several built-in algorithms of the Sacramento Soil Moisture Accounting Model (SAC-SMA) and empirical equations to quantify precipitation portioning to each soil layer, simulating infiltration, surface runoff, soil moisture storage, and subsurface and base flows (Burnash et al., 1973). Snowpack and melting processes are also simulated by the method by McCabe and Wolock (1999). The WaSSI model assumes that precipitation falling on impervious surfaces becomes direct runoff as a component of watershed water yield (Caldwell et al., 2012;Sun, Caldwell, et al., 2011), and ET from impervious surfaces is assumed to be negligible.
The main input data required by WaSSI (Sun, Caldwell, et al., 2011) included historical precipitation and air temperature , percentage of each of the 10 land cover types, and fraction of impervious surfaces within each land use for 2000, 2010, 2050, and 2100, mean monthly (2000-2012) LAI by land cover type, and 11 soil parameters derived from STATSGO-based soil properties (Table 1). The 10 LULC types included three forest classes (i.e., deciduous, evergreen, and mixed forest), shrubland, grassland, cropland, water, wetland, urban, and barren land. However, the ICLUS data sets (U.S. EPA, 2017) used for LULC inputs have only one land use class for forest land. Therefore, we equally divided the forest area by three to meet the data requirements of the WaSSI model. LAI values for each land use type were derived by overlaying MODIS LAI maps to ICLUS land use maps. Fractions of the impervious surface layer for each land use were derived by overlaying the impervious surface layer and land use layer. All gridded raster data were spatially aggregated to the HUC12 watershed level.
Water Resources Research LI ET AL.
The WaSSI has been extensively validated against ET data across CONUS using MODIS products (Sun, Caldwell, et al., 2011) and USGS measured streamflow data for selected undisturbed watersheds in different climatic zones and land uses (Caldwell et al., 2012). Overall, previous model performance comparison studies indicate that WaSSI is a reliable model and has advantageous over other watershed scale models for regional applications (Caldwell et al., 2015(Caldwell et al., , 2020)). The present study provides additional model validation using data from 717 watersheds located across the United States, the 2006 National Land Cover Database (NLCD), and data of impervious surface fraction from ICLUS V2.1 products, and LAI data products of 2006 (Zhao et al., 2005). Among these 717 watersheds, 608 watersheds represent USGS "reference" watersheds that are not influenced by human activities (e.g., interbasin water transfer, dams), and 109 watersheds are nonreference watersheds that have experienced rapid urbanization (Oudin et al., 2018) and possible hydrologic alterations (e.g., impoundment) found mostly in the Southeast (Wear, 2011). The impervious cover in these 608 "reference" watersheds ranges from 0% to 6.8% and urban land from 0% to 28%. These 109 "nonreference watersheds" had urban area fractions ranging from 0.1 to 1.0 and impervious surface fractions ranging from 0.01 to 0.67.
Because the size of a gauged USGS watershed may be greater (i.e., covering several HUC12 watersheds) or smaller than an HCU12 watershed, the simulation unit of WaSSI, modeled water yield was scaled to the gaged watersheds using an area weighted method. Validation was made for the 717 gaging watersheds using measured monthly streamflow from 1990 to 2009.
The WaSSI model was designed as a noncalibrated model (i.e., no adjustment of model parameters), and modeled water yield was directly compared to monthly and annual streamflow measurements (Caldwell et al., 2012;Sun, Caldwell, et al., 2011). Model performance statistics at both monthly and annual scales included Nash-Sutcliffe Efficiency (NSE; Nash & Sutcliffe, 1970), Coefficient of Determination of Linear Regression Model (R 2 ), and Root Mean Square Error (RMSE). NSE values that are >0.50, >0.65, and >0.75 for prediction of monthly streamflow have been viewed as indicative of satisfactory, good, and very good model performance, respectively (Moriasi et al., 2007).
The U.S. HUC system of watersheds consist of several hierarchy levels (Seaber et al., 1987). The WaSSI model simulations were conducted at the HUC12 level with approximately 88,000 watersheds (size from 0.2 to 9,238 km 2 , Mean ± SD 95 ± 66.7 km 2 ), but were summarized to an HUC8 level with approximately 2,100 watersheds (size from 184 to 22,965 km 2 , Mean ± SD 3,732 ± 2,253 km 2 ), for attribution analyses to determine the key factors controlling water yield responses to urbanization. A few HUC12 watersheds near the coastline with missing land use data or were entirely covered with water were excluded in this analysis.
As a result, a total of 81,900 watersheds were used for final analysis. Future LULC projections suggested that, among the 81,900 HUC12 watersheds, 30%-50% of watersheds were projected to have no changes in urban area for the three future study periods, 2010, 2050, and 2100. So, we focused our analysis on watersheds projected to increase in urban areas over time: 48,368 watersheds for year 2010, 51,640 watersheds for year 2050, and 54,705 watersheds for year 2100.
Based on previous studies (Ahmed et al., 2017;Kumar et al., 2018;Oudin et al., 2018;Weng, 2001), three groups of influential factors that control the water yield response were identified for in-depth attribution analysis. These factors included: (1) historical climatic variables including temperature (TEMP) and precipitation (P), ( 2) LULC characteristics at the baseline (year 2000) expressed as percentage of forest (For00), shrubland (Shru00), grassland (Gras00), cropland (Crop00), water (Wat00), wetland (Wet00), and urban (Urb00) covers, (3) change in LULC expressed as the absolute or relative change of a certain land cove type during 2000-2010, 2000-2050, and 2000-2100 periods. All the previous land cover characteristics and LULC changes were denoted by the first three or four letters of the land cover type with the source year or time period attached. For example, forest fraction in 2000 and its changes from 2000 to 2050 were denoted by For00 and For0050, respectively, and (4) change in impervious surface fraction during 2000-2010, 2000-2050, and 2000-2100, expressed by IMP0010, IMP0050, and IMP00100, respectively. All variables were standardized with a zero mean and standardized deviation of 1.0 for attribution analysis.
To test Hypothesis 1 (H1; i.e., the increase in water yield is caused by an increase in impervious surface, and loss of vegetation and ET), we used Standardized Stepwise regression to explore the relationships between absolute change in water yield (ΔQ) and absolute change in impervious surface, and the role of ET. Prior to stepwise regression analysis, independent variables with significant multilinearity (i.e., Variance Inflation Factor [VIF] >5 or tolerance <0.02) were removed. Because the independent variables were standardized, the stepwise regression coefficients were directly compared for determining the relative influences among the independent variables.
To test Hypothesis #2 (H2; i.e., spatially varied hydrologic responses to urbanization), ordinary Least Squared Regression (OLS), and Geographically Weighted Regression (GWR) were conducted (Li et al., 2017). OLS, as a global linear regression model, assumes spatial stationarity relationships between dependent and independent variables. In contrast, GWR, considered a local regression model, assumes spatially nonstationary relationship across variables and fits a regression model with a focus on neighboring observations around a watershed in this study. We used an adaptive bandwidth by golden section search and Gaussian function weighting methods to improve the goodness of fit of the GWR model with the minimum corrected Akaike Information Criterion (AICc). To evaluate the GWR against the OLS method, we used the same independent variables selected by the standardized stepwise regression model discussed above (Table S1 in the supporting information). The GWR analysis was conducted using the software of GWR 4.0 (National University of Ireland, Ireland and Ritsumeikan University, Japan). We used the F test, a built-in geographical variability test in GWR 4.0 software, to determine whether there is a spatial variable relationship between variables and ΔQ.
The global Moran's I index was adopted to test the spatial autocorrelation of the residuals for both the GWR and the OLS models using GeoDa0.9.5-I (Beta; The Regents of the University of Illinois, Urbana, Illinois, 10.1029/2019WR026574
Water Resources Research USA). Theoretical and algorithm descriptions of GWR method are found in Li et al. (2017). Because the large sample size (54,705 watersheds) exceeded the maximum computing capacity of the GWR software, this study only focused on data analysis at the HUC8 watershed scale that included 2,100 watersheds. The GWR analysis was conducted to demonstrate the advantages of GWR method over OLS in understanding the spatial variability of controlling factors explaining the hydrologic effects of urbanization across the CONUS.
WaSSI model validation results were analyzed for each of the 717 watersheds for a 20-year time period (1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009) at both monthly and annual scale. These watersheds covered a large gradient of climatic regime with annual average precipitation ranging from 226 to 3,019 mm, estimated annual PET from 332 to 1,321 mm, and measured annual streamflow (Q) varying from none to 2,500 mm. The modeled annual Q values significantly correlated with those from USGS measurements at both annual (adjusted R 2 = 0.88, p < 0.05, Figure 1a) and monthly (adjusted R 2 = 0.74, p < 0.05, Figure 1b) scales. Overall, the modeled annual Q values (Mean ± SD 472 ± 283 mm) were 5% higher than measurements (448 ± 342 mm) across the 717 watersheds for the 20-year study period (1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009).
Model performance as quantified by selected evaluation criteria varied greatly across space (Figures 2a-2d).
For example, about 422 or 59% of the 717 watersheds had an adjusted R 2 value higher than 0.8 and 5% or 35 watersheds had adjusted R 2 less than 0.4 at the annual scale (Figure 2a). The watersheds with low R 2 (<0.5) were located in Middle West regions and Texas where measured Q ranged 1-658 mm (mean = 200 mm) and modeled Q ranging 17-988 mm (mean = 261 mm). The NSE varied from negative values mostly in the Middle West regions (about 203 watersheds or 28%) to greater than 0.5 (377 watersheds or 53%) found in other regions (Figure 2c). Overall, 426 watersheds or 59% of the watersheds had NSE > 0.4 at the annual scale while 529 watersheds or 74% of the watersheds had NSE > 0.4 at the monthly scale.
Both NSE and R 2 varied greatly in space, and they did not correlate significantly. However, in general, watersheds that had low NSE values (<0.2) had wider range of R 2 (0.1-0.9) than watersheds having high NSE. For example, watershed that had high NSE (>0.5) had a high R 2 with a narrow range (0.5-0.9). Similarly, RSME (mean = 116 mm) varied greatly corresponding to the spatial pattern of NSE and adjusted R 2 , ranging from 15 to 603 mm at the annual scale. The spatial patterns of monthly-scale adjusted R 2 , NSE, and RSME were similar to those found at the annual scale (Figures 2b and2d).
The urban area and impervious areas increased rapidly from 2000 to 2100 in both relative and absolute terms (Figures 3, S1, S2, S3, and S4). For example, among the 54,705 HUC12 watersheds examined, the mean urban area fraction was 0.17, 0.21, 0.25, and 0.30 for 2000, 2010, 2050, and 2100, respectively. The number of watersheds with urban areas greater than 0.50 increased from 6,066 in 2000 to 7,984 in 2010, 10,398 in 2050, to 13,696 in 2100 (Figure 3a). This represents a relative change in urban area of 195%, 443%, and 870%, for the three periods (2010, 2050, and 2100), respectively (Figure 3b). Similarly, the number of watersheds with impervious surface fraction greater than 0.25 increased from 722 in 2000 to 1,770 in 2100 (Figure 3c), representing a relative increase of 20%, 84%, and 269% for the 2010, 2050, and 2100 time periods, respectively (Figure 3d).
Overall, urban growth from 2010 to 2100 was most apparent in the eastern United States and some western regions such as New Mexico, Arizona, Nevada, and California (Figures S1 andS3). However, the western United States is expected to see higher relative change in urban area and impervious areas than the eastern region (Figures S2 andS4). Urbanization occurred most rapidly in cropland and urban areas that had higher increase rates in impervious surfaces (Figure S5).
The modeled mean annual Q varied from less than 15 mm to over 4,600 mm (Figure 4a) across CONUS in 2000. The CONUS-level future mean absolute change in water yield was 2.8-11.7 mm representing relative change of 1.1-9.5% for the urbanized watersheds (a total of 48,368-54,705 out of 81, 900 HUC12 watersheds)
10.1029/2019WR026574
for 2010, 2050, and 2100 periods (Figures 4b-4d andS6). The mean ΔQ was estimated as 2.8 ± 5.7 mm, 6.2 ± 12.6 mm, and 11.7 ± 22.9 mm for 2010, 2050 and 2100 periods, respectively. For a few watersheds, ΔQ was as high as 254 mm or a 10 folds in relative change from 2000 to 2100 (Figure 4d). Similar to the spatial distribution of urbanization, ΔQ is most obvious in the eastern United States (Figures 4b-4d). However, the relative change in water yield was most pronounced in western United States where baseline water yield was low (Figure S6). Water yield increased by more than 50% in some watersheds in western regions such as New Mexico, Arizona, Nevada, and California (Figure S6).
Overall, watershed water yield increased with the increase in impervious surface (Figures 5, S7, S8, and S9), but not as obvious with urban area fraction (Figure S10). The increase in impervious area explained 80%-85% of the variance of water yield rise. In addition, climate apparently greatly influenced hydrologic responses (Figures 5, S7, S8, and S9). Wetter watersheds (wetness index P/PET≥1) generally had a higher ΔQ response to the increase in imperviopus area and urban area. Drier watersheds (P/PET<1) displayed a more varied response of urban expansion to water yield (Figures 5, S7, and S8). In some extreme cases,
the annual water yield increased by 50 mm with less than a 2% increase (absolute change) in impervious surface fraction (Figures 5 andS7). In other extreme cases, water yield was low (Figures 5 andS7) even under an increase of more than 20% in impervious surface, presumably due to the low precipitation and associated low runoff in these regions. Further analysis showed that ΔQ was also infuenced by the monthly variance of precipitation (Figure S9). Watersheds with higher precipitation variances, generally found in wet regions, had higher ΔQ.
Not surprisingly, the increase in impervious surface area had a negative relationship with the change in evapotranspiration (ΔET), mirroring the relationship between water yield and impervious surface at both the HUC12 and HUC8 levels (Figure S11). It appears that the variability of ΔQ and ΔET becomes larger with the increase in change in impervious area (Figures 5, S7, and S11) reflecting the influences of other factors (e.g., climate and original LULC). The number of watersheds with an annual ΔQ >50 mm increased from 50 in 2010, to 1,046 in 2050 and to 3,747 in 2100. A change in flow of 50 mm represents a great relative change for a large number of watersheds even for many of the "Water Rich" regions such as the coastal plain and piedmont of the Southeast where annual streamflow in forested watershed are often less than 250 mm (Sun, Caldwell, et al., 2011).
In addition to climate, the water yield responses were obviously different among watersheds grouped by land cover type as defined by a single land cover exceeding 50% the total area of a watershed (Figures 6 andS12).
The ΔQ was generally higher in watersheds that were previously dominated by urban land or wetlands (Figures 6 andS12) than other land uses. The relationships between change in impervious surface and ΔQ for urban, forest, wetland watersheds were much tighter, as indicated by a higher R 2 and/or a steeper slope, than most other land cover classes (Figures 7 andS13). The slope of the regression model for forested watersheds was the highest, suggesting a small change in impervious surface would result in a large change in runoff in forested watersheds that were often found under a wet climatic condition (P/PET>1).
Two examples (Figure 8) were provided to further illustrate the watershed water balances under baseline ( 2000) and future urbanization conditions. Both background climate as characterized by wetness index and temporal variance of LULC (forest vs urban) influenced the effects of urbanization. In both cases, ET is a large component, exceeding 50% of precipitation.
Standardized stepwise regression analysis provided further information to determine factors (e.g., magnitude of urbanization, previous land cover types, local climate) that might better explain ΔQ in future periods (Figure S14). For example, ΔQ had significantly positive correlations with change in impervious and precipitation, and baseline coverages of wetland, water, and urban (except 2010). In contrast, ΔQ had significantly negative correlations and change in land cover of forest, wetland, and baseline coverage for shrubland, cropland and forest. The coefficients of standardized stepwise regression models indicated that impervious surface and the precipitation were the most influential factors defining water yield response to urbanization (Figure S14).
We applied GWR to determine the spatial differences in terms of factors that explained the ΔQ at the HUC8 scale. The higher adjusted R 2 and lower AICc, residual sum of squares (SS) and spatial autocorrelations of residuals indicated a better model performance by the GWR than the OLS model (Table S2). The F tests showed that there were andS16). Strong positive correlations were observed between the ΔQ and changes in impervious surface for all time periods and the historical precipitation. Significant negative correlations between changes in wetland for all time periods and forest for 2050 and 2100 and the ΔQ. We also found the pattern of factors affecting ΔQ might be complex across space. For example, there is a significant negative relationship between ΔQ and baseline urban land area in the eastern United States, while insignificant correlations or significant positive correlations were found in the western United States (Figures 9, S15, and S16). In addition, the magnitude of local coefficients determined by GWR differed among influencing variables (Figures 10 andS17). Generally, the coefficients of urbanization represented by change in impervious surface and historical precipitation (P) were found to be the largest, suggesting they are the most important variables in explaining the variations of ΔQ.
In contrast to previous empirical studies on the effects of urbanization on streamflow in the United States (Boggs & Sun, 2011;Oudin et al., 2018;Wang & Hejazi, 2011), the present process-based study represents the first wall-to-wall assessment on the potential hydrologic responses to future urbanization across CONUS. Such a large scale study offers insights to a spectrum of hydrological responses to urbanization and identifies model strength and weakness under various conditions.
Extensive model validation with streamflow measurements at 717 gaging stations that included both references and nonreferenced watershed offered a few insights on large scale hydrologic modeling. First, spatial patterns of the accuracy of the uncalibrated WaSSI model was comparable to other calibrated, physically based models that require more climate and parameter data such as the Variable Infiltration Capacity (VIC) model (Yang et al., 2019). WaSSI model tended to overestimate water yield in the Midwest dry regions in general, but performed better in the wet southeastern United States (precipitation > 1,200 mm; Q > 500 mm) than in dry regions (Q < 500 mm) as judged by R 2 and NSE (Figure 2). Similar to McCabe and Wolock (2011), model bias, when expressed as a percentage of the mean-monthly runoff, can be very large in arid regions where runoff magnitudes are low. The WaSSI modeling results were consistent with findings in VIC for the United States (Yang et al., 2019) and globally (Lin et al., 2019). The relatively poor performance in arid and semiarid Middle West regions by VIC was attributed to both model structural and forcing deficiencies (Yang et al., 2019). Model calibration by adjusting soil parameters (e.g., thickness of soils) affecting infiltration and baseflow slightly improved model performance (Yang et al., 2019).
Similarly, McCabe and Wolock (2011) applied a monthly USGS water balance model across 735 USGS gauges over the conterminous United States, with a similar distribution of correlation coefficient between predicted and measured Q (i.e., median 0.78, 25th percentile 0.61, and 75th percentile 0.87) to that of this study (median 0.83, 25th percentile 0.72, and 75th percentile 0.88 for the referenced watersheds), and a similar spatial pattern of model performance at the annual scale. Other popular models applied to the United States also tended to overestimate runoff in the Great Plains and parts of the Southwest. Performance of SWAT-HUMUS (Arnold et al., 1999), the USGS model (Hay & McCabe, 2002), and the "abcd" model (Martinez & Gupta, 2010) exhibited similar regional patterns. Poor model performance in the west has been primarily attributed to the coarse model spatial resolution relative to precipitation distribution and in the mid-west to inadequate representation of irrigation (Arnold et al., 1999) and a lack of simulation of groundwater exchange processes (Nijssen et al., 1997). In the Northeast, Southeast, eastern Midwest, and Northwest where the WaSSI model performed well, these models also performed well. The monthly NSE reported for the USGS and "abcd" models were generally higher than those for the WaSSI model in the regions where all the models perform well. But that is to be expected due to the extensive calibration process used to parameterize these models, and the precipitation bias correction applied to the weather input data in the case of the USGS water balance model. The performance of the WaSSI model appears to be equal to or slightly better than the "abcd" model performance during the independent evaluation period.
The comparisons above among model performances suggested that human activities such as groundwater withdrawal for crop irrigation and methods of streamflow measurements at the USGS gaging stations might explain most of the modeling errors. In addition, water yield from uplands could be lost to groundwater through river bed recharge in ephemeral streams (McCabe & Wolock, 2011). This process was not considered in WaSSI and is not typically considered in large-scale hydrologic models in general. One generalized hydrological model may not fit all watersheds, even for undisturbed watersheds (i.e., "losing streams"). The WaSSI models were developed using generalized algorithms for ET, soil water routing, and simple treatments of groundwater and subsurface flow at a monthly scale. Similar to VIC and other models mentioned above, such a model structure appeared to work well for humid regions, but further model improvements and soil parameter calibrations are warranted for better describing watershed water balances in the Middle West region (Yang et al., 2019). Fortunately, model deficiencies are not likely to severely affect modeling results of the present study because most of the projected urbanization (Figures S1-S4) and its effects were found in the humid regions (Figure 4).
While ΔQ was small at the CONUS scale, it was as high as 250 mm/year for some watersheds that had previously experienced urbanization (Figures 4 andS6). Indeed, hydrologic effects of urbanization were rather local and were extremely variable across the CONUS in terms of both absolute and relative changes. Several factors emerged to best explain the variability of hydrologic effects of urbanization in the United States.
1. Impervious surface. As expected, water yield responded direct and positively to the increases in impervious surfaces area (Figure 5) and negatively to the reduction of ET over time. The increase in water yield is a direct result of increased partitioning of precipitation to overland flow and a reduction in ET. Both effects were associated with an increase in impervious land and removal of vegetated surfaces. These findings are consistent with previous empirical (e.g., Oudin et al., 2018;Shooshtari et al., 2017) and modeling studies (Anand et al., 2018;Kundu et al., 2017aKundu et al., , 2017b;;Marhaento et al., 2017). It is generally believed that urbanization increases in imperviousness, decreases in green areas, decreases soil infiltration capacity (Price, 2011), reduces ET (Boggs & Sun, 2011;Hao, Sun, Liu, Wan, et al., 2015), and thus elevates stormflow volume (Gwenzi & Nyamadzawo, 2014;Kundu et al., 2017b). Our study suggests that the increase in impervious surfaces and associated hydrologic change will be most pronounced in urban watersheds under future urban sprawl. In other words, existing urban watersheds will become more urbanized in the future and the hydrologic change is expected to be most obvious in these watersheds as impervious surface fraction rates continue to rise. 2. Local climate. Background climate, precipitation in particular, was identified to significantly influence the watershed hydrologic response to urbanization. While the absolute change in water yield in response to urbanization (i.e., increase in impervious surface) was found to be more pronounced in eastern United States (Figures 4 and5) where humid climate, large forest coverage, and high runoff ratio (i.e., Q/P) dominate the landscapes (Petersen et al., 2012), the relative change in water yield was more obvious in western United States. These results were consistent with global experimental studies on the effects of forest vegetation removal on streamflow and ET (Evaristo & McDonnell, 2019;Zhou et al., 2015). Regions with higher precipitation had higher change in direct runoff in response to the increase in impervious surfaces than drier regions. The response of watershed ET to vegetation conversion from forests with higher biomass and deeper roots to grass with lower biomass and shallower roots is not linear to wetness. ET for watersheds with the dryness index (P/PET) being close to unity is most sensitive to land cover change (Zhang et al., 2004). 3. LULC prior to urbanization. Previous LULC turned out to be an important factor explaining the variability of hydrologic response to urbanization. For example, watersheds dominated by forests or wetlands were most sensitive to change in impervious surface among all LULC (Figure 7). Because forest watersheds are located in wet region and forests and wetlands have higher ET, any change to impervious areas (ET reduced to zero) will have higher change in water yield. However, the magnitude of water yield change in a watershed depends on the total change impervious surface. As indicated by Figure S10, urban watersheds generally have higher or more change of impervious surface than forest watersheds as a result of urban sprawls, that is, urban watersheds are becoming more urbanized. Consequently, urban watersheds had the greatest response among all types of watersheds (Figure 6). Previous studies also found that the changes in water yield in urban dominated watersheds seemed to be more sensitive to a greater level of urban expansion than the nonurban dominated watershed (Kumar et al., 2018;Putro et al., 2016;Rouge & Cai, 2014).
Different watersheds have various processes in partitioning precipitation into ET, streamflow, and soil water storage depending on vegetation covers. Forested watersheds with high leaf area, deep roots, and high soil permeability generally have higher ET rates (G. Sun, Domec, et al., 2016) and thus lower water yield than highly urbanized watersheds (Boggs & Sun, 2011;Ekness & Randhir, 2015). Similarly, wetland watersheds have little soil water stress and thus ET are close PET (Sun, Alstad, et al., 2011), and when forests or wetlands are converted to "dry" impervious surfaces or lawns, ET is dramatically reduced (Hao, Sun, Liu, & Qian, 2015). In fact, this study assumes that ET is reduced to zero when all lands in a watershed are converted to impervious surfaces. Thus, change in ET or ΔQ is the highest in watersheds previously having highest ET such as wetlands or forests.
In summary, although the dominant factors controlling hydrologic responses varied across the CONUS and through time, the continued increase in impervious surface especially in previously urbanized areas, and background precipitation patterns contributed most to future ΔQ. Water yields in watersheds that are dominated by forests, wetlands, and urban lands are most responsive to further increase in impervious surfaces, or vulnerable to urban sprawls.
Our study found that increasing impervious surface areas resulted in elevated water yield through increased direct runoff and reduced water loss by ET. This finding is not new, but the spatial variabilities of hydrologic responses across CONUS quantified by this study provide insights about mechanisms of how future urbanization affects watershed hydrology.
The previous conceptual illustration by Livingston and McCarron (1992) has been widely cited in the literature (Arnold & Gibbons, 1996;Paul & Meyer, 2001) and used by U.S. EPA ( 2003) as a guide for stormwater management. However, the reported ET/P ratio of 40% for natural watersheds in those literture was much lower than what we found in the present study as demostrated in Figure 8 and previous studies (Boggs & Sun, 2011;G. Sun, Domec, et al., 2016;Sun, Caldwell, et al., 2011). Simairly, a USGS study on national level ET (Sanford & Selnick, 2012) indicated that ET/P is much higher than 40% in majority of lower 48 states of the United States. Thus, we argue that the "generic approximation" model developed by Livingston and McCarron (1992) might have sustantially underestimated watershed ET rates and the impacts of vegetation removal on stromwater runoff (ΔQ). Our study suggests that the role of vegetation in regulating water cycle (i.e., ET and water yield) in urban watersheds might have been underestimated previously.
Our findings have important implications to watershed management that aims at hydrologic impacts of urbanization. First, maintaining ET, the "biological drainage," is important in controlling urban stormflow (Hao, Sun, Liu, Wan, et al., 2015). Vegetated lands such as forested patches help to reduce frequent flooding risk (Palmer & Montagna, 2015) as well as urban non-point source water pollution (Li et al., 2016;Sun & Lockaby, 2012) due to the high ET rates as well as great water and nutrient cycling capacity of forests. Land use planners that aim at reducing storm runoff in urbanized watersheds should direct resources to urban green infrastructure and low impact development practices to maximize both ET and infiltration rates (Ekness & Randhir, 2015). Second, the hydrologic impacts of urbanization are highly variable in space as a result of climatic differences in the United States. To offset the negative hydrologic impacts of urban intensification across the humid southeastern United States, one of the most vulnerable regions identified by this study, watershed managers may consider practices that increase vegetation coverage, and create, restore, and protect existing wetlands (Sun & Lockaby, 2012). In contrast, planting trees or other greening efforts in dry and water-stressed regions (Gwenzi & Nyamadzawo, 2014) should take caution because city greening might bear high cost including irrigation and may aggravate water scarcity downstream and in groundwater (Lang et al., 2017;Wang et al., 2009). Thus, local and national planning and resource management agencies must consider local watershed and background climate conditions. In addition, the trade-off between runoff reduction and costs borne by landowners for building green infrastructure (Ekness & Randhir, 2015) or food security in populated areas (Bieger et al., 2015) should be considered.
This study integrated projected trends of LULCC, historical climate, vegetation and soil characteristics, and key watershed hydrological processes under a modeling framework. Using a set of consistent databases and a single validated model offered spatial comparisons of the likely range of magnitude of water yield response to urbanization at a middle (2050) and a long-term (2100) time horizons across the CONUS. The GWR model provides insights on the factors affecting hydrologic responses to future projected urbanization in difference regions in the United States. In spite of the advantages of this comprehensive approach, uncertainties exist in model structure and input data, and future studies are needed.
The hydrology of urbanizing watersheds with mixed LULC is complex, and many processes coexist simultaneously. For example, the WaSSI model assumes that the runoff from impervious surfaces goes directly to a 10.1029/2019WR026574
stream without having the opportunity to infiltrate the soil downslope of an area of impervious surface or to be retained in some sort of storm water control structure (e.g., detention ponds). This assumption might result in underestimates of ET by 1-5% (Lull & Sopper, 1969), and thus somewhat overestimate water yield, especially across dry regions of the CONUS. LAI is a major biophysical variable that control ecosystem ET (Sun, Alstad, et al., 2011). However, LAI products for urban lands are rare. This study estimated LAI values for urban core areas using MODIS LAI means of grid cells surrounding urban areas. This approximation might cause an overestimate of LAI for urban lands, thus overestimate ET, resulting in an underestimate of associated impacts on water yield. For future projections of LAI, because the MODIS LAI data set was independent and had a different spatial resolution from the ICLUS data, LAI of urban land could end up higher than previous land cover for nearly 3,000 watersheds. The direct effect was that future areas might have a higher ET and lower water yield during future periods. However, such scenarios (i.e., increase in ET under urbanization) could occur in certain urban areas where trees are planted or a significant amount of irrigation is used to maintain vegetation covers. In addition to vegetation and impervious surfaces, soil properties such as infiltration capacity, porosity, and hydraulic conductivity affect infiltration rates and timing of subsurface flows (Price, 2011). Change in soil properties was not considered in WaSSI, and this model deficiency might have caused underestimation of hydrologic response to urbanization, especially at the monthly scale.
To separate the effects of urbanization from climate change and variability, this study assumed that a static climate represented by a reference period of 1961-2010 would hold for future year 2050 and year 2100. However, climate change impacts, including increases in atmospheric CO 2 concentration, air temperature, and a higher frequency in extreme events, are expected in the 21st century (Wuebbles et al., 2017). These changes will no doubt affect watershed water balances (Martin et al., 2017;Vose et al., 2016) and water use and demand by humans (Sanchez et al., 2018). Thus, the hydrologic effects of urbanization are not likely to occur in isolation but act together with climate change. The effects of climate and urbanization can be additive or offsetting (Kundu et al., 2017a;Putro et al., 2016;Todd et al., 2007). Under multiple future stressors such as land use change, water demand, and climate change, projecting local water resources can be extremely complex and challenging (Sun et al., 2008). We recognize that climate is a major driver of hydrologic response to urbanization; therefore, climate change is essential for future comprehensive realistic assessments of urbanization impacts on water quantity and quality, and other emerging issues such as UHI and Urban Dry Island (Hao et al., 2018;Luo & Lau, 2019) and ecosystem productivity (Li et al., 2020).
We conducted the first of its kind urbanization impact study on watershed water balances at a national scale. We found that spatially varied hydrologic changes were closely associated to urban intensification patterns, LULC, and background climate. The hydrologic response was most pronounced in the southeastern United States, a region with generally higher precipitation amount and variances, forest coverage, and wetlands than in western United States. The increase in water yield was mainly due to the increase in impervious surfaces and decrease in ET associated with vegetation losses.
Our study confirms the hypothesis that "hydrologic impacts of urbanization are not created equal" across both time and space. Our study suggests that cost-effective environmental management measures and strategies must be designed to fit local watershed conditions. To reduce environmental impacts from urbanization, maintaining ecosystem ET capacity or "biological drainage" in urbanizing watersheds through conserving forests and wetlands or developing other "green infrastructure" is important in addition to other methods of minimizing impervious surfaces. Our study results support the idea of "Keeping forest lands as forests" in an urbanizing world to maintain watershed functions and many benefits that they provide to human-dominated urban ecosystems.