Landscape Changes in the Southern Coalfields of West Virginia: Multi-Level Intensity Analysis and Surface Mining Transitions in the Headwaters of the Coal River from 1976 to 2016

This study analyzes land-cover transitions in the headwaters of the Big Coal River in the Central Appalachian Region of the US, from 1976 to 2016, where surface mining was found as the major driver of landscape change. The land-change analysis combined Multi-Level Intensity Analysis for two-time intervals (1976–1996, 1996–2016) with Difference Components, to differentiate suspected misclassification errors from actual changes. Two land cover classifications were obtained with segmentation analysis and machine learning algorithms from historical high-resolution aerial images and ancillary data. Intensity Analysis allowed for the inspection of transitions across five land cover (LC) classes and measure the degree of non-stationarity of land change patterns. Results found surface mining-related classes and their transitions, including the effects of reclamation processes on areas mined before the enactment of the Surface Mining Control and Reclamation Act (SMCRA, 1977). Results included changes in settlement distribution, low vegetation, water bodies, and forest class transitions. The findings can be applied to infer similar land-change processes in the more extensive Appalachian region where Mountain Top Removal (MTR) operations are widespread. The overall method can be used to address similar problems and inform landscape managers with detailed data to support land use alternatives and conservation in regions that experienced intense changes and are characterized by anthropogenic disturbances and novel ecosystems.


Introduction
This study presents the findings of a land-change analysis conducted in the headwaters of the Coal river, during a 40 year interval, from 1976 to 1996. This is a subwatershed located in the central Appalachian coalfields region of the US [1,2], within the southern coalfields of West Virginia (WV) [3] (Figure 3). The land change processes in this area are of particular interest for the land-change transitions associated with surface mining and Mountain Top Removal (MTR). The overall knowledge derived by the analysis of the land change processes can be relevant for characterizing the region's landscape with information concerning the biophysical structure and function of ecosystems [4] through the identification of specific land cover (LC) transitions and classes.
The study's main objective is to address landscape changes in the study area and test the stationarity of land transitions applying a post classification comparison (PCC) method based on five relevant LC classes. To achieve this goal, the Multi-Level Intensity Land-change Analysis method [5] was applied to examine land-change processes in the Land 2021, 10, 748 3 of 32 reclaiming the land with vegetative covers. Both conditions will bring somehow contradictory outcomes in the region; the AOC to the leveled landscapes of MTR sites, while the fast-growing grasses will be used by coal companies to quickly meet the reclamation requirements [35,36]. Nevertheless, these reclamation procedures represented a response to the need of controlling and managing the unreclaimed landscape, the severe soil erosion, rockfalls, and landslides caused by strip mining operations before 1977 ( Figure 1) [32,[37][38][39].  The Clean Water Act of 1972 (CWA, 33 U.S.C. §1252) and the Surface Mining Control and Reclamation Act of 1977 (SMCRA, 30 U.S.C. §1201) represented a turning point for the US surface coal mining history. The SMCRA consistently readdressed, from individual state regulation to the federal level, the policies that drove the socio-technical processes of coal extraction, mine reclamation, and environmental regulation.
The SMCRA established the restoration of mining areas introducing the principle of "approximate original contour" (AOC), to eliminate the highwalls and the mining spoil for stabilizing potentially unstable terrain slopes [32]. Another requirement consisted in reclaiming the land with vegetative covers. Both conditions will bring somehow contradictory outcomes in the region; the AOC to the leveled landscapes of MTR sites, while the fast-growing grasses will be used by coal companies to quickly meet the reclamation requirements [35,36]. Nevertheless, these reclamation procedures represented a response to the need of controlling and managing the unreclaimed landscape, the severe soil erosion, rockfalls, and landslides caused by strip mining operations before 1977 ( Figure 1) [32,37,38].

Monitoring Surface Mining in the Region
Since the mid-1960s, remote sensing studies have assumed a critical role in monitoring the magnitude and the effects of surface mining in the region [32]. In the early applications, remote sensing methods were mainly based on the visual interpretation of aerial images, used for environmental applications, to assess surface mining's extents and their effects on the landscape [39]. By 1973 there was an increasing use of Landsat imagery. With the new satellites, the spatial, spectral, and radiometric resolution of new sensors progressively contributed to better results and applications. A comprehensive review of remote sensing studies applied to Mountaintop Mining (MTM) was published in 2001 [40].
In 2009, a land change study based on Landsat time series [7] examined land cover transitions in four watersheds in the Central Appalachian mountains, monitoring three LC classes (bare/urban, forest, and grassland/pasture/crop) across four time points (1976,1987,1999,2006). The study's focus was to assess changes in surface mining areas to link them to the losses of ecosystems' hydrologic functions [41]. The authors reported the inability to use the National Land Cover Classification (NLCD) products due to their low accuracies in classifying mine-relevant LC classes. To avoid the limitations derived by the low-resolution of Landsat images, the authors combined a set of ancillary data (mining permits, NLCD maps, data derived by manual photo interpretation) within a decision tree process to derive transitions based on class membership from time intervals comparisons. Results established the expansion of surface mining started in the 1950s and reached its peak in 1976. After the SMCRA adoption in 1977, there were two effects, mined areas tended to stabilize, and reclaimed areas progressively increased. The study highlighted that these results should be cautiously extended to the MTM region of southern West Virginia [7].
In 2018, a new land change study [17] proceeded with the revision of the surface mining extent in 74 counties in Central Appalachia using a dataset based on Landsat data built on a yearly imagery basis using the Google Earth Engine platform. Surface mines were mapped using two models corresponding to two different datasets. The main dataset involved the 1985-2015 time interval; for the interval 1976-1984, the authors reported some limitations in deriving a complete mapping of mines due to the lack of a comprehensive dataset. They also pointed to the opportunity of improving future studies before 1984 using data from earlier Landsat satellites. The results established that the higher value of land converted to new surface mines during each year was reached in 1999 (116 km 2 yr −1 ), and the lower in 2015 (31 km 2 yr −1 ). The study established that after 2010, there was a consistent decrease in new mining areas; starting from 2010, coal companies had to mine up to three times of land to obtain the same amount of coal they obtained in the 1990s. This trend pointed to a remarkable decrease in coal resources accompanied by a proportional increase in environmental costs [17].
Another set of studies were focused on land cover classification (LCC) problems in the Appalachian mountain region and on the specific problems derived by surface miningrelated LC classes, like the differentiation of mine-reclaimed grasslands from non-mining grasslands due to spectral similarities between barren land and developed land. These studies discussed the use of high-resolution images (0.6-2 m), the integration of ancillary data (e.g., topographic models, vector data), and the use of segmentation analysis combined with machine learning algorithms [42][43][44]. Table S1 reports a list of remote sensing and land change studies focused on surface mining in the Central Appalachian Region.

