Remote Sensing of Forest Structural Changes Due to the Recent Boom of Unconventional Shale Gas Extraction Activities in Appalachian Ohio

Abstract

Dense unconventional shale gas extraction activities have occurred in Appalachian Ohio since 2010 and they have caused various landcover changes and forest fragmentation issues. This research investigated the most recent boom of unconventional shale gas extraction activities and their impacts on the landcover changes and forest structural changes in the Muskingum River Watershed in Appalachian Ohio. Triple-temporal high-resolution natural-color aerial images from 2006 to 2017 and a group of ancillary geographic information system (GIS) data were first used to digitize the landcover changes due to the recent boom of these unconventional shale gas extraction activities. Geographic object-based image analysis (GEOBIA) was then employed to form forest patches as image objects and to accurately quantify the forest connectivity. Lastly, the initial and updated forest image objects were used to quantify the loss of core forest as the two-dimensional (2D) forest structural changes, and initial and updated canopy height models (CHMs) derived from airborne light detection and ranging (LiDAR) point clouds were used to quantify the loss of forest volume as three-dimensional (3D) forest structural changes. The results indicate a consistent format but uneven spatiotemporal development of these unconventional shale gas extraction activities. Dense unconventional shale gas extraction activities formed two apparent hotspots. Two-thirds of the well pad facilities and half of the pipeline right-of-way (ROW) corridors were constructed during the raising phase of the boom. At the end of the boom, significant forest fragmentation already occurred in both hotspots of these active unconventional shale gas extraction activities, and the areal loss of core forest reached up to 14.60% in the densest concentrated regions of these activities. These results call for attention to the ecological studies targeted on the forest fragmentation in the Muskingum River Watershed and the broader Appalachian Ohio regions.

1. Introduction

The active global market has created a huge demand for international fossil fuels and has promoted the new generation of extraction technologies. The advent of horizontal hydraulic-fracturing drilling technology led to the boom of shale oil and gas extraction in the United States [1]. Compared to the first generation of oil/gas extraction technology (pump-jacks) and the second generation (class II enhanced oil recovery wells [2,3]), the third generation of oil/gas extraction technology, which is generally referred to as the unconventional shale oil and gas wells [4], requires a more comprehensive extraction system (Figure 1). The unconventional shale oil and gas extraction system comprises a flat and open well pad with a size of about 100 m by 150 m, access roads connecting the well pad to the local transportation network, and a gathering pipeline connecting the well pads to compressor stations. Under a well pad, there is a containment vault in which the wells vertically drill approximately 1 mile deeper than freshwater aquifers and then conduct hydraulic fracturing with 1–2 miles horizontal drilling distance. The extracted gas is transmitted to compressor stations through underground pipelines between the well pads. Other facilities such as gathering and boosting stations, gas processing plants, and transmission compressor stations are constructed along the gathering and transmission pipeline network [5]. The construction of the unconventional shale oil and gas extraction system imposes direct impacts on the landcover and local environment (Figure 2). The well pad directly converts the vegetated landcover to the impervious surface, and retention ponds and other landcover disturbances can be observed surrounding the well pads [6]. The access roads commonly have a gravel surface and the landcover immediately adjacent to the access roads is changed. A buffer zone is determined with a width of 20 to 50 m along each pipeline between well pads, known as pipeline right-of-way (ROW) corridors [7,8]. Trees in pipeline ROW corridors are cut and cleared for the safety of the gathering pipelines.

Marcellus and Utica Shale are two of the largest and the most prolific formations in the northern Appalachian region in the continental United States. These two formations have experienced intensive extraction of shale oil and gas over the past decades [1,9]. By the end of 2018, 236,272 permits for oil or gas wells were issued in Ohio, and they were concentrated in the Appalachian counties in the east of Ohio [10]. Approximately 102,794 oil or gas wells (43.51%) are located in the Muskingum River Watershed. The construction of horizontal drilling pipes, hydraulic fracturing infrastructures, class II enhanced oil recovery wells, well pads, and other associated facilities have induced considerable landcover changes and forest fragmentation issues [11,12], water contamination issues [13,14,15,16], air pollution issues [17], ecosystem problems [18,19,20,21], and health problems [22].

This research focuses on the impacts of the recent boom of the unconventional shale gas extraction activities on forest structure in the Muskingum River Watershed in Appalachian Ohio. Given the data and statistics from Ohio Department of Natural Resources (ODNR), more than 3000 horizontal drilling shale gas wells were permitted in Ohio as of 2018, with 2115 of them actively producing [23]. These wells are aggregated within 801 well pads, and 356 of them are located in the Muskingum River Watershed. These well pads and their associated facilities, access roads, and gathering pipelines have left a bunch of footprints on the land surface and have opened many large, continuous core forests that are critical habitats for multiple species of forest-dependent wildlife [24]. The relevant deforestation activities directly cause the degradation of several key ecosystem services, including carbon storage in the format of biomass and in soils, as well as the regulation of water balance and river flow [25]. A systematic and quantitative investigation of these landcover changes and forest structural changes is necessary to inspect and model the magnitude of degradation in these key ecosystem services. For instance, landcover changes are critical to model the changes in drainage discharge and the suspended sediments carried by streamflow using distributed hydrological models [26,27], two-dimensional (2D) forest structural changes evaluate the capability changes of the forest as core habitats for wild species, and three-dimensional (3D) forest structural changes are useful to estimate the dynamics of forest biomass and the net release of carbon dioxide (CO₂) from forest to the atmosphere [28]. Great efforts have been devoted to elaborately investigate the landcover changes caused by these unconventional shale gas extraction activities, and most of them targeted the landcover change patterns on a county level [12,29]. However, regarding the increasing concerns related to the impacted ecological services, such as the habitat fragmentation and surface water safety, the investigation of landcover changes calls for more attention to natural geographical units, for instance, watersheds and catchments.

The objectives of this paper were to (1) quantify the spatial footprints of these unconventional shale gas extraction facilities in the Muskingum River Watershed and associated direct loss of forest area, (2) quantify the forest fragmentation conditions and their changes due to the unconventional shale gas extraction activities, and (3) estimate forest structural changes caused by the relevant deforestation, including both two-dimensional (2D) structural changes (the areal loss of core forest) and three-dimensional (3D) structural changes (the loss of forest volume). The rest of this paper first introduces the study area and the datasets, and then presents the method to digitize the footprints of these unconventional shale gas extraction activities and direct loss of forest area. Object-oriented classification was applied to the first-time node aerial images to capture the initial conditions of landcover, while the activity footprints were used to update the landcover conditions after changes. A set of landscape ecology metrics derived from forest image objects were applied to test the occurrence of forest fragmentation. Once significant forest fragmentation was verified, the initial and updated forest image objects were used to quantify the areal loss of core forest as the 2D forest structural changes. Similarly, the initial and updated canopy height models (CHMs) derived from airborne light detection and ranging (LiDAR) point-cloud data were used to quantify the loss of forest volume as 3D forest structural changes. Next, the reliability of the analysis and the possible ecological impacts of these activities are discussed. The final section draws the main conclusions.

2. Materials

2.1. Study Area

