Assessing the Effect of Land-Use and Land-Cover Changes on Discharge and Sediment Yield in a Rural Coal-Mine Dominated Watershed in Kentucky, USA

: The Appalachian Mountain region of eastern Kentucky is unique and contains high proportions of forestland along with coal and natural gas depositaries. Landscape changes due to extreme mining activities can eventually threaten the downstream ecosystems, including soil and water quality, resulting in excessive runoff and sedimentation. The purpose of this study is to assess the impacts of land-use and land-cover (LULC) changes in streamﬂow and sediment yield in Yellow Creek Watershed, Kentucky, USA, between 1992 and 2016. LULC, digital elevation model, soil, and weather data were inputted into the Soil and Water Assessment Tool (SWAT) to simulate discharge and sediment yield. The model output was evaluated on several statistical parameters, such as the Nash-Sutcliffe efﬁciency coefﬁcient ( NSE ), RMSE-observations standard deviation ratio (RSR), percent bias ( PBIAS ), and the coefﬁcient of determination ( R 2 ). In addition, two indices, P -factor and R -factor, were used to measure the prediction uncertainty. The calibrated model showed an increase in surface runoff and sediment yield due to changes in LULC in the Yellow Creek Watershed. The results provided important insights for studying water management strategies to make more informed land management decisions and adaptive practices.


Introduction
Land-use and land-cover (hereafter LULC) change is a phenomenon that has been taking place since the beginning of humankind. It has long been considered a factor affecting the global environment [1,2]. Numerous studies over the past decade have highlighted the role of LULC change on various aspects of the environment, such as hydrological response [3][4][5], climate change [6][7][8], ecosystem services [9,10], and food security [11,12]. Among the studied impacts of LULC change, hydrological responses and accompanying factors, such as water quality and quantity, have received a great deal of attention. The complexities in the LULC change process and the associated hydrological cycle vary across geographic regions [1]. Changes in LULC in the watershed alter hydrological responses, including increased runoff, high streamflow, soil erosion, and inferior water quality [4,5,13,14]. Moreover, studies have shown wide variations in hydrological responses due to variations

Description of the Study Region
The Yellow Creek Watershed falls within the Upper Cumberland River basin. The watershed is primarily located in Bell County, Kentucky, but extends into Claiborne County, Tennessee. The outlet of this watershed corresponds to the point where a monitoring station (USGS 03402000) is installed, which is at Yellow Creek near Middlesboro (latitude: 36 • 40 05 , longitude: 83 • 41 19 ) (Figure 1). From the outlet, it derives a catchment area of 157.18 m 2 . Except for the broad, alleviated valley of Yellow Creek at Middlesboro, most of the watershed is located in the rugged mountainous region with an elevation of 335 m to 960 m above sea level. The nearest United States Weather Bureau station is located at Middlesboro, at 358 m.a.s.l elevation. During the 1987-2016 period, the mean annual precipitation was 1283.64 mm. The mean annual temperature was 13.4 • C, with a mean winter minimum of −3 • C and a mean summer maximum of 24 • C. According to the Cropscape-Cropland data of 2016, the watershed was majority forested (76.14%), followed by urban (15.37%), pasture (6.11%), and barren land (2.0%).

