Mapping Urbanization and Evaluating Its Possible Impacts on Stream Water Quality in Chattanooga, Tennessee, Using GIS and Remote Sensing

Impervious surfaces (IS) produced by urbanization can facilitate pollutants’ movement to nearby water bodies through stormwater. This study mapped and estimated the IS changes in Chattanooga, Tennessee, using satellite imagery acquired in 1986 and 2016. A model was developed utilizing the Normalized Difference Vegetation Index coupled with density slicing to detect and map urbanization through IS growth. Urban growth was quantified at USGS HUC12 watershed level including stream riparian areas. The obtained results show a net growth of 45.12 km2 of IS with a heterogeneous distribution. About 9.96 km2 of this growth is within 90 m of streams, about 6% of the study site’s land cover. The Lower South Chickamauga Creek watershed experienced the largest urban growth with a change from 24.2 to 48.5 km2. Using the riparian zone percent imperviousness, a stream risk assessment model was developed to evaluate potential stream impairment due to this growth. Approximately 87, 131, and 203 km lengths of streams identified as potentially at high, very high, and extreme risks, respectively, to be impaired due to urban growth from the last 30 years. These findings would benefit to proactively implement sustainable management plans for the streams near rapidly urbanizing areas in the study site.


Introduction
Urban areas are often defined by anthropogenically created impervious surfaces such as concrete, asphalt, and metal [1,2]. Global urbanized areas, as of 2014, contain 54% of the global population and are estimated to increase to 66% by 2050. The rates of urbanization are not globally uniform. For example, in North America, 82% of the population resides in urban areas compared to other continents such as Africa and Asia with 40% and 48%, respectively. On the other hand, increases in urbanization are not consistent. The fastest growing areas are in Africa and Asia and are estimated to increase to 56% and 64%, respectively, with more urbanized areas such as Northern America growing at slower rates [3]. While slower than other areas globally, North America, and specifically the continental United States, is still experiencing significant growth both in urban population and urban development.

Urbanization and Surface Water Quality
As urban areas continue to densify and expand into their surrounding environment, their new development can alter the nearby landscape dynamics. Often these expanding urban areas encompass or encroach nearby surface water resources [18]. This spatial relationship between urban areas and water resources has long had the attention of the scientific community, whose work has shown that urban development can impact local surface water resources [19][20][21][22], especially the quality of water, as we can see in Table 1. All parameters listed in Table 1 have been found to be adversely affected by the listed land use types; with measured parameter levels all being higher than natural levels, (except dissolved oxygen which is found to be lower). Overall, urbanized watersheds have been found to have significant impacts on surface water resources including degraded water quality, habitat and biota. It is, therefore, crucial to monitor and investigate the quality of surrounding water resources near rapidly urbanizing areas to proactively implement water resource management plans.  [24,29,43] Note: The relationship is given by the general impact of the land cover on the parameter. A positive (+) or negative (-) symbol shows the parameter concentration response to increases in specified land cover. IS = Impervious Surfaces.

