Past, current, and projected landscape configurational effects on streamflow within the Rocky River HUC-8 watershed of the Charlotte metropolitan region

Recent studies within the Atlanta, Georgia (GA) metropolitan region (AMR) of the southeastern (SE) United States (US) have indicated that landscape configuration (LC) (pattern) was significantly correlated to watershed streamflow at various scales (Roberts 2016, 2017; Debbage & Shepherd 2018), including the Hydrologic Unit Code (HUC)-12 and HUC-10 (Roberts 2016, 2017). This finding is not surprising, as land cover/land use (LC/LU) change has been found to bring a range of physical changes to hydrologic systems and processes (Bian et al. 2017). LC/LU change-induced LC affects precipitation distribution via direct runoff and subsurface flows that are ultimately evidenced in the watershed runoff record as either discharge increases or decreases (Gao & Yu 2017). In general, LC/LU change alters arrangement and results in fragmentation that has also been found to be a main trend of LC/LU change in many global areas (Gao & Yu 2017). Thus, outside of the AMR, the entire SE US has been one of these global areas subjected to massive LC/LU changes that have induced permanent LC alterations.

Over the past 60 years, the SE US has undergone tremendous LC/LU change resulting in large amounts of natural land loss (Poudyal et al. 2016). Once a major agricultural producer over 100 years ago, the region on average has lost 45% of its agricultural land since 1920 (Ellenburg et al. 2016). Additionally, over the past 25 years, LC/LU change has been attributed to forest harvesting and regrowth, agricultural abandonment, and urbanization (Ellenburg et al. 2016). Over this time, one SE US state that has seen some of the highest LC/LU change rates resulting from larger bands of natural land loss is the State of North Carolina (NC) (Poudyal et al. 2016).

As of 2014, NC was ranked as the ninth most populous US state and experienced the sixth highest net population growth (Mo & Zhang 2016). Although unevenly distributed throughout the state, nearly 66% of the population is concentrated within the Southern Piedmont physiographic province (Mo & Zhang 2016). Thus, within this province is where the largest developed area known as the Charlotte, NC metropolitan region (CMR) is located and this region should be further evaluated for its LC-based runoff impacts.

The core of the CMR is located within Mecklenburg County, NC and has experienced substantial urbanization since 1976 (Pickard et al. 2017a). CMR developed land has increased from 2% in 1976 to 30% of the total land area by 2016 and has some of the densest SE US road networks (Pickard et al. 2017a). The CMR population has also increased 35.8% from 695,454 in 2000 to 919,628 in 2010 (Delmelle et al. 2014). Furthermore, a projected LC/LU change SE US study has found that urbanization will greatly increase over the next 50 years with the largest absolute change occurring in the CMR (Terando et al. 2014). Thus, with CMR available datasets of past, current, and future LC/LU, effects of LC/LU change-induced LC on streamflow could be quantified.

Although previous AMR studies have quantified LC effects on streamflow at the HUC-12 and HUC-10 watersheds scales (Roberts 2016, 2017), other SE US areas of significant development that include the CMR have only been examined recently to quantify these impacts (Debbage & Shepherd 2018). Although this study thoroughly examined these impacts on smaller watersheds via finite time scales, larger watersheds at extended temporal scales have also not been examined to decipher LC-based streamflow impacts. A finding that is important with previous research within the field of hydrology advocating a multi-scale approach in which the impacts of landscape pattern are characterized and compared at various spatial scales with the inclusion of time (Xiao et al. 2016). Furthermore, processes such as precipitation events, runoff, infiltration, and channel flow are known to act at various spatial and temporal scales (Eddy et al. 2017). Thus, spatiotemporal information on landscape pattern is of vital importance to provide insight into prospective relationships between LC and runoff (Xiao et al. 2016).