The Muskingum River Watershed is located in the east of Ohio. This watershed is defined as a hydrologic unit at level 4 in the United States Geological Survey (USGS) Watershed Boundary Dataset (WBD), and it has a total area of 20,849.43 km² (about 20% of the state) [30,31]. The Muskingum River Watershed completely or partially covers 27 counties of Ohio (Figure 3). Its population was about 2.73 million in 2018, mostly concentrated in small cities, e.g., Akron in Summit County, Canton in Stark County, and New Philadelphia in Tuscarawas County [32].

The Muskingum River Watershed sits at the west edge of the Appalachian Mountains, and it is a part of the Allegheny Plateau [31]. The Muskingum River Watershed has a humid subtropical climate type [33,34]. The average daily temperature is about 10.3 °C while the annual total precipitation is about 1075 mm [30]. The precipitation and temperature change significantly throughout the year, forming a wet warm summer and a dry cold winter. The topography in the watershed is low-lying in the south and mountainous on the other three sides. The Muskingum River is the dominant river in the watershed, and it flows into the Ohio River through the outlet at Marietta in Washington County. The Muskingum River has four main tributary rivers: the Tuscarawas River that originates from the northeast of the watershed, the Walhonding River that originates from the northwest of the watershed, the Linking River that originates from the west of the watershed, and the Wills Creek that originates from the east of the watershed [30,35]. The Tuscarawas River and the Walhonding River coalesce into the Muskingum River at the confluence near Coshocton, the Wills Creek flows in at the confluence 18.4 km downstream Coshocton, and the Linking River flows in at the confluence near North Zanesville (Figure 3a). The mean annual discharge of the Muskingum River was approximately 305 m³/s from 2013 to 2019, with its peak in April and minimum flow in October [36]. Most of the counties have an aquatic habitat areal coverage of less than 1.5% (

Table 1. Appalachian counties in the Muskingum River Watershed and available aerial imagery.

County Name	County Area (km2) 1	Forest Area (%) 2	Aquatic Habitat Area (%) 2	Capture Year of Available Aerial Images
				Start Year	Middle Year	End Year
Ashland 3	1111.36	32.83	0.87	2006	2013	2017
Athens	1315.09	73.87	0.62	No drilling	-	-
Belmont	1397.58	58.42	1.06	2007	2014	2017
Carroll	1033.41	55.37	1.34	2006	2014	2017
Columbiana	1386.58	46.16	0.76	2006	2012	2017
Coshocton	1467.65	58.40	1.00	2006	2014	2017
Crawford	1041.30	9.86	0.37	No drilling	-	-
Fairfield	1317.20	25.30	0.87	No drilling	-	-
Guernsey	1370.71	62.98	1.42	2007	2014	2017
Harrison	1069.47	62.48	2.03	2006	2014	2017
Holmes	1093.86	39.81	0.38	2006	2014	2017
Knox	1379.75	34.46	0.81	2006	2013	2017
Licking	1777.85	36.09	0.54	2006	2013	2017
Medina	1094.43	31.44	0.82	2006	2011	2017
Monroe	1186.73	77.08	0.40	2007	2014	2017
Morgan	1096.67	71.54	0.97	2007	2014	2017
Morrow	1042.35	25.14	0.45	No drilling	-	-
Muskingum	1745.62	57.65	1.41	2007	2014	2017
Noble	1042.99	69.68	1.23	2007	2014	2017
Perry 3	1062.06	59.92	0.69	2007	2013	2017
Portage	1308.73	38.33	3.00	No drilling	-	-
Richland	1293.71	32.97	0.86	2006	2013	2017
Stark	1504.99	24.72	0.98	2006	2012	2017
Summit	1088.90	29.58	1.72	2006	2011	2017
Tuscarawas	1485.49	53.17	1.12	2006	2014	2017
Washington	1662.45	66.72	1.01	2007	2014	2017
Wayne	1440.21	17.75	0.58	2006	2011	2017
Total	34,817.14	46.52	1.01	-	-	-

¹ Data are from the United States Census Bureau. ² Data are derived from National Land Cover Database (NLCD) 2016. ³ There is no drilling activities reported in Ashland and Perry, but several pipeline right-of-way (ROW) corridors go through them.

).

The majority of the Muskingum River Watershed is rural, and forest dominates the landcover (Figure 3b and Table 1). About 45% of the watershed area is covered by forest [30], comprising 38.50% deciduous forest, 0.69% evergreen forest, and 6.24% mixed forest. The main species of trees are deciduous, e.g., red oak, white oak, hickory, black cherry, yellow poplar, sugar maple, and red maple [37]. Following the seasonal changes in local climate, the annual forest growth in Ohio is divided into leaf-on seasons and leaf-off seasons. The considerable core forests deep inside the Ohio rolling mountains provide critical habitats for various wildlife species [38], for instance, songbirds, reptilians, and insects.

2.2. Datasets

The boom of these unconventional shale gas extraction activities in the Muskingum River Watershed started in 2010, reached a peak near 2014, and recently reached a trough in 2018. To trace the development pattern, the boom period was divided into two phases, the first rising phase and the latter declining phase. Three time nodes were set to capture these landcover conditions at the start and the end of the two phases: before the boom (before 2010), during the peak of the boom (2014), and the end of the boom (about 2018).

High-resolution remote sensing data are essential for the detection of landcover changes caused by the unconventional shale gas extraction activities on a large watershed scale. Field surveys indicated that the size of a well pad is approximately 150 m × 100 m, the width of a pipeline ROW corridor ranges from 20 to 50 m, and the width of an access road is approximately 3 to 5 m. Moderate-resolution satellite imagery, such as Landsat-7 and Landsat-8 imagery, is not adequate to precisely digitize and quantify the landcover changes induced by the unconventional shale gas extraction activities. New sources of satellite imagery, such as 10 m spatial resolution Sentinel-2 imagery, may be adequate for pipeline ROW corridors but still challenging for roads. Commercial sources of satellite imagery (e.g., SPOT-6/7, WorldView-2/3, and Planet imagery) provide adequate spatial resolution, but the cost for a large-scope investigation is still considerable. In this research, available public aerial images with the highest spatial resolution from the Ohio Statewide Imagery Program (OSIP) [39] and National Agriculture Imagery Program (NAIP) [40] for the entire Muskingum River Watershed were collected to digitize these shale gas extraction facilities and their impacts on landcover changes. The OSIP I dataset contains red/green/blue (RGB) natural-color aerial images acquired in either 2006 or 2007; since all of the unconventional shale gas extraction activities occurred after 2007 in the study area, these data were adequate to capture the initial landcover conditions before the boom. Most counties (14) had aerial images acquired in 2006 (Table 1). Eight counties did not have aerial images acquired in 2006 and, therefore, aerial images acquired in 2007 were used instead. The spatial resolution of these aerial images was 0.3 m (1 foot), which was the finest spatial resolution among all available public data was is sufficient to digitize and quantify the landcover changes induced by the unconventional shale gas extraction activities (Figure 2). To maintain the consistency of the high spatial resolution and geo-reference quality, the OSIP II dataset for the entire Muskingum River Watershed was selected to capture the landcover condition at the peak of the boom. The available aerial images from OSIP were acquired in 2014 for 12 counties, in 2013 for five counties, in 2012 for two counties, and in 2011 for three counties (Table 1). The spatial resolution of the OSIP II natural-color aerial images was also 0.3 m (1 foot). The RGB natural-color aerial images acquired in 2017 from the NAIP dataset were selected to capture the landcover condition at the end of the boom. NAIP 2017 aerial images were the most recent available public data for the entire Muskingum River Watershed with 1 m spatial resolution. Although coarser than OSIP data, NAIP 2017 aerial images were still adequate to precisely digitize the shale gas extraction activities during the declining phase and to capture the landcover conditions at the end of the boom.

