Assessing Ecosystem Services from the Forestry-Based Reclamation of Surface Mined Areas in the North Fork of the Kentucky River Watershed

Surface mining is a major driver of land use land cover (LULC) change in many mountainous areas such as the Appalachian region. Typical reclamation practices often result in land cover dominated by grass and shrubs. Assessing ecosystem services that can be obtained from a forest landscape may help policy-makers and other stakeholders fully understand the benefits of forestry-based reclamation (FRA). The objectives of this study are to (1) identify how surface mining and reclamation changed the LULC of a watershed encompassing the north fork of the Kentucky River, (2) assess the biophysical value of four major ecosystem services under the contemporary LULC condition, and (3) assess the benefits of the FRA scenario in the provision of ecosystem services. Geographic Information System (GIS) was used to study the LULC change and InVEST software models for ecosystem services assessment. The results indicate that watershed’s forest area has decreased by 7751 hectares from 2001 to 2011 and mining/reclamation activities may have contributed 65% of the overall changes in LULC. Barren and grassland land covers provide less carbon storage, yield more water, and export more sediments and nutrients than forests. At the watershed level, the FRA scenario increased carbon storage (13%) and reduced water yield (5%), sediment export (40%), and nutrient export (7%). The provision of these ecosystem services varies at the subwatershed level, and such spatial heterogeneity is primarily driven by land cover composition, precipitation, and topography. This study provides critical information regarding the ecological benefits of restoring mined land to assist policy and decision making at landscape scales.


Introduction
Land cover and land use (LULC) change is a major force of global environmental change [1], with as much as 50% of the earth's ice-free land surface having been transformed as a direct consequence of land uses [2,3].The changes of LULC help humans extract natural resources to meet society's growing needs of food, fiber, water, and shelter, but often at the expense of degrading environmental conditions [4,5], such as alteration of the climate system and atmospheric composition [6], land degradation [7,8], changes in hydrology [9], and loss of biodiversity [4].Understanding the pattern of LULC change and its impacts for human well-being is a key in land change science as a coupled system of human and environment [3].Although studying the cause and consequences of LULC change can be conducted at any spatial scales, it has been increasingly recognized that landscapes represent a pivotal scale domain for the research of LULC change because landscape scales are more directly linked with the land management, environmental policy, and application of sustainability [10].
Among all LULC change types, few are as intensive as surface mining, which extracts minerals near the Earth's surface [11,12].Surface coal mining generally involves a sequence of operations including vegetation clearing, topsoil removal, drilling, and blasting the hard strata over the coal seam, and then the subsequent extracting and transporting of coals [11,13].It has been predominantly conducted in the central Appalachian Mountains of the US that are mainly located in southern West Virginia, eastern Kentucky, southwestern Virginia, and northeastern Tennessee [14].Various studies have documented persistent negative impacts of mining on ecological integrity of Appalachia [14][15][16].One of the most prominent environmental impacts is the degradation of water quality in Appalachia as the overburden from surface mining has buried more than 2000 km of stream channels [15] and permanently altered drainage networks and topography [17,18].Forest clearing, soil loss, and subsequent substrate compaction induced by reclamation further contribute to water quality deterioration and the increased flood risk in downstream communities [16,17].Surface mining can also convert an area that was a carbon sink into a carbon source through the depletion of forest carbon storage and the burning of coal in electric power plants [14,19].
The Surface Mining Control and Reclamation Act of 1977 (SMCRA) was initiated to mitigate the negative environmental effects of coal mining in the United States of America.SMCRA requires reclamation of mountaintop-mined sites to a state that provides an equal or better use than the pre-mining condition.However, the law is vague on what constitutes equal or better use.Often the reclamation approach has resulted in plant communities dominated by persistent herbaceous species, grasses sown during reclamation, and early successional woody species [20], which is not an adequate substitute for the original diverse forest [21].To address this issue, the Appalachian Regional Reforestation Initiative (ARRI) was formed in 2004 by a coalition of citizens, government officials, and coal industry representatives dedicated to restoring forests in the abandoned coal mines [20,22,23].ARRI advocates a technique known as the Forestry Reclamation Approach, or FRA, a series of recommendations to guide successful regeneration of native forest on surface mined sites [24].The five-step guideline of FRA includes: creating a suitable growth medium, grading the top soil to get a noncompacted growth medium, planting less competitive ground covers that are compatible with trees, planting a mix of early successional woody species for wildlife and soil stability and commercially valuable crops, and using proper tree planting techniques to accommodate the seedling's root system [20,23].Experimental reclamation trials utilizing FRA techniques have been successful in Kentucky [24] and West Virginia [25].However, the reclamation approach is not in widespread use in the Appalachian region because its implementation is difficult and expensive for mining companies [23].In addition, local residents and the public are nonchalant toward forestry-based reclamation practices.Such situations might have arisen because the value of ecosystem services from forests in this landscape are not adequately valued.The public, mining companies, and policy-makers are not fully aware of the extent of ecosystem degradation induced by surface mining and the benefits to human wellbeing brought by the FRA approach compared to traditional reclamation practices.
Ecosystem services evaluation can help the decision making and the implementation of FRA techniques in reestablishing forest patches in the Appalachian Mountains.Evaluation of ecosystem services such as carbon storage, flood control, sediment, and nutrient retention at a local watershed or regional level can help establish sound ecological restoration policies [26] because it can help individuals and stakeholders appreciate and capture the value of ecosystem services to produce different outcomes [27].Ecological restorations can then in return help in the increased provision of biodiversity and ecosystem services [28].
Ecosystem services can be evaluated either in biophysical terms or monetary terms [29].In this study, a quantitative biophysical evaluation is used for multiple ecosystem services.The overarching study objective is to assess the value of major carbon-and water-related ecosystem services in the North Fork Kentucky River (NFKR) watershed of Kentucky, which is a watershed in Central Appalachia that has been severely impacted by surface coal mining [14,30], and evaluate the potential ecological service benefits of FRA reclamation at the landscape level.This study has three specific objectives: (1) Examining how LULC has changed in NFKR watershed from 2000s to 2010s and the contribution of mining and reclamation to the overall landscape change, (2) quantifying the value of major regulating ecosystem services (carbon storage, flood control, sediment retention, and nutrient retention) under contemporary LULC conditions, and (3) assessing the changes in the provision of these ecosystem services under a FRA scenario.Based on prior knowledge, it is hypothesized that (H1) there are considerable transitions from forest to barren and grasslands attributed to mining and reclamation in the NFKR watershed; (H2) the barren and grassland land covers provide the lesser amount of these ecosystem services than the forest land cover and consequently the subwatersheds with less mining/reclamation activities have higher mean ecosystem service provision; and (H3) FRA will produce greater ecosystem services compared to current reclamation practices, but there will still be considerable spatial variability of ecosystem service provision at the subwatershed level driven by LULC composition, topography and climate.