Post Classification Comparison Using Intensity Analysis and Difference Components
The land-change analysis in this study was obtained using a post-classification comparison (PCC) method. PCC is the most frequently used of several change detection techniques [46,47]. Among the main advantages, PCC provides a complete matrix of change directions, while major drawbacks generally originate from difficulties in deriving LCCs from historical datasets and assessing their accuracy [48]. PCC produces a crosstabulation matrix that summarizes the number of cells or pixels occurring in each potential pre-and post-change category combination to compare two timepoints, which some authors define as a transition matrix [49]. Errors in change map are multiplicative in combination, while the accuracy of the change map might be similar to the product of each map's accuracy [49,50]. A central complexity of PCC is rigorous error assessment, which is essential to differentiate misclassification error from actual landscape change [51,52]. Through aggregation across two time points, errors in the LCCs may produce mismatches that may be incorrectly identified as change. Indeed, LC products with levels of overall accuracy generally considered as satisfactory or adequate (≥85%), when combined, may be inadequate for accurate change detection [6]; an issue of particular concern for classes with low accuracy [53].
Intensity Analysis is a quantitative method that allows for the inclusion of multiple time points in PCC to test the stationarity of land transitions [5]. Intensity Analysis allows for the assessment of land-change components of transition and permanence [5,54] and the differentiation of systematic land-change processes from uniform processes [55,56]. Systematic processes of change are not due to a merely random or coincidental event or change; instead, they are characterized by specific processes, like a category's tendency to gain or lose land area from other specific categories [57].
Intensity Analysis tests the stationarity of change on three different analysis levels: the interval level, the category level, and the transition level. In the first level, the interval

Post Classification Comparison Using Intensity Analysis and Difference Components
The land-change analysis in this study was obtained using a post-classification comparison (PCC) method. PCC is the most frequently used of several change detection techniques [45,46]. Among the main advantages, PCC provides a complete matrix of change directions, while major drawbacks generally originate from difficulties in deriving LCCs from historical datasets and assessing their accuracy [47]. PCC produces a crosstabulation matrix that summarizes the number of cells or pixels occurring in each potential pre-and post-change category combination to compare two timepoints, which some authors define as a transition matrix [48]. Errors in change map are multiplicative in combination, while the accuracy of the change map might be similar to the product of each map's accuracy [48,49]. A central complexity of PCC is rigorous error assessment, which is essential to differentiate misclassification error from actual landscape change [50,51]. Through aggregation across two time points, errors in the LCCs may produce mismatches that may be incorrectly identified as change. Indeed, LC products with levels of overall accuracy generally considered as satisfactory or adequate (≥85%), when combined, may be inadequate for accurate change detection [6]; an issue of particular concern for classes with low accuracy [52].
Intensity Analysis is a quantitative method that allows for the inclusion of multiple time points in PCC to test the stationarity of land transitions [5]. Intensity Analysis allows for the assessment of land-change components of transition and permanence [5,53] and the differentiation of systematic land-change processes from uniform processes [54,55]. Systematic processes of change are not due to a merely random or coincidental event or change; instead, they are characterized by specific processes, like a category's tendency to gain or lose land area from other specific categories [56]. Intensity Analysis tests the stationarity of change on three different analysis levels: the interval level, the category level, and the transition level. In the first level, the interval level defines the speed and size of change relative to the time intervals; the second level describes the relative activity of gross loss and gross changes at the categories level, defining dormant versus active categories; the third level, the transition level, identifies if the gain of each category targets or avoids the other categories. The uniform line's position defines the stationarity of the land-change processes across the three levels of analysis. Identifying non-stationary land-change transitions across the time intervals may unveil some drivers, or causal mechanisms [57], that might be responsible for the specific patterns of change due to not merely random events [5].
The Intensity Analysis framework has been recently updated [58] to include the Difference Components method and the Quantity and Allocation disagreement metrics [59]. Xie et al. (2020) [6] integrated Intensity Analysis with Difference Components and error analysis in a unique framework to enable the evaluation of changes among LC classes and distinguish possible errors among land transitions.
The transition matrices allowed to examine the components of land-change, including gross Gains and gross Loss, Net change, Total change, and Swap. It has been noted that studies with discussions limited only to Net change quantities may offer inaccurate evaluations because the lack of Net change does not necessarily represent a lack of real changes on the landscape [60]. Indeed, the Swap component of change among LC classes in a landscape can be substantial without interfering on Net Change size. In the transition matrix, the absolute value of Net change represents the difference between the total cover of a specific class for the two time points, while the Total change can be expressed for each category as the sum of the gains and losses or the sum of Net change and Swap. Swap can be computed as the difference between Total change and the absolute Net change [60].
More recent updates have defined Net change as Quantity difference (or Quantity disagreement) and Swap as Allocation difference (or Allocation disagreement) [61,62]. Quantity difference measures the amount of difference between two maps in terms of unequal proportions of the categories; allocation difference analyzes differences in terms of spatial allocations of the categories [62]. The metrics can be used to compare both the differences between a reference map and a classified map in the context of error analysis in LCCs or assess differences between two temporal maps to measure land change. Allocation difference has been afterward expressed as the sum of two components: Exchange and Shift [59]. Exchange represents the components of difference that occurs between two categories without modifying the size of either category; while Shift represents the components of difference that is neither Quantity nor Exchange [58].

Study Area
The study area consists of 58,661 hectares of the Big Coal River's headwaters in WV (USA), and includes the two main tributary branches: Clear Fork and Marsh Fork ( Figure 3). The area is part of the Appalachian Plateau physiographic province, which is characterized by drainage networks with dendritic shapes, deeply incised and dissected valleys, steep slopes and a high degree of local relief, and narrow winding ridges [63]. The northern portion of the study area consists of the Clear Fork catchment and the Mash Fork basin's lower and middle catchments. Land cover in this area is mainly characterized by forest and surface mining-related classes and MTRVF sites. In the valley floors, built-up areas, infrastructures, and industrial facilities used for coal processing, storage, and transport are frequent.

Data Collection and Processing
All the LCCs were based on a pre-processing geographic object-based image analysis (GEOBIA) phase [65,66]. The variables derived from the segmentation analyses, which included ancillary data from Landsat satellites and digital elevation model (DEM) derivatives, were processed using machine learning (ML) supervised classification algorithms in R [67]. The use of historical images in this study introduced further difficulties due to the lack of spectral information in the primary databases and the limitations of building a proper validation sample. Other difficulties arose from classifying specific LC classes, like the scattered settlements (developed areas) and the fine-grain, linear patterns characteris- The southern part of the study area, consisting of the Upper Marsh Fork and Stephens Lake catchments, is a plateau with a noticeably different morphology than the rest of the study area. It has higher valley-bottom elevations but less local relief. Land cover in this part of the study area is dominated by open areas with grasslands and widespread settlements that densify to the East near the town of Beckley.

Data Collection and Processing
All the LCCs were based on a pre-processing geographic object-based image analysis (GEOBIA) phase [64,65]. The variables derived from the segmentation analyses, which included ancillary data from Landsat satellites and digital elevation model (DEM) derivatives, were processed using machine learning (ML) supervised classification algorithms in R [66]. The use of historical images in this study introduced further difficulties due to the lack of spectral information in the primary databases and the limitations of building a proper validation sample. Other difficulties arose from classifying specific LC classes, like the scattered settlements (developed areas) and the fine-grain, linear patterns characteristic of the strip mining operations (barren land). The 2016 LCC was obtained from an existing dataset of WV [44]. Each LCC included the error matrices and the number of samples derived from the validation polygons used during the classification stages. The study also included an intermediate post-processing phase conducted after the LCCs to reduce the amount of misclassified areas from the maps before the Intensity Analysis stage.

