Urban Runo ﬀ Simulation: How Do Land Use / Cover Change Patterning and Geospatial Data Quality Impact Model Outcome?

: With the increase in global urbanization, satellite imagery and other types of geospatial data have been extensively used in urban landscape change research, which includes environmental modeling in order to assess the change impact on urban watersheds. For urban hydrological modeling, as a focus of this study, several related research questions are raised: (1) How sensitive are runo ﬀ simulation to land use and land cover change patterning? (2) How will input data quality impact the simulation outcome? (3) How e ﬀ ective is integrating and synthesizing various forms of geospatial data for runo ﬀ modeling? These issues were not fully or adequately addressed in previous related studies. With the aim of answering these questions as research objectives, we conducted a spatial land use and land cover (LULC) change analysis and an urban runo ﬀ simulation in the Blue River watershed in the Kansas City metropolitan area between 2003 and 2017. In this study, approaches were developed to incorporate the Hydrologic Engineering Center Hydrologic Modeling System (HEC-HMS) model with remote sensing, geographic information systems (GIS), and radar rainfall data. The impact of data quality on the model simulation outcome was also analyzed. The results indicate that there are no signiﬁcant di ﬀ erences between simulated runo ﬀ responses in the two study years (2003 and 2017) due to spatial and temporal heterogeneity of urbanization processes in the region. While the metropolitan area has been experiencing remarkable urban development in the past few decades, the gain in built-up land in the study watershed during the study period is insigniﬁcant. On the other hand, the gain in vegetated land caused by forestation activities is o ﬀ set by a decrease in farmland and grassland. The results show that increasing spatial data resolution does not necessarily or noticeably improve the HEC-HMS model performance or outcomes. Under these conditions, using Next Generation Weather Radar (NEXRAD) rainfall data in the simulation provides a satisfactory ﬁt in hydrographs’ shapes, peak discharge amounts and time after calibration e ﬀ orts, while they may overestimate the amount of rainfall as compared with gauge data. This study shows that the developed approach of synthesizing satellite, GIS, and radar rainfall data in hydrological modeling is e ﬀ ective and useful for incorporating urban landscape and precipitation change data in dynamic ﬂood risk assessment at a watershed level.


Introduction
Flooding is one of the most severe and damaging natural hazards. As revealed by the Emergency Events Database (EM-DAT) Center for Research on the Epidemiology of Disasters (CRED) [1], from 2005 to 2014, floods accounted for 46% of all natural disasters and affected about 85 million people. During the 20th century, floods were the number-one natural disaster in the United States in terms of and effectiveness of integrating various forms of geospatial data in a unified runoff simulation. Further, our review indicates that the impact of the quality of geospatial data, specifically different spatial resolutions, on the model outcome has not been fully examined in previous studies, particularly in an urban watershed such as our case study area.
The response of LULC to precipitation varies spatially; thus, spatially distributed rainfall data, such as Next Generation Weather Radar (NEXRAD), are critical in distributed hydrological models, such as HEC-HMS for accurate runoff computation, which requires calculating the mean areal precipitation (MAP) for the watershed. Using this dataset, computing MAP explicitly considers the spatial variability of rainfall compared with ground-based gauge rainfall [24]. Some previous studies have tested the performance of different NEXRAD precipitation products against the ground-based gauge data in hydrological modeling, and results have shown the superiority of NEXRAD data because of their ability to capture the rainfall spatial variations for better outcomes [25][26][27]. Only a few studies have attempted to couple NEXRAD Level 3 rainfall with the HEC-HMS model for the assessment of LULC change impact on urban flooding. For example, Knebl et al. [17] used NEXRAD data in HEC-HMS simulations and found that the model tended to overestimate the runoff, and the calibration runs improved the overall results. On the other hand, McCormick [28] found that the NEXRAD data were effective for utilization with HEC-HMS, and the model produced reasonable results. Given these different findings, a further understanding of this issue is essential.
In summary, all the issues described above raised the research questions: (1) How sensitive are runoff simulation to land use and land cover change patterning? (2) How will input data quality impact the simulation outcome? (3) How effective are integrating and synthesizing various forms of geospatial data for runoff modeling? To address these questions, this study aims to (1) detect LULC changes and assess the impact of the change patterns on the watershed hydrological response (runoff); (2) examine the impact of data quality on simulation outcomes; (3) test the suitability of NEXRAD rainfall data when being used in integration with other geospatial datasets in runoff modeling.