The information of the shale gas well locations and their associated permission/completion dates, the extent of well pads, the centerline of access roads, and the horizontal pipe vector data were collected from the Ohio Department of Natural Resources (ODNR) Oil and Gas Well Locator [10]. The well pad data had an empty time attribute before 5 September 2017, and the horizontal pipe data did not have a time attribute. The LiDAR point-cloud topographical data acquired at the same time as aerial images before the boom were also collected from OSIP. The LiDAR data were used to construct the initial canopy height models (CHMs) before the boom of the unconventional shale gas extraction activities.

There are no relevant unconventional shale gas extraction activities in five counties in the west of the Muskingum River Watershed, namely, Athens, Crawford, Fairfield, Morrow, and Portage, which were excluded from the change detection analysis (Table 1).

3. Methods

The available geographic information system (GIS) data from ODNR were first used to analyze the temporal development and the structure of well pads, and then used to constrain the search area. The comprehensive well pad facilities and landcover changes in vector format were then digitized from the triple-temporal aerial images. The kernel density maps were estimated and used to uncover the spatial pattern of the unconventional shale gas extraction activities and deforestation. Object-oriented classification and statistical tests were employed in the hotspots of the densest activities to quantify the forest fragmentation conditions and to verify the occurrence of forest fragmentation. Lastly, the forest structural changes in 2D space (core forest areal loss) and 3D space (forest volume loss) were estimated through forest image objects and LiDAR-derived canopy height models (CHMs), respectively.

3.1. Spatial Footprints of Shale Gas Extraction Facilities

A group of spatial joins in ArcGIS 10.5 were used to estimate the time attributes of well pads since only the shale gas well locations had adequate time attributes. Spatial joins were first used to bond the well pads and the shale gas wells inside them, and then to bond the horizontal pipes and their adjacent shale gas wells. The earliest completion date of the shale gas wells within each well pad was assumed as the completion date of the well pad, while the completion date of a shale gas well was assumed as the completion date of its adjacent horizontal pipe. Other statistics were also derived from the spatial joins to analyze the structure of the shale gas extraction systems, including the well number in each well pad, the average and maximum well number in each well pad, and the average length of horizontal pipes in each well pad.

The comprehensive footprint GIS data of the unconventional shale gas extraction activities and landcover changes in vector format are critical for the analysis of forest structural changes (Figure 4), but available public GIS data are limited. Landcover change detection analysis and digitization using high-resolution aerial images represent the most effective and accurate approach since the direct changes are sensibly small. Because the study area is huge, the available well pad locations were first used as initial search sites to verify if there exists a well pad on the aerial image of each time node. For the detected and verified well pad locations, the landcover change detection analysis compared the later time node aerial image with the earlier time node aerial image. Through numerical analysis and on-screen interpretation, a group of well pad facilities were mapped out, e.g., the well pads, retention ponds, and access roads. The numerical attributes of each well pad, retention pond (e.g., centroid location and area), and access road (e.g., centerline, length, area, and width) were derived. From the verified well pads, the search space was expanded to identify and map the pipeline right-of-way (ROW) corridors between well pads and other facilities along the pipeline ROW corridors, e.g., the gathering and boosting stations and compressor stations. The centerlines and borders of pipeline ROW corridors were digitized when they went through forest regions. Outside the forest regions, straight-line segments were drawn to connect adjacent endpoints of pipeline ROW corridor centerlines as the centerlines of potential pipeline ROW corridors. The extent of the facilities along the pipeline ROW corridors and the centerlines and borders of their associated access roads were also mapped out, and their length and area attributes were derived. To visualize the spatial pattern of the spatial footprint from shale gas extraction facilities, a kernel density map was estimated from the well pad centroid locations using an identical weight by employing the ArcGIS 10.5 Kernel Density tool. The well pad number per km² was estimated using a 20 km searching radius.

Disturbance zones were also digitized from aerial images (Figure 4). Adjacent to the well pads and on the pipeline ROW corridors, irregular buffer zones were deforested or disturbed during the construction of the relevant facilities and then these zones quickly recovered to grassland or cropland. These buffer zones are called disturbance zones. The landcover change detection analysis also digitized the disturbance zones (location, area) immediately adjacent to well pads and on pipeline ROW corridors. Shown in Figure 5a–c, the disturbance zones adjacent to well pads were apparent and traceable. Because the well pads were always constructed in forest region or quarantined by a fence, the digitization traced the sharp edge of the deforestation or the fence, and determined the disturbance zones of well pads as the regions enclosed by these edges and beyond the well pad boundaries. In contrast, the disturbance zones of the pipeline ROW corridors were less traceable. A rare case of observable pipeline ROW corridor under construction is shown in Figure 5d–f. Pipeline ROW corridors left apparent footprints only when they went through forest where deforestation was required. However, when a pipeline ROW corridor went through a field between two forest patches, the disturbed land recovered so quickly and perfectly and, thus, was difficult to trace (Figure 5f). The disturbance zones of the pipeline ROW corridors in the forest regions were directly digitized from deforestation regions by comparing the image pairs on adjacent time nodes. Indicated by the footprints of pipeline ROW corridors in forest regions, a 30 m width buffer zone was created for each pipeline ROW corridor centerline outside the forest region as the potential disturbance zone.

The landcover change detection analysis further targeted direct deforestation caused by these unconventional shale gas extraction activities. Within the boundary of all the footprints, e.g., well pads, access roads, disturbance zones, and other facilities, the landcover change detection analysis mapped out the regions where the trees disappeared on the aerial images of the later time node. The centroid location and area attribute were derived from each deforestation region. Similarly, a kernel density map was estimated from deforestation region centroid locations while using the area attribute as the weight. The deforestation area (unit in ha) per km² was estimated using a 5 km searching radius.

The ecological serving capability was then used to simplify the complicated unconventional shale gas extraction system and to support the following ecological landscape analysis. The deforestation degraded several key ecosystem services, including the carbon storage in the format of biomass and in soils, as well as the regulation of water balance and river flow [25]. The landcovers were grouped into four categories: no ecological serving capability area, including impervious surfaces of the well pads, retention ponds, and the stations built on the concrete floor; low ecological serving capability area, including disturbance zones and access roads that are mainly grassland/cropland, bare soil, and gravel surface; unknown, mainly the under-construction regions that cause apparent deforestation but are not sure to cause how large impervious surface, grassland surface, or bare soil surface. In contrast, the forest area is considered as an area with high ecological serving capability and it is classified as the fourth category. The landcover conversion matrix was built given such a classification system to discuss the degradation of ecological serving capability in Muskingum River Watershed due to these unconventional shale gas extraction activities.

3.2. Forest Fragmentation Verification