The LCC of 1976
For the 1976 LCC, three datasets were used during the classification process. The primary dataset consisted of an orthomosaic derived from 31 single black and white aerial photographs obtained from the U.S. Geological Survey (USGS) collected in 1976 and 1977. The high-resolution digital images (2 m) were derived from the original black-and-white (BW) panchromatic film (V-F Panchromatic) [67]. The 3 m DEM for the 1976 classification was produced from digital line graph (DLG) contour lines derived from 1:24,000 scale, 7.5 × 7.5 Minute Map Series Historical Topographic Maps that range from 1965 to 1969 [68]. The orthomosaic was obtained using the Imagine Photogrammetry add-on in Erdas Imagine v16.5 [69]. Ground points were derived from 1996 Digital Orthophoto Quarter Quads (DOQQs) [67] as a reference for X and Y coordinates using points with recognizable locations, such as streets intersections and buildings; the heights (Z coordinates) were derived from the 1965-1969 DEM included in the photogrammetric model. The final Total Image Unit-Weight root mean square error (RMSE) for the orthomosaic was equal to 0.89 (meters). Tonal and brightness differences were corrected using MosaicPro in Erdas Imagine. A Landsat 1 Multispectral Scanner (MSS) image from 1976 (June 7) with four spectral bands was obtained [68]. The band files were stacked to a single raster using the Layer Stack tool in ERDAS Imagine.
The 1976 LC included six classes and was derived with a 2-m resolution ( Table 1). The general method used for the LCCs of 1976 and 1996 was based on a pre-processing GEOBIA phase undertaken using eCognition [70]. Variables derived for each object during the segmentation analysis were provided as predictor variables for a ML supervised classification in R using the caret package [71] with the ranger method [72]. To potentially improve the classification performance, additional variables were obtained from the ancillary data (satellite image and topographic derivatives). The 1976 LCC was based on an unbalanced training dataset with 108 variables obtained with the multiresolution segmentation algorithm in eCognition. The validation sample was derived using a two-stage stratified cluster random sampling design [46,73]. The sampling design process consisted of generating random points throughout the study area extent; then, image objects intersecting these points were selected to create the validation set. The objects were assigned to the closest class by visual inspection of the 1976-1977 orthomosaic. The number of objects used for the training and validation samples are reported in Table S5.

Class Name Description
Barren Non-vegetated areas generally associated with surface mine features, strip-mining roads, spoil banks, coal deposits, dump sites.

Deciduous forest
Areas dominated by broad-leaved trees and woodlands.
Evergreen forest Areas dominated by evergreen trees. It can include evergreen shrubs.

Low vegetation
Low vegetation areas such as grasslands, including pastureland, agricultural fields, pastures, and croplands.

Developed areas
Areas dominated by mixed development, residential areas and yards; areas characterized by impervious surface such as roads and parking lots. It can also include industrial complex and facilities used for coal processing.

Water
All waterbodies, including rivers, lakes, water ponds and coal sludge impoundments.
The overall accuracy for the 1976 LCC was estimated at 95.0%. Table 2 reports the error matrix for this classification with overall accuracy, user's accuracy, and producer's accuracy [52]. The error matrix obtained in Ranger (Table S2) was converted to an Estimated population Matrix using the PontiusMatrix42. Figure 4 reports the plots derived using PontiusMatrix42 [74], including the false alarm (commission or false positive) and miss (omission or false negative) components of error divided into the three parameters of Quantity, Exchange and Shift disagreement. PontiusMatrix42 allowed estimation and Comparison of the size of the errors to the size of the unbalanced samples used to validate the LCC [6].
The barren class obtained the lowest producer's accuracy of 89.7%; it presented the largest miss quantity and exchange components. Confusion exists between this class and the deciduous forest and developed areas. Deciduous forest obtained the largest false alarm quantity. The low vegetation class tended to be overestimated; those errors were present, especially in plateau areas characterized by low slope values where the importance of the topographic variables, used in the classifier, decrease. The water class only captured main waterbodies since rivers were generally too narrow to be distinguished.  The 1996 LC was based on seven classes with a 1-m spatial resolution (Table 3). It was obtained using the Ranger package in R with an unbalanced training dataset based on 139 variables derived with the multiresolution segmentation algorithm in eCognition, following the same method used in the LCC of 1976. The validation sample (Table S6) was obtained following the same procedures described in the classification of 1976, and by visual inspection of the aerial images of 1996. The overall accuracy was estimated at 87.3%. Table 4 reports the error matrix of the classification converted to an Estimated population Matrix using the PontiusMatrix42; Table S3 reports the error matrix obtained in Ranger. Figure 5 shows the miss and false alarm components of error derived with Pon-tiusMatrix42. The lowest producer's accuracy was obtained for the developed areas class, resulting from confusion with the barren clear class; overall, the class was underestimated, and misses were the prevalent source of error. Barren clear class presented false alarms, mostly due to confusion with developed areas (impervious surface) and deciduous forest (probably due to the presence of forest clear cuts). Another class that showed a relatively low producer's accuracy was the barren dark class, due to a similar spectral signature with the water class, which included coal sludge impoundments. In the case of barren dark, omission errors were prevalent. In contrast, commission errors characterized the water class. The classifier overestimated the water class at the expense of barren dark and developed areas classes, especially in the flatter plateau areas where slope and curvature variables derived from the DEM had negligible effects.

Class Name Description
Barren clear Non-vegetated areas generally not associated with impervious surface. This class includes surface mine features and non-impervious roads.
Barren dark Surface mining areas characterized by coal presence and large coal deposits.

The LCC of 1996
Three datasets were used to derive the 1996 LC. The main dataset used for the segmentation analysis was a mosaic derived from 35 USGS color infrared (CIR) Digital Orthophoto Quadrangles (DOQs) collected in 1996. Digital images have a 1m resolution and three bands (G, R, NIR). The images were collected as part of the National Aerial Photography Program (NAPP) [67]. The final mosaic was assembled in Erdas Imagine v16.5, where overlapping areas, tonal, and brightness differences were adjusted using the MosaicPro tools. Additional spectral information was added using seven spectral bands obtained with the Landsat 5 Thematic Mapper (TM) on 18 September 1995 [68]. To derive elevation measures and DEM derivatives, the 2003 WV Statewide Addressing and Mapping Board (SAMB) photogrammetrically derived DEM was used [75], since it had the best temporal alignment with the imagery.
The 1996 LC was based on seven classes with a 1-m spatial resolution (Table 3). It was obtained using the Ranger package in R with an unbalanced training dataset based on 139 variables derived with the multiresolution segmentation algorithm in eCognition, following the same method used in the LCC of 1976. The validation sample (Table S6) was obtained following the same procedures described in the classification of 1976, and by visual inspection of the aerial images of 1996. The overall accuracy was estimated at 87.3%. Table 4 reports the error matrix of the classification converted to an Estimated population Matrix using the PontiusMatrix42; Table S3 reports the error matrix obtained in Ranger. Figure 5 shows the miss and false alarm components of error derived with PontiusMatrix42. The lowest producer's accuracy was obtained for the developed areas class, resulting from confusion with the barren clear class; overall, the class was underestimated, and misses were the prevalent source of error. Barren clear class presented false alarms, mostly due to confusion with developed areas (impervious surface) and deciduous forest (probably due to the presence of forest clear cuts). Another class that showed a relatively low producer's accuracy was the barren dark class, due to a similar spectral signature with the water class, which included coal sludge impoundments. In the case of barren dark, omission errors were prevalent. In contrast, commission errors characterized the water class. The classifier overestimated the water class at the expense of barren dark and developed areas classes, especially in the flatter plateau areas where slope and curvature variables derived from the DEM had negligible effects.