Case Study Area: The Blue River Watershed
The case study area ( Figure 1) selected for the simulation is the Blue River watershed. It is a fifth-order stream basin that covers about one-half of the Kansas City metropolitan area south of the Missouri River. It extends about 686 square kilometers through two states (Missouri and Kansas), four counties (Johnson and Wyandotte in Kansas; Jackson and Cass in Missouri), and 11 municipalities [29], representing diverse LULC types. Most of the watershed consists of silt loams and silty clay loams, and more than a third of the watershed is under the urban land complex category. Approximately 60 to 70% of this category are impervious surfaces, which also contain silt loams and silty clay loam soils. These types of soil are associated with moderate to high runoff rates according to the hydrologic soil groups of the Natural Resources Conservation Service (NRCS) [30]. Generally, the Kansas City Metropolitan area has experienced remarkable urban development in the past few decades, which resulted in a significant increase in impervious surfaces (built-up land). However, this changing trend has varied spatially and temporally. As such, we selected this watershed for our case study, which has more diverse land use and development activities and a history of frequent floods, in order to better examine the LULC impact on the model outcome.

Methodology
Runoff-rainfall modeling requires a variety of spatial data inputs related to, for example, land use and land cover, terrain, soil, precipitation, and streamflow data records. These datasets are used as inputs for the HEC-HMS model. The method in this study can be described in several major stages ( Figure 2). It first derives the LULC maps with 6, 20, and 30 m spatial resolution using the Maximum Likelihood classification (MLC) and performs a change detection analysis, which is followed by delineating the drainage networks with 3 and 30 m spatial resolution DEMs. The process then prepares the soil data in order to derive the Soil Conservation Service (SCS) curve number (CN) values and then prepares NEXRAD rainfall data. The final step includes the model setup and simulation run.

Methodology
Runoff-rainfall modeling requires a variety of spatial data inputs related to, for example, land use and land cover, terrain, soil, precipitation, and streamflow data records. These datasets are used as inputs for the HEC-HMS model. The method in this study can be described in several major stages ( Figure 2). It first derives the LULC maps with 6, 20, and 30 m spatial resolution using the Maximum Likelihood classification (MLC) and performs a change detection analysis, which is followed by delineating the drainage networks with 3 and 30 m spatial resolution DEMs. The process then prepares the soil data in order to derive the Soil Conservation Service (SCS) curve number (CN) values and then prepares NEXRAD rainfall data. The final step includes the model setup and simulation run.
Water 2020, 12, x FOR PEER REVIEW 4 of 19 Figure 1. Location of the study area, the Blue River watershed, with major sub-watersheds and the locations of rain gauges, USGS streamflow gauges, and the Kansas City Radar rain station.

Methodology
Runoff-rainfall modeling requires a variety of spatial data inputs related to, for example, land use and land cover, terrain, soil, precipitation, and streamflow data records. These datasets are used as inputs for the HEC-HMS model. The method in this study can be described in several major stages ( Figure 2). It first derives the LULC maps with 6, 20, and 30 m spatial resolution using the Maximum Likelihood classification (MLC) and performs a change detection analysis, which is followed by delineating the drainage networks with 3 and 30 m spatial resolution DEMs. The process then prepares the soil data in order to derive the Soil Conservation Service (SCS) curve number (CN) values and then prepares NEXRAD rainfall data. The final step includes the model setup and simulation run.    (Table 1). ERDAS IMAGINE 2016, an image processing software of HEXAGON Geospatial that performs advanced remote sensing analysis and spatial modeling, was used for image preprocessing, classification, and accuracy assessment in this study. A satellite image preprocessing technique, such as geometric correction, was applied to all four images. Using the SPOT 2017 image as a georeferenced image, an image-to-image registration was conducted on the SPOT 2003 image and LANDSAT 2003 and 2017 images, resulting in an accuracy of below 0.5 root-mean-square error for each image. The selection of dates was based on the availability of images. The MLC classification was applied to derive land use and land cover maps, one of the required inputs for HEC-HMS. MLC is the most common classification used with remote sensing data [31]. It is based on the probability that a pixel belongs to a particular class and considers the variability of classes using a covariance matrix [32]. Spectral signature samples were carefully collected to ensure proper representation of the class spectral reflectance to obtain reasonable estimates of the conditional mean vector and covariance matrix. This process enables more accurate classification results. The error matrix method, the most effective and widely accepted measure to evaluate the accuracy of thematic maps, was used to evaluate the accuracy of the classification results. It compares information from reference sites with information on the classified map for a number of sample areas. Kappa is another measure of image classification accuracy that is based on the difference between the actual agreement in the error matrix (the agreement between the classification and the reference data) and the agreement expected by chance [33]. In this process, 250 reference points were randomly assigned to each of the four classified images. Some other available images from previous studies and Google Earth historical images were utilized for further verification of the reference points. The classification scheme was determined to better reflect major land use/cover types in the study area and also represent the common types of land use/cover used for the SCS CN method in modeling the runoff. They include water, forestland, farmland/grassland, and built-up areas. The water class includes rivers, lakes, and ponds. The built-up area class contains all impervious surfaces, such as buildings, driveways, paved parking lots, etc. The farmland/grassland category represents farms, grass, and bare soil land, whereas the forestland class represents woods and trees. For change detection analysis, post-classification change detection was performed on classified SPOT and LANDSAT images using ERDAS IMAGINE 2018. Resampling was conducted on the SPOT 2003 (20 m) image to obtain a uniform 6 m cell size for both images. There was no need to resample LANDSAT images because both have the same spatial resolution (30 m).