Study Area
The North Fork Kentucky River (NFKR) watershed is situated in the Eastern Kentucky Coal Field physiographic region, which is part of a larger physiographic region, the Cumberland Plateau [30].The watershed occupies most of Letcher, Perry, Knott, and Breathitt counties and a small portion of Lee County (Figure 1).These counties are part of the 65 counties in Central Appalachia where surface mining has been concentrated [14].The welfare of many residents in Central Appalachia have been affected by surface mining in the past [31].The NFKR watershed ranks among the groups with highest need for protection and restoration [30].
Forests 2018, 9, x FOR PEER REVIEW 3 of 23 potential ecological service benefits of FRA reclamation at the landscape level.This study has three specific objectives: (1) Examining how LULC has changed in NFKR watershed from 2000s to 2010s and the contribution of mining and reclamation to the overall landscape change, (2) quantifying the value of major regulating ecosystem services (carbon storage, flood control, sediment retention, and nutrient retention) under contemporary LULC conditions, and (3) assessing the changes in the provision of these ecosystem services under a FRA scenario.Based on prior knowledge, it is hypothesized that (H1) there are considerable transitions from forest to barren and grasslands attributed to mining and reclamation in the NFKR watershed; (H2) the barren and grassland land covers provide the lesser amount of these ecosystem services than the forest land cover and consequently the subwatersheds with less mining/reclamation activities have higher mean ecosystem service provision; and (H3) FRA will produce greater ecosystem services compared to current reclamation practices, but there will still be considerable spatial variability of ecosystem service provision at the subwatershed level driven by LULC composition, topography and climate.

Study Area
The North Fork Kentucky River (NFKR) watershed is situated in the Eastern Kentucky Coal Field physiographic region, which is part of a larger physiographic region, the Cumberland Plateau [30].The watershed occupies most of Letcher, Perry, Knott, and Breathitt counties and a small portion of Lee County (Figure 1).These counties are part of the 65 counties in Central Appalachia where surface mining has been concentrated [14].The welfare of many residents in Central Appalachia have been affected by surface mining in the past [31].The NFKR watershed ranks among the groups with highest need for protection and restoration [30].The watershed's geology is comprised of coals, sandstones, and shales [30].Land form is characterized by mountainous terrain, rapid surface runoff, and moderate rates of groundwater discharge [30].Elevation ranges from 193 m to 998 m.The NFKR headwaters are allocated in Letcher The watershed's geology is comprised of coals, sandstones, and shales [30].Land form is characterized by mountainous terrain, rapid surface runoff, and moderate rates of groundwater discharge [30].Elevation ranges from 193 m to 998 m.The NFKR headwaters are allocated in Letcher County.The main stem of the river is 270 km long and flows northwest through the communities of Whitesburg, Hazard, and Jackson to reach Beattyville where it joins with the South Fork to form the Kentucky River (Figure 1).
The NFKR watershed occupies a total area of approximately 3430 sq.km.There are 44 HUC-12 subwatersheds (Figure A1).HUC (Hydrologic Unit Code) is a scheme to describe hierarchy of watershed in the US, with HUC-2 representing the broadest region level and higher number of HUC digits representing spatially smaller levels.The size of the HUC-12 subwatersheds ranges from 38 sq.km to nearly 143 sq.km and the mean area is 76 sq.km (Table A1).