The LCC of 2016
The classification of 2016 was derived from a high-resolution state LCC produced by the Natural Resource Analysis Center (NRAC) of WV University in 2018. It is based on the 2016 US Department of Agriculture National Agriculture Imagery Program (NAIP) orthophotography and DEM derivatives obtained from the 10 m USGS National Elevation Dataset (NED). The 2016 LCC has a 1-m resolution and six LC classes ( Table 5). The methods used to obtain the LCC have been fully described in Maxwell et al. (2019) [44]. The sample size with the number of training and validation objects is reported in Table S7. The error matrix from Maxwell et al. (2019) [44] (Table S4) was converted to an Estimated population Matrix using the PontiusMatrix42 ( Table 6). The classes that presented lower accuracies were barren, impervious, and developed areas. Confusion was reported between the barren and impervious classes and among the mixed developed, low vegetation and impervious classes. The plots obtained with the PontiusMatrix42 highlighted the false alarm quantity component for the forest class, while for the low vegetation (grass) reported large omission errors; for the mixed developed areas and barren classes, the exchange components, relating to miss quantities, appeared to be dominant (Table 6, Figure 6). Table 3. The 1996 land cover classes definition.

Class Name Description
Barren clear Non-vegetated areas generally not associated with impervious surface. This class includes surface mine features and non-impervious roads.
Barren dark Surface mining areas characterized by coal presence and large coal deposits.

Deciduous forest
Areas dominated by broad-leaved trees and woodlands.
Evergreen forest Areas dominated by evergreen trees. It can include evergreen shrubs.

Low vegetation
Low vegetation areas such as grasslands, including pastureland, agricultural fields, pastures, and croplands.

Developed areas
It includes areas dominated by mixed development, residential areas and yards; areas characterized by impervious surface such as roads and parking lots. It can also include industrial complex generally used to process coal.

Water
All waterbodies, including rivers, lakes, water impoundments and coal sludge impoundments.

Class Name Description
Forest Areas dominated by tall, woody vegetation and mature forests. This class includes forest and woodlands.
Low vegetation Low vegetation such as grasslands, pastureland, agricultural fields, and croplands.

Barren
Non-vegetated areas not associated with impervious surface. This class includes bare soil, quarries, and surface mine features.

Water
All standing water, including rivers, streams, ponds, lakes, and impoundments.
Impervious surfaces All areas dominated by impervious surface, such as road surfaces, parking lots, airport runways, and buildings.

Mixed developed
Areas dominated by mixed development and mixed land cover, such as residential areas, yards, and development.

Study Area Land Cover Classifications Post-Processing Phase
A post-processing phase followed the LCCs [52]. It consisted of three separate steps: the first one concerned the optimization of the number of classes and resolution to allow their comparison during the change detection phase; the second phase focused on the LCCs of 1976 and 1996, where misclassified portions of the maps were extensively remapped; a final step involved overlaying some vector data to improve the quality of the LCCs.

Optimizing Land Cover Classes and Resolution
The three LCCs were aggregated into five classes (Table 7) to allow their comparison. In the 1976 LCC the number of classes was brought from 6 to 5, aggregating the deciduous forest and the evergreen forest classes. In the 1996 LCC, the number of classes was decreased from 7 to 5. The barren-clear class was combined with barren-dark, and the deciduous forest class was merged with the evergreen forest. In the 2016 LCC, the number of classes was decreased from 6 to 5; impervious and mixed developed classes were aggregated. All classes were aggregated and mapped following a consistent pattern to avoid inconsistencies among the portions of the maps.

Analyzing and Correcting the Misclassified Portions of the Maps
The error analysis derived by the category-level errors for the three LCCs (Tables 2,  4 and 6, and Figures 4-6) was followed by a visual inspection of the maps and correcting their obvious errors. The LCC maps were overlaid with their respective reference orthoimages to remove the most relevant misclassified areas; the areas misclassified by the ML algorithms were reallocated to the more appropriate LC class.
In the 1976 LCC, particular attention was given to removing the barren class omission. A large miss quantity for the barren class in the 1976 LCC corresponded to most of the false alarms reported for the deciduous forest class. Indeed, large exchange components may indicate systematic errors, like misclassification or misregistration, between a couple of categories [59]. Another category that was improved in the 1976 LCC was the

Study Area Land Cover Classifications Post-Processing Phase
A post-processing phase followed the LCCs [51]. It consisted of three separate steps: the first one concerned the optimization of the number of classes and resolution to allow their comparison during the change detection phase; the second phase focused on the LCCs of 1976 and 1996, where misclassified portions of the maps were extensively remapped; a final step involved overlaying some vector data to improve the quality of the LCCs.

Optimizing Land Cover Classes and Resolution
The three LCCs were aggregated into five classes (Table 7) to allow their comparison. In the 1976 LCC the number of classes was brought from 6 to 5, aggregating the deciduous forest and the evergreen forest classes. In the 1996 LCC, the number of classes was decreased from 7 to 5. The barren-clear class was combined with barren-dark, and the deciduous forest class was merged with the evergreen forest. In the 2016 LCC, the number of classes was decreased from 6 to 5; impervious and mixed developed classes were aggregated. All classes were aggregated and mapped following a consistent pattern to avoid inconsistencies among the portions of the maps.

Analyzing and Correcting the Misclassified Portions of the Maps
The error analysis derived by the category-level errors for the three LCCs (Tables 2, 4 and 6, and Figures 4-6) was followed by a visual inspection of the maps and correcting their obvious errors. The LCC maps were overlaid with their respective reference orthoimages to remove the most relevant misclassified areas; the areas misclassified by the ML algorithms were reallocated to the more appropriate LC class. Table 7. Description of the classes adopted in the three LCCs after the post-processing phase.

Class Name Description
Forest Areas dominated by trees and woodlands.

Low vegetation
Low vegetation areas such as grasslands, including pastureland, agricultural fields, pastures, and croplands.

Barren
Non-vegetated areas generally associated with surface mine features, strip-mining roads, spoil banks, coal deposits, dump sites.

Water
All waterbodies, including rivers, lakes, water impoundments and coal sludge impoundments.

Developed areas
Areas dominated by mixed development, residential areas, and yards; areas characterized by impervious surface such as roads and parking lots. Industrial facilities generally used for coal processing.
In the 1976 LCC, particular attention was given to removing the barren class omission. A large miss quantity for the barren class in the 1976 LCC corresponded to most of the false alarms reported for the deciduous forest class. Indeed, large exchange components may indicate systematic errors, like misclassification or misregistration, between a couple of categories [58]. Another category that was improved in the 1976 LCC was the lowvegetation class. It expressed a tendency to gain over the deciduous forest class; this happened mainly in the plateau's relatively flat areas, where slope variables had limited importance in the classification.
In the 1996 LCC, errors between industrial areas and barren-dark and between barrenclear and developed areas classes were generally removed. Moreover, the water class in this classification tended to be overestimated (commission errors), especially in the plateau's flat areas and in the presence of settlements in areas where the number of small impoundments was increased.
Water was consistently remapped in the three LCCs. As previous studies experienced, small water bodies in the region are difficult to map due to confusion with the shadows generated by the forest class [44].
The developed areas class in the three LCCs was subjected to limited editing. Due to the GEOBIA objects from which it was derived and the organization of the residential settlements in the area, this class presented an inherent confusion, particularly with the low-vegetation class. This inherent confusion was challenging to remove in the three LCCs. Moreover, in the study area, abandonment processes in which vegetation often gains on industrial settlement and built-up areas are not infrequent; they have also been documented by specific studies carried out in the area [16]. Therefore, we suggest that the transitions concerning the developed areas' quantities should be considered with caution. They may contain suspected amounts of error that are difficult to assess when combined in multiple classifications. We will focus on this class more carefully in the results phase, using the framework proposed by [6] and accompanying the results' discussion with the high-resolution aerial imagery's visual assessment.