Watershed Delineation
For the purpose of comparison, USGS DEMs with two resolutions were utilized to delineate the watershed drainage network in order to investigate the impact of different spatial resolution data inputs on the simulation results. The first is LiDAR-derived DEM data with a 1/9 arc-second (3 m) horizontal resolution with a vertical accuracy of 0.87 m, and the second dataset has one arc-second (30 m) horizontal resolution with a vertical accuracy of about 3 m. The HEC-GeoHMS toolkit was used to delineate the drainage network with both DEMs. HEC-GeoHMS is a geospatial hydrologic model extension for ArcMap released by the US Army Corps of Engineers [34]. The tool allows the modelers to analyze digital terrain data, derive stream and watershed drainage networks, and construct and prepare inputs for hydrological model software, such as HEC-HMS [35]. Using this tool, terrain preprocessing was performed in multiple steps, starting with filling the sinks and ending with delineating the catchment. A threshold value of 0.25% was chosen to define the streams in both DEMs. Terrain preprocessing outputs are composed of both raster and vector data, which are the primary inputs for the HEC-GeoHMS model setup. By defining the basin outlet point, the tool created the watershed boundaries and the stream networks and computed parameters. Two drainage watershed systems were extracted from both DEMs.

NEXRAD Rainfall Data Processing
For model calibration, hourly NEXRAD Level III rainfall data, obtained from the National Weather Service (NWS) Weather Surveillance Radar Doppler units (WSR-88D) [36], were utilized. The accuracy of NEXRAD precipitation data is based on the Z-R relationship model that is used to estimate rainfall from reflectivity. An automated GIS-based approach was applied to prepare and process the data. A multi-step procedure was implemented utilizing NCEL Radar Software, ArcGIS 10.7, and the standard UNIX Command. The decoding step produced about 2700 radar images for the selected events. Precipitation grids are required to be in the Data Storage System format (HEC-DSS) for use in the HEC-HMS model. HEC-DSS is a dataset system designed to effectively store and retrieve sequential data, such as textual or gridded time-series data [37]. Figure 3 shows some selected rainfall images for the 2017 flood event of the one-hour NEXRAD precipitation data after processing.
Water 2020, 12, x FOR PEER REVIEW 6 of 19 model extension for ArcMap released by the US Army Corps of Engineers [34]. The tool allows the modelers to analyze digital terrain data, derive stream and watershed drainage networks, and construct and prepare inputs for hydrological model software, such as HEC-HMS [35]. Using this tool, terrain preprocessing was performed in multiple steps, starting with filling the sinks and ending with delineating the catchment. A threshold value of 0.25% was chosen to define the streams in both DEMs. Terrain preprocessing outputs are composed of both raster and vector data, which are the primary inputs for the HEC-GeoHMS model setup. By defining the basin outlet point, the tool created the watershed boundaries and the stream networks and computed parameters. Two drainage watershed systems were extracted from both DEMs.

NEXRAD Rainfall Data Processing
For model calibration, hourly NEXRAD Level III rainfall data, obtained from the National Weather Service (NWS) Weather Surveillance Radar Doppler units (WSR-88D) [36], were utilized. The accuracy of NEXRAD precipitation data is based on the Z-R relationship model that is used to estimate rainfall from reflectivity. An automated GIS-based approach was applied to prepare and process the data. A multi-step procedure was implemented utilizing NCEL Radar Software, ArcGIS 10.7, and the standard UNIX Command. The decoding step produced about 2700 radar images for the selected events. Precipitation grids are required to be in the Data Storage System format (HEC-DSS) for use in the HEC-HMS model. HEC-DSS is a dataset system designed to effectively store and retrieve sequential data, such as textual or gridded time-series data [37]. Figure 3 shows some selected rainfall images for the 2017 flood event of the one-hour NEXRAD precipitation data after processing.

Urban Runoff Model: HEC-HMS
The watershed model prepared by HEC-GeoHMS is shown in Figure 4. The watershed is divided into the four major sub-watersheds to accurately calculate and simulate the runoff. For model calibration and validation, we used three flooding events, which were reported as flash flood events