Since the direct areal loss of forest was sensibly small, the precise verification of the occurrence of forest fragmentation was essential to justify the values to discuss the forest fragmentation issues. Two apparent hotspots were identified from the well density map; one covered Carroll–Harrison counties, and the other covered the Belmont–Guernsey–Monroe–Noble counties. An object-oriented approach through geographic object-based image analysis (GEOBIA) [41,42] was applied to two representative sites of the densest shale gas extraction activities (hydrologic units at the level of 12). Test Site 1 was the Dining Fork Catchment located in the Carroll–Harrison hotspot, and Test Site 2 was the Beaver Creek Catchment located in the Belmont–Guernsey–Monroe–Noble hotspot. Beyond the hotspots, the impacts of deforestation on forest structure were quite limited and, thus, are not discussed in this paper. Three landscape ecology metrics were used to indicate and quantify the forest fragmentation conditions on the image object level, including overall forest connectivity, anthropogenic fragmentation, and natural fragmentation [43,44]. Three landcover categories were mandatory to derive these three metrics, namely, forest (high ecological serving capability), non-forest natural landcover (low ecological serving capability), and anthropogenic landcover (no or degraded ecological serving capability). To increase the separability of image objects and to improve the classification accuracy, the non-forest natural landcover category was further divided into non-forest vegetation (e.g., grassland) and water body, while the anthropogenic landcover category was further divided into buildings and highways made of impervious surfaces. Meanwhile, regarding the similarity in ecological serving capability between grassland and cropland, the cropland was also included in the non-forest vegetation category, even though it is a production of human activities.

The initial landcover condition was mapped out using the first-time node aerial images, and the landcover changes were updated through the digitized shale gas extraction activity footprints. To form proper image objects, especially at the edge of the forest, and to obtain reliable topological information (adjacent landcover types), multisource data were avoided and only aerial images were used in the classification. For each test site, the multiresolution segmentation algorithm in eCognition software [45] was applied to the aerial images of the first-time node. The parameters of the segmentation were set to 150 for the scale, 0.3 for the shape, and 0.8 for the compactness. Then, 18 relevant feature variables were used to classify the image objects into five landcover categories. Those feature variables included nine spectral variables (the mean pixel value of the red/green/blue band within an image object, the pixel value standard deviation of the red/green/blue band, hue, saturation, and brightness of an image object), seven geometric variables (border length, length, width, border index, shape index, compactness, and roundness of an image object), and two texture variables (gray-level co-occurrence matrix (GLCM) homogeneity and dissimilarity in all directions within an image object) [42]. A machine-learning algorithm known as support vector machine (SVM) [46] was trained in eCognition software to perform the classification and to produce the landcover maps. A total of 2067 image objects were randomly interpreted as a training sample set to train the SVM classifier, and the other group of 3124 independent image objects were used as a test sample set. Because of the different occlusion of tree canopies in leaf-on and leaf-off seasons, the misclassification of forested coverage on the aerial imagery of adjacent time nodes caused unavoidable false-positive forest landcover changes. The error area even surpassed the true landcover changes, whereby the low signal-to-noise ratio (change percentage to error percentage ratio) indicated that the direct comparison of two classification maps on adjacent time nodes were not suitable to detect the true changes of forest and other landcovers. Regarding the second and the third time nodes, the dissolved union of all polygonal shale gas extraction activity footprints was used as an input thematic layer to the eCognition software. Using the assign class algorithm in the eCognition software, these image objects out of the thematic layer polygons were assigned as their initial landcover categories on the basis of first-time node classification maps, while these image objects within the thematic layer polygons were assigned as a new category called “shale gas extraction activities” (SGEAs), which was assumed to have a degraded ecological serving capability.

Three metrics were defined on an object level to quantify the forest fragmentation conditions on the basis of a feature variable called the relative border to nearby class [47]. The formula to calculate the feature variable is shown below.

$$Rel.BDto\#C=\frac{theborderlengthsharedby\#C}{totalborderlength},$$

(1)

where #C is a landcover category on the classification map, and Rel.BD to #C is the abbreviation of relative border length to landcover category C. The habitat transition was assumed to be equal in terms of connectivity, which means that the transition from forest to non-forest vegetation was the same as the transition from the forest to the road, etc. Accordingly, the three metrics, overall forest connectivity ($P_{ff}$), anthropogenic fragmentation ($P_{fa}$), and natural fragmentation ($P_{fn}$) are defined below.

$$P_{ff}= Rel. BD to Forest,$$

(2)

$$P_{fa}= Rel. BD to Highway+Rel. BD to Building+Rel. BD to SGEA,$$

(3)

$$P_{fn}= Rel. BD to Waterbody+Rel. BD to Nonforest Vegetation.$$

(4)

Statistical tests were used to determine the significant changes in forest fragmentation conditions. These three landscape ecology metrics were derived from each image object at the three time nodes, before and after the update of landcover changes near the peak and at the end of the boom. The mean values and standard deviations were estimated from all the forest image objects in the test sites before and after the landcover changes. A two-sample test of the difference with 95% confidence (α = 0.05) was used to determine the occurrence of forest fragmentation. Two tests were designed in each site; the first one tested the difference in forest fragmentation conditions near the peak of the boom compared to that before the boom, and the other one tested the conditions at the end of the boom compared to that before the boom. If the test results were significantly different from 0, a decision was determined that the forest fragmentation conditions changed and forest fragmentation occurred; otherwise, if the test results were not significantly different from 0, a decision was determined that there were no changes even though the numbers varied.

3.3. Loss of Core Forest Area and Forest Volume

Wild species are sensitive to the 2D forest structure since core forest is more inhabitable than the forest edge. A pipeline ROW corridor that goes through a large forest may cause only a small portion of deforestation but may tremendously change the 2D structure of the forest and destroy the inhabitable core forest. To further evaluate the 2D forest structural change pattern and quantify the loss of core forest area, an object-oriented forest structural analysis was applied to the forest landcover in the test sites. In this analysis, the structural components of the forest were defined as core, edge, perforated, and patch [43,48]. The forest image objects with more than 80% of their border connecting to non-forest were defined as patch forests, the forest image objects that were 100 meters away from any other landcover types were defined as core forests, the forest image objects within 100 m of outside non-forest landcover types were defined as forest edges, and the forest image objects within 100 m of inside non-forest landcover types were defined as perforated forests. These parameters were determined because of their extensive use in relevant habitat fragmentation studies [8]. After the classification of forest structure, the area statistics and percentage changes of each forest component were derived to analyze the change pattern of the forest structure, especially the areal loss of core forest.

Another impact of deforestation on ecological service is the loss of forest biomass and the net release of CO₂ from forest to the atmosphere. Furthermore, the deforestation volume is the fundamental quantity to estimate the loss of forest biomass. LiDAR point-cloud data were employed to evaluate the 3D forest structural change pattern and to quantify the loss of forest volume caused by these unconventional shale gas extraction activities. All LiDAR points were used to interpolate the digital surface model (DSM), and only the ground points were selected to interpolate the digital terrain model (DTM) [49]. A triangulation interpolation method with a linear option in ArcGIS 10.5 was employed. The DSM and DTM were interpolated with a cell size of 0.762 m (2.5 feet), to be consistent with OSIP DTM data. The raw canopy height model (CHM) was constructed from the DSM by subtracting the DTM [50] and was resampled to a 0.3 m (1 foot) cell size through a bilinear resampling method, to be consistent with OSIP imagery data.