Overlaying Geospatial Data
A final phase consisted of overlaying the three LCCs with the same vector geospatial data representing rivers and lakes. This step was necessary to define, within the threetimepoints, the main water bodies characterized by no change. Indeed, while no relevant changes occurred in the main rivers during the 1976-2016 interval, considerable changes involved the water class, such as the numerous coal sludge impoundments realized for mining uses. Table 7 reports the five classes adopted in the three LCCs, while Figure 7 reports the three LCC maps after the post-processing phase and a histogram with LC class size comparisons. Table S8 reports the quantities of LC classes in the three LCCs after the post-processing phase. Table 7. Description of the classes adopted in the three LCCs after the post-processing phase.

Class Name Description Forest
Areas dominated by trees and woodlands.

Low vegetation
Low vegetation areas such as grasslands, including pastureland, agricultural fields, pastures, and croplands.

Barren
Non-vegetated areas generally associated with surface mine features, strip-mining roads, spoil banks, coal deposits, dump sites.

Water
All waterbodies, including rivers, lakes, water impoundments and coal sludge impoundments.

Developed areas
Areas dominated by mixed development, residential areas, and yards; areas characterized by impervious surface such as roads and parking lots. Industrial facilities generally used for coal processing.

Land-Change and Multi-Level Intensity Analysis Results
The

Multi-Level Intensity Analysis and Difference Components
The intensity analysis relied on 1976, 1996, and 2016 LCCs obtained after the postprocessing phase. The three LCCs identify two intervals of 20 years for a total range of 40 years. The two transition matrices [58], from 1976-1996 and 1996-2016, were obtained with the Land-Change Modeler (LCM) in TerrSet [76], the graphics used to describe land transitions were derived using the intensity.analysis R package [77]. The graphics showing the Difference Components of change (Quantity, Exchange, and Shift) were derived using PontiusMatrix42 and PontiusMatrix41 [78]. The maps showing gains and losses were derived using the crosstab raster generated with the LCM in TerrSet.