Urban Runoff Model: HEC-HMS
The watershed model prepared by HEC-GeoHMS is shown in Figure 4. The watershed is divided into the four major sub-watersheds to accurately calculate and simulate the runoff. For model calibration and validation, we used three flooding events, which were reported as flash flood events according to NOAA Storm Events Database [38]. The selected events for model calibration occurred in the periods of 2-6 June 2005 and 21-23 August 2017. The 2008 event was selected for model validation, which occurred from 29 to 31 July 2008 ( Table 2). The available hourly and daily gauge data from NOAA/National Climatic Center [39] were used to evaluate NEXRAD Level III precipitation. Observed streamflow data from two USGS discharge gauges were available for all events. The study simulated the runoff hydrographs using HEC-HMS. For the accuracy of the runoff simulation, four separate models were used in the simulation to represent each element of the runoff process, including runoff volume, direct runoff, channel flow, and baseflow. The SCS CN loss method was utilized to compute the runoff volume. It estimates precipitation excess according to cumulative precipitation, soil type, land use/cover types, and antecedent moisture [24]. Soils were classified into A, B, C, and D categories according to their minimum infiltration rates, based on the hydrological soil group (HSG) (more details in National Engineering Handbook and TR-55 report [40,41]). The Soil Survey Geographic (SSURGO) data, available from the United States Department of Agriculture (USDA) [42], were prepared based on both the SSURGO and GeoHMS user's manual instructions ( Figure 4). Using the four SPOT and LANDSAT imagery-generated LULC maps and the soil data, the CN value for each cell in the elevation grid and the average CN value for each sub-watershed were calculated.  Table 2). The available hourly and daily gauge data from NOAA/National Climatic Center [39] were used to evaluate NEXRAD Level III precipitation. Observed streamflow data from two USGS discharge gauges were available for all events. The study simulated the runoff hydrographs using HEC-HMS. For the accuracy of the runoff simulation, four separate models were used in the simulation to represent each element of the runoff process, including runoff volume, direct runoff, channel flow, and baseflow. The SCS CN loss method was utilized to compute the runoff volume. It estimates precipitation excess according to cumulative precipitation, soil type, land use/cover types, and antecedent moisture [24]. Soils were classified into A, B, C, and D categories according to their minimum infiltration rates, based on the hydrological soil group (HSG) (more details in National Engineering Handbook and TR-55 report [40,41]). The Soil Survey Geographic (SSURGO) data, available from the United States Department of Agriculture (USDA) [42], were prepared based on both the SSURGO and GeoHMS user's manual instructions ( Figure 4). Using the four SPOT and LANDSAT imagery-generated LULC maps and the soil data, the CN value for each cell in the elevation grid and the average CN value for each sub-watershed were calculated.  The ModClark model was selected to simulate the direct runoff of the excess precipitation in the watershed to be utilized with NEXRAD precipitation data. The ModClark algorithm is the modified version of the Clark unit hydrograph that is suitable for the use of spatially distributed precipitation data ( Figure 5). This distributed parameter model accounts explicitly for variation in runoff travel time and storage. In this method, a grid is overlaid on the watershed, which allows for calculating the distances from all regions of the watershed to the outlet to compute the inflow and outflow. The model then combines them to determine the direct runoff hydrograph [43,44].  The ModClark model was selected to simulate the direct runoff of the excess precipitation in the watershed to be utilized with NEXRAD precipitation data. The ModClark algorithm is the modified version of the Clark unit hydrograph that is suitable for the use of spatially distributed precipitation data ( Figure 5). This distributed parameter model accounts explicitly for variation in runoff travel time and storage. In this method, a grid is overlaid on the watershed, which allows for calculating the distances from all regions of the watershed to the outlet to compute the inflow and outflow. The model then combines them to determine the direct runoff hydrograph [43,44]. The application of the model requires Storage coefficient R and Time of concentration tc. R represents the temporary storage rainfall excess while the drainage is processing to the outlet. It can be calculated by dividing the flow at the inflection point on the hydrograph's falling limb by the time derived from the flow. Alternatively, it can be estimated as: Time of concentration tc was estimated according to the SCS TR-55 method using HEC-GeoHMS for each cell of the model and is derived as: where tc is the time of concentration for the subwatershed and a function of basin length and slope, dcell the travel distance from the cell to the outlet, and dmax the travel distance from the cell furthest from the outlet. For routing (channel flow) modeling, the Lag method was applied. It is widely used, mainly in urban drainage channels, and is calculated as: where Ot = outflow hydrograph ordinate at time t; It = inflow hydrograph ordinate at time t; lag = time by which the inflow ordinates are to be lagged. As with other parameters, it can be estimated with the availability of observed streamflow hydrographs as the elapsed time between the time of hydrograph peaks. Baseflow was modeled using the Constant Monthly Baseflow method that represents the baseflow as a constant flow, which may vary monthly. Monthly baseflow values from the available USGS streamflow gauges were calculated for each flood event [24].