Although the LiDAR system can acquire dense point-cloud data to reconstruct the 3D structure of the forest, the laser likely hits the “shoulder” of the trees rather than the exact treetops [50]. Thus, the tree height and the forest profile from the reconstruction of LiDAR point-cloud data are indeed underestimated. According to the in situ measurements from field surveys, a linear model was developed through ordinary least squares (OLS) fitting to remove the bias between the tree height in the raw CHM and in situ observations. The linear model was then applied to correct the raw underestimated CHMs. The spatial differential of deforestation volume is the product of the deforestation area size and the relevant corrected CHM height in each terrain cell, while the total deforestation volume is the integral of the spatial differential within the 2D spatial domain.

4. Results

4.1. Footprints of Shale Gas Extraction Facilities

The boom of the unconventional shale gas extraction activities was spatially concentrated in the east of the Muskingum River Watershed (Figure 6) and its temporal development showed a clear unimodal pattern. The first horizontal pipe was drilled in 1991 in the Muskingum River Watershed, while the first well pad was constructed in 2010. The relevant unconventional shale gas extraction activities then showed gradual development (Figure 7a). These activities reached their peak in 2014 and then slowed down until 2018. Given the well location data and pad extent data from ODNR, by the end of 2018, 1481 shale gas wells were permitted and 1189 shale gas wells were constructed in 356 well pads. In total, 340 well pads were verified on the aerial images during the boom period from 2006 to 2017; 222 of them were constructed during the raising phase, while the rest were constructed during the declining phase. During the declining phase, about 21% of the new wells are drilled in old well pads rather than new ones, which further contributed to the decrease in well pad construction. The verified pad density was approximately 1.06 pads per 100 km² near the peak of the boom and 1.63 at the end of the boom. The average and the maximum number of shale gas wells per pad kept increasing until 2014 and 2015, respectively (Figure 7b), while 87.63% of the well pads comprised no more than six shale gas wells (Figure 7c). Moreover, with the development of horizontal drilling technology, the average length of a horizontal pipe continued to increase after 2000 (Figure 7d).

The landcover change detection analysis discovered an uneven temporal pattern of these unconventional shale gas extraction activities. In total, 291.00 ha of well pad impervious surface was constructed in the Muskingum River Watershed near the peak of the boom, while the area expanded to 435.07 ha at the end of the boom. The total length of the access roads to all relevant facilities was 156.94 km near the peak and 249.25 km at the end; the area of these paved access roads was 117.73 ha near the peak and 194.64 ha at the end. These facilities disturbed 927.37 ha of area that was immediately adjacent to the well pads near the peak of the boom, while the area expanded to 1356.04 ha at the end of the boom. The average well pad area during the raising phase and the declining phase was 1.31 ha and 1.28 ha, respectively. The average disturbed area of these well pad facilities during the raising phase and the declining phase was 4.18 ha and 3.99 ha, respectively. Approximately two-thirds of all constructed well pad facilities and all relevant landcover changes were completed during the raising phase, whereas one-third occurred during the declining phase. This pattern indicates a quick raising of these unconventional shale gas extraction activities, as well as a quick decay after the peak. The average well pad area and average disturbed area were close during these two phases, indicating a consistent construction pattern throughout the boom period, and the reduced landcover changes during the declining phase of the boom were simply caused by the decay of shale gas extraction activities. In addition to these well pad facilities, 324.01 km pipeline ROW corridors were found before the boom, and their length increased to 1235.42 km near the peak of the boom and 2084.90 km by the end of the boom. The disturbed area of these pipeline ROW corridors, either through the forest regions or out of forest regions, was about 2623.01 ha during the raising phase and 2467.32 ha during the declining phase.

The disaggregated statistics for each county within the Muskingum River Watershed showed uneven spatial patterns of these unconventional shale gas extraction activities (Figure 8). These shale gas extraction activities were concentrated in six counties, and all were located in the east edge of the watershed. Near the peak of the boom, the top six counties constituted 89.19% of the well pads, 90.46% of the immediately adjacent disturbance zones to well pads, and 59.43% of the pipeline ROW corridors. At the end of the boom, these numbers slightly changed to 90.29%, 91.25%, and 58.98%, respectively. There were two apparent hotspots clustered by the six most concentrated counties; one involved the Carroll–Harrison counties, and the other involved the Belmont–Guernsey–Monroe–Noble counties (Figure 9). The pattern was the same for well pad density and deforestation, and it was consistent during both raising and declining phases except for the increase in density. The stretch of deforestation to the west of the watershed was mainly caused by the pipeline ROW corridors, while the outstretched deforestation was less severe than that within the hotspot counties.

The spatial pattern of deforestation was consistent with the spatial pattern of these unconventional shale gas extraction activities (Figure 10). Correspondingly, the dense shale gas extraction activities in the Muskingum River Watershed totally caused 2201.61 ha of deforestation area near the peak of the boom, while this number was 3626.80 ha at the end of the boom. Moreover, 80.09% of the deforestation area occurred in the six most concentrated counties near the peak of the boom, while this number was 78.89% at the end, which is consistent with the unconventional shale gas extraction activities. Additionally, 52.93% of the deforestation area was caused by the pipeline ROW corridors near the peak of the boom, while this percentage was 58.68% at the end of the boom. The remainder of the deforestation was caused by the well pad facilities and their immediately adjacent disturbance zones. The landcover change detection analysis was only aimed at the deforestation caused by unconventional extraction activities and, actually, little deforestation was observed due to other activities in this rural region.

More than half of these landcover changes showed a degradation of ecological serving capability (

Table 2. The landcover conversion matrix showing the degradation of land ecological serving capability near the peak and at the end of the boom in the Muskingum River Watershed due to unconventional shale gas extraction activities.

	Ecological Serving Capability Near the Peak (ha)				Ecological Serving Capability at the End (ha)
	Low	None	Unknown	Total	Low	None	Unknown	Total
High capability	1891.44	100.97	209.20	2201.61	3450.74	156.62	19.44	3626.80
Low capability	1517.38	354.39	0.00	1871.77	3031.37	526.85	0.00	3558.22
Total	3408.82	455.36	209.20	4073.38	6482.11	683.47	19.44	7185.02

). Near the peak of the boom, in total, 4073.38 ha of the land surface was affected by the unconventional shale gas extraction activities, which represents about 0.20% of the entire Muskingum River Watershed. Within these affected regions, 37.25% remained with low ecological serving capability after the disturbance, whereas 57.61% degraded to a lower level of ecological serving capability. At the end of the boom period, the affected land surface expanded to 7185.02 ha, which is about 0.34% of the entire Muskingum River Watershed. The area which maintained low ecological serving capability and the degraded regions represented 42.19% and 57.54%, respectively. The remainder constituted under-construction area whose ultimate landcover conversion was unknown. The unknown area was far less prevalent at the end of the boom than near the peak, and this reduction was simply because of fewer new shale gas extraction activities in the declining phase of the boom. Only a small part of the land surface was affected in the Muskingum River Watershed due to unconventional shale gas extraction activities, while the majority of these affected regions suffered a degradation of ecological serving capability.