Land-Change and Multi-Level Intensity Analysis Results
The results of the two transition matrices, from 1976-1996 and 1996-2015, generated from the LCM in TerrSet are reported in Tables 8 and 9. The rows of the matrices show the classes for the first time point while the columns show the classes for the second time point. The transition matrix diagonal shows the persistence components, representing the areas that did not experience any change. The values off the diagonal correspond to transitions from one category into another. A row with gross loss and a column with gross gain were added in both matrices; those values show the total of row or column minus the persistence component. More details about calculation and matrix interpretation have been given in [60]. Table 10 shows the values of Gross Loss (the sum of gross loss and gross gains), Net change and Swap for the two intervals. Maps (Figure 8) show areas with persistence, gains and losses for the two time intervals.     At the first level of the intensity analysis, the interval level, it is possible to observe that the 1976-1996 interval presents the largest change magnitude, with the percent of interval's domain equal to 17.2%, while 1996-2016 interval has a percentage of change equal to 14.5% of the domain (Figure 9). The change size represents the total amount of change for each interval, expressed as a percentage of the study area. The intensity analysis defines the annual change rate in 1976-1996 as relatively faster, while the annual change rate of 1996-2016 is relatively slower when compared to uniform (Figure 9). Overall processes of change at the interval level are not stationary since they have a different annual rate of change. At the category level, the results were derived using the intensity.analysis package in R. Gain Intensity represents the annual intensity of gross gain (equation 3 in [5]). Loss Intensity is the annual intensity of gross loss (Equation (4) in [5]). The Uniform Category Intensity is the annual intensity of change (Equation (1) in [5]). Graphics ( Figure 10) were obtained with the intensity.analysis package in R, except the graphics (Figure 11) that express Quantity, Exchange, and Shift components of change were obtained with the Pon-tiusMatrix42.
At the category level, during the 1976-1996 interval, the largest gain is reported by forest, while the low vegetation class has the most extensive gross loss. Even though the forest category has the largest gross gain (5559.4 ha), forest is dormant in both loss and gain due to the class's large domain dimension. The most intense change is expressed by the barren class, with loss size almost equal to the gain size. The substantial Exchange components in this class highlight swap processes among categories, pointing at the processes of reclamation of the strip-mining from barren to forest for the losses, and mostly to the new emerging MTR operations for the new gains. The developed areas class expresses an active behavior, with almost equal annual gains and losses. Even in this case, the Exchange component is the most extensive one, indicating that potential errors may characterize this class [59]. The low vegetation class is active, with the annual losses larger than the annual gains. The large Exchange quantities ( Figure 11) highlight consistent substitution processes, but the class is characterized by a significative loss Quantity component. Finally, the latter category that expresses an active behavior in this interval is the water class; it is active only in terms of gain. It is worth remembering that the water category includes coal sludge impoundments, which overall number and size in the region considerably increased after SMCRA (1977).
During the second time interval, from 1996 to 2016, barren class, developed areas, low vegetation, and water classes confirmed their active behavior. Forest remains a dormant class even if it has the most considerable annual change in annual loss and annual gain. Forest maintains its constant growth trend, although the size of losses is greatly increased compared to the previous time interval, mainly due to the new MTR operations that occurred in this time interval. Low vegetation registered a considerable increase in terms of size, while the intensity of annual gains is larger than the annual losses. This behavior is the opposite of the one registered in the previous phase. This process can be attributed to the large new MTR operations, where the artificial grasslands emerge as an effect of the SMCRA reclamation processes. The barren class has the absolute highest annual change intensity, not only in this time interval but also when compared to the previous one, highlighting that very intense changes characterize this class. Its absolute size is diminished due to larger gross losses. The large exchange component indicates processes At the category level, the results were derived using the intensity.analysis package in R. Gain Intensity represents the annual intensity of gross gain (equation 3 in [5]). Loss Intensity is the annual intensity of gross loss (Equation (4) in [5]). The Uniform Category Intensity is the annual intensity of change (Equation (1) in [5]). Graphics ( Figure 10) were obtained with the intensity.analysis package in R, except the graphics (Figure 11) that express Quantity, Exchange, and Shift components of change were obtained with the PontiusMatrix42.
Land 2021, 10, x FOR PEER REVIEW 20 of 32 quantity ( Figure 11). The annual change intensity of the developed areas highlights intense processes of change associated with this class. Compared to the previous time interval, this class expresses a slight contraction in terms of size, with the annual loss size relatively larger than the annual gain. Moreover, processes of gain and loss are now less intense and with an opposite behavior. Overall Net Change reports that developed areas now appear to be shrinking after the relative growth registered in the previous interval, indicating an inverse trend that may be suspicious. Water has almost an equal balance in terms of loss and gain, very close to a uniform behavior, with a slightly positive annual gain change intensity. On the category level, gain and loss processes are stationary during both time intervals for all classes except water. The stationary processes occur when the behavior of gains and losses intensity, compared to the uniform line, is constant across the time intervals. Water class behavior in terms of loss is not stationary; in the first time interval, it has a dormant behavior, while in the second interval, it has an active behavior. At the third level of the intensity analysis, the transition level, only the most noteworthy results have been reported to examine the gains of some classes and recognize classes targeted in both intervals. Some results were omitted since errors may affect the size of some categories and affect the analysis at this level.
The annual transition size for the forest gain in both intervals was mainly obtained At the category level, during the 1976-1996 interval, the largest gain is reported by forest, while the low vegetation class has the most extensive gross loss. Even though the forest category has the largest gross gain (5559.4 ha), forest is dormant in both loss and gain due to the class's large domain dimension. The most intense change is expressed by the barren class, with loss size almost equal to the gain size. The substantial Exchange components in this class highlight swap processes among categories, pointing at the processes of reclamation of the strip-mining from barren to forest for the losses, and mostly to the new emerging MTR operations for the new gains. The developed areas class expresses an active behavior, with almost equal annual gains and losses. Even in this case, the Exchange component is the most extensive one, indicating that potential errors may characterize this class [58]. The low vegetation class is active, with the annual losses larger than the annual gains. The large Exchange quantities ( Figure 11) highlight consistent substitution processes, but the class is characterized by a significative loss Quantity component. Finally, the latter category that expresses an active behavior in this interval is the water class; it is active only in terms of gain. It is worth remembering that the water category includes coal sludge impoundments, which overall number and size in the region considerably increased after SMCRA (1977).
During the second time interval, from 1996 to 2016, barren class, developed areas, low vegetation, and water classes confirmed their active behavior. Forest remains a dormant Land 2021, 10, 748 21 of 32 class even if it has the most considerable annual change in annual loss and annual gain. Forest maintains its constant growth trend, although the size of losses is greatly increased compared to the previous time interval, mainly due to the new MTR operations that occurred in this time interval. Low vegetation registered a considerable increase in terms of size, while the intensity of annual gains is larger than the annual losses. This behavior is the opposite of the one registered in the previous phase. This process can be attributed to the large new MTR operations, where the artificial grasslands emerge as an effect of the SMCRA reclamation processes. The barren class has the absolute highest annual change intensity, not only in this time interval but also when compared to the previous one, highlighting that very intense changes characterize this class. Its absolute size is diminished due to larger gross losses. The large exchange component indicates processes of substitution within this class are dominant in this time interval when compared to net quantity ( Figure 11). The annual change intensity of the developed areas highlights intense processes of change associated with this class. Compared to the previous time interval, this class expresses a slight contraction in terms of size, with the annual loss size relatively larger than the annual gain. Moreover, processes of gain and loss are now less intense and with an opposite behavior. Overall Net Change reports that developed areas now appear to be shrinking after the relative growth registered in the previous interval, indicating an inverse trend that may be suspicious. Water has almost an equal balance in terms of loss and gain, very close to a uniform behavior, with a slightly positive annual gain change intensity. On the category level, gain and loss processes are stationary during both time intervals for all classes except water. The stationary processes occur when the behavior of gains and losses intensity, compared to the uniform line, is constant across the time intervals. Water class behavior in terms of loss is not stationary; in the first time interval, it has a dormant behavior, while in the second interval, it has an active behavior.
At the third level of the intensity analysis, the transition level, only the most noteworthy results have been reported to examine the gains of some classes and recognize classes targeted in both intervals. Some results were omitted since errors may affect the size of some categories and affect the analysis at this level.
The annual transition size for the forest gain in both intervals was mainly obtained at the expenses of grass and barren land ( Figure 12). All the transitions are stationary. While forest targeted barren class in both time intervals, it avoided low vegetation, developed areas, and water. It is worthwhile to mention that even though the low vegetation class experienced large losses due to forest growth, it was not targeted by the forest class. On the contrary, barren land is targeted by forest in both time intervals. In the case of forest gaining from low vegetation, it appears as an avoiding pattern that can be attributable to an overall expansion of forest at the expense of the generally neighboring low vegetation class. While in the case of forest gaining from barren class, it can be attributed to the intentional reclamation processes that characterize the strip mining areas after 1977 (SMCRA) that were generally recolonized by forest.
In the case of the low vegetation class, during both time intervals, the size of the transitions from forest to low vegetation was the largest one. This process can be attributed to the new MTRVF operations concentrated in formerly forested mountain ridges, replaced by new artificial grasslands. Nevertheless, in both time intervals, the most targeted category by grass in terms of annual transition intensity was the developed areas ( Figure 13), highlighting potential processes of abandoning. These results are nonetheless suspicious since the developed area class presents an inherent error component. It is difficult to define how much of this class was targeted because of abandonment processes and settlement reshaping or because the intrinsic confusion derived from the LCCs.  In the case of the low vegetation class, during both time intervals, the size of the transitions from forest to low vegetation was the largest one. This process can be attributed to the new MTRVF operations concentrated in formerly forested mountain ridges, replaced by new artificial grasslands. Nevertheless, in both time intervals, the most targeted category by grass in terms of annual transition intensity was the developed areas ( Figure 13), highlighting potential processes of abandoning. These results are nonetheless suspicious since the developed area class presents an inherent error component. It is difficult to define how much of this class was targeted because of abandonment processes and settlement reshaping or because the intrinsic confusion derived from the LCCs. Further investigations are possible by conducting visual comparisons of the high-resolution images. In the case of the low vegetation class, during both time intervals, the size of the transitions from forest to low vegetation was the largest one. This process can be attributed to the new MTRVF operations concentrated in formerly forested mountain ridges, replaced by new artificial grasslands. Nevertheless, in both time intervals, the most targeted category by grass in terms of annual transition intensity was the developed areas ( Figure 13), highlighting potential processes of abandoning. These results are nonetheless suspicious since the developed area class presents an inherent error component. It is difficult to define how much of this class was targeted because of abandonment processes and settlement reshaping or because the intrinsic confusion derived from the LCCs. Further investigations are possible by conducting visual comparisons of the high-resolution images.  During the first time interval, the developed areas gained mainly from low vegetation while, during the second time interval, they gained mostly from forest. In both time intervals, low vegetation was the most intensively targeted class. Intensity processes are stationary in both time intervals (Figure 14). The growth of settlement occurred mostly during the first time interval in the southern plateau areas characterized by low vegetation. This process is easily recognizable from the visual inspection of the highresolution aerial images.
Land 2021, 10, x FOR PEER REVIEW 23 of 32 During the first time interval, the developed areas gained mainly from low vegetation while, during the second time interval, they gained mostly from forest. In both time intervals, low vegetation was the most intensively targeted class. Intensity processes are stationary in both time intervals (Figure 14). The growth of settlement occurred mostly during the first time interval in the southern plateau areas characterized by low vegetation. This process is easily recognizable from the visual inspection of the high-resolution aerial images.

Land Change Errors Estimation
In land-change studies, errors emerge as a combined effect of the maps' misclassified areas and may propagate when deriving the land change quantities [6,49,51,53]. In this study, due to the specific methods used, the LCCs' errors may propagate due to several factors. Errors may derive from misregistration of the imagery, by errors derived from the segmentation analysis, and finally, by specific omission and commission errors produced by the classifiers.
In the case of image misregistration errors, we reported an RMSE relative to the 1976 orthomosaic equal to 0.89 m. As reported in Pontius (2019) [59], errors due to image misregistration are common and are generally reported in the Exchange component of change. Due to the scale of analysis and the multiple time-points used in the PCC, it is possible to assume that a component of this error is included in this study, even if the large size of the Exchange component can be considered as an effect of the land-change dynamics and not merely derived by errors.
In the case of GEOBIA, the segmentation phases generally produce objects that contain mismatches along the boundaries. Moreover, even if geo-objects are considered homogeneous, the polygons may include fuzziness [80]. In this study, this problem is more evident in the mixed developed class due to the character of settlement composition in the region, where a component of vegetation with trees or grasses is always present in the