Land Use and Land Cover Change (LULC) and Urbanization
Urban growth is the change in LULC for and by anthropogenic activities such as residential, commercial and industrial development. It is commonly quantified through LULC detection and analysis of selected classes such as agricultural, open field, urban, forest, wetland, and barren areas. Quantitatively assessing urban growth requires defining which LULC categories should be considered urban. For the USGS National Land Cover Database [48][49][50], urban land cover is defined by spatial intensity of developed land, which is the percentage of impervious surface (IS) per unit area. Other efforts to map urban growth involve using the amount of detectable IS to measure urban extent [51][52][53][54]. Impervious surfaces are defined as being composed of materials such as asphalt, concrete, and metal that slow or inhibit water infiltration to topsoil. Therefore, these surfaces include building rooftops, sidewalks, parking lots, and roads, and the amount of land cover composed of IS can then be related to urban growth. Sustainability 2020, 12, 1980 3 of 46 Remote sensing has been established as an effective method for investigating urban growth due to the differing visible and spectral differences between LULC classes [55][56][57][58][59]. Urban growth detection using remote sensing has been used to investigate changes across a range of temporal and spatial scales due to the advantages offered by various remote sensing technologies. Through these investigations, numerous image processing techniques have been developed or utilized to improve LULC research precisions including: Tasseled Cap Transformation (TCT), Normalized Difference Built-Up Index (NDBI), Index-based Built-Up Index (IBI), Normalized Difference Vegetation Index (NDVI), and supervised classification schemes such as density slicing [56,[60][61][62][63][64][65][66]. Impervious surface growth has been investigated using remote sensing and is often paired with water quality investigations due to the established relationship between the two variables [67][68][69].

Urban Growth and Riparian Areas
The land immediately surrounding water resources is considered riparian areas. These zones can be described as a range of area adjacent to water resources including areas 100 m away and further [70]. The conservation and protection of forest riparian land is critical as they act as buffers between the flows of matter from upland areas into hydrologic resources and have a strong influence on the surrounding water resource quality [71]. They are attributed with mitigating the flow of sediment and nutrients from surface and groundwater into adjacent water resources [72], effects on fish assemblages [73], and other physiochemical water quality variables [74,75]. As with other critical habitats, LULC development poses a significant threat to riparian areas by compromising the buffer's ecological integrity and altering buffer areas' landscape [76,77]. Research has suggested that vegetated buffers of a minimum of 30 m is needed to act as a functioning non-point source pollutant control [78].
Acting as a non-point source pollutant, IS proximity to water resources and the respective riparian buffer have also been found to affect water quality [78][79][80]. Riparian buffer LULC can have a better ability to predict stream surface water quality leading to an increase in research and government attention to understanding and protecting these corridors [81,82]. Other research has found that the composition of riparian land cover has a slightly better ability to predict stream water chemistry than watershed land cover [83]. Threshold percent imperviousness values for stream riparian zones that have been found to begin indicating stream impairment begin at 10% cover [84].

Mapping Urban Growth
Traditional methods for investigating LULC change has been primarily through in situ surveying with measurements conducted by aerial photographs. These methods can provide accurate, detailed information but are limited temporally and spatially due to resources. Remote sensing and geographic information systems (GIS) provide useful techniques for mapping and analysis of historical and current LULC across larger areas than traditional methods [85]. Remote sensing data acquisition also have the ability for routine temporal data collection of the same area of interest, which allows for more efficient change detection analyses [86]. Since remote sensing and GIS acquisition and analysis of LULC data has been established, it presents a more temporal and fiscal affordable option for studying urbanization.
Acquisition of LULC data is primarily conducted by satellite-based optical sensors due to the ability to differentiate different types of LULC classes by their spectral responses in the EMS [87,88]. The partnered USGS and NASA Landsat mission is noteworthy for its frequent use in LULC studies with the Landsat mission starting in 1972 and continued sensor improvements [89]. Data from the Landsat missions have been used for the creation of a National Land Cover Database (NLCD) by the Multi-Resolution Land Characteristic Consortium (MRLC) to generate an accurate nationwide land cover map every five years.
Optical satellite sensors, including Landsat, have been well established for mapping urban growth, including impervious surfaces (IS) [19,51,[90][91][92]. While optical sensors are the primary data source used for LULC detection, there have been efforts to use active sensor technologies such as Synthetic Aperture RADAR (SAR) or light Detection and Ranging (LiDAR) for urbanization investigations due to the unique physical characteristics of many primary land cover features associated with urban areas such as buildings and roads [93,94]. Studies combining the two sensing technologies have found that the data fusion can increase model accuracy and information extraction [95,96].
Impervious surface (IS) mapping using multispectral remote sensing technologies has been used successfully for many years [19,68,97]. The use of NDVI is common for IS mapping by finding threshold values for water and vegetation in an attempt to isolate IS [98][99][100]. Supervised and unsupervised classification algorithms to extract spectral signature classes that can be assigned LULC values are also commonly used, however, traditional classification algorithms can confuse IS with other LULC classes such as bare or dry soils, shadowed areas caused by oblique image collection, and wetlands due to similar spectral signatures [98]. Image spectroscopy-based techniques such as linear spectral unmixing and multiple endmember spectral mixture analysis have been found successful for IS mapping but require significant data collection and processing for accurate results [101][102][103][104][105].

Analysis Techniques
The use of remote sensing and GIS can be highly beneficial; however, current image processing and geospatial analysis techniques may not be sufficient for accurately characterizing non-linear or complex urban growth relationships or predicting future urban growth [106][107][108][109][110]. Regression models and machine learning algorithms are therefore used in these applications to increase modeling and information extraction accuracy and offer better predictive capabilities. Regression models have been used successfully to relate historic urban growth to various socioeconomic variables [24,111], changes in water quality parameters [112,113], and predict future urban growth [114,115]. Machine learning is a type of artificial intelligence-based computer program that uses a pre-determined algorithm to be optimized on training datasets such as the target dataset, which it analyzes for algorithm parameter optimization. The result is a trained algorithm for a specific task that can be used for regression or classification analyses either by supervised or unsupervised methods [116][117][118]. Machine learning-based algorithms using support vector machines, decision tree, random decision forest (random forest)-based classification, and regression tree algorithms have been found to be very efficient and accurate when used for urban growth studies [108,119,120]. Although found successful for urban growth research, machine learning algorithms often require large datasets for training and validation, which can be very expensive and time-consuming for some research.

Statement of the Problem and Scope of the Study
The City of Chattanooga, Tennessee, has grown substantially during the last several decades and has become the center of a series of urbanized sub-watersheds. This significant recent growth necessitates a much-needed effort to be aware of the location of the growth and its relationship to the surrounding water resources since the environmental impacts, especially the quality of surface waters due to this growth, have become a major concern for the sustainable developments of the greater Chattanooga area [17,121].
Chattanooga's water quality has been the topic of a few studies in recent years. Schorr et al. [122] conducted an assessment of water quality and aquatic biota in Chattanooga area streams partnering with the municipal Stormwater Management office. Long and Schorr [35] investigated urban land use at the watershed scale in the Chattanooga area and its effects on selected streams. They found that urbanized watersheds were negatively correlated with fish species diversity and biotic integrity and that this negative correlation was related, in part, to high levels of sedimentation. Other than Long and Schorr [35] and Schorr et al. [122] to date, there have not been any published attempts to quantitatively map urban growth and analyze the development's relationship to water resources in the Chattanooga, Tennessee, area.
Therefore, this study was designed to investigate the potential impact of present urban growth on water quality in and around the City of Chattanooga, Tennessee, using satellite observed geospatial data. The investigation was conducted by accomplishing the following three sequential and related Sustainability 2020, 12, 1980 5 of 46 objectives: (1) mapping the net impervious surface (IS) change in the study site between 24 January 1986 and 26 November 2016 using satellite imagery from NASA's Landsat satellite missions, (2) performing a quantitative analysis of the obtained IS changes in relation to local water resources, and (3) developing a risk assessment model to identify the potential areas of concern for surface water quality in the study site due to the proximity and quantity of IS growth.
This study can be considered as the first attempt to quantitatively map urban growth in the Chattanooga, Tennessee, area and analyze the development's spatial relationship to local water resources. Through the utilization of remotely sensed data and GIS analyses, this research aims to better understand if/how urbanization has impacted the land use, land cover and landscape dynamics of the greater Chattanooga area.

Study Site
The City of Chattanooga is a growing metropolitan area on the southeastern border of Tennessee ( Figure 1). It is settled along the Tennessee River. The city has a strong relationship with its surrounding environment, relying on it for utilities, commerce, recreation and tourism. Downtown Chattanooga lies directly adjacent to the Tennessee River, and has several major parks and trails following along the riverfront. The parks, along with other attractions such as the 120-million-dollar downtown Riverwalk, have earned Chattanooga multiple national awards for the best outdoor city [123,124]. The city of Chattanooga reported that tourism alone generated 1 billion dollars for the economy in 2015, with the largest attraction, the city's Aquarium, being along the Tennessee River [125]. These awards and the revenue generated by tourism reflect the impact that the environment and outdoor attractions have on the suitably nicknamed "Scenic City". The study site of this research includes seven USGS Hydrologic Unit Code 12 (HUC-12) watersheds located in Hamilton County, Tennessee, as seen in Figure 1. These watersheds were selected as the stream networks within them feed into the Tennessee River. Thus, having the potential for sediments to move from within the HUC-12 watersheds to the Tennessee River. Hamilton County, located in the southeastern portion of Tennessee, occupies part of the Appalachian Valley and is bisected by the Tennessee River. The total areas covered by land and water within the county are of 1403.8 km 2 and of 86.3 km 2 , respectively [133]. The elevation of the areas within the study site ranges from 191 m to 655 m [134]. The study site lies in a complex geologic subregion that is primarily composed of carbonate rock such as limestone and dolomite; however, there is a considerable amount of shale, including Chattanooga shale, sandstone and other rock types [135].
Previous studies found that the watershed characteristics such as geology, topography, soil type, and land use and land cover (LULC) had noticeable effects on the quality of surface water in Hamilton County [136]. However, it is also reported that since the Tennessee Ridge and Valley ecoregion's The city's legacy with water resource issues is tied to its industrial history with several large industries being on or nearby local waterways and causing heavy, historic contamination of the local water resources [126]. The most significant example of Chattanooga's historic impact on local waterways is Chattanooga Creek. This waterway is currently listed as an EPA superfund site and has Sustainability 2020, 12,1980 6 of 46 been since 1995 due to extensive industrial dumping. Since listing, the EPA, along with regional and local agencies and potentially responsible parties, have conducted several projects to clean up and restore the creek. However, Chattanooga Creek is still currently listed as a superfund site and a fish consumption advisory is still in place [127][128][129]. Other local Chattanooga waterways (e.g., Stringers Branch) have also had issues with historic pollution [130].
In recent years, the city of Chattanooga has experienced a steady rate of population growth, primarily due to the rapid economic growth in its six-county metropolitan area in southeastern Tennessee and northwestern Georgia [131]. Chattanooga is home to several large organizations such as Volkswagen, Unum, Tennessee Valley Authority, Blue Cross Blue Shield, Wacker, and Amazon that are driving the area's economic growth. The economic growth of the city is heavily attributed to the implementation of the fiber optic internet by the City's Electric Power Board (EPB) having 1 gigabit and now 10 gigabytes of service. With multiple large corporations continuing to expand, a nationally ranked internet infrastructure, and a supported nickname as the "Scenic City", Chattanooga's economic and social environment is becoming increasingly attractive to startup businesses. The city has also designated a 140-acre portion of the downtown area for startups, small businesses, nonprofits and government offices called the "Innovation District" [132].
The study site of this research includes seven USGS Hydrologic Unit Code 12 (HUC-12) watersheds located in Hamilton County, Tennessee, as seen in Figure 1. These watersheds were selected as the stream networks within them feed into the Tennessee River. Thus, having the potential for sediments to move from within the HUC-12 watersheds to the Tennessee River.
Hamilton County, located in the southeastern portion of Tennessee, occupies part of the Appalachian Valley and is bisected by the Tennessee River. The total areas covered by land and water within the county are of 1403.8 km 2 and of 86.3 km 2 , respectively [133]. The elevation of the areas within the study site ranges from 191 m to 655 m [134]. The study site lies in a complex geologic sub-region that is primarily composed of carbonate rock such as limestone and dolomite; however, there is a considerable amount of shale, including Chattanooga shale, sandstone and other rock types [135].
Previous studies found that the watershed characteristics such as geology, topography, soil type, and land use and land cover (LULC) had noticeable effects on the quality of surface water in Hamilton County [136]. However, it is also reported that since the Tennessee Ridge and Valley ecoregion's geology, topography, and soil composition remained unchanged for the past 50 years, spatiotemporal changes in surface water quality should, therefore, be associated with changes in LULC. This is also supported by the USGS National Water Quality Assessment that concludes that sediment contamination of surface stream water quality is persistent, and nutrient loadings in sub-basins for the Upper Tennessee River Basin are primarily influenced by land use and streamflow [137].

Data Collection and Processing
A pair of satellite imagery acquired by NASA's Landsat missions was used in this study. The first scene was acquired by Landsat 5 Thematic Mapper (TM) on 24 January 1986 and was used to study the land use and land cover (LULC) of 1986. The second scene was acquired by Landsat 8 Operational Land Imager (OLI) on 26 November 2016 and was used to study the LULC of 2016. The processing and analysis of the imagery and the overall methodological scheme are explained in Figure 2.
Both satellite images in the Red Green Blue (RGB) and Near Infra-Red (NIR) spectrums have a spatial resolution of 30 m. The true color displays of Landsat 5 TM and Landsat 8 OLI images for the study site can be seen in Figures 3 and 4, respectively. Since the two images were collected by different Landsat sensors with a significant temporal gap between acquisitions, atmospheric correction was needed to convert both images into the same radiometric scale [138]. The conversion of the image values to represent the surface reflectance requires the correction of atmospheric effects on exiting electromagnetic radiation, which is an important image pre-processing step [139,140]. The Landsat 5 TM scene reflectance values were generated using NASA's Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) software [141]. For the Landsat 8 OLI scene, NASA's Landsat Surface Reflectance Code (LaSRC) was used to generate the values [142].
Surface Reflectance Code (LaSRC) was used to generate the values [142].
A one-meter resolution multispectral aerial photograph acquired by the United States National Agricultural Imagery Program (NAIP) in 2014 over the study site was obtained for the accuracy assessment part of the study. The obtained NAIP imagery included both RGB and NIR bands. The Tennessee county and state boundary data were collected from the Tennessee State government GIS data server and the U.S. Census Bureau's data server, respectively. The watershed data were collected from the USGS National Hydrography Dataset Plus (NHDplus) and the City of Chattanooga Water Quality Program. The stream location data were collected from NHDplus and the City of Chattanooga GIS Office.

Mapping Impervious Surfaces
The satellite data used for this research were pre-processed including geo-rectification and restoration [143]. Once the Landsat data and NAIP imagery were collected, necessary bands for each dataset were masked to the study site and stacked to create the Red Green Blue (RGB) and Near Infrared (NIR) composite imagery. The model designed for mapping impervious surfaces (IS) is based on  A one-meter resolution multispectral aerial photograph acquired by the United States National Agricultural Imagery Program (NAIP) in 2014 over the study site was obtained for the accuracy assessment part of the study. The obtained NAIP imagery included both RGB and NIR bands. The Tennessee county and state boundary data were collected from the Tennessee State government GIS data server and the U.S. Census Bureau's data server, respectively. The watershed data were collected from the USGS National Hydrography Dataset Plus (NHDplus) and the City of Chattanooga Water Quality Program. The stream location data were collected from NHDplus and the City of Chattanooga GIS Office.

Mapping Impervious Surfaces
The satellite data used for this research were pre-processed including geo-rectification and restoration [143]. Once the Landsat data and NAIP imagery were collected, necessary bands for each dataset were masked to the study site and stacked to create the Red Green Blue (RGB) and Near Infra-red (NIR) composite imagery. The model designed for mapping impervious surfaces (IS) is based on the application of Normalized Difference Vegetation Index (NDVI) due to the breaks in values for different land covers, which can be used to isolate impervious surface (IS) areas. To assist in land cover isolation, the NIR band was used to separately classify water resources, and the green band for shadows in both images. The model developed for detecting and mapping impervious surfaces is represented by Equation (1).
where IS NDVI represents a NDVI image that has been classified into IS and vegetated areas, Shadow Green represents shadows (detected by the green band of a Landsat sensor), and Water NIR represents surface waterbodies (detected by the NIR band of a Landsat sensor). All image processing and analyses were conducted using ArcGIS Pro GIS software [144]. The statistical tests applied were conducted in the Jupyter Notebook environment using Python scripting language [145].

Normalized Difference Vegetation Index (NDVI)
The primary component of the developed IS mapping model is NDVI. This index, seen in Equation (2), uses the NIR and red (R) portions of the EMS.
where ρ NIR is the surface reflectance value of NIR and ρ R is the surface reflectance value of R. The Normalized Difference Vegetation Index was originally developed to map green vegetation, which has a sharp contrast in the amount of absorption of red wavelengths compared to near infrared [62]. Green vegetation has a high absorption of light in the red portion of the Electromagnetic Spectrum (EMS) while reflecting a significant portion of the NIR portion; this response is driven primarily by the vegetation's photosynthetic activity [146]. This difference in EMS absorption causes areas of healthy, green vegetation to have large, positive values from NDVI. This index can also be used to map other non-vegetated land covers, primarily water and IS, due to each land cover's distinct interactions with red and NIR wavelengths [55]. Where IS have close to equal reflection of both red and NIR light, they generate NDVI values close to 0 [146]. Negative NDVI values indicate non-vegetated areas often consisting of water [146]. A NDVI was generated for each image, to separate the vegetated and non-vegetated areas [64]. This concept of mapping vegetation to subsequently extract exposed IS has been previously established as an alternative method for estimating IS [147]. (4).

The Landsat 5 TM NDVI equation can be seen in Equation (3) and the Landsat 8 OLI NDVI equation can be seen in Equation
where ρ B4 and ρ B3 represent the reflectance values of NIR and R, respectively.
where ρ B5 and ρ B4 represent the reflectance values of NIR and R, respectively. Figure 5 shows the NDVI image generated for 1986 and Figure 6 shows the NDVI image generated for 2016.
areas of healthy, green vegetation to have large, positive values from NDVI. This index can also be used to map other non-vegetated land covers, primarily water and IS, due to each land cover's distinct interactions with red and NIR wavelengths [55]. Where IS have close to equal reflection of both red and NIR light, they generate NDVI values close to 0 [146]. Negative NDVI values indicate non-vegetated areas often consisting of water [146]. A NDVI was generated for each image, to separate the vegetated and non-vegetated areas [64]. This concept of mapping vegetation to subsequently extract exposed IS has been previously established as an alternative method for estimating IS [147]. The Landsat 5 TM NDVI equation can be seen in Equation (3) and the Landsat 8 OLI NDVI equation can be seen in Equation (4).
where ρB4 and ρB3 represent the reflectance values of NIR and R, respectively.
where ρB5 and ρB4 represent the reflectance values of NIR and R, respectively. Figure 5 shows the NDVI image generated for 1986 and Figure 6 shows the NDVI image generated for 2016.

Image Classification
To classify desired land covers, density slicing classification was performed. This classification method requires the researcher to determine the spectral response range of desired variables from a single layer image. This is conducted by an initial manual or statistical-based classification of an image based on the distribution of the pixel values. Then, continued re-classification is conducted by adjusting pixel value threshold until the desired land cover has been effectively isolated. Reclassification and verification can be improved using reference location data, which can be colored imagery and/or pre-labeled reference areas. Density slicing technique was applied on the generated Landsat 5 and Landsat 8 Normalized Difference Vegetation Index (NDVI) images to classify them into impervious surface (IS) and vegetated areas. Density slicing was also used to categorize waterbodies using the NIR band and shadows with the green band for each date. To assist in classification for each date, the true color Landsat image was used as a visual reference for each date. Reference polygons were used to assist in determining IS land cover threshold values for both 1986 and 2016. Figures 7 and 8 show the density slicing reference locations used. Figure 9 shows examples of reference polygons for both 1986 and 2016 with the NDVI image and true color image for comparison. During density slicing, repeated inspection of these anchor points were conducted to gauge the effectiveness of the used pixel threshold values. Table 2 gives the threshold values determined for each class.

Image Classification
To classify desired land covers, density slicing classification was performed. This classification method requires the researcher to determine the spectral response range of desired variables from a single layer image. This is conducted by an initial manual or statistical-based classification of an image based on the distribution of the pixel values. Then, continued re-classification is conducted by adjusting pixel value threshold until the desired land cover has been effectively isolated. Reclassification and verification can be improved using reference location data, which can be colored imagery and/or pre-labeled reference areas. Density slicing technique was applied on the generated Landsat 5 and Landsat 8 Normalized Difference Vegetation Index (NDVI) images to classify them into impervious surface (IS) and vegetated areas. Density slicing was also used to categorize waterbodies using the NIR band and shadows with the green band for each date. To assist in classification for each date, the true color Landsat image was used as a visual reference for each date. Reference polygons were used to assist in determining IS land cover threshold values for both 1986 and 2016. Figures 7 and 8 show the density slicing reference locations used. Figure 9 shows examples of reference polygons for both 1986 and 2016 with the NDVI image and true color image for comparison. During density slicing, repeated inspection of these anchor points were conducted to gauge the effectiveness of the used pixel threshold values. Table 2 gives the threshold values determined for each class.

Data Integration
Using raster algebra, the final IS dataset for each date was created by combining the classified NDVI, green band and Near Infra-Red (NIR) band datasets using a nested conditional statement executed with the ArcGIS raster calculator tool. The conditional statement isolated each desired land covers from each dataset and returned a single dataset containing all the classified classes. The resulting datasets were then reclassified into two classes: (1) impervious surfaces and (2) pervious and other surfaces. This reclassified thematic image represents the final IS dataset. The IS dataset is seen in Figures 10 and 11 for 1986 and 2016, respectively.

Accuracy Assessment
The accuracy assessments of the resulting impervious surface (IS) datasets were conducted using two different techniques: (1) Zonal Statistics and (2) Kappa Statistics.

Zonal Statistics
Zonal statistics was used for the accuracy assessment of the (IS) models for both 2016 and 1986. For the 2016 model, this was conducted by digitizing 140 randomly located reference IS areas within the study site. Figure 12 shows the locations for the 2016 reference IS areas. Each HUC-12 watershed contains 20 of the digitized reference polygons. Every reference IS polygon covers an area equal to 15 m 2 . The NAIP dataset was used as the source of these reference IS polygons. A majority count zonal statistic was applied to the digitized reference polygons overlaying the impervious dataset. This was conducted to determine if the majority of the classified pixels covered by the reference polygons were impervious or pervious. Due to resource constraints, it was not possible to obtain high resolution aerial photography to generate reference polygons to evaluate the 1986 dataset. In this case, a true color Landsat 5 TM image was used to obtain the 140 reference polygons. The reference polygons digitized using the NAIP data that contained 100% IS on the Landsat 5 TM image were used for the accuracy assessment. Only 60 polygons digitized with NAIP imagery represented 100% IS cover in 1986. These polygons were then used to assess the accuracy of the 1986 IS dataset following the methods used for the 2016 dataset.

Accuracy Assessment
The accuracy assessments of the resulting impervious surface (IS) datasets were conducted using two different techniques: (1) Zonal Statistics and (2) Kappa Statistics.

Zonal Statistics
Zonal statistics was used for the accuracy assessment of the (IS) models for both 2016 and 1986. For the 2016 model, this was conducted by digitizing 140 randomly located reference IS areas within the study site. Figure 12 shows the locations for the 2016 reference IS areas. Each HUC-12 watershed contains 20 of the digitized reference polygons. Every reference IS polygon covers an area equal to 15 m 2 . The NAIP dataset was used as the source of these reference IS polygons. A majority count zonal statistic was applied to the digitized reference polygons overlaying the impervious dataset. This was conducted to determine if the majority of the classified pixels covered by the reference

Kappa Statistics
A confusion matrix and Cohen's Kappa-coefficient (k) were calculated for each date. The Kappa statistic has been shown to be an effective measure of a classification model in remote sensing and other scientific fields, due to its ability to evaluate the interclassifier agreement and remove the bias [148][149][150]. Kappa statistics evaluated the performance of the IS classification models by calculating the User's accuracy (Type I error or false positive) and Producer's accuracy (Type II error or false negative) for each class, the proportion of pixels correctly classified (PCC), and Kappa-coefficient of agreement. The Producer's accuracy reports the error of omission created by the model, quantifying the amount of each class that have not been correctly classified. The equation for the Producer's accuracy is shown by Equation (5).
where C T is the amount of a class C that is correctly classified, and O C is the sum of other classes that were classified as C. The proportion of pixels correctly classified (PCC) is calculated by dividing the total number of correctly classified pixels by the total number of pixels classified. The equation for calculating overall accuracy of a classification model is shown as Equation (7).
where nT is the sum of correctly classified subjects and n is total subject sample size.
The Kappa-coefficient of agreement (k) describes the accuracy of the classification model compared to random classification on the assumption that some of the pixels could have been classified correctly by chance. The equation for calculating the Kappa-coefficient is shown by Equation (8). where PO is the observed agreement of classification, PE is the chance agreement of classification and 1 represents maximum agreement. To calculate PO and PE Equations (9) and (10)  The User's accuracy reports the error of commission which is the amount of other classes that have been classified incorrectly. The equation for the User's accuracy is shown by Equation (6).
where C E is the sum of the class that are incorrectly classified. The proportion of pixels correctly classified (PCC) is calculated by dividing the total number of correctly classified pixels by the total number of pixels classified. The equation for calculating overall accuracy of a classification model is shown as Equation (7).
where n T is the sum of correctly classified subjects and n is total subject sample size. The Kappa-coefficient of agreement (k) describes the accuracy of the classification model compared to random classification on the assumption that some of the pixels could have been classified correctly by chance. The equation for calculating the Kappa-coefficient is shown by Equation (8).
where P O is the observed agreement of classification, P E is the chance agreement of classification and 1 represents maximum agreement. To calculate P O and P E Equations (9) and (10) were used.
where O D is the sum of the observed frequencies along the diagonal in the confusion matrix, E D is the sum of expected frequencies along the horizontal and n is the total number of subjects.
To generate the confusion matrix, 250 randomly stratified accuracy assessment points were generated for each date. These points were validated using the relevant Landsat image. For 2016, the NAIP imagery from 2014 was also used in this regard. The locations of the assessment points of the confusion matrix generation can be seen in Figures 13 and 14 for 1986 and 2016, respectively. where OD is the sum of the observed frequencies along the diagonal in the confusion matrix, ED is the sum of expected frequencies along the horizontal and n is the total number of subjects.
To generate the confusion matrix, 250 randomly stratified accuracy assessment points were generated for each date. These points were validated using the relevant Landsat image. For 2016, the NAIP imagery from 2014 was also used in this regard. The locations of the assessment points of the confusion matrix generation can be seen in Figures 13 and 14 for 1986 and 2016, respectively.

Water Resource Proximity Analysis
To understand the spatial relationship of the impervious surface (IS) growth regarding water resources within the study site, IS change per HUC-12 watershed (and stream riparian areas) was measured. The obtained IS datasets were clipped by each of the selected HUC-12 watershed boundaries for quantification of the IS change in each HUC-12 watershed. IS area was calculated by multiplying the total number of IS pixels with the area of a single pixel (30 m x 30 m). Finally, the obtained values were converted to km 2 .
For the stream analyses, streams within the study site were obtained from the NHDplus dataset and combined with auxiliary data from the City of Chattanooga Water Quality Program. The first stream analysis was performed by creating buffer distances of 30 m, 60 m, and 90 m on both sides of the selected streams. A map showing the stream buffers generated can be seen in Figure 15. After the buffer polygons were generated, they were clipped by each of the HUC-12 watersheds to create separate watershed riparian buffer boundaries for each watershed within the study site. The obtained polygon features were then used to extract impervious surface (IS) pixels from both IS datasets. IS area was then calculated in km 2 for each of the buffer features using the quantity of IS pixels and the area coverage for a single pixel. The proximity analysis was conducted in ESRI's ArcGIS Pro GIS environment and the bar charts for reporting the results were generated using the Python scripting language in the Jupyter Notebook environment.

Water Resource Proximity Analysis
To understand the spatial relationship of the impervious surface (IS) growth regarding water resources within the study site, IS change per HUC-12 watershed (and stream riparian areas) was measured. The obtained IS datasets were clipped by each of the selected HUC-12 watershed boundaries for quantification of the IS change in each HUC-12 watershed. IS area was calculated by multiplying the total number of IS pixels with the area of a single pixel (30 m × 30 m). Finally, the obtained values were converted to km 2 .
For the stream analyses, streams within the study site were obtained from the NHDplus dataset and combined with auxiliary data from the City of Chattanooga Water Quality Program. The first stream analysis was performed by creating buffer distances of 30 m, 60 m, and 90 m on both sides of the selected streams. A map showing the stream buffers generated can be seen in Figure 15. After the buffer polygons were generated, they were clipped by each of the HUC-12 watersheds to create separate watershed riparian buffer boundaries for each watershed within the study site. The obtained polygon features were then used to extract impervious surface (IS) pixels from both IS datasets. IS area was then calculated in km 2 for each of the buffer features using the quantity of IS pixels and the area coverage for a single pixel. The proximity analysis was conducted in ESRI's ArcGIS Pro GIS environment and the bar charts for reporting the results were generated using the Python scripting language in the Jupyter Notebook environment.

Assessing Stream Risk
To assess the probability of stream impairment due to impervious surface (IS) change, a model was developed to independently assess streams for their potential risks. The model evaluates stream segments to provide further detail regarding sections of streams that may have higher risk due to new or existing developments. More specifically, the model accounts for both the quantity of changes in IS and the proximity of those changes to the stream segments evaluated. Percent imperviousness was calculated for each buffered polygon and was incorporated with the model. The risk assessment was conducted in ESRI's ArcGIS Pro GIS environment. The bar charts were generated using the Python scripting language in the Jupyter Notebook environment [145].

Stream Segmentation
The segmentation of streams was done by generating points at every 90 m along the stream with points also placed at the ends of the stream and at the intersection of stream branches. These points were then used to splice the original stream datasets and to act as end points for new, unique stream segments. A total of 12,910 stream segments were generated within the study site. Some segments were shorter than 90 m due to stream branching and/or small stream segment lengths.

Riparian Zone Generation
Both IS datasets were converted from pixels to points, which represented the centers of the original pixels. The obtained IS point datasets were used to perform a spatial analysis technique called 'spatial join'. Spatial join summarizes and joins the attributes of features that meet a predetermined spatial criterion, to be performed for points within 30 m, 60 m, and 90 m of each stream

Assessing Stream Risk
To assess the probability of stream impairment due to impervious surface (IS) change, a model was developed to independently assess streams for their potential risks. The model evaluates stream segments to provide further detail regarding sections of streams that may have higher risk due to new or existing developments. More specifically, the model accounts for both the quantity of changes in IS and the proximity of those changes to the stream segments evaluated. Percent imperviousness was calculated for each buffered polygon and was incorporated with the model. The risk assessment was conducted in ESRI's ArcGIS Pro GIS environment. The bar charts were generated using the Python scripting language in the Jupyter Notebook environment [145].

Stream Segmentation
The segmentation of streams was done by generating points at every 90 m along the stream with points also placed at the ends of the stream and at the intersection of stream branches. These points were then used to splice the original stream datasets and to act as end points for new, unique stream segments. A total of 12,910 stream segments were generated within the study site. Some segments were shorter than 90 m due to stream branching and/or small stream segment lengths.

Riparian Zone Generation
Both IS datasets were converted from pixels to points, which represented the centers of the original pixels. The obtained IS point datasets were used to perform a spatial analysis technique called 'spatial join'. Spatial join summarizes and joins the attributes of features that meet a pre-determined spatial criterion, to be performed for points within 30 m, 60 m, and 90 m of each stream segment. This zone of influence considers all land cover within these distances from any point along the segment, including the start and end points. The attributes of each IS dataset were then summarized per segment of the stream in the zone.
These results allowed for quantifying and locating the IS growth around each stream segment within the study site. Quantification of the amount of IS development were normalized using the number of points within each zone. The normalization by area generated a percent imperviousness value per zone of influence surrounding each stream segment. By subtracting the 30 m zone results from the 60 m zone results and the 60 m zone from the 90 m zone, the area with each segment was partitioned into three evenly spaced, 30m wide riparian zones. This relationship between spatial join buffer sizes and riparian zone creation can be seen in Table 3. Figure 16 demonstrates how each segment considers the points within the zone of influence and how the percent imperviousness for the riparian area is generated. The model is based on the understandings that riparian areas up to 150 m from the stream can have significant impacts of stream water quality [151]. The proximity of the development is considered to have a stronger impact on the stream than the development which occurred father away [152,153].
summarized per segment of the stream in the zone.
These results allowed for quantifying and locating the IS growth around each stream segment within the study site. Quantification of the amount of IS development were normalized using the number of points within each zone. The normalization by area generated a percent imperviousness value per zone of influence surrounding each stream segment. By subtracting the 30 m zone results from the 60 m zone results and the 60 m zone from the 90 m zone, the area with each segment was partitioned into three evenly spaced, 30m wide riparian zones. This relationship between spatial join buffer sizes and riparian zone creation can be seen in Table 3. Figure 16 demonstrates how each segment considers the points within the zone of influence and how the percent imperviousness for the riparian area is generated. The model is based on the understandings that riparian areas up to 150 m from the stream can have significant impacts of stream water quality [151]. The proximity of the development is considered to have a stronger impact on the stream than the development which occurred father away [152,153].  Figure 16. Visual description of how stream segment percent imperviousness was calculated for the stream risk assessment. Assigned color and number ID are used to aid in identifying segment-buffer pairs. An example chart is given to describe the attributes calculated from the segment-buffer pair and the resulting percent imperviousness. Figure 16. Visual description of how stream segment percent imperviousness was calculated for the stream risk assessment. Assigned color and number ID are used to aid in identifying segment-buffer pairs. An example chart is given to describe the attributes calculated from the segment-buffer pair and the resulting percent imperviousness. A linear weighted model was then developed and applied to each segment to calculate the probability of surface water impairment due to the distance and magnitude of percent imperviousness in each riparian zone. This model can be seen in Equation (10).

Relating IS Development Proximity and Quantity to Risk
where X 1 , X 2 , and X 3 , represent the percent imperviousness (normalized) within the three separate 30 m riparian zones, and W 1 , W 2 , and W 3 , represent the IS-stream influence weight. The values for the IS-stream influence weight are shown as the subscript numeral. The weights were incorporated arbitrarily using the values of 1 to 3, where the value of 3 represents the most influential, 2 represents the moderately influential and 1 represents the least influential. These values reflect that usually the influence of the growth of IS on the stream will have more influence from the nearby development.
The output of this model indicates the portions of the streams, within the study site, that have an increased probability of being impaired due to the surrounding IS development.

Impervious Surface (IS) Mapping
According to the results of zonal statistics, the accuracy for the IS classification of 1986 was 88.3%. That is, 53 of the 60 reference polygons were correctly classified as IS. The accuracy for the IS classification of 2016 was 90%. That is, 126 of the 140 reference polygons were correctly classified as IS. The confusion matrix prepared for 1986 is shown in Table 4. The overall accuracy reported by the confusion matrix for 1986 is 90.0%. The Kappa-coefficient calculated is 0.624, showing that the model for 1986 classified 62.4% better than a random classification of the data. The confusion matrix prepared for 2016 is shown in Table 5. The overall accuracy reported by the confusion matrix for 2016 is 84.8%. The Kappa-coefficient calculated is 0.545 showing that the model for 2016 performed 54.5% better than a random classification. The values of Kappa statistic ranging from 0.41 to 0.6, and from 0.61 to 0.80 can be interpreted to represent moderate and substantial agreement for a classification, respectively [154]. According to this interpretation, the accuracy of the classification model for 1986 dataset can be described as substantial. Similarly, on the other hand, the accuracy of the classification model for 2016 dataset can be described as moderate.  The total HUC-12 IS area calculations performed on this study's model show that there has been a significant growth in the study site as shown in Figure 17 and described in Table 6. The net growth within the study site was 45.12 km 2 . This growth was not spatially equal in its distribution and occurred heavily in the Lower South Chickamauga Creek watershed with an increase of 24.3 km 2 of IS. The Lower South Chickamauga Creek watershed had the largest percent increase in impervious development with a change from 24.2 to 48.5 km 2 , equaling slightly more than a 100% increase in impervious area. All but the Chattanooga Creek watershed showed an increase in IS with development being less than 10 km 2 in each area. Finally, the Chattanooga Creek watershed had a small decrease in IS area, decreasing by 0.01 km 2 .

Water Resource Proximity Analysis
The impervious surface (IS) growth estimated shows that there has been significant IS growth near selected streams in the study site. The results of this analysis can be seen in Figure 18 and summarized in Table 7

Water Resource Proximity Analysis
The impervious surface (IS) growth estimated shows that there has been significant IS growth near selected streams in the study site. The results of this analysis can be seen in Figure 18 and summarized in Table 7. In 1986, IS cover accounted for 14.9% of the total area within 30 m of streams, increased to 18.1% by 2016. Compared to the change within 90 m of streams, where in 1986, IS cover occupied 15.2%, increased to 20.9% by 2016. The net growth of IS within 90 m of streams was 9.96 km 2 . The total increase of IS within the first 30 m was 2.04 km 2 , with an increase of 3.43 km 2 between 30 m and 60 m, and 4.49 km 2 between 60 m and 90 m of streams. The average change in IS area for the 30 m, 60 m and 90 m buffers in each HUC-12 watershed was 0.29, 0.49, and 0.64 km 2 , respectively. IS development within 90 m of streams accounted for 22% of the total IS development detected. Changes in IS extent within the stream riparian areas within the study site are shown in Figure 18, which shows the coverage of IS in riparian areas in green and the added extent in 2016 in red.    Creek watershed experienced a decrease in IS within 30 m and 60 m of streams but did experience a small increase within 90 m. There was growth for the remaining five watersheds at each buffer distance with similar trends found in the total HUC-12 IS area analysis. The Lower South Chickamauga Creek watershed showed the largest growth of IS development in each of the three distance zones. It experienced an increase of 1.35 km 2 within 30 m of streams, 3.29 km 2 within 60 m, and 5.63 km 2 within 90 m. This equated to being 66.3%, 60.1% and 56.4% of the total IS development regarding each buffer distance. A portion of the Lower South Chickamauga Creek IS change in relation to stream buffers is shown in Figure 19.

Riparian Percent Imperviousness
For the 1986 dataset, the first two 30 m riparian interval zones on average were 15.2% covered by IS with the third interval having an average of 15.5%. In the 2016 dataset, the average percent impervious cover for all the 30 m riparian zones had increased to 17.9%, 20.1%, and 21.6% respectively, with the amount of growth increasing for further riparian zones. Descriptive statistics for the riparian development analysis can be seen in Table 8. A visual description of the distributions of percent imperviousness for each zone in 1986 and 2016 can be seen below in Figure 20 with outliers being visualized as plus symbols above the top whisker. The distributions also show that the two further riparian zones show increases in the second and third quartiles compared to 1986. Maps showing the percent imperviousness for each riparian zone can be found in Appendix A (Figures A1-A6).

Riparian Percent Imperviousness
For the 1986 dataset, the first two 30 m riparian interval zones on average were 15.2% covered by IS with the third interval having an average of 15.5%. In the 2016 dataset, the average percent impervious cover for all the 30 m riparian zones had increased to 17.9%, 20.1%, and 21.6% respectively, with the amount of growth increasing for further riparian zones. Descriptive statistics for the riparian development analysis can be seen in Table 8. A visual description of the distributions of percent imperviousness for each zone in 1986 and 2016 can be seen below in Figure 20 with outliers being visualized as plus symbols above the top whisker. The distributions also show that the two further riparian zones show increases in the second and third quartiles compared to 1986. Maps showing the percent imperviousness for each riparian zone can be found in Appendix A ( Figures  A1-A6).

Stream Risk Assessment
The range in values for the stream impairment risk model was 0 through 6, allowing for different levels of probable risk to be assigned. Table 9 shows the levels of risk related to the equivalent range in risk model values. The basis for these levels originated from the amount of development that must be present to achieve this score. For a stream segment to be assigned the lowest value (0.0) there must

Stream Risk Assessment
The range in values for the stream impairment risk model was 0 through 6, allowing for different levels of probable risk to be assigned. Table 9 shows the levels of risk related to the equivalent range in risk model values. The basis for these levels originated from the amount of development that must be present to achieve this score. For a stream segment to be assigned the lowest value (0.0) there must be 0% imperviousness in all three riparian zones while the maximum value (6.0) requires 100% imperviousness in all three riparian zones. Therefore, a score between 0 and 1 shows a stream segment with little to no risk of impairment related to IS development. A score between 5 and 6 shows areas that have significant impervious surfaces in their riparian zones, and thereby should have higher risk of impairment related to IS development. Figures 21 and 22 show the possible risk of stream impairments due to IS development for 1986 and 2016, respectively.   Table 9. Model scores for potential stream risk assessments due to riparian imperviousness.    The average score per stream segment was 0.92 and 1.06 for 1986 and 2016, respectively. Statistics describing the distribution of the model's results are shown below in Table 10. For both dates, 75% of the model scores were below 2, a classification of low risk of impairment. Table 11 gives the quantity of segments with each level of risk, the total length of segments for each level of risk, and the change in quantities between 1986 and 2016. The distributions of risk scores for both dates can be visualized in boxplots shown in Figure 23.

Discussions
This study aimed to: (1) estimate and map the net impervious surface (IS) change in and around the City of Chattanooga, Tennessee, between 1986 and 2016 using satellite imagery from NASA's Landsat satellite missions, (2) perform a quantitative analysis of the obtained IS changes in relation to local water resources, and (3) develop a risk assessment model to identify the potential areas of concern for surface water quality in the study site due to the proximity and quantity of IS growth. All three tasks were successfully completed in a reasonable sense. However, there were challenges associated with each task that should be taken into consideration to evaluate the outcome of this study. The important aspects of the obtained results for each task along with the limitations and challenges encountered are discussed in the following sections.

Discussions
This study aimed to: (1) estimate and map the net impervious surface (IS) change in and around the City of Chattanooga, Tennessee, between 1986 and 2016 using satellite imagery from NASA's Landsat satellite missions, (2) perform a quantitative analysis of the obtained IS changes in relation to local water resources, and (3) develop a risk assessment model to identify the potential areas of concern for surface water quality in the study site due to the proximity and quantity of IS growth. All three tasks were successfully completed in a reasonable sense. However, there were challenges associated with each task that should be taken into consideration to evaluate the outcome of this study. The important aspects of the obtained results for each task along with the limitations and challenges encountered are discussed in the following sections.

Impervious Surface (IS) Mapping
This study shows that there has been an overall increase in the net IS development between 1986 and 2016 within the study site. This amount of IS growth fits the substantial growth in population, housing units, and economic activity in the area, as other research has found that increases in population can be related to proportionally larger increases in IS [155,156]. The larger proportional increase in IS area could be viewed as the result of residential and commercial development driven by the population growth in the greater Chattanooga areas. The contribution of both population and economic growth to IS growth are considered for the EPA's IS Growth Model [157]. The results are further supported by the EPA's IS Growth Model since the three primary contributors to IS cover are buildings, roads, and parking lots, listed in order of contribution [157].

Impervious Surface (IS) Classification Accuracy
The results of the confusion matrix accuracy assessments show that the classification of impervious and non-impervious surfaces using the model developed in this study was reasonably successful. The model's classification accuracy for pervious surfaces for both dates are much higher than impervious surfaces, which could result from the assignment of points using the random stratification method, which placed a much larger number of accuracy assessment points for assessing the pervious surfaces. The lower quantity of assessment points in impervious surfaces could, therefore, be more heavily affected by outlying classification errors. Differences in the User's Accuracy and Producer's Accuracy show the variation of different types of errors that the classification model experienced for each class. Having a higher Producer's Accuracy shows that the model was more likely to have more false positives and less false negatives. On the other hand, a higher User's Accuracy indicates the presence of the opposites [154]. In this study, for IS classification, both dates showed a higher Producer's Accuracy. This result indicates that the model had a difficulty to differentiate between subject's spectral responses effectively and consistently and, as a result, classified some non-urban areas as urban IS. The main source of confusion is believed to be derived from areas covered with open dry soils.
The zonal statistics accuracy assessments show that both dates had a similar performance in IS classification. This indicates the ability of the model to correctly classify moderately sized urban areas using both Landsat 5 TM and Landsat 8 OLI imagery.

Water Resource Proximity to Impervious Surfaces
The accuracy of the proximity of the impervious surface (IS) development to streams accuracy is dependent on the accuracy of the detected IS and the positional accuracy of the stream data, since the locations of the streams are assumed to be accurate due to the U.S. National Map accuracy standards met by the NHDplus dataset [68]. Errors from the IS detection portion of this study are believed to be the cause of the detected decrease in IS with 30 m of streams in the Chattanooga Creek and North Chickamauga Creek Lower watersheds and within 60 m of streams for the Chattanooga Creek watersheds. The results of the IS growth near streams also show that the growth is unequal across HUC-12 watersheds, thereby signifying that some may have increased impairment. This study shows that the Lower South Chickamauga watershed experienced the largest amount of urban growth between 1986 and 2016. This finding parallels with the finding that the Lower South Chickamauga watershed includes the largest amount of overall IS growth. Having the largest portion of IS development within the watershed and near the streams suggests that the Lower South Chickamauga watershed surface water quality could have the most noticeable and significant impairment relative to the other watersheds within the study site.
Differences in percent imperviousness between the three riparian zones show an average increase for all three zones with the increase being larger in the distant riparian zones. Shifting in the second and the third quartile percent imperviousness in the 2016 dataset compared to the 1986 dataset supported the conclusion that in 2016, stream riparian zones had experienced increases in IS development. The largest increases in the farthest zone suggests that the growth has begun to develop closer to streams. Increases in anthropogenic development and activity closer to streams increases the risk of potential impairment of stream surface water quality and riparian habitat [152]. Large portions of stream riparian zones are still present and not heavily disturbed within the study site. However, the findings of this research clearly shows the need for local policy action for riparian zone protection that could help mitigate the potential increase in stream impairment.

Impairment Risk Due to IS Development
The potential risk of stream impairment increased from 1986 to 2016. The increase in the total number of stream segments at risk being impaired suggested that stream surface water quality health in 2016 was at a higher risk than it was in 1986. Segments with extreme risk of impairment experienced the third largest increase between the two datasets. The increase in average risk can be attributed to the increase of risk for 7.64% of stream segments, where 987 stream segments have experienced increases in risk. This is coupled with the 27.75% of segments that have maintained risk of impairment since 1986. These findings show that there is a moderate proportion of streams within the study site that could have significant impairment due to long-term exposure or rapid introduction to urban development.
The risk model developed in this research showed that the locations of the potentially impaired stream segments due to impervious surface development are within the immediate riparian areas.

Limitations and Challenges
The accuracy assessment results for both datasets are at an acceptable level with regards to the effectiveness of the data and environmental conditions at the time of data acquisition. Moderate scale resolution sensors, such as the Landsat Thematic Mapper (TM) and Operational Land Imager (OLI), have been utilized for mapping impervious surfaces (IS), however, a consensus has been found that with data of this scale, it is not an effective ability to detect smaller areas of IS [158]. This study found that pixels determined as being IS are those where the majority to entirety of the land cover is impervious.
Environmental factors at the time of image acquisition are believed to have been an influence for both IS dataset generation. Water, shadow, and dry soil have been found to be difficult for many classification models due to the spectral confusion [159]. Shadows, which are present in both images, can be caused by the angle of the sun and/or sensor. The effect of shadowed areas was accounted for in the model and is not believed to have had any significant effect on the results of the study. For the 2016 image, a strong drought affected the study site and produced an increased amount of dry soils, river recession, and wildfires [160]. For the 1986 image, several clear-cut areas were seen in the study site exposing large areas of dry soil. The evaluation of the obtained IS data shows the presence of some dry and compacted soil, which were classified as IS. This probably stems from the conflicting spectral response of dry soil and IS (in some cases). However, it is believed that detection of dry soils as IS, in some cases, did not contribute to a major portion of detected areas [161,162].
The effects of noise and errors from the data, model, and environmental systems are also believed to have an influence on the results. The estimated area of IS for both dates is believed to be conservative, even with the detection of dry soil as it is inferred through the increase in housing units and population in the study site. There has been a net growth of suburban areas between 1986 and 2016, much of which is believed to have not been detected by this model based on visual inspection of the true color images. Suburban land cover can include IS such as roofs, drive-ways, and sidewalks. Although attempts were made to detect IS in suburban areas with the aid of true color images, it was not possible to accurately detect them due to sensor's coarser spatial resolution and erroneous detection of dry soil.

Conclusions
This study successfully utilized remote sensing and GIS technologies to estimate and map historical and current impervious surfaces (IS) of the greater Chattanooga area to determine its net spatial growth across seven HUC-12 watersheds and their relationship to the streams. IS cover was used to represent areas of anthropogenic development, specifically urban and suburban areas. Utilizing the multispectral satellite imagery acquired by the USGS Landsat 5 and 8 programs, IS detection was conducted using the Normalized Difference Vegetation Index (NDVI) and density slicing technique. This study found that there had been a net increase in IS within the study site with significant growth occurring near many of the streams. The IS change estimation showed that the overall growth was not equally distributed. Most of the IS development occurred in the Lower South Chickamauga Creek watershed. It was found that dry soil from transitioning land cover or drought caused erroneous detection of IS in some cases. Areas of suburban development are not believed to be fully mapped due to the spectral dominance of the pervious surfaces such as tree cover.
The results of the stream -IS proximity analysis show that there is an overall increase in percent imperviousness surface growth in the first three 30 m stream riparian zones in the study site. The largest increase in percent imperviousness of stream riparian zone occurred between 60 m to 90 m from streams, indicating that urban growth is beginning to encroach on critical, immediate riparian zones. The model for potential risk of stream surface water quality impairment reflects that there is an increase in risk for some portions of the streams due to riparian IS development. The HUC-12 watersheds in the study site directly feed into the Tennessee River, thereby increasing the possible impact of the land use and cover change for areas downstream of the study site. This study also found that decreases in pervious surfaces in stream riparian zones signaled a decrease in buffer capacity for filtering impaired surface and ground water. IS development detected within the watersheds could be the sources of potential new or continued surface water quality degradation. Finally, this study can conclude that between 1986 and 2016, there was an increase of at least 45.12 km 2 of IS cover in the seven HUC-12 watersheds within the greater Chattanooga, Tennessee, area and 9.96 km 2 of this growth were within close proximity of the streams, which should cause concerns for the local stream and river surface water quality.
Author Contributions: J.H. processed all the GIS and satellite imagery, analyzed the results, wrote major portion of the manuscript. A.K.M.A.H. provided expertise in remote sensing and GIS to carry out the data processing and analysis. He also wrote part of the manuscript and supervised the project. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded in part by the University of Tennessee at Chattanooga, Department of Biology, Geology and Environmental Science.

Acknowledgments:
The authors thankfully acknowledge the contribution of Mark Schorr and Hong Qin for their valuable suggestions and comments during the entire length of the project. The authors also would like to acknowledge the useful advice given by Kyle Jones, Charlie Mix and the rest of the IGT Lab GIS technicians.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A
their valuable suggestions and comments during the entire length of the project. The authors also would like to acknowledge the useful advice given by Kyle Jones, Charlie Mix and the rest of the IGT Lab GIS technicians.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. Appendix A Figure A1. Percent impervious cover within 30 meters of streams in 1986. Streams are visualized in segments generated for the stream impairment risk model. Percent imperviousness was calculated by