4.2. The Occurrence of Forest Fragmentation

The classification strategy developed in this research produced classification maps with sufficient accuracy. An SVM classifier in eCognition software was trained using a stratified sample set containing 2067 interpreted image objects, and the accuracy of the SVM classifier was validated through an absolute independent stratified sample set containing 3124 interpreted image objects. The confusion matrix (

Table 3. Confusion matrix and the classification accuracy of the support vector machine (SVM) classifier trained for this research. The assessment is conducted on the image object level.

		Interpreted Landcover Type
		Forest	NFV 1	Water Body	Building	Highway	Total
Classification	Forest	1525	29	1	7	0	1562
	NFV 1	3	780	0	0	1	784
	Water body	0	42	177	1	0	220
	Building	0	3	0	219	1	223
	Highway	0	9	0	46	280	335
Total		1528	863	178	273	282	3124
Producer accuracy		99.8%	90.4%	99.4%	80.2%	99.3%	-
User accuracy		97.6%	99.5%	80.5%	98.2%	83.6%	-
KIA 2 per class		0.996	0.872	0.994	0.787	0.992	-
Overall Accuracy		95.4%					-
KIA 2		0.931					-

¹ NFV: non-forest vegetation. ² KIA: Kappa index of agreement.

) indicated the high classification accuracy of forest, whose producer accuracy was 99.8%, user accuracy was 97.6%, and Kappa index of agreement (KIA) per class was 0.996. Two sorts of misclassification can be observed; one was between non-forest vegetation and water bodies, and the other was between buildings and highways. However, when combining non-forest vegetation and water body categories, the non-forest natural landcover category had an aggregated producer accuracy of 96.0% and an aggregated user accuracy of 99.5%. When combining building and highway categories, the anthropogenic landcover category had an aggregated producer accuracy of 98.4% and an aggregated user accuracy of 97.8%. Since building and highway image objects served together to derive anthropogenic fragmentation ($P_{fa}$, Equation (3)), while non-forest vegetation and water body image objects served together to derive natural fragmentation ($P_{fn}$, Equation (4)), the aggregated accuracy still indicates reliable estimation of these critical landscape ecology metrics.

The statistical tests on forest fragmentation conditions were further conducted in two test sites within the two apparent hotspots (Figure 9a). The initial and updated object-oriented classification maps were employed to derive critical metrics. The overall forest connectivity ($P_{ff}$), anthropogenic fragmentation ($P_{fa}$), and natural fragmentation ($P_{fn}$) were derived through Equations (2)–(4) for each forest image object. Mean values and standard deviations of these metrics were estimated for the two-sample test of difference.

The statistical tests indicated different forest fragmentation processes in these two hotspots (

Table 4. The changes in forest fragmentation conditions in Carroll–Harrison counties.

	Before the Boom		Near the Peak		Z of the Difference 1	p-Value 2	Significance
	Mean	SD	Mean	SD
$P_{ff}$	0.799917	0.267146	0.794027	0.268770	1.742	0.041	Significant
$P_{fa}$	0.045157	0.119444	0.054618	0.129225	5.172	<0.001	Significant
$P_{fn}$	0.154926	0.230302	0.151355	0.227833	1.180	0.119	Insignificant
			At the End		Z of the Difference 3	p-Value 2	Significance
			Mean	SD
$P_{ff}$			0.792344	0.268851	2.235	0.013	Significant
$P_{fa}$			0.056763	0.131100	6.287	<0.001	Significant
$P_{fn}$			0.150893	0.227374	1.332	0.091	Insignificant

¹ Sample size: N_before = 28,549, N_peak = 28,081. ² α = 0.05. ³ Sample size: N_end = 28,016.

and

Table 5. The changes in forest fragmentation conditions in Belmont–Guernsey–Monroe–Noble counties.

	Before the Boom		Near tde Peak		Z of the Difference 1	p-Value 2	Significance
	Mean	SD	Mean	SD
$P_{ff}$	0.830812	0.255493	0.827336	0.255726	1.218	0.112	Insignificant
$P_{fa}$	0.053855	0.135275	0.059433	0.139807	3.184	0.001	Significant
$P_{fn}$	0.115333	0.208435	0.113231	0.206689	0.872	0.192	Insignificant
			At the End		Z of the Difference 3	p-Value 2	Significance
			Mean	SD
$P_{ff}$			0.825775	0.255953	1.757	0.039	Significant
$P_{fa}$			0.063294	0.143385	5.289	<0.001	Significant
$P_{fn}$			0.110931	0.204968	1.828	0.034	Significant

¹ Sample size: N_before = 34,962, N_peak = 34,652. ² α = 0.05. ³ Sample size: N_end = 34,264.

). Near the peak of the boom, the overall forest connectivity in the test site of Carroll–Harrison counties significantly decreased, while the overall forest connectivity in the test site of Belmont–Guernsey–Monroe–Noble counties also decreased, but the change was not significantly different from 0. Furthermore, at the end of the boom, the decrease in overall forest connectivity at both sites was statistically significant. These results are consistent with the denser hotspot in Carroll–Harrison counties and the less dense hotspot in Belmont–Guernsey–Monroe–Noble counties. The unconventional shale gas extraction activities first quickly caused forest fragmentation in Carroll–Harrison counties near the peak of the boom and advanced the forest fragmentation level until the end of the boom. In contrast, the forest fragmentation conditions slowly changed in Belmont–Guernsey–Monroe–Noble counties, and significant forest fragmentation occurred only toward the end of the boom, representing a delayed response compared to the Carroll–Harrison counties. The anthropogenic fragmentation level significantly increased in both hotspots and at the end of both phases. The significant increases were apparent because the landcover changes were caused by the anthropogenic activities of the recent boom of the unconventional shale gas extraction, and these changed landcovers adjacent to forest showed a sensible degradation of ecological serving capability. The change in natural fragmentation was considered insignificant because natural fragmentation is less likely to occur in such a short period and no relevant natural fragmentation procedure was included in the analysis. However, the natural fragmentation significantly changed in Belmont–Guernsey–Monroe–Noble counties at the end of the boom, indicating that the forest structure changed severely, whereby the contribution of the natural fragmentation to the overall forest fragmentation conditions also significantly changed.

The test results also suggested lower value in discussing the forest fragmentation issues beyond the two hotspots of densest unconventional shale gas extraction activities. The overall forest connectivity insignificantly decreased in Belmont–Guernsey–Monroe–Noble counties near the peak of the boom, highlighting the forest fragmentation issues as negligible in this region near the peak of the boom since the change was not significantly different from 0. The well pad density, in this case, was not denser than 9.0 pad per km². Using the density of 9.0 pad per km² as a threshold, the counties beyond the hotspots either near the peak or at the end of the boom had a well pad density lower than the threshold. These counties had a lower probability for the occurrence of forest fragmentation and, hence, it may be meaningless to discuss the issues of forest fragmentation. Thus, the discussion of these counties was excluded in the sections below.

4.3. Loss of Core Forest and Forest Volume

