Land Cover Change Effects on Stormﬂow Characteristics across Broad Hydroclimate Representative Urban Watersheds in the United States

: Urban development alters stormﬂow characteristics and is associated with increasing ﬂood risks. The long-term evaluation of stormﬂow characteristics that exacerbate ﬂoods, such as peak stormﬂow and time-to-peak stormﬂow at varying levels of urbanization across different hydroclimates, is limited. This study investigated the long-term (1980s to 2010s) effects of increasing urbanization on key stormﬂow characteristics using observed 15 min streamﬂow data across six broad hydroclimate representative urban watersheds in the conterminous United States. The results indicate upward trends in peak stormﬂow and downward trends in time-to-peak stormﬂow at four out of six watersheds. The watershed in the Great Plains region had the largest annual increasing (decreasing) percent change in peak stormﬂow (time-to-peak stormﬂow). With the current change rates, peak stormﬂow in the Great Plains region watershed is expected to increase by 55.4% and have a 2.71 h faster time-to-peak stormﬂow in the next decade. series of The Two in and had decreasing in despite implementation This study provided an approach to monitor quantify the effects of developing areas on key stormﬂow characteristics using the observed streamﬂow data. Similar studies could be extended to identify watersheds with large changes in stormﬂow characteristics and for implementing stormwater management practices to minimize future ﬂood risks and associated