Land Change Errors Estimation
In land-change studies, errors emerge as a combined effect of the maps' misclassified areas and may propagate when deriving the land change quantities [6,48,50,52]. In this study, due to the specific methods used, the LCCs' errors may propagate due to several factors. Errors may derive from misregistration of the imagery, by errors derived from the segmentation analysis, and finally, by specific omission and commission errors produced by the classifiers.
In the case of image misregistration errors, we reported an RMSE relative to the 1976 orthomosaic equal to 0.89 m. As reported in Pontius (2019) [58], errors due to image misregistration are common and are generally reported in the Exchange component of change. Due to the scale of analysis and the multiple time-points used in the PCC, it is possible to assume that a component of this error is included in this study, even if the large size of the Exchange component can be considered as an effect of the land-change dynamics and not merely derived by errors.
In the case of GEOBIA, the segmentation phases generally produce objects that contain mismatches along the boundaries. Moreover, even if geo-objects are considered homogeneous, the polygons may include fuzziness [79]. In this study, this problem is more evident in the mixed developed class due to the character of settlement composition in the region, where a component of vegetation with trees or grasses is always present in the geo-objects. It was difficult to isolate and remove this error, even in the post-processing phase, due to its intrinsic characteristics.
Finally, other errors were derived from the supervised classification phase, where the classifier tended to reproduce commission and omission errors depending on the classification problem and the specific dataset. The errors derived from the classifiers were reported in the LCCs sections. Those errors were generally consistently removed in the post-processing phase, but still, some errors may persist.
The framework proposed by Xie et al. (2020) [6] compares the error size of the classifications to the size of change derived from the land-change study. We utilized this approach to estimate the developed areas class's errors by comparing the land change quantities derived from the LCCs of 1976 and 1996 to the estimation derived from the land change results of 1976-1996. Indeed, we suspect that this class may include the largest quantity of errors concerning its size for the reasons above mentioned.
The size of temporal difference (total amount of change) between 1976 and 1996 is equal to 17.2% of the study area (Table 10). The 1976 LCC reported an overall error of around 5.0% (Table 2), while for the LCC of 1996, the overall error was equal to 11.8% (Table 4). Therefore, not all of the difference during the time interval 1976-1996 is real change. In the LCC of 1976, the developed areas class had user's accuracy (commission) equal to 93.9% (commission error equal to 6.1%) and producer's accuracy (omission) equal to 96.4% (omission error equal to 3.6%); the commission error was larger than the omission error, then the map overestimated this category. On the contrary, in the classification of the 1996 omission errors (23.1%) of the developed areas class was larger than commission error (14.9%). Therefore, the category was underestimated.
The gross gains in 1996 for the developed areas class were equal to 973 ha. The exchange component in the developed areas class during 1976-1996 (Figure 11d), which generally highlights most of the errors derived from misclassifications and misregistration [58], was equal to 882 ha. The difference between these two quantities indicates that even in the worst scenario where all the Exchange quantities were errors, the developed areas class during the interval 1976-1996 was gaining. Figure 15a-c reports the LCCs' category level errors intensity before the post-processing phase and estimated as a percentage of the validation objects. Graphics were obtained with PontiusMatrix41; the uniform line in the figures represents the overall LCC error [6].
Nonetheless, it should be noted that the gaining of developed areas expresses the overall tendency of the class at the scale of the study area. It does not mean that, during this temporal interval, there was a uniform tendency of built-up areas in gaining across all the catchments of our study area; instead, we should expect opposite trends of contraction and expansion, but still, overall, the class experienced gain. In the LCC of 2016, developed areas were classified with relatively low accuracies (58.8% user's accuracy and 55.5% of producer's accuracy); therefore, we should expect that during the second time interval , there were more errors in land transitions that include this class. During 1996-2016, the developed areas reported gross loss larger than gross gains. It indicates built-up areas were shrinking; this inversion may point at potential abandonment processes that can be supported by the fact that low-vegetation targeted the developed class during this interval, as shown by the intensity analysis (Figure 13d). Nevertheless, the transitions are suspicious, even because of the large Exchange component, equal to 1268 ha (Figure 11d). More details about changes in this class are reported in the discussion sections, where inferences were derived from the observation of historical maps, LCCs, and high-resolution imagery.

Methods
The method proposed in this study presented a complete framework based on a PCC approach to quantify relevant components of temporal land transitions. The method combined the use of LCCs based on high-resolution historical images with the framework of Intensity Analysis and Difference Components [6,59,81]. Except for limited examples [82], previous studies based on Intensity Analysis have generally grounded on low-resolution satellite images with few or no data about the LCCs error matrices [83][84][85][86][87][88][89].
The high-resolution images allowed us to collect the validation samples from the geospatial polygons derived from the segmentation analysis and assign them the best membership to the specific classes. This approach allowed us to obtain the error matrices for the LLCs of 1976 and 1996. The inclusion of the errors derived from the historical LCCs achieved in-depth information about errors propagation across the LC classes identifying the suspicious transitions within the land change analysis. Moreover, the PontiusMatrix42 allowed estimating false alarms and missing components of errors within each LCC. The inclusion of further ancillary datasets, like vector polygons, may be implemented within this method to improve the LCCs [45]. In the US, surface mining databases are being completed [90]; similar information can be obtained in other countries from historical maps or to apply this method to other LCC problems.

Methods
The method proposed in this study presented a complete framework based on a PCC approach to quantify relevant components of temporal land transitions. The method combined the use of LCCs based on high-resolution historical images with the framework of Intensity Analysis and Difference Components [6,58,80]. Except for limited examples [81], previous studies based on Intensity Analysis have generally grounded on low-resolution satellite images with few or no data about the LCCs error matrices [82][83][84][85][86][87][88].
The high-resolution images allowed us to collect the validation samples from the geospatial polygons derived from the segmentation analysis and assign them the best membership to the specific classes. This approach allowed us to obtain the error matrices for the LLCs of 1976 and 1996. The inclusion of the errors derived from the historical LCCs achieved in-depth information about errors propagation across the LC classes identifying the suspicious transitions within the land change analysis. Moreover, the PontiusMatrix42 allowed estimating false alarms and missing components of errors within each LCC. The inclusion of further ancillary datasets, like vector polygons, may be implemented within this method to improve the LCCs [44]. In the US, surface mining databases are being completed [89]; similar information can be obtained in other countries from historical maps or to apply this method to other LCC problems.
The integrated framework of Intensity Analysis and Difference Components, through a group of free and accessible tools, enables researchers to consider the components of errors and assess whether those errors may affect the results of the temporal analyses. The method makes it possible to exercise control over the error in land-change analysis based on PCC, encouraging the recognition of the more intense processes of change.