The changes in 2D forest structure components (core, edge, perforated, and patch) were much more severe than the direct areal loss of forest (

Table 6. Forest structural changes at two representative test sites.

	Before the Boom (ha)		Near the Peak (ha)		Change (ha)		Percentage
	Site 1 1	Site 2 2	Site 1	Site 2	Site 1	Site 2	Site 1	Site 2
Core	535.86	909.65	482.84	859.53	−53.02	−50.12	−9.89%	−5.51%
Edge	2836.73	2788.04	2841.62	2851.49	4.88	63.45	0.17%	2.28%
Perforated	450.63	721.68	404.83	647.03	−45.80	−74.65	−10.16%	−10.34%
Patch	108.65	85.51	112.46	85.36	3.81	−0.15	3.51%	−0.18%
Total	3931.88	4504.89	3841.75	4443.42	−90.12	−61.47	−2.29%	−1.36%
			At the End (ha)		Change (ha)		Percentage
			Site 1	Site 2	Site 1	Site 2	Site 1	Site 2
Core			457.61	830.55	−78.25	−79.10	−14.60%	−8.70%
Edge			2859.28	2841.43	22.55	53.39	0.79%	1.91%
Perforated			391.42	624.49	−59.21	−97.19	−13.14%	−13.47%
Patch			112.84	85.04	4.19	−0.47	3.86%	−0.56%
Total			3821.15	4381.51	−110.73	−123.38	−2.82%	−2.74%

¹ Site 1 was the Dining Fork Catchment in Carroll–Harrison counties. ² Site 2 was the Beaver Creek Catchment in Belmont–Guernsey–Monroe–Noble counties. Test sites are extended to include the complete forest patches on the edge of the catchments.

). Near the peak of the boom, the total direct areal loss of forest was approximately 2.29% in the test site of Carroll–Harrison counties and 1.36% in the test site of Belmont–Guernsey–Monroe–Noble counties. The areal loss of core forest was approximately 9.89% in Carroll–Harrison counties and 5.51% in Belmont–Guernsey–Monroe–Noble counties. At the end of the boom, these numbers increased to 2.82%, 2.74%, 14.60%, and 8.70%, respectively.

The direct areal loss of forest represented only a small part of the entire forest regions in the two hotspots of the most concentrated unconventional shale gas extraction activities, but the areal loss was quadruple to quintuple for the core forest that served as key habitats for local ecological communities. These lost core forests were caused by the facilities of these unconventional shale gas extraction activities; together with their inner perforated forests, they were converted to outside forest edges, resulting in an increase in forest edges. About two-thirds of the core forest areal loss occurred during the rising phase of the boom, while one-third occurred during the declining phase. This temporal pattern is consistent with the development of the unconventional shale gas extraction activities. Moreover, the total areal loss of forest increased only slightly during the declining phase in Carroll–Harrison counties, while the total areal loss of forest doubled during the declining phase in Belmont–Guernsey–Monroe–Noble counties. Considering the larger areal loss of forest and quicker speed, the forest structure showed a more drastic change in Belmont–Guernsey–Monroe–Noble counties during the declining phase than during its previous raising phases. The change during the declining phase was also more severe than that in Carroll–Harrison counties during both phases. The results are consistent with the statistical tests of forest fragmentation conditions in Section 4.2, where the changes in the three landscape ecology metrics in Belmont–Guernsey–Monroe–Noble counties were all significant at the end of the boom.

The raw CHM was separately derived from LiDAR point-cloud data in each county. A linear regression model through OLS fitting was developed to correct the raw CHMs through 44 in situ measurements of tree height collected in Carroll County (Figure 11a). The plot shows an apparent underestimation of the tree heights from LiDAR-derived raw CHM, while the R² of 0.9204 of the model indicates an overall good fit to correct the heights.

The deforestation volume was also spatially uneven, and the majority was concentrated in the two hotspots. The total volume of deforestation was 1.09 × 10¹⁰ cubic feet near the peak of the boom, while this number was 1.81 × 10¹⁰ cubic feet at the end of the boom. The top counties with the largest volume of deforestation (Figure 12) were the same as those with the largest area of deforestation (Figure 10) and were ranked exactly the same. The average tree height of deforested trees was 45.85 feet near the peak and 46.20 feet at the end of the boom. The average tree heights in each county and for the entire Muskingum River Watershed only slightly varied near the peak and at the end of the boom (Figure 11b), indicating a consistent deforestation pattern during the two phases of the boom. The average tree heights of deforestation varied from county to county, ranging from 38.00 feet to 65.30 feet (Figure 11b), indicating the different forest targets for the deforestation in different counties. In particular, regarding the two hotspots of unconventional shale gas extraction activities, the deforestation in Carroll–Harrison counties targeted higher forest, with a height of about 50 feet on average. In contrast, the deforestation in Belmont–Guernsey–Monroe–Noble counties targeted lower forest, with a height of about 42 feet. Other counties are not discussed here because of their overall minor contribution to the deforestation in the Muskingum River Watershed.

5. Discussion

5.1. The Reliability of Analysis

Limited biases were included in the three time nodes when capturing the recent boom of the unconventional shale gas extraction activities in the Muskingum River Watershed. Since all well pads were constructed after 2010, the available OSIP high-resolution aerial images acquired by 2006/2007 were sufficient to capture the initial landcover condition before the boom. On the other hand, the most recent available NAIP aerial images obtained by 2017 were used to capture the landcover condition at the end of the boom. Although not exactly at the trough of the boom, these images are close enough and a few of the lost well pads (16 well pads in 2018 and fewer later) had trivial impacts on analyzing the pattern of these unconventional shale gas extraction activities. The peak of these unconventional shale gas extraction activities occurred in 2014; restricted by the available OSIP data, not all the counties captured the landcover conditions exactly at the peak of the boom. However, the top six counties constituted about 90% of well pads and disturbance zones, about 80% of the newly constructed pipeline ROW corridors, and about 80% of the deforestation during the raising phase. The aerial images of these six counties were obtained by 2014, which was exactly the peak of the boom. Thus, the bias of analysis was sensibly negligible near the peak of the boom in this research.

The current available highest-resolution remote sensing data for the entire Muskingum River Watershed enabled the accurate quantification of the spatial footprints of these unconventional shale gas extraction activities and landcover changes. The spatial resolution of the aerial imagery used in this research was 0.3 m for the time node before the boom and near the peak of the boom, and 1.0 m for the time node at the end of the boom. The high resolution contributed to the accurate digitization of the footprints of the unconventional shale gas extraction activities. The well pad extent data from ODNR were used as a reference to check the areal accuracy and the location accuracy of the digitized footprints. The statistical analysis of 34 well pads indicated that the areal accuracy of the digitized well pads was 96.49%. The average location displacements of well pad centroids were 3.22 m along the X-axis and 0.91 m along the Y-axis, and the bias was within the spatial accuracy of the aerial imagery (±5 m). These high-quality data accurately and reliably quantified the spatial distribution of these unconventional shale gas extraction activities. Moreover, the 0.762 m (2.5 feet) spatial resolution CHMs derived from dense OSIP LiDAR data also contributed to the accurate estimation of deforestation volume and tree heights. These results helped to elaborately understand the spatiotemporal development pattern of the unconventional shale gas extraction activities and their impacts on landcover, especially in forest regions.