The SWAT Model
The SWAT model is a physically based, watershed-scale simulation model jointly developed by the USDA Agricultural Research Service and Texas A&M AgriLife Research, part of the Texas A&M University System [22,23]. In SWAT, a watershed is divided into a number of sub-watersheds or sub-basins, which are further divided into hydrologic response units (HRUs) [23,47]. HRUs are lumped areas within the sub-basin that are comprised of unique land cover, soil, and slope. The model estimates relevant hydrologic components such as evapotranspiration, surface runoff, groundwater flow, and sediment yield for each HRU [23]. The hydrologic cycle simulated by SWAT is based on the water balance equation [23,47]: where SW t is the final soil water content on day i, SW 0 is the initial soil water content on day i, t is the time in days, R day is the amount of rainfall on day i, Q surf is the amount of surface runoff on day i, E a is the amount of evapotranspiration on day i, W seep is the amount of water entering the vadose zone from the soil profile on day i, and Q gw is the amount of base flow or return flow on day i. The surface runoff was estimated by using the United States Department of Agriculture (USDA) Soil Conservation Services (SCS) runoff curve number (CN) approach [22,23,48]. The CN is the function of land use, soil permeability, and antecedent soil moisture condition. The SCS equation [48] is: (R day +I a +S) (2) where Q surf is the accumulated runoff or rainfall access, R day is the rainfall depth of the day (mm H 2 O), I a is the initial abstractions, which includes surface storage, interception and infiltration prior to runoff (mm H 2 O), and S is the retention parameter (mm H 2 O). The retention parameter (S) is dependent spatially to changes in soils, land use, management, and slope, and temporarily to changes in soil water content [23]. The retention parameter (S) [23] is given by the equation where CN is the curve number for the day. The initial abstraction, I a , is usually approximated as 0.2S, so equation 2 can be rewritten as The SWAT offers three methods to estimate potential evapotranspiration (PET): the Penman-Monteith method, the Priestley-Taylor method, and the Hargreaves method [23]. In our study, PET was estimated using the Penman/Monteith method, which requires solar radiation, air temperature, relative humidity, and wind speed as inputs. The Penman-Monteith equation [49] is where λE is the latent heat flux density (MJ m −2 d −1 ), E is the depth rate evaporation (mm d −1 ), ∆ is the slope of the saturation vapor pressure-temperature curve, de/dT (kPa • C -1 ), H net is the net radiation (MJ m −2 d −1 ), G is the heat flux density to the ground (MJ m −2 d −1 ), ρ air is the air density (kg m −3 ), c p is the specific heat at constant pressure (MJ kg −1 • C −1 ), e 0 z is the saturation vapor pressure of air at height z (kPa), e z is the water vapor pressure of air at height z (kPa), γ is the psychometric constant (kPa • C −1 ), r c is the plant canopy resistance (s m −1 ), and r a is the diffusion resistance of the air layer (aerodynamic resistance) (s m −1 ) [23]. Upon determination of PET, the actual evapotranspiration (ET) was calculated with the method similar to that suggested by Richtie [50].
The sediment component for each HRU was calculated using the modified universal soil loss equation (MUSLE). The MUSLE [51] is sed =11.88 ·(Q surf · q peak · area hru ) 0.56 · K USLE · C USLE · P USLE ·LS USLE · CFRG (5) where sed is the sediment yield on a given day (metric tons), Q surf is the surface runoff volume (mm H 2 O/ha), q peak is the peak runoff rate (m 3 /s), area hru is the area of HRU (ha), K USLE is the USLE soil erodibility factor (0.013 metric ton m 2 hr/m 3 -metric ton cm), C USLE is the USLE cover and management factor, P USLE is the USLE support practice factor, LS USLE is the USLE topographic factor, and CFRG is the coarse fragment factor.   Table A1 for a complete description of SSURGO soil classes). Table 1 highlights the source and description of the data used in the SWAT model. The weather data was downloaded from PRISM Climate Group, which includes maximum and minimum temperature and precipitation. Wind speed, relative humidity, and solar radiation were simulated for the nearest weather station using the weather generator in SWAT. Monthly discharge (1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004) and sediment data (1987)(1988)(1989)(1990)(1991)(1992) were obtained from the gauging station (USGS 03402000) at the outlet of the watershed. Since sediment data was only available for parts of the simulation periods (until 1992), we used a LOAD Estimator (LOADEST) to simulate sediment data to match the period with discharge (1987-2004). The LOADEST is a FORTRAN-based model developed by the U.S. Geological Survey [52] to estimate the constituent loads in streams and rivers given a time series of stream flows, additional data variables, and constituent concentrations. It uses a regression model to estimate the constituent loads. Streamflow, the time factor, and the user-specified variables are the independent variables of the regression model which is used to estimate load over the user's specified time interval. The regression model is used to estimate mean load and the statistical properties: standard errors and 95 percent confidence interval on a monthly basis.

Hydro-Climatic Data
LOADEST uses three statistical estimation methods for the calibration [52]. The first two methods, adjusted maximum likelihood estimation (AMLE) and maximum likelihood estimation (MLE), are useful when a model error (residuals) is normally distributed, whereas the third method, least absolute deviation (LAD), is useful when residuals are not normally distributed. The valuable information inside the LOADEST output, like diagnostic tests and warnings, helps users select the appropriate estimation method and interpret the estimated loads. The stream flows and their constituent are used to calibrate LOADEST to determine the best-preloaded models inside the LOADEST for load estimating.