Land Cover Change Analysis
The Kentucky portion of the NLCD for 2001 and 2011 were obtained from the Kentucky Geography Network.The NLCD is a comprehensive Landsat-based, 30-m resolution land cover product.The NLCD has 16 land cover classifications applied across the continental United States.The first NLCD dataset was published for 1992, but its classification scheme does not match with that for the subsequent years of 2001 and forward, so the land covers of 2001 and 2011 (the latest available with the 2016 data are scheduled for release on 28 December 2018) have been chosen to study the land cover change over a decade.
The NLCD classification system is based on the Anderson Land Cover Classification System Level II [32].There are generally a few Level-II LULC classes for each Level-I LULC type.For instance, LULC type Forest has three sub categories: Deciduous, Mixed, and Evergreen forest.Developed Land as a broad land cover has four sub categories with varying intensities (i.e., percentage) of impervious surfaces.These sub classes were aggregated into a single land cover type for simplicity in this research.In addition, the wetlands were aggregated with water class, and the cultivated crops were reclassified into pasture because there were negligible cultivated areas and visual comparison to the crop landscape dataset of CropScape [33] showed that most crops classified in the NLCD were pastures in this watershed.
The GIS data of watershed boundary was obtained from the National Hydrography Dataset (NHD).This watershed boundary data was used in the GIS operation Clip to extract NLCD data only for the NFKR watershed (Figure 1).The Not Equal operation in ArcGIS-Spatial Analyst toolbox was used to determine where in the watershed land cover has changed from 2001 to 2011.The output from this operation was a new raster layer in which cells/pixels valued 1 indicated land cover change and 0 indicated no change.To identify the magnitude of LULC changes among the 44 HUC-12 subwatersheds, a Zonal Statistics operation was performed on the output of Not Equal tool to compute the total number of cells/pixels that changed land cover within each subwatershed.Then the Join Field tool was used to join the contents of the zonal statistics table to the watershed shapefile for mapping the amount of area with LULC changes at the subwatershed level.The Tabulate Area operation was applied to calculate the LULC transiting matrix tabulating amount of area changing from one land cover type in 2001 to another in 2011.The transition matrix was then used to assess overall LULC changes in the entire watershed.To quantify the contribution of mining/reclamation to overall LULCC in this watershed, a separate LULC transition matrix was computed for only the areas where active mining was identified to have occurred at least one year from 2001 to 2011.This cumulatively mined area mask was derived from the annual mining extent dataset of Pericak et al. (2018) [34].The contribution of mining and reclamation to the overall LULC change was based on the percentage of transitions from vegetated land covers to barren land and the transition from barren land to vegetated land cover types occurring within the cumulative mining area mask.

Ecosystem Services Assessment
The InVEST (Integrated Valuation of Ecosystem Services and Tradeoffs) software (InVEST 3.3.3,Natural Capital Project) was used to quantify the provision of four critical ecosystem services provided by the NFKR watershed under the most recently available LULC data (2011).InVEST software was developed by the Natural Capital Project, Stanford University [35].It uses ecological production functions to generate spatially explicit predictions of ecosystem service supply with inputs of LULC maps, corresponding biophysical attributes, and additional GIS data representing environmental conditions such as climate, soil, and topography.For this study, InVEST's Carbon Storage and Sequestration, Water Yield: Reservoir Hydropower Production, Sediment Delivery Ratio, and Nutrient Delivery Ratio (nitrogen and phosphorus) modules were used to estimate and map the annual delivery of the corresponding regulating services: carbon storage, flood control, soil retention, and surface water quality (Table 1).We chose to assess these four ecosystem services because they are characteristic regulating ecosystem services provided by the forests in Appalachian region that are important to human wellbeing but the degradation of these services due to deforestation are less obvious without a quantitative assessment.

• Carbon Storage and Sequestration
The InVEST carbon model aggregates the amount of carbon stored in four carbon pools: aboveground biomass, belowground biomass, soil, and dead organic matter to produce total amount of carbon storage.The input includes a user-defined biophysical attribute table (Table A2) that quantitatively describes each of these four biomass pools for each land use land cover type, and a LULC map.The model generates a map of total carbon storage by summing these four carbon pools for each grid cell based on its corresponding LULC type in million grams per hectare (mg per ha).In this study, the reclassified LULC map of NFKR watershed derived from NLCD 2011 was used as the input.The coefficient values used for the four carbon pools were from published data [36] and the InVEST manual (Table A2).