Object-oriented classification maps were used to test the occurrence of forest fragmentation in the hotspots of the densest unconventional shale gas extraction activities. Given the test set of more than 3000 interpreted image objects, the user accuracy of forest was 97.6%, the user accuracy of non-forest natural landcover was 99.5%, and the user accuracy of anthropogenic landcover was 97.8%. Compared to other generally used landcover products from Landsat data [51,52], the accuracy of the classification maps in this research was sufficient for the forest fragmentation analysis. On the other hand, the commission error (1—user accuracy) of forest in terms of area was about 2.5%, which is comparable to the maximum forest areal change (about 2.8%) at the end of the boom. The ratio of areal change percentage to error percentage (signal-to-noise ratio) was about 1. This excessively low ratio indicates that it would not be proper to directly compare the classification maps to detect the forest landcover change because the detected changes were probably dominated by the classification errors rather than the true changes. The strategy used in this research, whereby the initial classification maps were updated through digitized footprints, avoided such issues, thus promising accurate maps of forest changes.

The forest image objects also enabled the accurate quantification of forest fragmentation conditions. Conventional pixel-based forest fragmentation analysis through overall forest connectivity, anthropogenic fragmentation, and natural fragmentation is generally based on moderate-resolution landcover and land use maps derived from Landsat TM data [43,44]. Pixel-based analysis through a moderate-resolution map overestimates the forest fragmentation conditions because the coarse representation of forest landcover destroys the critical long narrow forest components called “bridges” [53], correspondingly creating more tiny individual isolated forest patches. The following analysis is then “unfavorably inadequate” for the in-demand accurate evaluation of the habitat conservation status [53]. On the contrary, because of the redundant spatial information and spatial connectivity, the analysis causes a dramatic underestimation of forest fragmentation conditions if the pixel-based approach is directly applied to the extreme high-resolution aerial imagery [42]. The applied approach in this research integrated high-resolution aerial imagery and object-oriented image analysis to address these issues, providing two advantages. On one hand, image objects adequately maintained the border of forest on the high-resolution aerial imagery; on the other hand, image objects reduced the redundant connectivity of extremely high-resolution pixels through clustering the adjacent homogeneous pixels [54].

5.2. Landcover Change Pattern and the Degradation of Ecological Serving Capability

Similar to relevant studies, a fairly small magnitude of direct landcover changes and direct areal loss of forest were found due to these unconventional shale gas extraction activities, while the indirect impacts, e.g., the forest fragmentation and the areal loss of core forest, were much greater [8,12,55,56]. The majority of these landcover changes and disturbed zones suffered from the degradation of ecological serving capability. The volume of deforestation and the carbon storage inside them caused the net release of carbon dioxide to the atmosphere, which is a stimulator of the greenhouse effect of the Earth [57]. Moreover, without the interception of tree canopies, surface runoff has a high probability of increasing within these dense deforestation regions [58]. The extra surface runoff increases the risks of flooding in local regions and elevates the streamflow in local rivers, especially the Muskingum River to the Ohio River. Excessive suspended sediment may also occur due to the stronger surface runoff, consequently resulting in a higher risk to surface water quality [59]. Although overall small regarding the entire Muskingum River Watershed, the landcover changes may still cause significant degradation of ecological serving capability and higher risks of flooding and contaminated surface water in these six most concentrated counties, namely, Belmont, Carroll, Guernsey, Harrison, Monroe, and Noble. The densest deforestation in these six counties reached up to 2.8%, which calls for more attention.

5.3. Forest Fragmentation Issues

These shale gas extraction facilities, especially the pipeline ROW corridors, break the forest canopy and convert large continuous forest zones into small isolated forest zones. The deforestation caused by these pipeline ROW corridors leads to a degradation of ecosystem serving capability from forest to grassland or shrub. Such degradation is difficult to recover from with respect to the safety of the buried gathering pipelines, and these corridors result in irreversible changes in microclimate conditions [60,61]. Moreover, the deforestation regions that go through the deep core forests act as isolation belts for the natural species, e.g., area-sensitive or forest-interior songbirds and amphibians [24], and these isolation belts can fragment wildlife populations [62] and alter the movement of wildlife [63]. The unconventional shale gas extraction activities and corresponding deforestation in the Muskingum River Watershed were mainly concentrated in two hotspots; one was in Carroll–Harrison counties and the other one was in Belmont–Guernsey–Monroe–Noble counties. The change in forest fragmentation conditions was quicker and worse in Carroll–Harrison than that in Belmont–Guernsey–Monroe–Noble counties, and both hot-spots exhibited significant forest fragmentation at the end of the boom. Similar ecological issues may already exist in these regions. The statistical tests in this research indicate that the forest fragmentation may have occurred in the regions with a well pad density larger than 9.0 pad per km². Such a value could be used as a threshold to search a broader region for forest fragmentation issues caused by these unconventional shale gas extraction activities.

The loss of core forest and the forest fragmentation had strong impacts on area-sensitive species and ecosystem service in central Appalachian, and it is of particular importance to limit them [8,64,65]. Regarding the recent boom of these unconventional shale gas extraction activities in the Muskingum River Watershed, fewer well pad facilities were constructed during the declining phase of the boom but the construction of pipeline ROW corridors remained active throughout the boom period. Due to their strong impacts on the local ecological communities, continuous attention to and ecological investigations of these existing and ongoing pipeline ROW corridors are necessary.

6. Conclusions

This research investigated the most recent boom of unconventional shale gas extraction activities and their impacts on the landcover changes and forest structural changes in the Muskingum River Watershed in Appalachian Ohio. High-resolution aerial images were used to precisely capture the landcover conditions of the entire Muskingum River Watershed before the boom, near the peak of the boom, and at the end of the boom. Geographic object-based image analysis (GEOBIA) was applied to high-resolution aerial images to accurately quantify the forest fragmentation conditions, to verify the occurrence of forest fragmentation, and to derive the 2D forest structural changes. High-resolution canopy height models (CHMs) derived from dense LiDAR point-cloud data were used to estimate the 3D forest structural changes.

The results indicated decrease in construction activities of well pad facilities but the same level of active construction of pipeline ROW corridors before and after the peak of the boom. About two-thirds of the well pad facilities and half of the pipeline ROW corridors were constructed during the raising phase of the boom, while the remaining construction occurred during the declining phase of the boom. About 90% of the well pad facilities and their immediately adjacent disturbance zones, about 60% of the pipeline ROW corridors, and about 80% of deforestation were concentrated in six counties. These six counties clustered into two hotspots; one was in the Carroll–Harrison counties and the other was in the Belmont–Guernsey–Monroe–Noble counties.

These unconventional shale gas extraction activities caused limited landcover changes and areal loss of forest. In these changed and disturbed areas, the majority suffered from a degradation of ecological serving capability. These unconventional shale gas extraction activities had a stronger impact on forest structural changes. Forest fragmentation occurred in both hotspots at the end of the boom, and the tree targets of deforestation varied in space.

The comprehensive footprint GIS data of the unconventional shale gas extraction activities are essential for subsequent investigations on local ecological communities in the Muskingum River Watershed. The GEOBIA approach applied in this research is transferable and applicable to habitat fragmentation investigations in the broader Appalachian Ohio regions.