Model Performance
The model performance was evaluated by assigning the following indicators: Nash-Sutcliffe Efficiency (E) E indicates how well the result of observed flow versus the simulated flow fits the 1:1 line. The value of E ranges between 1.0 (perfect fit) and −∞. An efficiency value below zero means that the value of the observed data would have been a better predictor than the model [45]. The application of the model requires Storage coefficient R and Time of concentration t c . R represents the temporary storage rainfall excess while the drainage is processing to the outlet. It can be calculated by dividing the flow at the inflection point on the hydrograph's falling limb by the time derived from the flow. Alternatively, it can be estimated as: Time of concentration t c was estimated according to the SCS TR-55 method using HEC-GeoHMS for each cell of the model and is derived as: where t c is the time of concentration for the subwatershed and a function of basin length and slope, d cell the travel distance from the cell to the outlet, and d max the travel distance from the cell furthest from the outlet. For routing (channel flow) modeling, the Lag method was applied. It is widely used, mainly in urban drainage channels, and is calculated as: where O t = outflow hydrograph ordinate at time t; I t = inflow hydrograph ordinate at time t; lag = time by which the inflow ordinates are to be lagged. As with other parameters, it can be estimated with the availability of observed streamflow hydrographs as the elapsed time between the time of hydrograph peaks. Baseflow was modeled using the Constant Monthly Baseflow method that represents the baseflow as a constant flow, which may vary monthly. Monthly baseflow values from the available USGS streamflow gauges were calculated for each flood event [24].

Model Performance
The model performance was evaluated by assigning the following indicators: Nash-Sutcliffe Efficiency (E) E indicates how well the result of observed flow versus the simulated flow fits the 1:1 line. The value of E ranges between 1.0 (perfect fit) and −∞. An efficiency value below zero means that the value of the observed data would have been a better predictor than the model [45].

Percent Bias (EPBIAS)
EPBIAS measures the average tendency of the simulated data to be larger or smaller than their observed counterparts. A low value indicates a more accurate simulation. Positive results mean underestimation, whereas negative values mean overestimation [46].
The Root-Mean-Square Error Standard Deviation (RMSE Std Dev) The resulting value of the RMSE Std Dev varies from 0, as an optimal fit, to a larger positive value. The lower the value, the better the simulation performance [46]. EPBIAS measures the average tendency of the simulated data to be larger or smaller than their observed counterparts. A low value indicates a more accurate simulation. Positive results mean underestimation, whereas negative values mean overestimation [46].

LULC Maps, Classification Accuracy, and Change Detection
The Root-Mean-Square Error Standard Deviation (RMSE Std Dev) The resulting value of the RMSE Std Dev varies from 0, as an optimal fit, to a larger positive value. The lower the value, the better the simulation performance [46].   (Table 5). As a result, during the study period, farmland and grassland remained the largest land use type, followed by forestland, built-up area, and water bodies. satisfactory results for the purpose of this study, in which LULC maps with different resolutions are produced as inputs for the runoff modeling. SPOT and LANDSAT imagery classification results show an increase in water bodies and forestland and a decrease in farmland/grassland. There is a marginal gain in built-up area ( Table 5). As a result, during the study period, farmland and grassland remained the largest land use type, followed by forestland, built-up area, and water bodies. Change detection statistics were calculated and are shown in Tables 6 and 7. Both classifications reveal an increase in the water class (0.19% from SPOT imagery and 0.15% from LANDSAT imagery). Forestland increased by a higher percentage (11% from SPOT imagery and 4% from LANDSAT imagery). These change patterns might have been caused by efforts to restore the Blue River watershed's ecosystem, such as a partnership project, called Renew the Blue, that restores riparian forests, upland habitats, and wetlands in the study watershed [47]. The increase in forestland through forestation activities was mostly at the cost of losing farmland and grassland, according to image classification results. A small decrease in built-up areas might also be associated with the gain in forestland. The built-up area class shows, in SPOT classification results, a decrease of 0.02%, which can be considered as no change, while LANDSAT classification results indicate that the decrease in built-up areas was about 3.67%. These change patterns are consistent with the similar findings of Ji et al.'s [48] study, which found that, between 1992 and 2010, forestland increased by about 4%; farmland and grassland decreased by about 16%; built-up areas increased by about 10%. Tables 6 and 7 represent the total land cover transformation between 2003 and 2017. For instance, according to SPOT imagery classification results, forestland increased at the expense of farmland/grassland (89.5 km 2 ) and built-up areas (22.25 km 2 ). The total loss of farmland and grassland was about 167.61 km 2 ; 89.5 km 2 was converted to forest or trees, and the rest (76.15 km 2 ) was converted to built-up areas. In general, forestland has increased the most, whereas the farmland/grassland class has decreased the most. According to LANDSAT imagery classification results, the farmland/grassland class was the larger contributor to the increase in forestland, and both forestland and built-up area classes contributed to the loss of farmland/grassland. The farmland/grassland class was the main land use type in both SPOT and LANDSAT results after 2003. In 2003, the built-up area was the second-largest land use; however, by 2017, forestland became the second-largest land cover in the watershed. Both results show an increase in water bodies, and in general, all classes contributed to this increase. This change detection analysis is essential for understanding the effect of past and current LULC regimes and the effect of LULC change that has been taking place in the watershed and its response to flood events. The study of Ji et al. [48] indicated that all watersheds in the entire metropolitan area gained about 7% in built-up areas from 1992 to 2010. The above change statistics suggest that our study area (Blue River watershed) had a relatively low rate of urban development and more complex vegetation change patterning during the study period from 2003 to 2017. Naturally, further examination of whether the runoff simulation outcome can reflect such LULC change patterning was carried out, as described in Section 3.5.