• Water Yield and Reservoir Hydropower Production
The InVEST water yield and reservoir hydropower production model estimates the total annual water yield (Y) for each pixel of the watershed as total annual precipitation minus total annual actual evapotranspiration [37].The input of the water yield model includes five biophysical parameters as georeferenced raster inputs (Table A3).These inputs are precipitation (mm), average annual potential evapotranspiration (mm), depth to root restricting layer (mm), plant available water content (AWC, as a proportion), and LULC [38].The precipitation data were obtained from the PRISM climate group of Oregon State University [39].The precipitation data are 30-year Normal (1981-2010) dataset with a resolution of 800 m.The precipitation layer was resampled to a spatial resolution of 30 m by an interpolation method in ArcGIS.The PET was obtained from CGIAR-CSI website [40] (http://www.cgiar-csi.org).The depth to root restricting layer and the AWC were extracted from the SSURGO (Soil Survey Geographic database) [41] in the Soil Map viewer of the Web Soil Survey (USDA Natural Resources Conservation Service (NRCS).The LULC 2011 was obtained from the Kentucky Geoportal Network and its coordinate system was projected to meters.
The InVEST Water Yield model also requires several tabular values for each LULC class (Table A2).These values include an attribute indicating whether the land cover class is vegetated or not (1 being vegetated and 0 being not vegetated), maximum rooting depth for vegetated LULC and plant evapotranspiration coefficient (Kc).Kc is used to modify the reference evapotranspiration to obtain potential evapotranspiration.The reference ET is based on a 15 cm tall surface of actively growing, well-watered alfalfa.The plant ET coefficient (Kc) is a decimal number between 0 and 1.5 based on plant physiological characteristics.These tabular values (Table A2) were obtained from the biophysical attribute table compiled by Sharp et al. (2015) [37].

•
Sediment Delivery Ratio (SDR) The InVEST Sediment Delivery Ratio model quantifies average annual sediment delivery per subwatershed and produces a map representing per-pixel contribution to yield.For each pixel, the model first computes the amount of eroded sediment or soil loss based on precipitation pattern, soil properties, and topographic conditions.The model then estimates the sediment delivery ratio (proportion of soil loss reaching the subwatershed outlet) based on the pixel's hydrologic connectivity [42] and further estimates sediment export based on the product of soil loss and sediment delivery ratio.The raster inputs required for the InVEST SDR model are a Digital Elevation Model (DEM), Rainfall Erosivity Index, Soil Erodibility, and LULC (Table A3).The DEM and LULC were derived from Kentucky Geoportal Network.The Rainfall Erosivity Index was obtained from European Soil Data Centre.The soil erodibility raster layer was acquired from the Soil Map Viewer program in GIS using SSURGO database.The tabular data needed for the SDR model includes biophysical parameters of cover management factor and support practice factor for the Universal Soil Loss Equation (USLE) (Table A2) [42].These factors are floating-point values between 0 and 1 for each land cover.These biophysical parameters are from Sharp et al. (2015) [37].

•
Nutrient Delivery Ratio (NDR) The InVEST Nutrient Delivery Ratio (NDR) model aims to quantify nutrient export (nitrogen and phosphorus) and maps the transport of nutrients from watershed sources to the stream network.The model uses a mass balance approach, describing the long-term, steady flow of nutrients based on nutrient sources (nitrogen and phosphorus) associated with different LULC and the retention properties of pixels belonging to the same path [37].The input raster layers required for the InVEST NDR model are DEM, LULC, and precipitation (Table A3).The DEM and LULC used are same as the SDR model inputs, obtained from the Kentucky Geoportal Network.The precipitation raster is from the PRISM Climatic Group of Oregon State University [39].The tabular coefficient values for nitrogen and phosphorus loading for each land use category required (Table A2) are from Sharp et al. (2015) [37] and Line et al. (2002) [43], respectively.

FRA Scenario
To assess the ecological benefits brought by the FRA reclamation scenario, a new LULC map was created and used in InVEST models, while all other GIS input data and biophysical parameters were kept the same as the ecosystem service assessment of the contemporary LULC conditions.This scenario reflects a best-case scenario to provides a bench mark for us to assess the potential benefits gained by the FRA at watershed scales.The new LULC map was derived from NLCD 2011 map using the reclassification tool in ArcGIS in which all the barren, grassland, and shrubland were reclassified into forests.The ecosystem services indicators assessed in this scenario were carbon storage, water yield, sediment export, and nutrient export.The output of the various InVEST models were then analyzed in ArcGIS to examine differences in the biophysical indicators between this FRA scenario and the business as usual (BAE) scenario at the watershed level.

Statistical Analysis
To quantify the spatial heterogeneity of ecosystem services provision at the HUC-12 subwatershed level, a GIS zonal statistics operation was performed to compute mean value of the biophysical indicators within each subwatershed under both the contemporary and FRA scenarios.Zonal statistics was also applied with the 2011 land cover map to compute mean ecosystem services production of different land cover types.In addition, mean and standard deviation of elevation, slope steepness, precipitation, and percentage of forest land cover were computed at the subwatershed level, and were used as the explanatory variables in the subsequent driving factor analysis.Multiple linear regression (MLR) model was applied to evaluate the effect size and relative importance of LULC composition, topography, and climate in determining the spatial variability of the five biophysical indicators of ecosystem services under the FRA scenario.Since explanatory variables represent different biophysical properties, a scaling approach was used to standardize all variables into z-scores, which represents the distance between the raw value and the mean value in units of the standard deviation.After standardization, the effect size of each explanatory variable was quantified as the estimated slope from the MLR model.A model averaging method was used to estimate the slope of each explanatory variables to account for model selection uncertainty [44].We used the "MuMIn" package in R 3.4.1 [45] to identify a set of alternative models whose delta Akaike Information Criterion (AIC) is less than 3.The relative importance of a given explanatory variable was quantified as the sum of Akaike weights over all of the selected models in which the term appeared [44,46].

Land Cover Change
Six of the seven land use land cover categories increased their areal coverage between 2001 and 2011, only forest was reduced in its area (Table 2).The forest area reduction is notable, with a decline of 7751 hectares from 2001 (266,256 ha) to 2011 (258,505 ha).The loss in forests is accompanied with an increase in barren lands: 3844 hectare gain in 2011.Similarly, 3522 hectares of grasslands have been created.The area occupied by shrub land covers increased slightly by 53 hectares.
In terms of spatial distribution of land cover change, a distinct variability can be observed in the landscape.The subwatersheds at the central locations of the NFKR watershed generally experienced the greatest magnitude (>600 ha) of LULC changes (Figure 2).In addition, the size of these land cover change patches is quite large.The rest of the NFKR watershed is not free of land cover changes, with the subwatersheds located at the upper portions of the watershed (i.e., downstream of the North Fork Kentucky River) having the lowest magnitude of LULC changes.In addition, the land cover change patches are smaller in size and they are more spread out throughout the landscape (Figure 2).Of the total area of 342,895 hectares, 15,414 hectares of land covers were converted to other land cover categories between 2001 and 2011, which accounts for 4.5% of the watershed.The greatest transition is observed from forest in 2001 to barren land covers in 2011 (4840 hectares).This is followed by conversion of forest to grassland, where 4594 hectares of forest in 2001 are changed to grassland in 2011.The change of barren land covers to grassland cover is 2189 hectares in size (Table 3).The transition matrix shows the area of LULC that has transformed from one category to another (off-diagonal).
The diagonal of the matrix shows size of area that did not change LULC between two time periods.

Ecosystem Service Assessment under Contemporary LULC Conditions
Ecosystem services production varied across the watershed, with the central locations of the watershed having lower ecosystem services production than other areas (Figure 3).Different land cover types produce different amount of ecosystem services.Forests generally produce greater ecosystem services than any other land cover types; they are associated with highest carbon storage It is assumed that the transformation of forest to barren lands and grasslands, and barren to grasslands occurring in the cumulatively mined areas is due to mining and reclamation activities.The area of such LULC changes are 4125 ha, 2521 ha, and 2059 ha respectively (Table A4).The other noteworthy transition within the cumulatively mined areas is the conversion of grasslands to barren (1298 ha).This might be due to the fact that previously reclaimed mined areas are re-mined (Figure A2).These four types of land cover transitions within the mined areas make up a total of 9994 ha, which contribute approximately 65% of total land use land cover transitions that have occurred in the entire watershed.The results thus show mining and reclamation as a major driver of overall land use land cover change in the watershed.

Ecosystem Service Assessment under Contemporary LULC Conditions
Ecosystem services production varied across the watershed, with the central locations of the watershed having lower ecosystem services production than other areas (Figure 3).Different land cover types produce different amount of ecosystem services.Forests generally produce greater ecosystem services than any other land cover types; they are associated with highest carbon storage and lowest water yield, sediment export, and nutrient export.In contrast, barren lands provide the lowest amount of ecosystem services among all non-urban land cover types; they produce lowest carbon storage and highest water yield, sediment export, and nutrient export (Table 4).
The total modeled carbon storage for the entire watershed is 71,343,168 mg.Mean carbon storage at the subwatershed level ranges from 138 mg per ha to 231 mg per ha (Table A5), with the lowest ranking subwatersheds occupying central locations of the watershed where surface mining is the most (Figure 3a).Carbon storage is highest in the forested lands with a mean of 250 mg per ha.The shrub land cover ranks second in storing carbon (100 mg per ha).Pasture ranks third with a carbon storage of 51 mg per ha.Grassland ranks fourth with a mean carbon storage of 45 mg per ha.Barren lands produce zero carbon storage (Table 4).
The total estimated water yield for the watershed is 2,219,528,435 mm.Mean water yield at the subwatershed level ranges 522 mm to 702 mm (Table A5), with the subwatersheds of highest mean water yield (hence lowest flood control capacity) occupying the central location in the watershed where the barren lands are concentrated (Figure 3b).Forested areas produce the least water yield (hence the greatest flood control capacity), with an average of 534 mm and grasslands have higher average water yield (665 mm) than the forests (Table 4).
Sediment export for the entire watershed is 52,848,288 kg.Mean sediment export ranges at the subwatershed level ranges from 101 kg per ha to 313 kg per ha (Table A5).The spatial distribution shows hotspots of sediment exports at the central locations of the watershed (Figure 3c).Among the land cover types, the barren areas have the highest mean sediment export (971 kg per ha).Pasture lands rank second with a mean sediment export of 628 kg per ha.The grasslands have a mean sediment export of 404 kg per ha.Forest land covers have the least mean sediment export (110 kg per ha) (Table 4).
The total nitrogen export of the entire watershed is 320,525 kg.Mean nitrogen export at the subwatershed level ranges from 0.67 kg per ha to 1.40 kg per ha (Table A5).The hotspots of nitrogen export are located around urban land covers (Figure 3d).Among all the five non-urban terrestrial land covers, nitrogen export is highest in barren areas (2 kg per ha), followed by pasture (1.6 kg per ha) and grassland (1.6 kg per ha) (Table 4).The forest and shrub lands both have least nitrogen export with a mean value of 0.5 kg per ha (Table 4).The total phosphorus export for the watershed is 18,384 kg.Mean phosphorous export at the subwatershed ranges 0.027 kg per ha to 0.14 kg per ha (Table A5).Phosphorus export is highest in pasture (0.051 kg per ha) and then in grassland (0.02 kg per ha) (Table 4).Forest and shrub lands show a similar pattern, with a mean phosphorus export of 0.003 kg per ha.The barren lands have the least phosphorus export (0.001 kg per ha) (Table 4).

Ecosystem Service Assessment under FRA Scenario
Although more than 91% of the watershed was covered by forest under the forestry reclamation approach scenario, there are still conspicuous spatial variability of ecosystem service production (Figure 4), and the spatial patterns of ecosystem services differ from those under the contemporary LULC condition (Figure 4 vs. Figure 3).The spatial patterns under the FRA scenario largely reflect spatial configuration of LULC, particularly the juxtaposition between forest and developed land (i.e., urban and road networks).However, other driving factors such as topography and climate were also important contributors.Influences of forest land cover proportion, slope steepness, elevation, and precipitation on five biophysical indicators of ecosystem services at subwatershed level differed substantially in terms of both effect size and relative importance (Table 5).Forest land cover proportion was the sole important variable in explaining the variability of mean carbon storage among the 44 subwatersheds.Mean precipitation and slope, and standard deviation of slope steepness were the top three important variables and they all exerted positive influences (i.e., with positive slope coefficients in the model) on water yield.Sediment export was mainly determined by

Ecosystem Service Assessment under FRA Scenario
Although more than 91% of the watershed was covered by forest under the forestry reclamation approach scenario, there are still conspicuous spatial variability of ecosystem service production (Figure 4), and the spatial patterns of ecosystem services differ from those under the contemporary LULC condition (Figure 4 vs. Figure 3).The spatial patterns under the FRA scenario largely reflect spatial configuration of LULC, particularly the juxtaposition between forest and developed land (i.e., urban and road networks).However, other driving factors such as topography and climate were also important contributors.Influences of forest land cover proportion, slope steepness, elevation, and precipitation on five biophysical indicators of ecosystem services at subwatershed level differed substantially in terms of both effect size and relative importance (Table 5).Forest land cover proportion was the sole important variable in explaining the variability of mean carbon storage among the 44 subwatersheds.Mean precipitation and slope, and standard deviation of slope steepness were the top three important variables and they all exerted positive influences (i.e., with positive slope coefficients in the model) on water yield.Sediment export was mainly determined by slope steepness, with steeper topographic slopes associated with higher amount of soil sediment.Forest land cover proportion was found to be the most important variable in explaining subwatershed level variations of nutrient export.The subwatersheds with higher forest land cover proportion were associated with lower export of nitrogen and phosphorous as the corresponding slope coefficients were negative values.Slope steepness and elevation had positive effects on the export of these two nutrients, but stronger effects were found on phosphorus than nitrogen (Table 5).The benefits of forest reclamation approach, in comparison to the contemporary 2011 LULC scenario, are evident in terms of the production of all the ecosystem services assessed in this study.The total carbon storage of the watershed was estimated to be 80,633,377 mg of carbon under the FRA scenario.The difference is 9,290,209 mg of carbon which makes up an increase of 13% (Table 6)  The benefits of forest reclamation approach, in comparison to the contemporary 2011 LULC scenario, are evident in terms of the production of all the ecosystem services assessed in this study.The total carbon storage of the watershed was estimated to be 80,633,377 mg of carbon under the FRA scenario.The difference is 9,290,209 mg of carbon which makes up an increase of 13% (Table 6) compared to the 2011 LULC scenario.The result suggests the climate regulation ecosystem service would be improved under FRA scenario.In addition, total water yield decreased by 5.2%, sediment export decreased by 40.5%, nitrogen export decreased by 22.6%, and phosphorous export decreased by 7.1% (Table 6), suggesting the watershed's capacity in regulating flood, retaining soil, and maintain surface water quality would also be enhanced under the FRA scenario.

Discussion
This study investigated how LULC has changed in the NFKR watershed and the contribution of mining and reclamation to the overall land cover land use change, followed by the valuation of major ecosystem services under contemporary LULC conditions and the assessment of the benefits of the Forestry Reclamation Approach (FRA).The results show that even during a coal industry decline period, mining and reclamation are still prominent drivers shaping the structure and functions of many landscapes in Central Appalachia.Forests provide greater regulating ecosystem services than other land cover types, but the enhancement of ecosystem service production induced by forestry-based reclamation would differ for different types of ecosystem services and at different subwatersheds.

Surface Mining as a Prominent Driver of LULC Change
Rapid economic and social development since the industrial revolution has greatly geared up human demand for mineral resources.However, mining for these resources can cause extensive changes in LULC at landscape scales [47,48].Coal is a primary energy source, and the three states from Appalachian area (Pennsylvania, West Virginia, and Kentucky) account for over 50% of the total coal production in the US [12].It is estimated that more than 500,000 hectares of land have been mined for coals in this region [21].Although coal production from Appalachian surface mining has declined by nearly 50% since its peak in 2008 due to depletion of high-quality coal seams, strict environmental regulation policy during the Obama administration era, and the collapse of coal mining industry [12].A recent study that reconstructed mining activities across 30 years shows that cumulative mining area is continuing to expand in Central Appalachia, albeit at a much lower expansion rate [16].The continuing expansion of cumulative mining area suggests that mining may still be a major driving force of the LULC change in some of the Appalachian landscapes, and this study demonstrates it is the case in the NFKR watershed.Based on the GIS analysis of the NLCD, the forested area was reduced by 7751 hectares (2.3% of the watershed) from 2001 to 2011 and barren and grassland LULCs increased 3844 and 3522 hectares, respectively (Table 2).A total of 15,414 hectares of land covers were converted to other land cover categories between 2001 and 2011, which accounts for 4.5% of the watershed.In contrast, within the cumulatively mined area of 29,189 ha, a total of 10,785 hectares of land with changing LULC, which accounts for 37% of mined area (Table A4).The forest land cover within the mining area decreased from 9675 ha in 2001 to 3539 ha in 2011, with a total reduction of 6136 ha, which accounts for 79% of forested area reduction at the watershed scale.Overall LULC changes and forest clearing were far more extensive within the cumulatively mined area than the surrounding area, suggesting mining is a major driver of LULC change in the NFKR watershed.
Our LULC change analyses show that areas changed from barren lands to grassland is much greater than the conversion of barren to forest.Only 179 hectares of barren lands in 2001 have converted to forest by 2011, whereas the conversion of barren lands to grasslands is 2189 hectares (Table 3).Tree plantation was prevalent in the abandoned mining sites prior to the SMCRA implementation, but reclamation to grassland became the dominant approach post-SMCRA [49].Grasslands provide less regulating ecosystem services (e.g., climate regulation, flood control, and sediment retention) in comparison to forests.The InVEST modeling showed that grasslands have a mean carbon storage of 45 mg per ha, which is much less than what forests can store (250 mg per ha).Comparing to grassland, forests have lower mean water yield (534 vs. 665 mm), sediment export (110 vs. 404 kg per ha), nitrogen (0.5 vs. 1.6 kg per ha), and phosphorus export (0.001 vs. 0.02 kg per ha) than grasslands (Table 4).Although SMCRA mandates restoring the post-mining land to a condition capable of supporting the uses to the level similar to or higher than that prior to any mining, the majority of reclamation has failed to meet such standards when ecosystem services are considered as the evaluation criteria.

Ecosystem Services Assessment
Although ecosystem services are considered important, there is a lack of studies in Appalachia to evaluate how different land use and land cover might affect the provision of various ecosystem services.However, there are studies that are consistent with our findings that LULC changes driven by mining can significantly reduce the potential of a watershed to provide ecosystem services.According to Burkhard et al. (2009) [35], highly modified land cover types, such as mine sites, have very low or no relevant capacities to provide ecosystem services.Some studies have also stressed that unique impacts are brought about by mining, particularly by mountain top mining such as altering landform shape and structure and burying headwater streams [50].All these changes adversely affect the functioning of ecosystems and results in reduced capacity of the landscape to regulate climate and flood, to retain sediments and nutrients, and to conserve and purify water.
The Millennium Ecosystem Assessment classifies ecosystem services into four categories, namely supporting, regulating, provisioning, and cultural ecosystem services.Our study only assessed the regulating ecosystem services because, unlike provisioning services such as food and timber production that are readily appreciated by the public, regulating ecosystem services are often overlooked in the decision-making as they do not have a market value.These regulating ecosystem services are critical in determining the production of many other ecosystem services and decline of these services foreshadows the decline of other ecosystem services [51].Understanding the impacts of mining on the provision of critical regulating ecosystem services at landscape scales provide a better or understandable language to contribute positively to the formulation and evaluation of the relevant environmental policies [52,53].
The ecosystem assessment conducted with the InVEST modeling showed that even under a best-case reclamation scenario in which all the barren and grasslands were to be restored to forests, the enhancement of the regulating ecosystem services would still vary by ecosystem service type and subwatershed.This suggests that environmental benefits of forestry-based reclamation are landscape-specific and a thorough analysis regarding which ecosystem services need to be enhanced in which locations is warranted.In addition, forest based reclamation is not the only alternative or the best reclamation option in every mining site.An abandoned mining site can be reclaimed to grassland, pasture, agricultural land, etc.So trade-offs are bound to occur when the mining sites are reclaimed to various land uses.These trade-offs which arise from management choices can change the type, magnitude, and relative mix of ecosystem services provided by different ecosystems [54].Reclamation would also require to consider land suitability.For example, mining sites that are close to developed areas may be more suitable to be developed to residential, commercial or industrial area and mining sites away from human settlement can be restored to forest [55].The proper reclamation planning should consider both land suitability and ecosystem services that are needed to restore for the people living in the surrounding communities.

Limitations and Uncertainty
An important limitation of this study is the land cover dataset.The NLCD is a broad dataset and analysis for this study is for a relatively small spatial area of the continental scale dataset, so compromises are inherent in the classification of land use land covers.The NLCD classifications were based on the information from multiple years prior to 2001 and 2011, meaning classification used in the NLCD may not truly represent ground truth in all pixels for a given year.There are also known inaccuracies because of the techniques used to collect and classify the remotely sensed data.Some pixels identified as LULC change may be artifacts due to classification inaccuracies, not real changes.This is proven a major issue when conducting LULC change analysis with 1991 NLCD data, but a less serious issue with 2001 and 2011 NLCD data [56].Generally, the patches misclassified as LULC changes are small in size, and should have negligible impacts when assessing major LULC transitions at landscape scales.
The InVEST software has its own modeling limitations.For the InVEST Carbon Storage and Sequestration Model, the model only estimates the temporally average carbon storage for each LULC hence assumes that none of the LULC types in the landscape are gaining or losing carbon over time.Changes in carbon storage simulated in this model can only be induced by the changes from one land cover type to another.The InVEST Water Yield and Reservoir Hydropower Production Model is based on annual averages, which neglects extremes and does not consider the temporal dimensions of water supply.It does not consider complex land use patterns or underlying geology, which may induce complex water balances.The main limitation of the InVEST Sediment Delivery Ration Model is its reliance on the Universal Soil Loss Equation.Even though this equation is widely used, it only represents rill erosion process, which is the removal of soil by concentrated water running through little streamlets.The InVEST Nutrient Delivery Ration Model is highly sensitive to inputs, so small errors in the empirical load parameter, will have a large effect on predictions of nutrient delivery.The tabular values used are not entirely from the study area because of data limitation in the study area; they have been acquired from published sources and the master table of the InVEST manual.
This study assumed that the FRA forests on the abandoned mining sites can provide regulating ecosystem services at a similar level to the natural forests on the unmined sites.This assumption was made to assess the ecological benefits at a best-case scenario.Although site conditions on those FRA forests are altered via mining, field studies of FRA effects show that forests reclaimed on the mining sites can quickly recover many important ecological processes characteristic of the native forest ecosystems [20].This is particularly the case for the processes related to regulating ecosystem services.For example, because the FRA emphasizes replacement of loose soil materials, FRA reclamation practices can help the mine site "evolve" hydrologically, developing infiltration/runoff patterns more like unmined landscapes with time [24,57].
For the ecosystem service assessment conducted in this study, a major implication is for the authorities and the general public to appreciate the value of ecosystem services; to gain knowledge about the loss of ecosystem services due to land conversions like mining, and the potential improvement in the delivery of ecosystem services when forest-oriented reclamation practices are applied.This study adds a block to the ecosystem services' study of Appalachia as such studies are scarce in this region.The outcomes presented in tables and maps illustrate different land cover types have different capacity in providing ecosystem services, and the spatial characteristics of LULC, topography, and climate are important drivers of the provision of ecosystem services at watershed and subwatershed levels.The maps produced in this study provide important spatially explicit information to support managers to identify areas where the ecosystems are produced in larger quantities and where they are not.Despite model and data limitations, this study can be helpful in enforcing and popularizing reclamation strategies like the forestry reclamation approach and prioritizing locations for restoration at landscape scales.

Conclusions
Assessing ecosystem services under different land use land cover conditions can provide critical information regarding the benefits of the ecological restoration.Such information is helpful in the decision-making process for achieving landscape sustainability.In mountainous areas such as the Appalachian region, the change of LULC is often extensive and intensive due to the large-scale extraction of natural resources, and this study demonstrated the extent of this land degradation problem and its cascading effects on ecosystem services provision in an Appalachian watershed.Comparing to current reclamation strategies in which most abandoned mined sites are reclaimed to grassland, forestry-based reclamation can greatly enhance the provision of many regulating ecosystem services, particularly sediment retention and water purification.It further showed that the different types of ecosystem services can have different response to the forestry-based reclamation strategy, and the spatial heterogeneity of ecosystem services provision is largely driven by LULC composition, topographic variability and climate.The results showcase the importance of considering landscape ecology and land change science in ecosystem service assessment and policy-making of mined land restoration.The carbon storage estimates for the four carbon pools is measured in mg per ha.Kc is the plant evapotranspiration coefficient.Root depth is the maximum root depth (mm) for vegetated land use land covers.usle_c and usle_p are the cover management factor and support practice factor for the USLE respectively.Load_n and load_p are the nutrient loading for each land use (kg per ha per year).

Figure 1 .
Figure 1.The North Fork Kentucky River Watershed showing the counties occupied and their towns.

Figure 1 .
Figure 1.The North Fork Kentucky River Watershed showing the counties occupied and their towns.

Forests 2018, 9 , 23 Figure 2 .
Figure 2. Area of land use land cover (LULC) changes by HUC-12 subwatershed in the NFKR watershed between 2001 and 2011, overlaid with locations of where LULC changes occurred.The two insert maps show the extent of active mining area in 2001 and 2011, respectively based on the published annual mining extent data of Pericak et al. (2018) [34].

Figure 2 .
Figure 2. Area of land use land cover (LULC) changes by HUC-12 subwatershed in the NFKR watershed between 2001 and 2011, overlaid with locations of where LULC changes occurred.The two insert maps show the extent of active mining area in 2001 and 2011, respectively based on the published annual mining extent data of Pericak et al. (2018) [34].

Figure 4 .
Figure 4. Spatial distribution of modeled (a) carbon storage (mg per ha), (b) water yield (mm), (c) sediment export (kg per ha), (d) nitrogen export (kg per ha), and (e) phosphorus export (kg per ha) across the NFKR Watershed under the FRA scenario.

Figure 4 .
Figure 4. Spatial distribution of modeled (a) carbon storage (mg per ha), (b) water yield (mm), (c) sediment export (kg per ha), (d) nitrogen export (kg per ha), and (e) phosphorus export (kg per ha) across the NFKR Watershed under the FRA scenario.

Figure A2 .
Figure A2.An example of land conversion from grassland in 2001 to barren in 2011.Figure A2.An example of land conversion from grassland in 2001 to barren in 2011.

Figure A2 .
Figure A2.An example of land conversion from grassland in 2001 to barren in 2011.Figure A2.An example of land conversion from grassland in 2001 to barren in 2011.

Table 1 .
Ecosystem services assessed in this study and their biophysical indicators.

Table 2 .
Area of each land use land cover (LULC) type in 2001 and 2011 (unit, hectares) and its changes (unit, percent).

Table 3 .
Size of area that has experienced transition from one land cover category to another between 2001 and 2011 (unit, hectares) in the North Fork Kentucky River (NFKR) watershed.

Table 4 .
Mean production values of the five biophysical indicators of the ecosystem service by each land cover type under the 2011 LULC scenario.The standard deviation values are presented in the parenthesis after the mean production values.
Note: There is no variability in carbon storage within a land cover type calculated by the InVEST model.The sample size (i.e., number of cells) of each LULC class is the division of total area (as shown in Table2) by cell size (0.01 ha).

Table 5 .
Slope coefficient (Slope) and relative importance (RIMP) of the predictor variables estimated from a model average approach for explaining variability of the five biophysical indicators of ecosystem services at the subwatershed level.
Note: Forest, STD_Forest are mean and standard deviation of forest land cover proportion within a subwatershed.Elevation and STD_Elevation are mean and standard deviation of elevation within a subwatershed.Slope and STD_Slope are mean and standard deviation of slope steepness.Precipitation and STD_Prep are mean and standard deviation of precipitation.The significance level was denoted with code ** p < 0.01, * p < 0.05, and ^ p < 0.1.

Table 5 .
Slope coefficient (Slope) and relative importance (RIMP) of the predictor variables estimated from a model average approach for explaining variability of the five biophysical indicators of ecosystem services at the subwatershed level.
Note: Forest, STD_Forest are mean and standard deviation of forest land cover proportion within a subwatershed.Elevation and STD_Elevation are mean and standard deviation of elevation within a subwatershed.Slope and STD_Slope are mean and standard deviation of slope steepness.Precipitation and STD_Prep are mean and standard deviation of precipitation.The significance level was denoted with code ** p < 0.01, * p < 0.05, and ˆp < 0.1.

Table 6 .
The modeled ecosystem service productions under the contemporary 2011 LULC condition and the FRA scenario and their differences.

Table A1 .
The ID, name, area, and HUC of the 44 Subwatersheds in the NFKR watershed.

Table A2 .
Carbon pool estimates and other biophysical attributes used for the InVEST models.

Table A4 .
Size of area that has experienced transition from one land cover category to another between 2001 and 2011 (unit, hectares) within the mined sites in the NFKR watershed.

Table A4 .
Size of area that has experienced transition from one land cover category to another between 2001 and 2011 (unit, hectares) within the mined sites in the NFKR watershed.

Table A5 .
Mean ecosystem service production per subwatershed under 2011 LULC conditions.