One watershed spatial scale that should be examined in this context is that of the HUC-8. A HUC-8 watershed is known as a cataloguing unit within the US Geological Survey (USGS)'s hierarchy of drainage area characterizing watersheds at the fourth level and have been utilized within NC for seasonal streamflow forecasts (WRRI 2017). It is also the next largest scale in watershed hierarchy after HUC-10 and HUC-12 watersheds that have shown significant AMR LC streamflow effects (Roberts 2016, 2017). Temporally, past and projected LC/LU should also be examined to decipher previous and future LC runoff impacts. From a projected LC/LU change standpoint, climatic data should also be forecasted to quantify its LC impacts, as climatic changes will have significant impacts on watershed hydrology, including extreme events such as drought (Andrew et al. 2016; Tong et al. 2016). NC HUC-8 watersheds have also been used to quantify climate impacts on streamflow (WRRI 2017). Thus, in this paper, the effects of past and current LC on streamflow in a selected CMR HUC-8 watershed will be addressed. If LC over the past and current time periods is found to have significant streamflow impacts at this scale, projected LC streamflow effects on the selected watershed will also be addressed.

The CMR consists of 11 counties within the State of NC and a HUC-8 watershed representing the northern and eastern urbanized core and suburbs was chosen for analysis (Figure 1(a)). This HUC-8 watershed, known as the Rocky River (RR), was synthesized from 43 watershed catchments based upon the USGS's SPAtially Referenced Regression On Watershed Attributes (SPARROW) SE US's water quality model of the South Atlantic, eastern Gulf of Mexico, and the Tennessee River Basin (SAGT) (Hoos et al. 2008). The 43 catchments were analogous to the USGS's HUC-12 watersheds and based upon the US Environmental Protection Agency's (USEPA) 1:500,000-scale Reach File 1 (RF1), a national dataset of more than 60,000 stream segments (Hoos et al. 2008). Additionally, the 43 HUC-12 watersheds and seven HUC-10 watersheds that are nested within the RR HUC-8 watershed were verified from data obtained from the NC Department of Environmental Quality Online GIS (NCDEQ 2017). The RR HUC-8 watershed had a total drainage area of 3,553 square kilometers (km²).