SWAT Model Setup
ArcSWAT 2012, an ArcGIS extension and interface for SWAT, was used to simulate the SWAT model. ArcSWAT interface is a public domain software, and as such may be used freely [47]. ArcSWAT extension for ArcGIS and documentation are available as downloads (https://swat.tamu.edu/software/arcswat/ accessed on 20 December 2020). Based on ArcSWAT, the SWAT model data preparation is characterized by three modules: (1) the SWAT Watershed Delineator, which allows the users to discretize the watershed and subwatersheds using the data derived from the digital elevation model (DEM); (2) the SWAT HRU Analysis Tool, which combines data derived from LULC, soil characteristics, and slope variations to discretize the HRUs; and (3) the SWAT Input Editor, which allows the user to create an input database and modify various model parameters.
A threshold critical source area of 350 ha (between the suggested range by ArcGIS watershed delineation interface) was used for this study, which divided the watershed into 21 sub-basins ( Figure 3a). It was based on the understanding that a smaller threshold generates a denser stream network in the watershed [53]. The LULC maps of 1992 (NLCD) and 2016 (CropScape-Cropland Data) were aggregated into seven LULC types: Water, Urban, Barren, Forest, Pasture, Agriculture, and Wetland. SSURGO, a higher resolution soil map, defined 38 different soil types in the watershed. The slopes were classified into four classes: 0-15%, 15-30%, 30-45% and >45%. We opted out threshold (LULC/Soil/Slope: 0/0/0[%]) in HRU definition, which resulted in the further partition of sub-basins into 1847 HRUs (Figure 3b).

SWAT-CUP Premium and SWAT Parameter Estimator (SPE) Algorithm
We used SWAT-CUP Premium (SWAT-CUPP), a computer program developed for calibration of the SWAT model (https://www.2w2e.com/home/SwatCupPremium; accessed on 17 February 2021). SWAT-CUPP is an improved version of SWAT-CUP, which allows for behavioral and multi-objective calibration [54].
SWAT-CUPP offers two algorithms, SWAT Parameter Estimator (SPE) and Particle Swarm Optimization (PSO). We used the SPE algorithm (previously Sequential Uncertainty Fitting (SUFI-2)) for model sensitivity analysis, calibration, uncertainty analysis, and validation [5,13]. In SPE, the algorithm maps all uncertainties (parameter, conceptual model, input, etc.) on the parameters (expressed as uniform distributions or ranges) and tries to capture most of the measured data within the 95% prediction uncertainty (95PPU) of the model in an iterative process determined at the 2.5% and 97.5% levels of cumulative distribution of output variables obtained through Latin hypercube sampling [54,55]. Users are provided with several choices of objective function (11 functions including the Nash-Sutcliffe Efficiency coefficient (NSE), the coefficient of determination (R 2 ), and percent bias (PBIAS)). We selected NSE as our objective function in this study because it is recommended to be the best objective function for reflecting the overall fit of a hydrograph [56,57]. NSE is a normalized statistic that determines the relative magnitude of the residual variance compared to the measured data variance [58]. NSE is computed as where Y obs i is the ith observation for the constituent being evaluated, Y sim i is the ith simulated value for the constituent being evaluated, Y mean i is the mean of observed data for the constituent being evaluated, and n is the total number of observations [59]. The details of the SWAT-CUPP and SPE algorithms can be found in the SWAT-CUPP user manual [54].

Model Sensitivity Analysis
Sensitivity analysis is the foremost step in the calibration and validation process in SWAT-CUPP, which determines the most sensitive parameters for a given watershed [55]. SPE is embedded with two systems for sensitivity analysis (i.e., one-at-a-time and global). This study utilizes both: one-at-a-time sensitivity analysis to identify the sensitive parameters and global sensitivity analysis to define the rank of the model parameters.
One-at-a-time sensitivity analysis was conducted for each parameter while keeping other parameter values constant [54]. In global sensitivity analysis, a multiple regression system regresses the Latin hypercube-generated parameters against objective function values to determine the sensitivity of the parameters [54]. The sensitivity is calculated according to where g denotes the objective function, b is the parameter, α is the regression constant, β corresponds to the technical coefficient attached to the variable b, and m is the number of parameters [60].
In addition, statistical measurements t-stat and p-value were used to identify the sensitive parameters, with larger t-stats and smaller p-values representing greater sensitivity. Based on one-at-a-time sensitivity analysis, Abbaspour et al. [60], and the relevant literature [45,47,61], a total of 30 parameters (see Table A2) were selected for global sensitivity analysis. To perform global sensitivity analysis, we ran our first iteration among 23 selected parameters for discharge. Our first iteration was run with 600 simulations. In our analysis, the larger the value of t-stat (in absolute value) and the smaller the p-value (p < 0.05), the more sensitive the parameters were [54]. Once the satisfactory calibration performance was obtained for discharge, sensitivity analysis was carried out for the sediment parameters with a similar approach. The initial parameter ranges were set according to the SWAT manual and available guidelines [23,47,55] (see Table A2).

Model Calibration and Uncertainty Analysis
The second step is model calibration, in which parameter value ranges are adjusted to improve the fit between model predictions and real-world observations [47]. SWAT-CUPP allows users to select specific model parameters for auto-calibration within their acceptable ranges (see Table A2) and execute hundreds of SWAT simulations to find the optimal set of parameter values that minimizes the error between model predictions and observed data [54,55]. During calibration of parameters, only sensitive parameters ( Table 2) were calibrated based on the results of the sensitivity analysis in SWAT-CUPP.
Calibration is characteristically subjective and, therefore, closely linked to model output uncertainty [24,55]. The term uncertainty analysis refers to the propagation of all model input errors and uncertainties to model outputs. In SWAT-CUPP, two indices referred to as "P-factor" and "R-factor" are used to quantify the fit between simulation results [24,47,55,60]. The P-factor is the fraction of measured data (plus its error) bracketed by the 95PPU band, and the average distance d between the upper and the lower 95PPU (or the degree of uncertainty) estimated from where k is the number of observed data points, l is a counter and X L (2.5th) and X U (97.5th) represent the lower and upper boundaries of the 95PPU [54,55,60]. The P-factor varies from 0 to 1, where 1 indicates 100% bracketing of the measured data within model prediction uncertainty (i.e., a perfect model simulation considering the uncertainty). The R-factor is the ratio of the average width of the 95PPU band and the standard deviation of the measured variable. The R-factor is estimated by where σ X is the standard deviation of the measured variable X. The P-factor value of >0.7 (at least 70%) and R-factor value of <1.5 (at most 150%) are desirable and recommended to be adequate [60,62]. However, depending on the project scale and availability of input and calibration data, these requirements can be adaptable [24,60,62]. In general, a larger P-factor can be achieved at the cost of a larger R-factor. A balance between the two must, therefore, always be considered. When appropriate values of R-factor and P-factor are achieved, the parameter ranges are taken as the calibration parameters in the final iteration. However, it is important to understand that there is no single good solution representing model output but rather an envelope of good solutions expressed by the 95PPU, with certain parameter ranges [24].

Model Validation
After calibration, the model is validated with the calibrated parameter ranges and by comparing predictions to additional observed data, which is not used for calibration. Based on the level of agreement between predictions and the additional observations, the model is validated for further use, or model inputs and parameters are revisited for further calibration [47,56]. Similar to calibration, validation results are also quantified by the P-factor and R-factor.
The initial model obtained from ArcSWAT was divided into two periods, 1990-1997 for calibration and 1998-2004 for validation in SWAT-CUPP. Discharge was calibrated first since it is the primary controlling variable [54]. After running the one-at-a-time sensitivity analysis and following the literature [5,24,45,61], the model was parameterized, and ranges were assigned. The model was run for three iterations (600 simulations each) for calibration. After each iteration, new parameter ranges were suggested, which were used for another round of iteration. Finally, the model was run for another iteration (600 simulations) with the calibrated parameter ranges for validation. A similar procedure was applied for the calibration and validation of sediment.
Parallel processing was utilized to speed up the calibration process. The parallel processing module utilized all eight CPUs where, for each iteration, 600 simulations were divided into eight simultaneous runs of 75 each per CPU. This substantially improved the runtime of the calibration and validation process [54].

Model Performance Indices
In addition to two indices used for prediction uncertainty analysis, P-factor and R-factor, multiple indices are made available to check the performance of the SWAT model [59]. In this study, the Nash-Sutcliffe efficiency coefficient (NSE), RMSE-observations standard deviation ratio (RSR), percent bias (PBIAS), and the coefficient of determination (R 2 ) were used to evaluate calibration and validation performance of the SWAT model. These are among some of the most widely used indices. NSE indicates how well the observed versus simulated data plot fits within the 1:1 line [59]. NSE ranges between −∞ and 1, with NSE = 1 being the optimal value. RSR is the ratio of RMSE and the standard deviation of observed value, with a lower value being the better model simulation performance [59]. PBIAS, with an optimal value of 0, evaluates the average tendency of the simulated data to be greater or lesser than their observed counterparts [63]. The low values of PBIAS imply accurate model simulation, positive values indicate model underestimation, and negative values indicate model overestimation bias [59]. The proportion of the variance in measured data explained by the model is described by R 2 [63]. The value of R 2 ranges from 0 to 1, with higher values indicating less error variance [59,64]. More details on the thresholds for evaluation of model performance can be found in Moriasi et al. [59,65] (see Table A3).

Model Calibration, Uncertainty Analysis, and Validation
After performing three iterations (600 simulations in each) with modifications to the parameters, the model was observed to have a good fit between the observed and simulated discharge and between the LOADEST-simulated and SWAT-simulated sediment. The statistical results of the model performance for discharge and sediment during both calibration and validation periods are summarized in Table 3. The calibration and validation of discharge and sediment are illustrated as a 95PPU graph in Figures 4 and 5, respectively. Graphical results during both calibration and validation ( Figure 4) indicated adequate performance over the range of discharge. The calibrated model achieved both the P-factor and R-factor of 0.58. This explains that 58% of the measured discharge was bracketed by 95PPU within model prediction uncertainty. The 95PPU bracketed about 62% of the measurement in the validation period, which is higher than the calibration period. Similarly, the R-factor had the desired values of 0.58 and 0.69 (R-factor < 1.5). Overall, the NSE values (0.87 and 0.84) indicate a good match between measured and simulated monthly discharge for the calibration and validation period. The RSR values were 0.36 and 0.40 during the calibration and validation period. These values indicated a very good model performance for residual discharge variation. The R 2 (0.89 and 0.86) showed a strong linear correlation between measured and simulated monthly discharge. PBIAS values were 2.8 and -0.4 for calibration and validation, respectively. These values were also within the "very good" range (PBIAS < ±5) according to the guidelines [59,65].  The 95PPU bracketed about 55% of the observations in the case of sediment load calibration ( Figure 5). Additionally, NSE (0.86), RSR (0.38), and R 2 (0.88) were excellent statistics. However, the PBIAS value (−18.8) could only be of a satisfactory level for the calibration period. In comparison to a calibration period, the validation period ( Figure 5) reached a lower P-factor (0.46). However, considerably impressive results were recorded for NSE (0.81), RSR (0.44), and R 2 (0.84). Similar to calibration, the validation period also achieved a satisfactory result with PBIAS (−15.9). In sediment load calibration and validation, the model achieved R-factor values of 3.15 and 4.19, respectively. Overall, it can be evaluated as a good model based on the criteria by Moriasi et al. [59,65].

Impact of LULC Change on the Water Balance Components in Yellow Creek Watershed
The annual summaries of water balance components between 1992 and 2016 are listed in Table 4. Results show that total yearly precipitation, lateral flow, and groundwater experienced a decrease of 2.76%, 45.3%, and 22.43%, respectively, whereas surface runoff increased by 66.85% from 1992 to 2016. Loss of water due to percolation accounted for about a 15.43% decrease from 1992 to 2016. There was a slight increase in evapotranspiration (0.32%), while the potential evapotranspiration witnessed a 20.68% increment from 1992 to 2016. Overall, the total annual water yield decreased from 462.37 mm in 1992 to 424.19 mm in 2016. Between 1992 and 2016, sediment yield increased from 5.34 t/ha to 14.66 t/ha.

Sediment Yield and LULC Change in Yellow Creek Watershed
The spatial distribution of sediment yield in 1992 and 2016 within the Yellow Creek Watershed is illustrated in Figure 8. In 1992, the total sediment yield ranged from 0.14 t/ha to 15.63 t/ha, where sub-basin 21 recorded the highest sediment yield, followed by subbasin 18. In 2016, the sediment yield increased significantly, ranging from 0.09 t/ha to 42.29 t/ha. Sub-basins 8,9,18,20, and 21 yielded the highest sediment (>15.64 t/ha). Except for sub-basins 3, 4, 5, 11, 12, 14, and 15, all sub-basins experienced an increase in sediment yield.
As illustrated in Table 5, a comparison between 1992 and 2016 reveals that all LULC classes have undergone some changes. In the last 25 years, the changes in each LULC ranged from 0.05% to 9.37% of the total watershed area. The highest percentage of change occurred in forested areas, by a total of 14.74%. The change can be attributed to the conversion of forest area into urban development, pasture, and barren land, which increased by 9.43%, 5.02%, and 1.01%, respectively. Similarly, agriculture and water also declined by 0.45% and 0.18%, respectively. The wetland area declined by 0.079%, which only covered 0.02 ha in 2016.  Furthermore, the spatial distribution of LULC classes in 1992 and 2016 is illustrated in Figure 9. The northwestern region of the watershed experienced the most intensive changes in forest cover. Sub-basins heavily forested in 1992 converted into barren lands, pastures, and urban developments in 2016. Overall, the change in forest cover is more or less found to be distributed in all sub-basins. Conversion of LULC in the form of barren land in 2016 is primarily attributed to the development of coal-mining lands, which can be seen in sub-basins 7, 8, 9, and 10 ( Figure 9). As evident in Figure 8, these sub-basins contributed the highest sediment yield in 2016. Moreover, the increase in sediment yield seems to be mainly associated with sediment transport from barren lands (mining fields), urban development activities, and pasturelands [9,[66][67][68]. This implies that the sediment yield increased due to the direct and indirect consequences of LULC changes in the watershed. Figure 9 also shows barren land areas in 1992 (Sub-basins 1, 3, & 19) converted into other classes. This may be the result of various post-mining reclamation practices that have been implemented in the watershed [69,70]. Additionally, some forest-dominated regions (sub-basin 20) witnessed an increase in sediment in 2016. This might be attributed to comparatively steeper slopes in some regions. As seen in this study, sub-basins with the greatest sediment yield had a slope of more than 30% (Figures 2 and 8). Furthermore, studies have demonstrated a significant correlation between sediment yield and slope gradient [71,72].

Discussion
The resulting statistics of our study were acceptable based on the guidelines [59,65], and were similar to values found in other studies conducted in the United States [5,13,45,[73][74][75][76][77] and around the world [4,[45][46][47][48][49][50]. Spruill et al. [45] reported the SWAT model as an effective model for simulating monthly runoff in a small watershed in central Kentucky, with NSE values of 0.89 (calibration) and 0.58 (validation). Coffey et al. [78] stated similar results for the same watershed. Shrestha et al. [5] used SWAT to study industrial wood pellet impact on streamflow in Georgia. In the latter study, predicted values compared well with observed values, with NSE and R 2 values above 0.75, and predicted streamflow 3.57-7.28% lower than the observed streamflow. In their study in North Carolina, Ayivi and Jha [79] reported a good agreement for discharge, with both NSE and R 2 values greater than 0.70. In another study conducted by Tadesse et al. [73] in a watershed located in Tennessee and Alabama, the model adequately simulated sediment with R 2 and NSE values of 0.71 and 0.70, respectively. Compared to these studies, our model achieved adequate statistics, which indicates a much better correlation with NSE values higher than 0.80 for both calibration and validation periods.
The results from this study imply that LULC modification might have a significant impact on the annual water balance in the watershed. The decline in percolation, lateral flow, and groundwater might be mainly due to changes in LULC and altered soil properties due to coal-mining activities in the watershed. Similarly, increased runoff indicated the lower infiltration capacity of the surface. Evidence from previous studies also suggests the historical LULC change impact in basin-scale water balance, such as with the Raccoon River Watershed in Iowa [13], the Little River Watershed in Tennessee [3], the Kentucky River Basin [46], and the Reedy Fork-Buffalo Creek Watershed in North Carolina [79].
Our findings are consistent with previous studies that suggest increasing sediment yields, water quality degradation, and increasing flood events pertaining to LULC changes [9,73] caused by the clear-cutting of forests. The increase in sediment yield indicates an impact of LULC changes in the Yellow Creek Watershed. As seen, sub-basins dominant with barren lands and pastures contributed most in sediment yield. It may be due to mining activities and emerging pastures in the reclaimed mine lands [9,69,70]. According to the report prepared by the Appalachian Regional Commission, coal production has fallen by more than 45% in Appalachia between 2005 and 2015 due to the depletion of high-grade coal seams [80]. However, a recent study showed that the cumulative mining area is still increasing in central Appalachia, though at a significantly lower growth rate [9,81]. The continuous growth of cumulative mines indicates that mining might still be a key driving factor in the LULC transformation in some Appalachian Mountains areas [9,69,70].
Several studies show higher peak and total storm runoff from mined lands compared to forested lands in the Appalachian region [82,83]. These results suggest that surface mining may increase the risks of more flooding hazards and recommend more studies for quantification and evaluation of the effects of mined land conversions in the region. Knowledge of the extent of mining is critical to managing or mitigating the potential impacts of surface mining on sedimentation.
SWAT is an effective and widely accepted tool in predicting the impact of LULC changes [23,47]; however, success of SWAT depends on the quality of input data, such as the resolution of the DEM, LULC data, and soil [23,84]. In our study, given the availability of data, we were limited to moderate resolution data (30 m × 30 m) such as NLCD LULC data. Precipitation is a crucial factor affecting runoff events and sediment export. In SWAT, precipitation is simulated in each sub-basin according to the nearest gauge to the sub-basin centroid [85]. Having a sufficient number of gauges in a watershed is, therefore, incredibly beneficial for hydrological simulations in SWAT. However, our analysis was restricted by the absence of ground observation data, and substitute methods such as PRISM climate data were utilized.

Conclusions
This study was successful in simulating sediment yield and quantifying the impact of LULC change on the hydrology of the watershed. The results showed a very good agreement between observed and predicted discharge and sediment at the outlet. The results showed that the Yellow Creek Watershed experienced LULC changes over 25 years between 1992 and 2016. It consisted primarily of a decrease in forest cover which mostly converted into coal-mine lands, urban development, and growth in pastureland. Similarly, the comparison of the distribution of water balance components, sediment yield, and LULC changes in the watershed implied that LULC changes, specifically ongoing surface mining activities, increasing pastureland, and urban development, have contributed to augmenting the sediment yield in the watershed.
This study contributes to studying land management practices in watersheds impaired by coal-mine operations. Addressing impacts of LULC change issues strategically and timely is essential for the effective management of coal-mine lands and water resources. Identifying critical areas and selecting best management practices (BMPs) in watershed scale are necessary for reducing sedimentation. The BMPs in watershed scale, for example, should focus on approaches for reducing the environmental footprints of mining by implementing reclamation practices that are feasible in the watershed.
Despite model and data limitations, this study fills the gap that exists in modeling LULC impacts on water and land management in the coal-mine lands of the Appalachian region in Kentucky. The findings presented provide a plethora of information on LULC changes and their impact on hydrology. The maps and tables produced in this study provide essential spatial information on the sub-basin scale to researchers in identifying and implementing appropriate land management practices.
There is a dearth of systematic research and data to quantify the effects of coal-mine lands on increased surface runoff and sediment yield. This research is a case study modeling impacts of LULC on runoff and sediment yield. We recommend continuing similar studies in the other watersheds to discern the role of change in LULC on runoff and sediment yield using high-resolution land cover data and integrating such data with human populations and infrastructure changes such as access roads to coal mine areas. Future studies should include local rainfall characteristics and evapotranspiration data. Rainfall patterns and their impacts on coal-field areas will be important to study in this region.
In future studies, this approach can further be utilized in other similar watersheds to explore the impacts of surface-mining activities on the hydrological properties and to evaluate the effectiveness of best management practices (BMPs) in managing sediment yields for environmental stewardship at the watershed level. We will expand this research to include other watersheds from the Appalachian region to explore the impacts of land-use change on water discharge, as well as the effects of current agricultural and forestry-related BMPs in sediment control and management.

Data Availability Statement:
The data used to support the findings of this study are available from the corresponding author upon request.