1. INTRODUCTION
In the Southern Coalfields of West Virginia, MTR/VF is currently the leading cause of land cover change [1-4]. Multiple watersheds in West Virginia have more than 10% of their surface area disturbed by surface mining [5], which results in the loss of forest and a conversion to barren land cover [6]. It has been estimated that all surface mining in Appalachia has resulted in a net loss of 420,000 ha of forest [3], and the interior character of the forest is threatened by the introduction of non-forest edges [7]. During the MTR/VF process, forests are cleared, top soil is removed, overburden material is blasted away to uncover coal seams, and overburden material is placed in adjacent valleys, filling stream segments and creating valley fills [5,8]. In Kentucky, it has been estimated that greater than 660 km of headwater streams were buried between 1985 and 1999 [9]. Later reclamation produces grasslands or shrub/scrub land cover; however, productivity is often limited due to poor soil conditions [6]. Attempts are often made to preserve soil; however, it becomes homogenized and soil horizons are not maintained [10].
In addition to the land use and land cover (LULC) change, the landscape and terrain are recontoured with modified watershed ridges, altered vegetation conditions, and modified soil character [8]. References [1,11] estimated that surface mining was responsible for displacing more material in the Southern Coalfields of West Virginia than river systems and natural geomorphic processes. Furthermore, MTR/VF methods have resulted in more moved material and faster landscape alterations as compared to traditional surface mining methods, such as auger, contour, and highwall mining [12].
Elevation differencing of Digital Elevation Models (DEMs) representing different temporal conditions has been explored to describe topographic change. For example, [13] utilized NASA Shuttle Radar Topographic Mission (SRTM), US Geologic Survey National Elevation Dataset (NED), and Terrain Resource Information Management Program (TRIM) data to describe and map alpine glacier changes in southeast Alaska and northwest British Columbia. Reference [14] utilized LiDAR-derived DEMs to calculate dune volume changes over a 1-year period at sites along the Cape Hatteras National Seashore.
Landscape and geomorphic change resulting from MTR/VF disturbance were specifically analyzed in Perry County, Kentucky, USA [15] utilizing NED (pre-mining) and SRTM (post-mining) data. The study highlighted the complexity of such an analysis when timestamps of the NED data are variable since they were created relative to best available data and 1:24000 scale topographic quadrangle maps. In a similar study, the West Virginia Department of Environmental Protection (WVDEP) utilized interferometric synthetic aperture radar (IFSAR) and elevation raster data produced from digital line graph (DLG) hypsography to map valley fill extents throughout nine counties in southern West Virginia. Because the radar data did not adequately penetrate the tree canopy, it was necessary to remove forested areas from the analysis. It was found that a complete inventory of the fills required additional visual classification [16].
To the best of our knowledge no previous work has quantified the post-mining landscape in terms of changes in terrain characteristics from pre-mining conditions using landscape-scale categorical terrain data. An understanding of the terrain alteration using landform data is appropriate to assess the impact of the topographically altered mountaintops.
This paper expands upon earlier differencing and topographic change work in the MTR/VF region by implementing a methodology relying on a categorical representation of the landscape as landforms. This data differentiates the landscape into the following classes: cliff, steep slope, slope crest, upper slope, flat summit, sideslope, cove, dry flat, moist flat, wet flat, and slope bottom. This method was adopted after [17] because such features adequately represent the Southern Coalfields at the landscape scale.
Attempts have been made to link landform data to habitats using predictive modeling. For example, [17] attempted to link the ecological community types developed by the Nature Conservancy to landforms, elevation, and lithology. These categories represent landforms of ecological significance at the landscape scale. Our goal was to investigate how the distribution of these categories was impacted by surface mining. Our results may help to reestablish a post-mining landscape that can benefit terrestrial habitat, ecosystem health, and biodiversity since these have been shown to be linked to landform clasess [17].
2. METHODOLOGY
2.1. Study Area
The Coal River Watershed is a 230,755 ha Hydrologic Unit Code (HUC) 8 headwater watershed completely within West Virginia, USA. The Coal River Watershed exists within the Appalachian Plateau physiogeographic province, a dissected, westward-tilted plateau dominated by Pennsylvanian bedrock. Pennsylvanian stratigraphy is characterized as cyclic sequences of sandstone, shale, clay, coal, and limestone [18]. The terrain is dissected by a dendritic stream network and shows fine texture with moderate to strong local relief. In comparison to the northern Appalachian Plateau, the Southern Coalfields is generally more rugged due to resistant strata [19]. The terrain analysis using landforms suggests that this terrain is naturally dominated by steep slopes.
According to LULC estimates derived from aerial photography, 8.8% of this watershed was disturbed by active surface mining or mine reclamation in 2009. It should be noted that this estimate does not take into account historical mining areas which have since been reforested. Figure 1 shows the watershed location within West Virginia. This watershed was selected as a case study because it is heavily impacted by mining, there has been continued mining and reclamation between 2003 and 2010, and because LiDAR data were available for the extent.
2.2. Digital Elevation Data Utilized
The LiDAR data, representing recent topographic conditions, were collected between April 9th, 2010 and April 18th, 2010 at a flight height of 1524 meters above ground level (AGL), a pulse frequency of 70 kHz, a scan frequency of 35 Hz, and a scan angle (full field of view) of 36˚. The swaths were flown with a 30% overlap, at an average speed of 250 km/hr, and with an average width of 979 m. An Optech ALTM 3100 C sensor was used to collect the data. The scan and flight specifications were selected to support a 0.7 m contour interval and a 1 m nominal ground post spacing. Ground points were filtered from the all return data to produce ground point files in LAS 1.2 format.
An interpolated raster surface was created from point data using ArcMap 10, and average point spacing of 0.01 m. Inverse distance weighting (IDW) was used to interpolate a raster surface since sample points were in an evenly distributed pattern. This process resulted in a 1 m resolution, floating point elevation dataset.
2.3. Assessing Landform Changes between 2003 and 2010
Landform changes were assessed in filled and cut or excavated areas resulting from mining and reclamation between 2003 and 2010. The 2010 LiDAR-derived DEM was compared to a 2003 photogrammetrically-derived DEM representing 2003 conditions. The 2003 DEM data were provided by the West Virginia Statewide Addressing and Mapping Board (SAMB) who contracted BAE Systems ADR to create a stereo photogrammetric-derived digital terrain model (DTM) from statewide, spring, leaf-off, 1:4800 scale aerial photography. The DEM was created in compliance with National Dataset standards (1/9th arc second) and produced from mass points and breaklines. This dataset supports a vertical accuracy of +/− 3.048 m, root mean square error (RMSE) and has a cell size of 3 m [20].
In order to match the resolution of the SAMB DEM, we resampled the LiDAR-derived raster to 3 m resolution using bilinear interpolation and snapped the grid to match the extent of the SAMB DEM. This process assured that the raster grids were completely aligned. This process resulted in two 3 m DEMs representing 2003 and 2010 conditions that could be compared. A vertical data transformation was unnecessary since both DEMs had a vertical reference of NAVD88 orthometric.
In order to detect systematic error between the 2003 and 2010 DEM data, a total of 281 points were collected in the field throughout the Coal River Watershed on flat, paved surfaces using Pacific Crest realtime kinematic (RTK) survey equipment. The elevation measurements from the 2003 and 2010 elevation raster data were obtained at these locations. The measurements from each DEM were then compared. The mean difference was −0.4 m with a maximum difference of 0.8 m and a minimum difference of −2.0 m. Based on this analysis we did not correct for systematic difference between DEMs.
Once DEMs were obtained, processed, and prepared, they were subtracted using Raster Calculator within the Spatial Analyst Extension of ArcMap [21]. This produced a grid of elevation differences throughout the watershed. Negative values indicated potential cuts or excavations while positive values indicated potential fills. However, difference could simply have resulted from errors in the DEMs or in methodology. As a result, once this elevation difference model was calculated, it was necessary to determine a tolerance or threshold that would constitute true change and not simply error or noise between the digital elevation datasets. The photogrammetric DEM had an error tolerance (RMSE) of +/− 3.048 m while the LiDAR DEM had and error tolerance of only 15 cm (0.15 m). An equation suggested by [22] was used to estimate this threshold. The RMSE of each dataset was squared, the results were summed, and the square root of the sum was then taken. The result was then multiplied by three to determine a cut off representing values that were greater than three standard deviations from the mean. This method assumes a Gaussian distribution. Elevation change measurements outside of this range were considered true elevation differences and not a result of error or noise. Reference [15] used the same methodology to derive an error tolerance for the analysis conducted in Perry County, Kentucky, USA.
According to this method, it is not certain whether differences less than +/− 9.2 m represented true topographic change. As a result, the elevation difference grid was reclassified so that values between −10 m and +10 m were considered no change, error between the DEMs, or noise. Pixels with values less than −10 m were considered potential cuts while pixels with values greater than 10 m were considered potential fills. This process produced a cut and fill mask.
A LULC dataset was created for the region from 2009 National Agriculture Imagery Program (NAIP) orthophotography. Optimally, imagery collected at the time of the LiDAR collection would have been utilized; however, such data were not available. The imagery was classified using an object-based feature extraction methodology, augmented with GIS decision rules and manual digitizing. Forest, grass, and barren land cover were extracted from the raw imagery.
Combining these data with mine permit boundaries from the West Virginia Department of Environmental Protection (WVDEP) made it possible to delineate barren and grassland land cover in mine permits. Grasslands in permits represent potential reclamation while barren areas represent areas of active mining or areas that have not yet been reclaimed. The extents of valley fill faces and slurry impoundments were digitized. According to an error assessment based on manual aerial photograph interpretation at 100 randomly selected points, the LULC dataset had an overall accuracy of 94% (KHat of 93%).
It was possible to further remove erroneous pixels as potential cut or fills by utilizing additional data, such as the high resolution land cover. For example, potential cut or fill pixels could be removed if they were outside of disturbance resulting from surface mining. The potential cut or fill pixels that existed within areas of mine disturbance (barren, grasslands, or valley fill faces) were considered. Slurry impoundments were excluded because of errors associated with water level. To complete this analysis, it was necessary to convert the raster data to polygons that represented contiguous areas of cut or fill. No smoothing or simplification was applied.
Small areas of topographic change could also be removed as error or noise. Thresholds were determined from the hand digitized valley fill faces. The average size of these faces was 73,257 m2 (7.3 ha) with a standard deviation of 67,013 m2 (6.7 ha). Based from this distribution, any contiguous area of cut or fill that was smaller than 60,000 m2 (6 ha) was considered noise and removed. Although this threshold was somewhat arbitrary, it was selected since this study focused on detecting large expanses of excavation and fills.
The result of this analysis was a vector layer of potential cut or fill extents in which topographic change had occurred between 2003 and 2010. The extents were used to define areas of potential terrain change for analysis and comparison of the distributions of landforms. Although there was error associated with this methodology, it was adequate for this analysis since we were not interested in conducting an assessment of absolute elevation change or volumetric change but in defining areas in which to compare categorical terrain data and their impact on forest communities.
2.4. Creation of Landforms
The DEM data were used to classify the terrain into the following landforms: cliff, steep slope, slope crest, upper slope, flat summit, sideslope, cove, dry flat, moist flat, wet flat, and slope bottom. The methodology described by [17] provides a means to utilize DEM data to classify the landscape into different units of ecological significance. DEM derivatives including slope in degrees, a hydrologically filled DEM, flow direction, and flow accumulation were created using ArcMap 10. The slope and flow accumulation grids were utilized to calculate a moisture index. The moisture index is a relative measure of the moisture of a specific cell, and it assumes that the moisture level is a function of how much water flows into the cell, predicted by flow accumulation, and how fast the water can flow out, described by slope [20-25].
Landscape position was calculated following [26] to divide the landscape into the following categories: ridge, wide ridge, slope/flat, slope/cove. The approach is based on a local neighborhood analysis of elevation values in which a cell is compared to the mean elevation of all values within 3 by 3 windows. Landscape position was combined with slope and moisture index data to derive the final landform classes. This procedure was conducted for both the 2003 and the 2010 DEMs. The comparisons of the derived landforms included:
• Areas filled between 2003 and 2010 (mining related);
• Areas cut or excavated between 2003 and 2010 (mining related);
• Areas of mine reclamation in 2009, areas of mining disturbance or reclamation in 2009, and areas in WVDEP mine permits but not disturbed (forested in 2009).
This enabled the comparison of post-mining landform categories with pre-mining conditions.
3. RESULTS
Figure 2 shows the distribution of landforms in a MTR/VF site representing 2010 conditions. Comparing the mine site with the surrounding landscape, the distribution of landforms is greatly altered by mining even after reclamation has occurred. Figure 3 shows a different MTR/VF site and offers a comparison of landforms within areas cut or filled between 2003 and 2010, and