Using previous methodologies discussed by Roberts (2016, 2017) in AMR watersheds, the newest available CMR LC/LU (the 2011 National Land Cover Dataset (NLCD) (Multi-Resolution Land Characteristics Consortium (MRLC) 2017) was also incorporated to quantify the effects of current LC within the RR HUC-8 watershed. For determination of the past LC impacts on the RR HUC-8 watershed, the two previous time snapshots of the 2001 and 2006 NLCD were used. Consisting of 15 finer LC/LUs captured at the 30-meter (m) resolution, the 2001, 2006, and 2011 NLCD products were reorganized into four coarser LC/LU groups (Table 1), as detailed by Merwade (2012). Merwade's (2012) method allowed for NLCD reorganization into courser categories to create a curve number (CN) grid using the HEC-GeoHMS model, a hydrological simulation also utilized for monthly runoff computations which has similarities to the Thornthwaite Water Balance (TWB) model implemented within this study (Rao et al. 2014). The four regrouped LC/LUs representing the two previous and one current time snapshots were then overlain by the RR HUC-8 and clipped to consist only of watershed area (Figure 1(b)–1(d)).

Due to previous work in the AMR indicating that current LC represented as the contagion (CON) metric was shown to have significant influences on streamflow at the HUC-12 and HUC-10 scales (Roberts 2016, 2017), CON was also utilized to quantify the effects of past and current landscape arrangement on the RR HUC-8 watershed. CON refers to the degree to which mapped LC/LU classes are clumped into patches of the same LC/LU classes (Leitao et al. 2006). The metric takes into account if LC/LU classes of the same type are grouped into larger distributions or fragmented into smaller patches throughout the landscape. The metric takes into account all LC/LUs within the landscape and is normalized, meaning that watersheds of varying size (area) will not affect its calculation (Ritters et al. 1996). However, the distribution or spatial arrangement of LC/LU within the watershed will have a significant impact. CON is defined as 1 minus the sum of the proportional abundance of each LC/LU patch type multiplied by the proportion of adjacencies between cells of that LC/LU patch type and another LC/LU patch type, multiplied by the logarithm of the same quantity, summed over each unique LC/LU adjacency type and each LC/LU patch type, divided by two times the logarithm of the number of LC/LU patch types and multiplied by 100 to convert to a percentage (Leitao et al. 2006). Thus, CON ranges in value from 0 to 100% and for a complete definition of the full equation utilized in its calculation, it is suggested the reader consult Leitao et al. (2006). In regards to the RR HUC-8 watershed for the two previous (2001, 2006) and one current (2011) NLCD time snapshots, CON was tabulated via FragStatsBatch and ArcGIS 9.3 on the regrouped LC/LU classes.

The Thornthwaite Water Balance (TWB) model is widely used to estimate streamflow with data of monthly climate, LC/LU, and soil type to estimate hydrologic inflows, storages, and outflows (Tong et al. 2016). At the HUC-12 and HUC 10 scale within the AMR, the TWB easily implemented CON into variables utilizing LC/LU to aid in the deciphering of impacts of LC on surface runoff and could be extrapolated to the RR HUC-8 watershed (Roberts 2016, 2017). Thus, a USGS TWB model was used for this evaluation and required inputs of only the direct runoff factor (DRF) (%), runoff factor (RF) (%), soil moisture storage capacity (SMSC) (millimeters (mm)), latitude of location (LOL) (degrees (°)); rain temperature threshold (RTT) (degrees Celsius (°C)), snow temperature threshold (STT) (°C), and maximum melt rate (MMR) (%) (McCabe & Markstrom 2007).

The RR HUC-8 CNs for the two previous (2001, 2006) and one current (2011) NLCD time snapshots were tabulated via SHGs downloaded as 10-m resolution rasters from the US Department of Agriculture's (USDA) Natural Resources Conservation Service (NRCS) Geospatial Data Gateway's gridded Soil Survey Geographic Database (gSSURGO) for NC (NRCS 2017). Ranging from A (well drained soils) to D (poorly drained/saturated soils), SHGs were extracted to the RR HUC-8 watershed and combined with the regrouped LC/LUs to create a land/soil file for CN calculation for each past and current NLCD time snapshot using ArcGIS 10.5.1 and the ArcCN ArcScript (Zhan & Huang 2004). Computed CNs are based on CN IIs (Zhan & Huang 2004).

Some may question why barren land was included within the agricultural regrouped class using the methodology outlined by Merwade (2012) (Table 2). Using the NLCD and typical root zone depth values described for LC/LU types and SHGs within Table 9 of the USGS publication entitled A modified Thornthwaite–Mather Soil-Water-Balance code for estimating groundwater recharge (Westenbroek et al. 2010), barren land was found to exhibit CN properties under each SHG closest to hay/pasture and cultivated crops LC/LU types and justified its inclusion within the agricultural regrouped class. The values from this publication were utilized based upon the premise that each LC/LU class listed was derived verbatim from the original NLCD LC/LU classification and was also directly correlated to specific root zone depths based upon each SHG that ranged from A to D. Thus, this led to excellent correlation between all regrouped NLCD classes and CNs within this analysis, even in the case of barren land that was determined to be a small amount of the watershed LC/LU utilized in this study, as 2001, 2006, and 2011 NLCD percentages were found to be 0.10, 0.24, and 0.20%, respectively. Through this excellent correlation between regrouped NLCD classes and CNs, impacts of aggregated LC/LU on original (ORG) and CON-adjusted CNs should be minimized.

Please note that since blended CON-adjusted CN models will have much lower CNs ranging from 0 to 100% of their blended ORG values based upon the multiplied CON value, these simulations will have much higher volume of their estimated streamflow being attributed to baseflow due to the assumption of greater infiltration and permeability within the watershed soil profile. On the other hand, please also note that since blended ORG CN models will typically have higher CNs, these simulations will have much higher volumes of their forecasted runoff being synthesized from overland (direct) runoff due to the assumption of reduced hydraulic conductivity throughout the watershed soil profile. Thus, blended ORG models and blended CON-adjusted models can have very similar overall predicted total streamflow volumes over the entire same time period of their simulations based upon varying proportions of these contrasting components of runoff generation.

The selection of the 3,553 km² RR HUC-8 watershed was based upon having no other drainage area adding to streamflow above it and a robust 15-year (180-month) USGS stream gage period beginning in January 1999 and ending in December 2013. Thus, this gage period corresponded to the time snapshot connected with the CON and CN values tabulated from the two previous (2001, 2006) and one current (2011) NLCD. The selected RR HUC-8 watershed was located north and east of the CMR core within: Anson, Cabarrus, Iredell, Mecklenburg, Rowan, Stanly, and Union counties (Figure 1(a)). Over the time period from 2001 to 2011, NLCD LC/LU shows that the RR HUC-8 watershed consisted predominantly of forest and agricultural cover types (Figure 1(b)–1(d)). The RR HUC-8 watershed is an easterly flowing tributary of the Yadkin/Pee River Basin and experiences a typical SE US humid subtropical climate (Mo & Zhang 2016).

From the previous two (2001, 2006) to one current (2011) NLCD time snapshots, ORG CN and CON values varied slightly from 74.1 to 74.3 and 40.0 to 41.6%, respectively, with mean readings of 74.2 and 40.6%, respectively (Table 2). The CON-adjusted CN also indicated a small range from 29.8 to 30.8 over these three time snapshots (Table 3). The locations and characteristics of the USGS stream gaging site (USGS 2017) and National Weather Service (NWS) meteorological stations (National Oceanic and Atmospheric Administration National Centers for Environmental Information (NOAA NCEI) (NOAA NCEI 2017) incorporated for simulation are given in Figure 1(b)–1(d) and Table 4. Due to the larger drainage area of the RR HUC-8 watershed, seven meteorological stations providing watershed-based temperature (T) in °C and precipitation (P) data in mm were averaged to run the monthly TWB models.

The latitude of location (LOL) was determined to be 35° North (N) for the RR HUC-8 watershed. Additionally, rain temperature threshold (RTT) and snow temperature threshold (STT) were kept constant at 0 °C and the maximum snowmelt rate (MMR) was assumed 100% in all months, as no month during the period having a mean T below freezing. In conclusion, for the RR HUC-8 watershed, a TWB simulation consisted of analyzing the blended ORG CN and SMSC, as compared to the blended CON-adjusted CN and SMSC each month from January 1999 to December 2013. January 1999 to December 2003 corresponded to the past NLCD 2001 time snapshot, whereas January 2004 to December 2008 referenced the past NLCD 2006 time snapshot, and January 2009 to December 2013 represented the current NLCD 2011 snapshot.

As discussed in depth within the initial AMR studies on this topic, competing models were tested for performance using Nash–Sutcliffe efficiency (NSE) (dimensionless), root mean square error (RMSE)-standard deviation ratio (RSR) (dimensionless), percent bias (PBIAS) (mm), and Akaike Information Criteria (AIC) (Roberts 2016, 2017). For more information on how these indices quantify hydrologic model performance, please consult the research by Bian et al. (2017) and Eddy et al. (2017).

In addition to the aforementioned indices, two other indices can be quantified to indicate statistical significance of improvement between competing hydrological models and were utilized here. These indices are the mean squared error (MSE) and the relative correlation coefficient (RCC) (Gupta et al. 2009; Gupta & Kling 2011; Hwang et al. 2012). The MSE is the average of the squares of the forecasting error and is very useful as an unbiased estimate of the variance of the random component (Hwang et al. 2012). MSE values depend on the units of the predicted variable and can range from 0 to infinity (∞) (Gupta et al. 2009). The RCC is a new special measure that evaluates the relative forecasting efficiency or relative improvement of forecasting performance within data-driven hydrologic models (Hwang et al. 2012). In essence, the RCC is a good measure to quantify the performance of estimated datasets by comparing the forecasting accuracy rate in correlation with persistence of the dataset (Hwang et al. 2012). Persistence of the data refers to any random process that indicates a strong autocorrelation structure that increases the forecast accuracy, whereas efficiency in a forecasting model is the accuracy improvement rate in comparison with the naïve model or unchanged condition without other knowledge of the data except the calibration-period data (Hwang et al. 2012). The RCC is measured in units of % and also ranges from 0 to ∞. For more information on how the MSE and RCC measure hydrologic model performance, please consult the research by Gupta et al. (2009), Gupta & Kling (2011), and Hwang et al. (2012).

As conveyed in recent research by Ayele et al. (2017) using model performance rating classifications for several of the indices measured here, both RR HUC-8 TWB models incorporating either the blended CON-adjusted CN or blended ORG CN were found to exhibit a very good (0.65–0.75) NSE, a satisfactory PBIAS (±15≤ to <± 25), and a good (0.50< to ≤0.60) RSR. The RR HUC-8 TWB model run incorporating the blended CON-adjusted CN simulated streamflow over the 15-year (180-month) period slightly more accurately than the blended ORG CN model, as indicated by NSE and RSR values of 0.733 and 0.516 and 0.708 and 0.540, respectively (Table 5). NSE indicated an increase of 0.025 using the blended CON-adjusted CN model, whereas RSR decreased 0.024 when implementing the blended CON-adjusted CN model (Table 5). The PBIAS decreased by over 4 mm when utilizing the blended CON-adjusted CN model, as opposed to the blended ORG CN model (Table 5). Thus, this decrease in PBAIS conveyed that the blended CON-adjusted CN model performed more admirably in regards to the overall overestimation of streamflow by including landscape arrangement as a variable. Although differences in the NSE, RSR, and PBIAS values between the two models may seem slight, choosing the optimal statistically significant model based upon these small incremental differences is a valid approach and has been incorporated in recent hydrological simulation comparison research by Tegegne et al. (2017). In reference to AIC, the value for the blended CON-adjusted CN model was slightly smaller (6.977) than the value (7.067) of the blended ORG CN model and could be seen as the model to select for streamflow analysis at the HUC-8 scale within this watershed (Table 5).

In examining the MSE values for the two competing models, the blended CON-adjusted CN model displayed a much smaller value, by nearly 100 mm, of 1,031.73 mm, as compared to the blended ORG model value of 1,128.46 (Table 5). A smaller MSE value between competing models indicates a better overall simulation that yields a greater unbiased forecasting accuracy (Gupta & Kling 2011; Hwang et al. 2012). In evaluating the RCC values, the blended CON-adjusted CN model was found to have a value of 121.52%, whereas the blended ORG CN model had a tabulated value of 113.42%. For RCC values, the greater the number is over 100%, the more improved the model is as an RCC value of 100 indicates that the model is no better than one in which the estimate is the lag-k measured record and this represents an unchanged situation (Hwang et al. 2012). A RCC value of less than 100% conveys that the model is performing worse than a simulation in which the estimate is the lag-k measured record. Thus, according to this statistically significant index, the blended CON-adjusted was shown to have the best forecasting accuracy.

In regards to the RF for both of these simulations, it was determined that they were almost identical. The blended CON-adjusted CN model for the RR HUC-8 watershed predicted an MRO value of 61% with an MCO value of 39%. The blended ORG CN model indicated MRO and MCO values of 60 and 40%, respectively (Table 5). A graphical depiction of the total comparison between measured streamflow, estimated streamflow incorporating a blended ORG CN, and estimated streamflow using a blended CON-adjusted CN during the 15-year (180-month) time period is given in Figure 2. In conclusion and as discussed in previous studies of the AMR watersheds, larger spreads between measured monthly streamflow and estimated values incorporated the blended ORG CN and blended CON-adjusted CNs witnessed in Figure 2 may be correlated to wastewater treatment discharges, groundwater withdrawals, and other additions or subtractions, such as irrigation not easily measurable within the TWB framework.

After the initial calibration period and based on the hydrographs indicated in Figure 2, it was evidenced that the largest differences between the blended ORG model and blended CON-adjusted model occurred during the initial year (1999) of the evaluated time period. Thus, to compare how the blended ORG model fared against the blended CON-adjusted model over the remaining 14-year (168-month) time period, Figure 3 was synthesized. Measured runoff for the time period of January 2000 to December 2013 in the RR HUC-8 watershed was approximately 8,950 mm and both simulations were found to overestimate runoff by similar volumes. The blended ORG model overestimated runoff by 19.1% with a value of 10,661 mm, whereas the blended CON-adjusted model over-predicted observed streamflow by 18.7% with a value of 10,622 mm. Within both simulations, differences ranged from –151 and –157 mm, respectively, for the blended CON-adjusted model and blended ORG model in March 2003 to +64 and +69 mm, respectively, for the blended CON-adjusted model and blended ORG model in September 2008. For the blended CON-adjusted model, perfect agreement (±0 mm) between estimated and measured runoff occurred on several different dates that included June 2002, July 2002, September 2006, December 2009 and August 2011. Regarding the blended ORG model, perfect agreement (±0 mm) between simulated and observed streamflow was also recorded during various monthly time steps that included June 2002, June 2007, June 2009, November 2010, July 2011 and July 2012. In conclusion, although both simulations were found to overestimate runoff within the RR HUC-8 watershed over the remaining 14-year time period, these results also indicate that the CON-adjusted model was slightly better equipped to project overall total streamflow in the remaining 168-month model comparison snapshot.

Although the blended CON-adjusted CN model outperformed the blended ORG CN model using the hydrologic model performance evaluators over the total 15-year (180-month) period, the examination of LC on streamflow simulation for the two previous (2001, 2006) and one current (2011) NLCD time snapshots should be discussed. With that in mind, each of the time snapshots just mentioned were evaluated separately within each of these five-year (60-month) windows and compared between blended ORG CN and blended CON-adjusted simulated runs (Table 6).

During the initial five-year period of January 1999 to December 2003 corresponding to the past 2001 NLCD time snapshot, both the blended ORG CN model and the blended CON-adjusted model streamflow estimates were the closest to runoff measured at the USGS gage in Norwood, NC. In fact, the blended CON-adjusted model was able to capture 98.3% of the measured runoff over this time period, just underestimating the measured value by 60 mm (Table 6). On the other hand, the blended ORG CN model value was found to overestimate the measured streamflow by over five times (318 mm) the underestimation seen within the blended CON-adjusted CN simulation.

In regards to the middle five-year period of January 2004 to December 2008 correlated to the past 2006 NLCD time shot and the last five-year period of January 2009 to December 2013 representing the current 2011 NLCD time shot, both blended CON-adjusted CN models and blended ORG CN models acted quite similarly. For each of these time snapshots, both types of models overestimated measured streamflow in values that ranged from 727 to 776 mm (Table 6). Furthermore, during the current 2011 NLCD time snapshot representing the period of January 2009 to December 2013, the overestimation was nearly identical for both the blended ORG CN model (727 mm) and blended CON-adjusted CN model (728 mm) (Table 6). The spread in streamflow overestimation for the past 2006 NLCD time snapshot for dates including January 2004 to December 2008 was somewhat higher at 776 mm for the blended ORG CN model compared to 755 mm estimated for the CON-adjusted CN model (Table 6). Although on the larger overestimation side, the January 2009 to December 2013 blended CON-adjusted CN model corresponding to the current 2011 NLCD time snapshot had the highest NSE (0.750) and lowest RSR (0.500), AIC (5.541), and MSE (228.56 mm) values of any of the five-year model runs, but indicated a middle of the pack RCC (127.42%) value (Table 6). However, the largest take away from this comparison is that the least quantities of underestimation and overestimation were achieved during the past 2001 NLCD time snapshot correlated to January 1999 to December 2003 and that this period had the highest measured runoff (3,544 mm) of any of the three five-year time periods by nearly 18%. A finding that points towards LC influences on streamflow simulation being potentially more quantifiable during periods of greater precipitation within the watershed.

As previously mentioned, if LC over the past and current time periods was found to have significant streamflow impacts at RR HUC-8 scale, projected LC streamflow effects on the selected watershed was to be addressed. Fortunately, and within the RR HUC-8 watersheds, data products corresponding to 2030 LC/LU change are available and will be used to find projected impacts. This is important, as one of the best methods to evaluate watershed response to development is to apply a hydrologic model and examine runoff characteristic changes under varying urbanization scenarios (Bian et al. 2017).

Projected LC/LU data in the RR HUC-8 watershed was furnished by the FUTure Urban-Regional Environmental Simulation (FUTURES) model. The model is an open source simulation developed to project regional-scale urban change at the 30-m pixel level and was initially calibrated and validated using several studies in the CMR and the state of NC (Meentemeyer et al. 2013; Dorning et al. 2015; Pickard et al. 2017a, 2017b). With this in mind, projected (Year 2030) FUTURES outputs for the RR HUC-8 watershed that represented five scenarios including: (1) status quo (SQ) or business as usual growth on its current trajectory; (2) high per capita land conversion (HPCLC); (3) low per capital land conversion (LPCLC); (4) sprawl (SPR); and (5) infill (INF), were obtained. Under each scenario, urban LC/LU replaced non-urban LC/LU using the current (2011) NLCD time snapshot and regrouped LC/LU classes. Thus, any regrouped non-urban LC/LU class that was not converted to urban LC/LU within each of the projected (2030) FUTURES scenarios was kept at its current (2011) NLCD class in 2030 within the RR HUC-8 watershed (Figure 4(a)–4(e)).

For the five 2030 FUTURES scenarios, ORG CNs also displayed slight variation ranging from 73.3 to 73.8 and were nearly identical to the CNs tabulated utilizing the past and current NLCD values with a mean of 73.6 (Table 7). On the other hand, CON values increased by 2030 under all FUTURES scenarios, as compared to the past and current NLCD values, with values ranging from 43.0 to 46.8% and a mean value of 45.4 (Table 7). This finding points towards the clumping of the regrouped LC/LU classes, especially as it pertains to urban LC/LU. For the five 2030 FUTURES scenarios, the CON-adjusted CNs also displayed a slight increase over the past and current NLCD values with new values ranging from 31.7 to 34.5 (Table 8). Again, indicating that the proposed doubling of urban LC/LU from the current (2011) NLCD value of 19.2% (Table 2) to a projected mean of 37.3 (Table 7) under all 2030 FUTURES scenarios would cause an increase in direct (overland) runoff.

Finally, utilizing projected meteorological data from January 2028 to December 2032 would give an indication of how climate change could impact the RR HUC-8 watershed over the projected (2030) FUTURES scenarios. Thus, projected climatic data from each of the seven meteorological stations for T (°C) and (P) (mm) were derived from the National Center for Atmospheric Research (NCAR)'s Climate Change Scenario (2017) using mean climate change data for the five-year (60-month) time period above to run the projected TWB model with a blended CON-adjusted CN.

A graphical depiction of the estimated streamflow using a blended CON-adjusted CN during the projected (2030) FUTURES scenario time period of January 2028 to December 3032 (five-year (60-months)) is shown in Figure 5. Furthermore, the projected (2030) FUTURES scenarios time snapshot using a blended CON-adjusted CN model estimates that streamflow at the USGS gage in Norwood, NC over the five-year period of January 2028 to December 2032 will be 3,839 mm (Table 9), an interesting finding as this value is about 2.6% higher in volume than what was estimated for the current (2011) NLCD corresponding to January 2009 to December 2013 using a blended CON-adjusted TWB model.

This research has indicated that a hydrologic simulation incorporating LC may be utilized to more accurately simulate past and current streamflow in a larger urbanizing watershed, as opposed to a hydrologic model devoid of implementing landscape pattern over time, especially over time periods of greater precipitation and runoff. Furthermore, within a calibrated hydrologic simulation incorporating LC, using projected climatic data in conjunction with forecasted LC/LU can potentially provide watershed planners, resource managers, and other stakeholders with a hypothesis of how streamflow may respond due to these estimated changes in the future. In addition, past and current hydrologic simulations incorporating LC may be parsed into smaller time snapshots to allow stakeholders the ability to receive evaluation data and information as to how streamflow has been responding to these conditions and how changes and practices in watershed LC/LU can be better managed in the near-term and future at various scales.

This is because when information is available at a scale commensurate to a management objective and validated against available observed data, a model becomes a tool that aids in the understanding of the underlying hydrological processes and the impacts of human modifications (such as LC), climatic condition, and management options (Eddy et al. 2017). Thus, this research has attempted such a desperately needed endeavor as there is a need for models designed to be applied at various catchment management scales that can balance hydrologic modelling complexity with available observation-based data inputs (Eddy et al. 2017).