Introduction
The quantification of streamflow responses to rapidly changing land cover is critical, especially in developing urban areas where floods result in devastating consequences to lives and the economy around the world [1 -3]. Urbanization lowers soil infiltration rates and channels more water from precipitation to surface runoff, creating more favorable conditions for flooding and subsequently altering streamflow responses [4,5]. Floods were ranked as the third deadliest hazard worldwide over the past 50 years (1970-2019) in a report by the World Meteorological Organization [1]. With increasing urban developments and an increasing proportion of population exposure to floods [6], close monitoring of the urban growth and hydrological flux that exacerbates floods is important.
Streamflow from a watershed includes baseflow and stormflow. Baseflow is underground water flow through saturated soil or groundwater with a relatively slow velocity, whereas stormflow is typically quick overland flow with higher volume, potentially leading to devastating floods following major precipitation events. Several studies have reported varying magnitudes of streamflow changes depending on different levels of urban development [5,[7][8][9][10][11][12]. The hydrologic impacts of increasing urbanization are known for increasing stormwater runoff volume, increasing peak flows, and decreasing flow times [13][14][15]. A study at three different levels of urban development reported a larger annual stormflow volume in a high urbanized watershed compared to that in a low urbanized watershed [16]. Increases in peak stormflow and decreases in time-to-peak stormflow associated with increasing urbanization have been reported in event scale studies [13,14,17]. However, the rate of change in peak and time-to-peak is unique to each watershed, and the rates of change derived from event scale studies could vary over a long period of time. Long-term interannual evaluations at several watersheds with varying levels of urbanization across different hydroclimate are limited. A long-term observation-based study in the eastern United States (Washington, DC) reported higher annual stormflow discharge from an urban watershed than those from non-urban watersheds [18]. Extension of such studies at several watersheds with varying levels of urbanization provides detailed information, including variations in overall trends and interannual change rates of key stormflow characteristics (peak and time-to-peak) across different hydroclimate regions. The current study tries to explore these with observed streamflow data as a real case study.
Historical land cover data for the United States was scarce until the recent publication of annual land cover maps from the Land Change Monitoring, Assessment, and Projection (LCMAP) [19,20]. The continuous long-term (1985-2020) land cover data from LCMAP provides an opportunity to evaluate the long-term interannual effects of urban development on stormflow characteristics. Additionally, the near-instantaneous (15 min) streamflow data of numerous streamgages by the U.S. Geological Survey (USGS) are unparalleled for monitoring the changes in stormflow characteristics at a finer temporal scale. Taking advantage of these datasets, quantification of changes in stormflow characteristics, specifically those that are associated with urban floods, such as peak stormflow and time-to-peak stormflow, provides important information for understanding the linkage with changes in urban developments and for practical applications, including designing urban stormwater management infrastructures.
With increasing urban developments, the stormflow characteristics from urban watersheds are expected to have increasing peak stormflows and decreasing time-to-peak stormflows [13,14,17]. The expected long-term effects and the extent of these changes, however, may vary with the levels of urbanization of watersheds, implementation of urban stormwater management infrastructures, and hydroclimate. The main aim of the study was to investigate the long-term annual scale effects of urbanization on key stormflow characteristics (peak stormflow and time-to-peak stormflow) by generating a series of unit hydrographs (UHs) and quantifying the changes in UHs (described in Section 2.3) of six representative urban watersheds across six broad hydroclimate regions across the conterminous United States.

Selection of Urban Watersheds
Six representative urban watersheds were selected across six broad hydroclimate regions [21] of the conterminous United States (Figure 1). At each hydroclimate region, a watershed was selected which met the following criteria: (i) largest increase in urban area percent between 1985 and 2019, (ii) long-term (>15 years) near-instantaneous (15 min) streamflow data between 1985 and 2019, (iii) relatively small (<200 km 2 ) watersheds assuming a more uniform rainfall distribution for generating unit hydrographs (UHs), and (iv) less than 1% watershed area as open water to exclude potential large lakes and reservoirs that could alter natural stormflow characteristics. The watersheds meeting these criteria are generally located in the suburbs of the major cities ( Figure 1). The boundaries of the watersheds were delineated by applying digital elevation model 30 m maps [22] in ArcMap version 10.7.1 (Esri Inc., West Redlands, CA, USA). The areas of the watersheds range from 54.5 km 2 to 197.6 km 2 , with an average of 135.8 km 2 (Table 1).

Streamflow and Land Cover Data
For generation of UHs, near-instantaneous (15 min) streamflow data were obtained from the USGS National Water Information System [23]. The selected watersheds have streamflow data starting as early as the late 1980s (Midwest region watershed), although most of them start during or after the 1990s (Table 1). Land cover data were obtained from the LCMAP Collection 1.1 [19,20]. The LCMAP has a suite of 10 science products including the primary land cover with 8 land cover types (Figure 1d) for the conterminous United States. The primary land cover maps were used to generate annual time series of urban (developed) area percent and increase in urban areas during the study period at six watersheds. Although LCMAP data were available since 1985, streamflow data were only available from the late 1990s or early 2000s for most of the watersheds (Table 1). Thus, the start study years and study period (18 to 30 years) varied across study watersheds based on the availability of streamflow data.

Generation of Unit Hydrographs
The unit hydrograph of a watershed is defined as the direct runoff (total streamflow minus baseflow) resulting from a unit (e.g., 1 mm) of excess rainfall generated uniformly over the watershed area during a specified time. The event scale UHs generated from observed data, such as the USGS streamgage data applied in this study, reflect the characteristics of both watershed and rainfall. To reduce the influence of rainfall characteristics, several event scale UHs for a given year were averaged to represent the specific year at each watershed [24]. Thus, the interannual changes in average UH characteristics (peak and time-to-peak) can be attributed to changes in land cover, considering that other watershed physical characteristics (e.g., topography, river/stream network, slope, and soil type) remain unchanged.
The first step for the generation of event scale UHs requires the separation of stormflow from streamflow for those events with reasonably clean rising and falling limbs, i.e., not contaminated by previous or following flow events [25]. Upon visual inspection of such clean storm events, the constant slope method [26] was applied to separate stormflow from streamflow by using the Hydrograph-py version 1.0.1 Python package [27]. The slope recommended by the constant slope method's original authors [26] was too steep for the watersheds in this study. Thus, the slope was reduced to 1.83 L s −1 km −2 d −1 , as suggested and implemented in previous studies [28,29]. To include the stormflows generated from more uniform precipitation over the entire area of watersheds, only stormflow events that lasted for more than four hours were considered. The stormflow (m 3 s −1 ) was divided by excess runoff (total stormflow volume from each event divided by watershed area) to compute the ordinates of UHs in the units of m 3 s −1 mm −1 . For further refinement for reducing the influence of rainfall characteristics, the UHs with peak (m 3 s −1 mm −1 ) and time-to-peak (hour) values above the 90th percentile and below the 10th percentile for each year (January-December) were excluded. The retained event scale UHs for each year were then averaged to obtain an average UH representing that specific year. This process was repeated for all study years at all watersheds to generate several average UHs to quantify the peak and time-to-peak of UHs representing the peak stormflow and time-to-peak stormflow, respectively. The minimum number of stormflow events for a given year was determined using a threshold of mean number of stormflow events for all years minus one standard deviation; years with fewer stormflow events were excluded from the analysis.

Analysis
Changes in streamflow characteristics and land cover were conducted by trend analysis using the ordinary least squares regression during the study years at each watershed. If the urban area (% of watershed area) changed significantly (p-value ≤ 0.05) during the study years, then any changes or trends in stormflow characteristics were considered to be affected by changes in land cover. For practical applications, the quantification of expected increase or decrease in peak stormflow and time-to-peak stormflow due to a unit increase in excess rainfall is useful for designing stormwater management infrastructures. Therefore, the annual percent changes in peak stormflow and time-to-peak stormflow were computed as a ratio of the slopes of regression lines divided by their respective average values during the study years. Based on the slopes of regression lines, the expected increase or decrease in peak stormflow and expected faster or slower time-to-peak stormflow in a decade (10 years) are computed. The supporting data for this study are publicly available as a USGS data release [30].

Changes in Urban Area
The long-term (1985-2019) changes in land cover at all watersheds are shown in

Change in Peak Stormflow and Time-to-Peak Stormflow
Interannual variation in peak stormflow and time-to-peak stormflow at six

Change in Peak Stormflow and Time-to-Peak Stormflow
Interannual variation in peak stormflow and time-to-peak stormflow at six watersheds is presented in Figure 3. Upward trends in peak stormflow were observed at four out of six watersheds (Table 2 and Figure 3). On average, the peak stormflows during the study years were 2.64, 2.21, 1.11, 1.60, 5.74, and 12.73 m 3 s −1 mm −1 at the Northeast, Southeast, Midwest, Great Plains, West, and Pacific Northwest region watersheds, respectively (results not shown). Based on these peak stormflows, the largest increasing percent change (slope of trend line divided by average peak) in peak stormflow is at the Great Plains watershed, which is 5.54% (Table 2). With this increasing percent change in peak stormflow, the Great Plains watershed expects to increase peak stormflow by 55.4% in the next decade. For this watershed, the urban area increase percentage was also the largest among all watersheds ( Table 2). The increasing percent changes of peak stormflow at other watersheds were between 0.11% (1.1% in a decade) at the Pacific Northwest watershed, and 1.10% (11.0% in a decade) at the Southeast watershed. Among six watersheds, two watersheds (Southeast and Great Plains) show significant upward trends in peak stormflow attributed to the significant upward trends in urban areas ( Table 2). Table 2. Trends in peak stormflow, time-to-peak stormflow, and urban (developed) area. Upward or downward trends are based on the positive or negative slopes of trend lines, respectively. The annual percent change of increase (positive values) or decrease (negative values) in peak stormflow (second column) and time-to-peak stormflow (third column) at each watershed is computed as a ratio of slopes of trend lines divided by their respective average values during the study years. The annual percent change in peak stormflow was multiplied by 10 to obtain the expected decadal percent change in peak stormflow (second column). The slopes of trend lines of time-to-peak stormflow were multiplied by 10 to compute expected changes in decadal time-to-peak stormflow (third column). Contrary to what was expected, downward peak stormflow trends were observed at two watersheds (Northeast and Midwest regions) with substantial (>9%) increases in urban areas. However, these decreasing peak stormflow percent changes are relatively small, with −0.04% (−0.4% in a decade) at the Northeast watershed and −0.16% (−1.6% in a decade) at the Midwest watershed ( Table 2). One of the potential reasons for the decreasing peak stormflows at these watersheds could be the effect of urban stormwater management infrastructures. Such infrastructures are usually smaller (e.g., rain gardens, swales, or retention ponds) for small watersheds and were difficult or even impossible (e.g., pervious pavements, cisterns) to detect in the LCMAP land cover maps with a 30 m spatial resolution. The Northeast region watershed has 1177 stormwater best management practices of 29 different types [31]. Similarly, the Midwest region watershed in Kendall County has a stormwater management plan [32]; however, public data are not available for further details. These stormwater management practices were likely implemented after 2004 when the U.S. Environmental Protection Agency required many cities and counties to integrate such practices, especially green stormwater infrastructure, to obtain compliance with their permit [33,34]. Reduction of peak flows due to stormwater management infrastructures has been reported in urban watersheds [35,36], similar to the observations in this study. A few years of data missing for three watersheds (e,g,k) due to limited streamflow data and/or due to limited stormfl events to generate unit hydrographs. Table 2. Trends in peak stormflow, time-to-peak stormflow, and urban (developed) area. Upwa or downward trends are based on the positive or negative slopes of trend lines, respectively. T annual percent change of increase (positive values) or decrease (negative values) in peak stormfl (second column) and time-to-peak stormflow (third column) at each watershed is computed a ratio of slopes of trend lines divided by their respective average values during the study years. T annual percent change in peak stormflow was multiplied by 10 to obtain the expected deca percent change in peak stormflow (second column). The slopes of trend lines of time-to-pe stormflow were multiplied by 10 to compute expected changes in decadal time-to-peak stormfl (third column).  A few years of data are missing for three watersheds (e,g,k) due to limited streamflow data and/or due to limited stormflow events to generate unit hydrographs.

Region
Additionally, the effectiveness of stormwater management may vary based on the types, number, and location/distribution of such infrastructures within a watershed. For example, distributed stormwater infrastructures are more efficient to mitigate runoff volumes and peak flows compared to centralized stormwater management [36]. All six watersheds selected in this study have implemented some forms of stormwater management practices [31,32,[37][38][39][40], started as early as the 1990s in the Pacific Northwest region watershed [40]; however, the effectiveness of such practices might not have outweighed the larger effects of increasing urbanization to reflect downward trends in peak stormflows, as observed in the Northeast and Midwest region watersheds. Although this study is unique in using historical observed data (>15 years) as opposed to modeling studies, it did not explore the effects of stormwater management on key stormflow characteristics directly. Previous studies have investigated the effects of stormwater management on several flow characteristics [35,36]. Further studies would benefit from exploring how the flow characteristics vary among watersheds using detailed temporal and distributed information on the types and extent of stormwater management implementations.
Time-to-peak stormflow had a downward trend at four out of six watersheds (Table 2). On average, the time-to-peak stormflows were 14.44, 7.78, 33.76, 6.67, 6.61, and 2.91 h at the Northeast, Southeast, Midwest, Great Plains, West, and Pacific Northwest region watersheds, respectively. The decreasing percent changes of time-to-peak stormflow are relatively higher (−0.49% to −4.06%) at watersheds with a larger urban area increase (>9%) than watersheds with a relatively smaller urban area increase (<6%) with lower increasing percent changes (0.06% to 0.21%) ( Table 2). As expected, the largest percent change in time-to-peak stormflow was at the Great Plains watershed at −4.1%, which had the largest increase (20.3%) in the urban area. With this decreasing percent change in time-to-peak stormflow, a 2.71 h faster time-to-peak stormflow is expected in the Great Plains watershed in the next decade ( Table 2). The lowest percent change was for the West watershed (0.06%), which had the smallest increase (0.7%) in the urban area. Statistically, a significant downward trend of time-to-peak stormflow was attributed to a significant upward trend in urban areas, only at the Great Plains watershed ( Table 2).

Potential Uncertainties
The results presented are based on annual scale and are not representative of seasonal, monthly, or specific precipitation event-derived changes in stormflow characteristics. The generation of UHs assumes uniform precipitation within a watershed; this assumption may not be valid for all precipitation events that generate stormflows. However, the average UH approach applied is expected to reduce such uncertainties. The observed streamflow data from the USGS are based on the stage-discharge relationship, and this method has an uncertainty of ±6% for ideal conditions [41], and may increase depending on the condition of flow control structure and channel stability [42]. The LCMAP land cover data have a pixel-level overall accuracy above 80% [43], and errors may exist at the pixel level. Potential changes in stormflow characteristics from urban stormwater management infrastructures [7,36] and non-urban areas, such as drainage systems in agricultural croplands [44] or silviculture in forest areas [45], are not quantified. Other static features of watersheds that can affect stormwater characteristics, such as drainage area, length of the longest watercourse, length to center of area, and mainstream slope [46], were considered to have similar effects on two stormflow characteristics during the study years.

Further Research Direction
Only six watersheds were selected due to the difficulty of finding relatively small (<200 km 2 ) watersheds with long-term 15 min streamflow data and increasing urban area. However, similar studies could be further extended to several watersheds with shorter temporal coverage (<15 years considered in this study) with frequent and severe flooding events. This study captured the overall pattern and trends on two stormflow characteristics at an annual scale based on observed data at the watershed outlets. Further studies can explore sub-watershed scale analysis to identify the location and patterns of land cover changes [47] and compute the relative contribution of sub-watershed areas to stormflows. Additionally, finer-scale studies, such as seasonal, monthly, and event scales can provide better information when needed [18]. Furthermore, additional filters to differentiate the stormflows (e.g., high vs. low flows) could be applied to extend the analysis. Several studies have reported increasing flood risks under growing urbanization [48][49][50], which signifies the increasing importance of similar future studies across flood-prone areas. Upon quantification of altered stormflow characteristics using the UH approach, as applied in this study, the related stakeholders could make watershed-specific decisions in the design and implementation of stormwater management plans and infrastructures.

Conclusions
This study investigated the effects of changes in land cover on key stormflow characteristics (peak streamflow and time-to-peak streamflow) of six representative urban watersheds across six broad hydroclimate regions in the conterminous United States. Longterm time series of peak stormflow and time-to-peak stormflow data were generated for the study period varying from 18 to 30 years (1980s to 2010s) at each watershed. The long-term upward trends in peak stormflow and downward trends in time-to-peak stormflow were observed at four out of six watersheds. The watershed in the Great Plains region with the largest percent increase in urban area (20.3%) had the largest increasing percent change in peak stormflow (5.54%) and the largest decreasing percent change in time-to-peak stormflow (−4.06%) among all study watersheds. With this increasing percent change in peak stormflow and decreasing percent change in time-to-peak stormflow, a peak stormflow by 55.4% and 2.71 h faster time-to-peak stormflow is expected in the Great Plains region watershed in the next decade. However, these expected percent changes in the next decade could be affected by legal (urban planning), natural, or physical (e.g., space for expanding urban areas) limits. Two watersheds in the Northeast and Midwest regions had decreasing percent changes in peak stormflow despite increasing urban areas, likely attributed to the implementation of stormwater management practices. This study provided an approach to monitor and quantify the effects of developing urban areas on key stormflow characteristics using the observed streamflow data. Similar studies could be extended to identify watersheds with large changes in stormflow characteristics and for implementing stormwater management practices to minimize future flood risks and associated impacts.