Watershed Extraction and CN Values
Total stream length and watershed catchment area are often used to measure the level of details extracted from DEMs [20,[49][50][51]. The watershed drainage networks extracted from 3 and 30 m DEMs are shown in Figure 7. The 3 m DEM represents a total stream length of 230 km and a catchment area of 685.56 km 2 , whereas the 30 m DEM represents a total stream length of 178.60 km and a catchment area of 685.59 km 2 . The 3 m DEM is sensitive to minor topographic variations, so it captures more details of topographic features. With a vertical resolution of less than 1 m (0.87 m), the 3 m DEM modeled about 51.4 km more streams. Even though the 3 m DEM modeled more stream length, the difference between both results can be considered insignificant. The 3 m DEM produced a more detailed stream network and slightly lower values of basin parameters compared with the 30 m DEM; however, the difference between both DEM results is not significant. Additionally, both DEMs produced the same watershed total area. This outcome agrees with the findings of other studies that have compared the impact of different DEM resolutions on hydrologic and hydraulic modeling results [51][52][53]. Thus, for hydrologic modeling, the use of a moderate-resolution DEM provides reasonable results. Estimated CN values using SPOT and LANDSAT LULC maps are shown in a color range from 30 to 100 in Figure 8. Utilizing LULC maps with different spatial resolutions to estimate the CN values produced different results. LANDSAT LULC maps produced slightly higher values. The watershed, in general, has a high CN value even though the built-up area occupies just about 30-35%, and much of the area is vegetated land, which is occupied by farmland, grassland, and forestland. The reason for this is that the majority of the soil covering the watershed is from categories C and D, based on HSG, which usually have low infiltration and high potential runoff.

NEXRAD Level III Validation
NEXRAD data were validated against gauge rainfall records from three ground-based gauges that were only available for Upper Blue River and Indiana Creek sub-watersheds. Scatter plots in Figure 9 show a significant correlation between NEXRAD and gauge precipitation data for 2008 (R 2 Estimated CN values using SPOT and LANDSAT LULC maps are shown in a color range from 30 to 100 in Figure 8. Utilizing LULC maps with different spatial resolutions to estimate the CN values produced different results. LANDSAT LULC maps produced slightly higher values. The watershed, in general, has a high CN value even though the built-up area occupies just about 30-35%, and much of the area is vegetated land, which is occupied by farmland, grassland, and forestland. The reason for this is that the majority of the soil covering the watershed is from categories C and D, based on HSG, which usually have low infiltration and high potential runoff. Estimated CN values using SPOT and LANDSAT LULC maps are shown in a color range from 30 to 100 in Figure 8. Utilizing LULC maps with different spatial resolutions to estimate the CN values produced different results. LANDSAT LULC maps produced slightly higher values. The watershed, in general, has a high CN value even though the built-up area occupies just about 30-35%, and much of the area is vegetated land, which is occupied by farmland, grassland, and forestland. The reason for this is that the majority of the soil covering the watershed is from categories C and D, based on HSG, which usually have low infiltration and high potential runoff.

NEXRAD Level III Validation
NEXRAD data were validated against gauge rainfall records from three ground-based gauges that were only available for Upper Blue River and Indiana Creek sub-watersheds. Scatter plots in Figure 9 show a significant correlation between NEXRAD and gauge precipitation data for 2008 (R 2

NEXRAD Level III Validation
NEXRAD data were validated against gauge rainfall records from three ground-based gauges that were only available for Upper Blue River and Indiana Creek sub-watersheds. Scatter plots in Figure 9 show a significant correlation between NEXRAD and gauge precipitation data for 2008 (R 2 = 0.8) and 2017 events (R 2 = 0.6), while there is a less significant correlation for the 2005 event (R 2 = 0.30). The Upper Blue River sub-watershed data for the 2008 event display the highest correlation, with an R 2 value of 0.8. In contrast, for the 2005 event, both Upper Blue River and Indiana Creek sub-watersheds show the lowest correlation (R 2 = 0.3). There were some missing data in the gauge records for this event, which affected the correlation results. Validation outcomes of NEXRAD and gauge precipitation usually do not reveal a good relationship because of the possible errors in both of them [26,54,55]. Figure 10 compares the total precipitation of NEXRAD with gauge data. NEXRAD precipitation amounts were overestimated for two flood events compared with gauge rainfall. There is no difference between gauge and NEXRAD rainfall data for 2005 event for both sub-watersheds; however, it is about 0.4 percent for the 2008 event. For the 2017 event, the difference is about 0.2 percent for Indiana Creek and 0.5 percent for Upper Blue River sub-watersheds. The overestimation is probably a result of area (radar grid cell) and point (gauge) errors [54]. On the other hand, the gauge precipitation might be underestimated due to issues with gauge funnels, which might have been temporarily blocked by grass, bird debris, or other objects, and as a result, the precipitation could be missing or delayed and underestimated [55]. There were some missing data in the gauge records for this event, which affected the correlation results. Validation outcomes of NEXRAD and gauge precipitation usually do not reveal a good relationship because of the possible errors in both of them [26,54,55]. Figure 10 compares the total precipitation of NEXRAD with gauge data. NEXRAD precipitation amounts were overestimated for two flood events compared with gauge rainfall. There is no difference between gauge and NEXRAD rainfall data for 2005 event for both sub-watersheds; however, it is about 0.4 percent for the 2008 event. For the 2017 event, the difference is about 0.2 percent for Indiana Creek and 0.5 percent for Upper Blue River sub-watersheds. The overestimation is probably a result of area (radar grid cell) and point (gauge) errors [54]. On the other hand, the gauge precipitation might be underestimated due to issues with gauge funnels, which might have been temporarily blocked by grass, bird debris, or other objects, and as a result, the precipitation could be missing or delayed and underestimated [55].

Model Calibration and Validation
The preliminary simulated runoff hydrographs show a reasonable overall fit with observed data for the watershed. Although the runoff was overestimated in the validation results, the calibration attempts provided better results. Calibration and validation results at the Blue River watershed outlet are shown in Figure 11.

Model Calibration and Validation
The preliminary simulated runoff hydrographs show a reasonable overall fit with observed data for the watershed. Although the runoff was overestimated in the validation results, the calibration attempts provided better results. Calibration and validation results at the Blue River watershed outlet are shown in Figure 11.

Model Calibration and Validation
The preliminary simulated runoff hydrographs show a reasonable overall fit with observed data for the watershed. Although the runoff was overestimated in the validation results, the calibration attempts provided better results. Calibration and validation results at the Blue River watershed outlet are shown in Figure 11. Model calibration is the derivation of a set of model parameter values that produces the best fit to the observed data; on the other hand, the validation process is running the model for different events without changing the parameters. Calibrated models of 2005 and 2017 events provide better simulation results compared to the 2008 validation models (uncalibrated). There is almost no difference between the observed and calibrated amounts of peak discharge, whereas there is a 0.16 to 0.35 percent difference in direct runoff volume. Validation models provide relatively satisfactory results, although the peak discharge and runoff might be overestimated due to the precipitation amount overestimation in the NEXDAD data. The difference between observed and simulated peak discharge is about 0.36 percent, and direct runoff 0.5 percent. Better estimation for the rainfall amount in NEXRAD data could mitigate the runoff volume and simulate the peak discharge more accurately. At a sub-watershed level, the model tended to overestimate the runoff for some sub-watersheds, but for the rest of the sub-watersheds, the model more accurately estimated the runoff. Compared with the observed data, the calibrated hydrographs indicate relatively good performance. The difference in statistical results between all models with different data input resolutions and LULC conditions is not significant. All models simulated the runoff in a similar fashion. For instance, the results for the model with lower resolution and 2003 LULC data are E = 0.6, EPBIAS 60%, and RMSE Std Dev = 0.6, and the results for the model with higher resolution and 2017 LULC data are E = 0.7, EPBIAS = 32%, and RMSE Std Dev = 0.5. Thus, using moderate-resolution data (e.g., 30 m) for modeling can still provide satisfactory simulation outcomes. The validation hydrographs of the 2008 event might have been overestimated due to the overestimation of the rainfall in the NEXRAD data. All validation models performed similarly.  There is almost no difference between the two discharge amounts. This insignificant variance demonstrates that the simulation reflects outcomes that are consistent with the LULC change patterns in the study area and period. As discussed before, the built-up area in the watershed did not change notably during the study period, while the increase in forestland was offset by decreased farmland/grassland. Comparing the results of model 1.2, which uses higher spatial resolution data, with model 2.2, which uses lower resolution data, reveals that for each flood event all simulated values are close. For instance, in model 1.2 (uses 6 m LULC and 3 m DEM data), the simulated peak discharge is 33,103.3 cms for the 2017 event, while simulation with model 2.2 (with 30 m LULC and DEM data) results in a peak discharge of 33,295.7 cms. This also shows consistency with the low variation in the CN value and DEM processing outputs derived from different spatial resolution data, which led to similar simulation outputs and performance statistics. In addition, the HEC-HMS model is a generalized model system and applies mathematical models to represent the watershed flow; hence, small differences between derived inputs from 3 and 30 m DEMs or 6 and 30 m images are not enough to notably affect the model outputs. In previous studies, the impact of spatial and temporal heterogeneity of urban development processes on runoff model outcomes was not fully and adequately addressed, particularly in an urban watershed, such as in our case study area. There is a lack of understanding of the impact of different spatial resolutions on the model outcome and the suitability of coupling NEXRAD rainfall with the HEC-HMS model to assess LULC change impact on urban flooding. Our study attempts to further address this issue. In our 14-year study period, we found that the LULC in the Blue River watershed did not significantly change. Since 2003, the buildup area occupied about 1/3 of the watershed, and the majority of the area is vegetated land (forestland, farmland, and grassland). This less dynamic regime of the LULC in the watershed led to similar responses to precipitation and storm events. For runoff simulation applications, it is useful to utilize the standard available data resolution (e.g., 30 m) for satellite images to generate LULC data and CN values, or for DEMs to extract watershed parameters, which may cost less compared with data with higher spatial resolution. As LULC, NEXRAD data vary spatially, which provides a more accurate representation of the rainfall over each area compared with the ground-based gauge data that represent the rainfall amount only in the location area. In our case study, the use of (4 km) NEXRAD precipitation data for the simulation provides reasonably accurate runoff hydrographs. Using the dataset may contribute to runoff overestimation or underestimation due to the uncertainty and error that is usually associated with radar precipitation. At a sub-watershed scale, we found that using HEC-HMS, NEXRAD data can represent the rainfall amount and simulate the runoff for small sub-watersheds more accurately than large ones. Tests of other types of NEXRAD Level III precipitation data, such as the products of "dual-pol," which clearly identify and detect rainfall, might reduce rainfall overestimation and provide better results.

Conclusions
In this study, urban runoff was simulated in an urban watershed that has diverse LULC and development activities, and its streams tend to be frequently flooded. The study aims to understand how LULC change patterning and the quality of input geospatial data affect the simulation outcomes as well as examine the effectiveness of NEXRAD rainfall data in such modeling settings. The study results reveal that the simulation with the HEC-HMS model sensitively responds to the spatial and temporal patterning of LULC dynamics indicated by the change rate of imperious surfaces and vegetated land-use processes. Specifically, the simulation outcomes reflect the slow-down period of urban development and associated ecological restoration efforts in the case study watershed during the study period. The study indicates that the developed simulation approach can better tolerate small variations in derived input parameters, such as CN values, watershed boundaries, parameters, and stream networks, suggesting that input data with a moderate spatial resolution (e.g., 30 m) are suitable for urban runoff simulation at a watershed scale. Further, this study illustrates that, in such modeling, applying spatially distributed precipitation data, such as the one-hour NEXRAD Level III data, provide reliable and satisfactory outcomes after calibration efforts, particularly in hydrograph shape, peak discharge amounts, and time. In addition, the rainfall amount overestimation in NEXRAD data results in higher peak discharge and runoff volume as compared with observed data. Finally, the study proves the feasibility and effectiveness of incorporating satellite imagery-based LULC maps with related geospatial data, including DEM and distributed radar precipitation, in hydrological simulation to assess the watershed's hydrological response to flooding events. Funding: SPOT 2017 image purchase was funded by the University of Missouri-Kansas City School of Graduate Studies Research Grant. The other SPOT images were acquired through a previous remote sensing study, which was supported by the US EPA Grant CD 97701501 and by the grant of the Friends of the Library of the University of Missouri-Kansas City. A.E. received financial support from Umm AL-Qura University.