Landscape Transitions in the Coal River Headwaters
The findings revealed that the land change processes derived through the analysis of five relevant LC classes were not limited to surface mining sites and industrial facilities. Although coal mining has been the dominant driver of change in the region [7], its legacies on the landscape have interacted with other dynamics that are the outcome of previous and enduring socio-ecological processes.
Intensity analysis allowed us to compare the rate of change of the two time intervals , showing us how, during the first interval, land-change processes were faster since they involved a larger percentage of change.
At the category level, forest class gained mainly during the first time interval with transitions characterized by different intensity levels. The barren class was continuously targeted by forest across the two intervals. A large part of this new forest grew over the former strip-mining barren areas; a process that is clearly distinguished by the LCCs and the high-resolution photographs. The intensity analysis highlights the non-uniform nature of this process, distinguishing it from a simple random event. The possible explanation identifies it as the outcome of intentional reclamation stages initiated before or after the SMCRA enactment. Limited information is available about the characteristic of these forested areas grown on barren mining land. In 1951 a USDA report [29] examined 105 plots in WV, describing the occurrence of natural tree and shrub species that generally spread out in the strip-mine spoils areas and above the overburdens. In those plots, black locusts were the most widespread trees, followed by red maples, black cherries, and sassafrasses. These reclamation processes were similar to those described by Townsend et al. (2009) [7] in areas where MTR never occurred.
During the first time interval, the low-vegetation class was the largest class from which forest gained, while during the second interval, it was the second larger one. Nevertheless, in both time intervals of the intensity analysis, the forest class did not target the lowvegetation class. Moreover, our results did not report any relevant forest growth within the surface mining permits of MTR operations [90]; these areas in the 2016 LCC are predominantly characterized by barren and low-vegetation classes.
During the first time interval , we observed an overall transition from barren to forest class that generally occurred on former strip-mining areas that were already mined in 1976-1977. During 1996, inside the permit boundaries, reclamation processes of former strip mining areas often occurred in conjunction with new surface minings. These new operations already presented the characters of the large surface mining operations typical of MTRVF, with dominant barren, low-vegetation, and water (sludge impoundments) classes. In the 2016 LCC, the areas characterized by MTRVF operations became much more considerable. The expansion of these areas occurred both on the predominantly undisturbed mountain ridges covered by forest and the former strip-mining areas characterized by ongoing reclamation processes and young trees. This study confirms that the areas characterized by former surface mines can be distinguished into at least two broad groups. The first group consists of the former strip-mining areas on which the forest class has apparently assumed stable characters; the second group is characterized mainly by MTRVF areas, which may overlap with previous mining operations. These operations are characterized by "remining" practices generally allowed because the reclamation was considered not completed [36]. Further considerations and in-depth studies could be based on more detailed information derived from the analysis of the mining permit boundaries since they generally present very articulated developments.
It is worth noting that the overall quantity of barren land from 1976 to 2016 diminished by approximately 37% in the study area. The loss of grasslands in the valley floors and the plateau was paradoxically compensated by the newly reclaimed areas, where novel ecosystems [91,92] characterized by poor topsoils spread with non-native and fast-growing grasses [93] emerged as outcome of the regulatory reclamation processes. Moreover, surface mines' distribution changed from a spatially sprawled and more superficial one to a more compacted and deeper pattern. Recent studies based on DEMs measured the geological magnitude of these processes on a regional scale [23].
With regard to the general contraction of the low vegetated areas, it is worth noting that this process should not be directly ascribed only to coal mining. Indeed, it should probably be seen in continuity with the overall agricultural land abandonment processes reported in the Appalachian region at least from 1910 [94]. Among the variety of factors reported, the stagnation of farm productivity in Appalachia [15] and the decline of the traditional forest farming system [95] should be seen as those responsible for facilitating the introduction of new competitive practices of land use and the acquisition of land by the extractive industries of wood and coal.
The changes observed in the civilian and industrial settlements (the developed areas LC class) were articulated and followed fluctuating dynamics within the different subcatchments of the study area. Even though this study presents some limitations due to the intrinsic fuzziness reported in classifying this category, it is worth reporting some patterns of changes in this class, because of the overall lack of studies conducted with similar methods in this region. The settlements' dynamics of contraction and expansion in the area were mainly related to the economic cycles of growth and recession followed by the coal industry. During the interval 1880-1930, the first phase of rapid industrialization required labor to sustain the growth of coal industries [96]. Coal companies created several new towns in the Southern Appalachia, called company towns or coal towns. Following the mid-1920s, the coal industries bust contributed to the company towns' progressive decline and, after the 1950s, to their definitive closure and dismantling [97]. During 1950-1970, WV lost 13% of the population [98]. The intensification of coal production regained strength after 1962, with the second phase of strip-mining intensification after 1970 [96]. Nevertheless, the population numbers in most old towns connected to the coal economy never regained the values reached during the first industrialization phase. Moreover, the intensification of surface mining operations led to consistent effects for the settlements located in the valley floors that experienced "the devastating externalities of coal mining" ( [96], p. 28). As emerged from the intensity analysis, during both intervals (1976-1996 and 1996-2016), developed areas were the most targeted class by low vegetation. We already reported that this effect could be attributed to the fuzziness between low-vegetation and the developed areas class. Nevertheless, processes in which settlements were partially substituted by low vegetation class due to their abandonment or dismantlement have been recently documented in the study area [16]. In the lower catchments of the study area characterized by surface mining activities, we observed a pattern of settlement densification in the proximity of the most downstream areas. Settlement abandonment processes characterized the areas near the stream's headwaters and closer to surface mining operations. On the contrary, in the plateau area, close to the town of Beckley, the developed areas class was characterized by a constant gain, with a more vigorous intensity during the first time interval. Here, the widespread phenomenon of suburbanization observed in the US and the more suitable landscape morphology may explain the processes of settlement intensification [99,100]. By way of example, the process of establishment and abandonment of the settlement of Montcoal during the interval 1912-2016 is reported in Figures S1 and S2 together with some detailed images derived from the LC maps. The settlement was established with the opening of a mine in Montcoal Mountain in 1916; it was one of the leading coal towns founded in the Coal River's headwaters [96,97].

Conclusions
This article combined Intensity Analysis, Difference Components, and error analysis to assess land change processes applying a previous framework [6]. Additionally, this study provided three LCCs with accuracies assessments. The results confirm the method's Land 2021, 10, 748 28 of 32 validity and reliability. Its application demonstrated that not all the temporal changes in the study area indicate actual landscape change. Moreover, the framework coupled with category accuracies from the historical classifications provided more in-depth details about how errors may propagate. The findings of this study allowed us to improve the characterization of the study area, deriving specific LC transitions and classes, like those described for the reclamation processes, and improving previous limited information about the land change processes before and after the SMCRA enactment in the region in 1977. Our study detected that the extent of disturbed land by surface mining in the study area increased by 38.9% compared to previous research.
The combined framework of Intensity Analysis and Difference Components allowed us to isolate and discuss suspicious transitions. A further step to the proposed approach, to produce a more complete framework, would be to integrate an additional error assessment phase following the LCCs post-processing. This would help to further detail all suspicious transitions and provide a more comprehensive quantitative assessment after the postprocessing phase, which, although improved the overall accuracy, introduced an additional degree of indeterminacy. In addition, we suggest that given the multiplicity of co-occurring land-change dynamics and drivers generally present in the landscape, further analyses should include the comparison of different regions, catchments, or watersheds, to derive specific trends of growth or shrinkage and aggregate statistics using different spatial zones; some studies used this method using Intensity Analysis [81,84].
Information derived by land change studies can be of great importance for natural resource managers and future planning phases. It allows the characterization of large areas and the definition of historic-based and spatially explicit datasets. These data can inform models that can assess or readdress the efficacy of specific policies and inform conservation and natural resource management. In addition, assessing the analysis of the response produced by the socio-technological [101,102] reclamation practices can help in analyzing possible institutional failures or the achievements of specific policies [103].
Our study confirms the importance of land-change studies based on high-resolution and spatially explicit datasets. Indeed, detailed spatial scales of analysis may disclose details about overall landscape dynamics that characterize large regions. In the case presented, anthropogenic processes profoundly influenced the environmental and ecological dynamics of the study area, including relevant changes in ecosystem services provision, carbon sequestration, soil degradation, hydrologic cycles, water quality, habitat integrity, and biodiversity loss [104,105]. The study demonstrated the importance of site-specific studies and land-change analyses when dealing with deeply altered and fragile landscapes where human-based drivers, and not only the environmental gradients, have been responsible for changes that profoundly altered the landscape, its morphology, and the ecosystem structure and composition.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/land10070748/s1, Table S1: Synoptic table of remote sensing and land change studies concerning surface mining in the Central Appalachian Region.