Spatiotemporal wave propagation of the shale oil and gas development in the Marcellus Shale in the past one century and a half

ABSTRACT Dense shale oil and gas extraction activities have a long history in the Marcellus Shale and the investigation of their evolution draws increasing interest. The spatiotemporal analysis of these shale oil and gas extraction activities is essential to understand their evolution in space-time, and important for the estimation of their historical impacts on local environments and communities. These shale oil and gas extraction activities formed numerous scattered hot spots, and it demands tremendous efforts to trace them over time. This particular issue then challenges the conventional spatiotemporal analysis approaches which are effective for single or limited spatiotemporal objects. This study models the spatiotemporal developments of the shale oil and gas extraction activities as wave propagation. By using the wave representation, the complicated transmission of these shale oil and gas extraction activities is automatically traced within a sub-space-time as free propagation. To implement the proposed wave model, density surfaces of active wells from 1868 to 2018 were estimated annually. A 15-year focal statistic filter was applied to the density curve at each location to identify temporal objects. The waves were formed from temporal objects using a spatiotemporal expansion algorithm while maximally allowing 1 gap year. Finally, comprehensive attributes were derived for each wave. The results indicate that chaotic and fragmentary waves propagated at the western and northern edges of the Marcellus Shale. In contrast, large and simple wave propagation occurred in the middle of the Marcellus Shale. The propagation of large waves has a southward general trend, the trend follows the Appalachian region’s mountainous topography. Different types of drilling activities show diversified wave propagation patterns and the conventional gas wells formed the most chaotic wave propagation. The wave representation is successfully applied to the evolution of the shale oil and gas extraction activities, and this study also develops a comprehensive set of attributes to depict the wave. The scattered hot spots of dense shale oil and gas extraction activities are automatically assembled in waves. These concrete and detailed attributes of the waves are effective in summarizing the general spatiotemporal evolution pattern of chaotic and fragmentary spatiotemporal objects, e.g., the dense shale oil and gas extraction activities in this case. The developed method is also transferable and applicable to the evolution of other geographic phenomena, which may comprise of either limited or numerous spatiotemporal objects.


Introduction
The drilling activities for shale oil and gas have a long history of about two centuries in the Marcellus Shale in the northeast of the United States (Harper 2008;Kargbo, Wilhelm, and Campbell 2010). These shale oil and gas extraction activities manipulate three generations of facilities. The first generation uses pump-jacks, which periodically pumps up oil from the shale along with tremendous noise and the leakage of gaseous and liquid fossil masses. The second generation starts to use Class II injection wells (Clark, Bonura, and Van Voorhees 2005). The injection wells improve production efficiency through hydraulic fracturing. They are more environmental-friendly by encapsulating the extraction in a closed system. The recent booming third generation is nowadays referred to as the unconventional wells (Vengosh et al. 2014), they further dramatically enhance production efficiency by incorporating hydraulic fracturing, well pad, and horizontal drilling. Great efforts have been devoted to descriptively investigate the shale oil and gas extraction activities in the Marcellus Shale and their impacts on local environments and human health, but the investigation on the spatiotemporal developments of these activities is still in urgent demand.
The shale oil and gas extraction activities have complicated impacts on local environments and human health on multiple spatial scales (Table 1). Evidence shows that the environmental and health impacts of shale oil and gas extraction activities range from local scale (about 1-mile distance) to regional scale (about several hundreds of km). Small-scale studies provide evidence of these environmental impacts and health impacts within a distance of several kilometers. Payne et al. have provided evidence that the residents within 1.6 km (1 mile) of a natural gas compressor in the unconventional shale gas extraction system may be exposed to rogue methane (Payne et al. 2017). Alawattegama et al. have provided evidence that communities near the unconventional shale gas extraction activities (within 4 km) suffered from degraded groundwater quality (Alawattegama et al. 2015). Moreover, Weinberger et al. have observed that multiple health symptoms are correlated to the exposure to unconventional natural gas development within 1 km (Weinberger et al. 2017). Middle-scale studies provide evidence of the environmental impacts of the shale oil and gas extraction activities in a distance from 10 to 20 km. Multiple sorts of evidence in diversified studies have been found on the environmental impacts caused by these shale oil and gas extraction activities. These environmental impacts are not limited to alkalization of soil (Pärn 2001), air quality issues (Liblik and Pensa 2001;Olaguer 2012), groundwater quality issues (Hildenbrand et al. 2016), surface water contamination (Burgos et al. 2017), coastal water contamination (Panigrahy et al. 2014;de Souza et al. 2014), etc. Accidents and failures in these shale oil and gas extraction activities may lead to hazardous environmental impacts on a much larger spatial scale in coastal and marine areas (Bayramov and Buchroithner 2015;Gomiero et al. 2015;Goto, Fan, and Sakai 2010), the diameter of the affected area reaches up to 2 arc-degrees (more than 200 km).
Besides the multitudinous and concrete impacts on local environments and human health, the spatiotemporal development of these shale oil and gas extraction activities has drawn increasing public attention. The average production life of these shale oil and gas wells is about 10 ~ 30 years (Dilmore et al. 2015). The short production life of a well in comparison to the entire long history of drilling activities indicates incessant boom and decay throughout the spatial extent of the Marcellus Shale. The drilling activities always form hot spots around the initial successful detection of oil and gas. The revisit of the hot spots is then unlikely after the dry out of the wells until the full recovery of the shale or until the new technologies improve the yield. This pattern of production forms the continuous transmission of clustered drilling activities from one place to another. Numerous spatiotemporal analysis approaches and models are available to describe similar dynamic geographic phenomena as sequential changing spatiotemporal objects, e.g., the view snapshot (Armstrong 1988), the unified spatiotemporal information model (Worboys 1994), the three-domain representation (Yuan 1999), the event representation (Yuan 2001;Peuquet and Duan 1995), the geo-atom representation (Goodchild, Yuan, and Cova 2007), the continuous spatiotemporal model ( Van de Weghe et al. 2014), etc. The general paradigm of these approaches first controls the time variable and detects  Gomiero et al. (2015)~2 arc-degrees Ocean gas platform exploration & exploitation Marine environmental pollution the spatial objects (or other similar elements) at each time node, then traces the evolution of the spatial objects over time. These approaches are intuitively understandable and particularly effective for limited spatiotemporal objects. However, most of them are not adequate to summarize the chaotic, fragmentary, and free-directional transmission of these shale oil and gas hot spots on diversified spatiotemporal scopes. The dense shale oil and gas extraction activities are easy to form scattered hot spots in space-time, and it may take tremendous efforts to trace them over time. Sadahiro (Sadahiro 2001) proposed a wave model through the spatiotemporal analysis of surface changes using the concepts of tessellation. This attempt provides valuable insight to modeling the transmission of dynamic geographic phenomena as wave propagation. Similar wave analogy and representations are also applied to urban growth (Dietzel et al. 2005), etc. The wave representation uses a different paradigm, the formation of a spatiotemporal wave is an automatic trace of transmission. A wave itself includes the spatiotemporal evolution of local activities and it is easier to interpret.
This study develops a concept of temporal object and a systematic wave representation to model and inspect the spatiotemporal evolution patterns of the shale oil and gas extraction activities in the Marcellus Shale. Particularly, a comprehensive set of attributes is developed for the spatiotemporal wave to summarize the general patterns of these chaotic and fragmentary activities on multiple spatiotemporal scopes. The rest of this paper is organized as follows: Section 2 introduces the study area and the data set, Section 3 details the identification of temporal objects, the formation of waves, and a group of attributes that are defined to quantify the waves. Section 4 shows the results of the modeling and the analysis, while Section 5 discusses the accuracy, driving force, and environmental impacts of these waves. Finally, conclusions are drawn in Section 6.

Study area
The Marcellus Shale, namely Marcellus Middle Devonian age organic-rich formation, sits in the northeast of the United States. This shale stretches in the subsurface from the Appalachian Mountains to central Ohio and extends from the south half of New York to the northeast edge of Kentucky ( Figure 1). The footprint of the Marcellus Shale covers over 2.46 × 10 5 km 2 (95,000 square miles), with a prospective areal coverage of about 1.86 × 10 5 km 2 (72,000 square miles). The depth of the oil window limits to the subsea depth range of 0.3-1.5 km while the depth of the gas window limits to the subsea depth range of 1-2 km (Popova et al. 2009;Eia 2017).
The majority of the area in the Marcellus Shale is the rural region. 64.55% of the land is covered by forest, including 49.86% deciduous forest, 2.37% evergreen forest, and 12.32% mixed forest. 5.09% of the land is covered by aquatic habitats (mainly Lake Erie) and 8.44% is the developed area. The terrain above the Marcellus Shale belongs to the northeast quarter of the Ohio River Basin. The upstream Ohio River is the dominant river. The streamflow of the Marcellus Shale region discharges at the outlet near Huntington, WV. The streamflow of the Ohio River at Huntington has an annual discharge of 2130 m 3 /s (75240 ft 3 /s) and an average gage height of 8.46 m (27.77 feet), its peak is in March and trough in September (USGS 2018). Several upstream tributary rivers discharge into the Ohio River, including the Allegheny River, the Monongahela River, the Muskingum River, the Little Kanawha River, the Hocking River, the Kanawha River, and the Guyandotte River. The climate categories in the Marcellus Shale are mainly humid subtropical and humid continental (Koppen 1936;Beck et al. 2005). Seasonal changes in this region are apparent, showing a warm wet summer and cold dry winter. The annual average wind speed is less than 5 m/s (Gunturu and Schlosser 2012), which is less efficient for the diffusion of air pollutants. Considering the geographic settings of the study site, the spatiotemporal analysis employed an inland mode to model the spatiotemporal wave propagation. The maximal influence distance is set to 20 km, it is the maximal influence distance for air quality, groundwater, and inland surface water (Table 1). The spatial resolution is set to 1 km, and it is sufficient to capture the spatial variance of the dynamics. In comparison, marine shale oil and gas extraction should employ a coastal mode. A coastal mode should have a longer influence distance and coarser spatial resolution. However, it is beyond the discussion of this study.
The Marcellus Shale spatially overlaps with mainly five states, namely Pennsylvania, Ohio, West Virginia, Virginia, and New York. It is partially or completely covered by 204 counties. The population of these 204 counties is about 17.89 million by 2018, and most of them are concentrated in big cities, including Pittsburgh and Scranton, PA, Cleveland and Akron, OH, Buffalo, NY, and Charleston, WV (U.S. Census Bureau 2018). The overlap of the Marcellus Shale with Tennessee, Kentucky, and Maryland is fairly small, and these three states are excluded from the following analysis.

Data sets
The spatial, temporal, and type information of the shale oil and gas extraction activities were collected from four data sources. The well data in Pennsylvania were downloaded from the Pennsylvania Spatial Data Access (PASDA https:// www.pasda.psu.edu/uci/DataSummary.aspx?data set=1088). The dataset is created and maintained by the Pennsylvania Department of Environmental Protection (PADEP), and it was updated to December 18th, 2018 when downloaded. The well data in Ohio were collected from the Ohio Oil & Gas Well Locator (https://gis.ohiodnr.gov/ MapViewer/?config=oilgaswells) of the Ohio Department of Natural Resources (ODNR), and it was updated to December 23rd, 2018 when downloaded. The well data in West Virginia were downloaded from the GIS Server (https://tagis.dep.wv. gov/site/) of the West Virginia Department of Environmental Protection (WVDEP), and it was updated to January 1st, 2019 when downloaded. The well data in New York were downloaded from the Downloadable Well Data (https://www.dec.ny. gov/energy/1603.html) of the New York State Department of Environmental Conservation (NYSDEC), and it was updated to January 2nd, 2019 when downloaded. The New York dataset is a table. The well locations in the New York dataset were displayed as point features in ArcGIS 10.5 in the NAD 1983 State Plane (FIPS 3102) coordinate system with reference to NAD 1983 datum. These point features were then exported as an individual shapefile. The datasets of the other three states all are shapefiles of point features. A small set of well information (about 8,600 wells) was collected from the Virginia Department of Mines Minerals and Energy (DMME). However, these data do not include critical temporal attributes, thus are not useful in this study. Adequate well information has been collected for the spatiotemporal modeling. In total, 606,911 wells were collected from the four data sources. 486,552 of them contain useful date attributes in the attribute table (e.g. permit date, complete date, and plug date), and the rest have void temporal attributes. Since the temporal attribute is critical in the spatiotemporal modeling, only these 486,552 wells were used in this study. 94% of them are conventional shale oil and gas wells, and the rest are nowadays referred to as unconventional shale oil and gas wells ( Table 2). The earliest dated conventional oil well was drilled in 1826, the earliest dated conventional gas well was drilled in 1815, and the earliest dated unconventional oil and gas well was drilled in 2005. The initial stage (1815-1867) of the shale oil and gas extraction activities has drilled only a few wells, which is of little value to model and discuss. Hence, the study period is determined as from 1868 to 2018, the amount adds up to one century and a half. The 486,552 wells from four datasets were merged into a single shapefile of point features in the NAD 1983 State Plane (FIPS 3402) coordinate system with reference to North American 1983 datum.
70% of these wells are four types of producing wells, including 83,364 (17%) conventional oil wells, 169,514 (35%) conventional gas wells, 62,629 (13%) conventional hybrid wells that produce both oil and gas, and 24,372 (5%) unconventional wells that mainly produce gas. Since different types of producing wells may have different levels of impacts on local environments, the spatiotemporal evolution of these four types of shale oil and gas wells are further modeled, respectively. Other wells have various types and different functions in the shale oil and gas industry. The amount of each type is small and less valuable to discuss. These non-producing wells are hence not modeled individually.

Methods
To model the spatiotemporal evolution of these shale oil and gas extraction activities in the Marcellus Shale, the density surface of the active producing wells was first estimated each year. A focal statistic filter was then employed to identify the temporal objects from the time-well density plot at each location (surface cell). These temporal objects were used to form the spatiotemporal waves of these active shale oil and gas extraction activities. Finally, a group of spatial and temporal attributes and events were derived for the spatiotemporal wave to analyze the general patterns of these shale oil and gas extraction activities.

Estimate the density surfaces of active producing wells
The density surface of active shale oil and gas wells in the Marcellus Shale was used to represent the general spatial pattern of these shale oil and gas extraction activities. For each year from 1868 to 2018 (one century and a half), the active wells were selected from the entire oil and gas well database as a subset. The environmental impacts (for instance, air pollutions, surface water contamination, and groundwater quality degradation) are assumed to stop soon when the production of oil and gas is over and the well is plugged. Thus, the active wells are defined as those wells whose life cycle from completion year to plug year covers the single time node of density surface estimation. The subset of wells was used to estimate a surface density through the ArcGIS 10.5 Kernel Density tool (Hart and Zandbergen 2014).
Considering the spatial scale of the environmental impacts caused by these shale oil and gas extraction activities, the searching radius was set to inland maximal influence distance, which is 20 km ( Table 1). The cell size of the surface was set to 1 km since the spatial variance of the environmental impacts within this distance is small enough to ignore. ArcGIS 10.5 Kernel Density tool used the quartic distribution while conducting the kernel density estimation, and the identical weight was used for each well location. In total, 151 time nodes were set and 151 density surfaces of active shale oil and gas wells were consequently estimated to model the spatiotemporal evolution of the shale oil and gas development in the Marcellus Shale. These density surfaces were sequentially stacked to form a space-time cube for the following identification of spatiotemporal waves. These waves were then used to represent the spatiotemporal propagation of these shale oil and gas extraction activities.

Identify the temporal objects as spatiotemporal wave footprints
When a spatiotemporal wave of dense shale oil and gas extraction activities passes through space, the local magnitude of well density shows a waveform curve over time as the footprint of the spatiotemporal wave. Inversely, to detect the spatiotemporal wave, it's critical to first identify the wave's footprint at each location within the domain. In this section, a concept of temporal object is defined to represent a wave's footprint.
Critical years were identified on the time-well density plot at each location to capture the wave footprint (temporal object). At each location (surface cell) in the study area, the well density of the one century and a half were extracted and plotted against time.
Since the shale oil and gas extraction activities are locally dynamic and noisy, numerous minor crests and troughs may hamper the detection of valuable crests and troughs in the general pattern ( Figure 2). A focal statistic filter was employed to skip the minor crests and troughs and to automatically find the critical crests and troughs in the general pattern. A period of P years was set as a temporal window to sequentially scan the time-well density plot. While scanning, the temporal window was mounted on acentered year, and P/2 years before and after the centered year were included in the temporal window. The well density value of the centered year is defined as the central well density value. The central well density value was compared with the statistics of all the well density values within this temporal window. If the central well density value equals the minimum value within the window, the centered year is defined as a trough year in this surface cell; otherwise, there is at least one year in which the well density is smaller than the centered year and the centered year is skipped. Similarly, if the central well density value equals the maximum value within the window, the centered year is defined as a crest year in this surface cell; otherwise, the centered year is skipped.
The temporal object is defined as a period that comprises of the following components: centered on the crest year, the previous adjacent trough year is defined as the start year of the wave in the cell, and the posterior adjacent trough year is defined as the end year of the wave in the cell. The phase bounded by the start year and the crest year is defined as the boom phase, while the phase bounded by the crest year and the end year is defined as the decay phase ( Figure 2). The well density value of the crest year is extracted as the crest well density of the temporal object. The temporal objects were identified at each location (surface cell) in the study area. A cell may have several temporal objects since multiple spatiotemporal waves of dense shale oil and gas extraction activities may revisit the same hot spots of these activities. A temporal object indicates the beginning, crest, end, and the maximum well density when a spatiotemporal wave passes through a location.

Form spatiotemporal waves to represent the propagation of shale oil and gas extraction activities
A space-time expansion algorithm has been developed in ENVI 5.3 and IDL 8.5 to form the spatiotemporal waves from discrete temporal objects. The proposed algorithm is a novel algorithm that is particularly designed to form spatiotemporal waves. The expansion procedure of the algorithm searches both space and time domains, and it automatically traces the free transmission of wave's footprints (temporal objects).
To implement the algorithm, a void space-time cube (ST-cube) was created in a finite difference method to contain these identified temporal objects. The void ST-cube fully covers the study area, starts from 1868, and ends by 2018. It has a uniform spatial resolution of 1 km and a uniform temporal resolution of 1 year, which are consistent with the surface density ST-cube. To contain a temporal object, the ST-cell at the same location as the temporal object and in the crest year of the temporal object is marked as the crest year cell of this temporal object (Figure 3(a)). Similarly, at the same location, the ST-cell in the start year is marked as the start year cell and the STcell in the end year is marked as the end year cell (Figure 3(a)). The ST-cells bounded by the start year cell and the end year cell compose the sub-ST of the temporal object in the ST-cube. The crest year cells were then used to represent the temporal objects in the formation of spatiotemporal waves (Figure 3(b)). The first step of the developed algorithm is to place a seed. From the initial year of the ST-cube, the algorithm scanned the study area to search the earliest crest year cell, which represents the temporal object that has the earliest crest year. If no crest year cell is found, the algorithm skips to the next year and repeats the scan (Figure 3(c)). If a crest year cell is found, the crest year cell is marked as a seed with a unique wave id (Figure 3(d)). Correspondingly, the temporal object behind the marked crest year cell is included in the wave.
The developed algorithm then searched the adjacent space-time of the seed to include spatiotemporally adjacent crest year cells. The adjacent spacetime is defined as the cells in the ST-cube that are 1) spatially adjacent to the seed cell, and 2) temporally within a period of G years before/after the year of the seed (Figure 3(e)). Once a crest year cell is within the adjacent space-time of the seed, the crest year cell is marked as a new seed with the wave id ( Figure 3(f)). Correspondingly, the temporal object behind the newly marked crest year cell is included in the wave. The temporal objects behind the newly marked crest year cells are considered as the spatiotemporal expansion of the identified spatiotemporal wave. The expansion was repeated until no crest year cells (temporal objects) were found in the adjacent space-time of all the seeds. All the included temporal objects were assembled as a single spatiotemporal wave. The wave was exported as a new object (Figure 3(h)), and all the critical cells of these included temporal objects (crest year cells, start/end year cells) were reassigned as the void cells in the ST-cube to avoid duplicated detection (Figure 3(g)). The algorithm restarted the scan with a new unique wave id. The scan-expansion-export procedure was repeated until there was no temporal object left in the ST-cube. The crest years of temporal objects in a wave were integrated and interpolated as a propagation surface, the start years and end years of temporal objects in a wave were integrated and interpolated as the start and end time boundary of the wave, respectively (Figure 3(i,j)).
A group of attributes was derived from the spatiotemporal wave in ArcGIS 10.5 (Figure 3(k)). The spatial extent of the wave is the minimum area that contains all the temporal objects of the wave. The emergence location is defined as the geometric center of the temporal objects that have the earliest crest year. Similarly, the disappearance location is defined as the geometric center of the temporal objects that have the latest crest year. The propagation direction is the orientation from the emergence location to the disappearance location. The propagation distance is the distance between the emergence location and the disappearance location. The emergence year of the propagation is the earliest crest year and the disappearance year of the propagation is the latest crest year. The duration of the wave is defined as the latest crest year subtracted by the earliest crest year. The maximum well density is defined as the maximum value of crest well densities.
Because of the spatial clustering of these shale oil and gas extraction activities and the relatively coarse temporal resolution (1 year), the footprints of the spatiotemporal wave may form small adjacent hot spots. The geometric centers of these small hot spots were derived, respectively. A path is then defined as the line segment connecting a geometric center to its nearest geometric center of the earlier/later adjacent year. All the paths of a wave together form the route of the propagation. Besides emergence and disappearance, three extra events are defined on each center of these small hot spots based on the propagation route. A split event is defined as the connection with multiple later paths, a merge event is defined as the connection with multiple earlier paths, and a spread event is defined as the connection between a single earlier path and a single later path. A spread event indicates the monodirectional transmission of the densest shale oil and gas extraction activities from one location to the next. A split event indicates that the densest activities develop into multiple separated groups and they propagate into different directions. In contrast, a merge event indicates multiple groups of densest activities coincidently find the same location to develop and cluster into one.
The spatiotemporal waves indicate the trend of the shale oil and gas extraction activities on different spatiotemporal scales. A single wave includes the transmission and events of these activities within itself, and one wave indicates a single trend. A larger and longer wave occupies a larger sub-ST in the STcube, and it indicates a more general trend in the spatiotemporal evolution of these activities. In contrast, small and short waves disappear very quickly, and they cannot show apparent propagation patterns. These small and short waves are like sparkles and they should be excluded from summarizing the general evolution patterns. In this study, all the waves are extracted, but only these long waves (longer than a decade) and large waves (larger than 10 3 km 2 ) are used to analyze the general patterns of these shale oil and gas extraction activities.

Waves to represent the propagation of extraction activities
Two critical parameters of the spatiotemporal wave representation control the extent of the waves in the space domain and the duration of them in the time domain. One is the period P of a temporal window to identify temporal objects and the other one is the gap year G to form the spatiotemporal waves. A larger P is more insensitive to ignorable lowcontrast momentary crests while keeping apparent crests. Consequently, a larger P is capable to identify a few longer-lasting significant temporal objects while avoiding temporally fragmentary objects. Considering that the average life length of an active producing well is about 30 years (Dilmore et al. 2015), a P = 15 was determined to detect the apparent crests for the formation of spatiotemporal waves. An example in Ohio (Figure 4) shows that the determined P is sufficient to detect two apparent crests, one is in 1948 and the other is in 1997, and numerous observable but trivial crests are skipped as noise. As a result, two temporal objects were identified. One started in 1898, crested in 1948, and ended in 1968. The other one started in 1968, crested in 1997, and ended in 2018 since it's the end of the time series. Meanwhile, the crests with well density lower than 0.3 well/km 2 were filtered off as background noise. Only these significant crests with a well density larger than and equal to 0.3 well/km 2 were used to identify temporal objects.
The gap year G is a parameter defined in the time domain but it has effects on both the extent and duration of the wave. A larger G allows the assemble of more crest year cells/temporal objects and consequently results in larger and longer waves. Since the available crest year cells in space-time are limited, G thus also controls the abundance of waves along the axes of duration and areal size (Tables 3 and 4). Using a larger G to form waves, the logarithmical decrease was observed in the total number of waves and the total number of spatially and temporally fragmentary waves (areal size 1 ~ 10 km 2 and duration 1-4 years). Meanwhile, extremely large and long waves occurred only in the case of using large G. However, a large G may lead to a particular wave that shows the revisit of the previous booming locations, and the revisit is difficult to visualize on 2D maps. In this study, a G = 1 was determined to form rigorously continuous waves. Despite these fragmentary waves, large waves (areal size > 10 3 km 2 ) and long waves (duration > a decade) are still abundant enough to disclose the general spatiotemporal evolution patterns of these shale oil and gas extraction activities in the Marcellus Shale.
The space-time expansion algorithm was used to form waves from identified temporal objects. A wave is constructed as a three-layer structure in space-time ( Figure 5(a)), which includes the main propagation surface interpolated from the rigorously continuous crest years, an earlier start time boundary, and a later end time boundary interpolated from start/end years. The space-time bounded by the three layers is the sub-space-time occupied by the wave. Mapping the propagation surface onto the 2D surface of the earth enables the delineation of the propagation progress in space. A group of features and attributes are derivable from the propagation surface ( Figure 5(b,c)).
The propagation direction from the emergence location to the disappearance location shows the general direction of the evolution process. The route map delineates the details of the propagation, whose temporally longest path is definitely consistent with the propagation direction. Route map contains many interesting facts of the local evolution process in a single wave, e.g., the circular propagation of the shale oil and gas extraction activities in Guernsey-Muskingum-Perry-Morgan-Noble-Monroe counties in Ohio from 1982 to 1993 ( Figure 5(c)). Spatial extent and propagation direction are more general for the comparison between waves. Hence, the spatial extents and propagation directions of the large and long waves are used in the following sections to disclose the general spatiotemporal evolution patterns of these shale oil and gas extraction activities in the Marcellus Shale.

General patterns of the waves
In total, 527 rigorously continuous waves were formed in the Marcellus Shale in the past one century and a half using the spatiotemporal wave representation (P = 15, G = 1). Five of them are temporally long waves whose duration is longer than a decade (Figure 6), and 39 of them are spatially large waves whose extent area size is larger than 10 3 km 2 (Table 5). A wave includes the spatiotemporal evolution of the shale oil and gas extraction activities on a local scale, and it shows a single trend in the chaotic dynamics of these activities. A larger and longer wave shows a more general trend in space-time since it occupies a broader sub-space-time. In contrast, short and small waves do not indicate apparent wave propagation of these shale oil and gas extraction activities. These sparkle-like waves are not significantly different from  noises and they are then not used to analyze the general patterns of the shale oil and gas extraction activities. This study uses the 39 largest waves to demonstrate the most general propagation trends of the shale oil and gas extraction activities in the Marcellus Shale. The average duration of large waves is 6.46 years, the average extent area size is 5128.54 km 2 , and the average propagation distance is 47.19 km. The propagation direction has an apparent contrast along the north-south direction and a weaker contrast along the west-east direction. Thirty-one out of the wave propagation directions have a southward The east of Ohio is the most chaotic region for dynamic propagation of shale oil and gas extraction activities, 27 out of the 39 large waves are involved with this region. The west of West Virginia is second to the east of Ohio, which involves 12 waves. The west of Pennsylvania and the west of New York each involves four large waves, respectively. Apparently, the western edge of the Marcellus Shale is beneath the west edge of the Appalachia Mountains and it has comparatively flatter topography (Figure 1). The flat topography allows the free propagation of these shale oil and gas extraction activities and makes the evolution patterns complicated in the east of Ohio and the west of West Virginia. In contrast, the propagations of these shale oil and gas extraction activities in the deep rolling mountainous Appalachian regions are dominated by the trend of the mountains and the corresponding waves have less freedom in the east/west directions. The topographical pattern leads to the simpler propagation patterns in Pennsylvania and New York.
The spatially largest wave is also the temporally longest wave (Figure 6, wave ID: 461) and it extends along the west side of the Appalachia Mountains. This wave partially covers all of the four states and exists from 2005 to 2018 (the end of the time series). The propagation originated from two emergence locations. The first was in the southwest of New York and it spread southward along with the mountain's trend. The other one originated from the east of Ohio and spread to Pennsylvania and West Virginia in three directions. An interesting revisit (overlap of wave spatial extents) is observed in West Virginia between wave ID 350 and wave ID 461. The wave ID 350 is from west to east through West Virginia, and it lasts from 1989 to 2000 when the shale oil and gas extraction activities were still dominated by the conventional Class II injection wells. The Wave ID 461 (2005-2018 lived through the booming of the unconventional shale oil and gas extraction activities that incorporate hydraulic fracturing, horizontal drilling, and well pad. The new techniques stimulated a new broader wave that revisited the entire wave ID 350 and followed a parallel direction.

Different propagation patterns of four types of extraction wells
The disaggregated spatiotemporal evolution patterns of the four types of shale oil and gas extraction activities have been modeled, respectively. The diversified patterns of propagation were observed (Figure 7). The propagation of conventional gas wells is the most chaotic, 21 corresponding large waves stretch out all the four states in the study area. The propagation of conventional hybrid wells (producing both oil and gas) is second, 11 corresponding large waves stretch out mainly in the east of Ohio and the north of Pennsylvania. The propagation of conventional oil wells is comparatively simple and local, only 9 corresponding large waves are found. The unconventional wells are newly emerged in the last two decades, the short booming period has not yet formed a chaotic propagation pattern. Only one large wave is found in the Marcellus Shale that is purely from the unconventional wells.
The conventional gas wells show the earliest large wave. The earliest large wave started in 1908 and ended in 1913, which is on the west edge of the Marcellus Shale (Figure 7 The propagation also shows observable diversified propagation directions for the four types of shale oil and gas extraction activities. The propagation directions of conventional gas wells have a significant southward trend and an eastward trend, with an average propagation distance of 53.71 km. Since a wave may have multiple propagation directions, 32 propagation directions were found in the 21 large waves of conventional gas wells. Twenty-two of them have a southward component and 20 of them have an eastward component. In contrast, the number of northward components and westward components both are 7, respectively. The propagation of conventional oil wells shows a southward trend, with an average propagation distance of 30.14 km. Among the 9 formed large waves, the numbers are 7, 6, 4, and 3 for the southward, westward, eastward, and northward components, respectively. The propagation of conventional hybrid wells shows a westward trend, with an average propagation distance of 44.12 km. Among the 11 formed large waves, the numbers are 10, 8, 8, and 6 for the westward, northward, southward, and eastward components, respectively. The westward propagation of conventional hybrid wells in northeast Ohio is the most apparent. The unconventional wells have only two propagation directions in a single large wave, propagating to the east at an average distance of 68.99 km. Revisits are common in the propagation of the four types of shale oil and gas extraction activities, especially in these chaotic regions along the west edge of the Marcellus Shale. Particularly, a complete revisit is observed for conventional oil wells in the west edge of the Marcellus shale. The first wave (Figure 7(a), wave ID: 57) started in 1953 and ended in 1960, propagating northward. 24 years later, the second wave (Figure 7(a), wave ID: 107) revisited this region, from 1984 to 1987, propagating westward. The two waves overlapped well with each other. This revisit is after the first (1973) and second (1979) oil crisis (Venn 2016), and it is possibly triggered by the first booming of the shale oil and gas extraction activities in the Marcellus Shale.
Another thing to notice is that the largest wave in Figure 6 (wave ID: 461) does not show up in any wave map of the four types of shale oil and gas extraction activities. This absence indicates that this largest wave is not purely contributed by any single type of extraction activities. The same interpretation is applicable for most of the waves. In contrast, the wave in Ohio (Figure 6, wave ID: 283, 1980ID: 283, -1993 shows a similar spatial extent and duration as the wave formed by conventional hybrid wells (Figure 7 (c), wave ID: 21, 1982ID: 21, -1993. This coincidence indicates that the wave in Ohio from 1980 to 1993 may be primarily contributed by the hybrid shale oil and gas extraction activates. However, the wave of conventional hybrid wells (Figure 7(c), wave ID: 21) altered the propagation direction. Its propagation direction is consistent with the circular propagation in the aggregated wave ( Figure 6, wave ID: 283), which has a roundabout path.

Spatial and temporal accuracy of the wave representation
Bi-temporal high-resolution natural-color aerial images collected from the Ohio Statewide Imagery Program (OSIP) (http://gis5.oit.ohio.gov/geodata download/) are employed to assess the spatial and temporal accuracy of the wave representation. The footprints of historical conventional wells were small and might fully recover after they stopped producing and were plugged. In contrast, the recently emerged unconventional wells require a large well pad, and they are verifiable from high-resolution remote sensing data. Two images were collected for Carroll County in Ohio from OSIP, one was captured in 2006 and the other in 2014. These images have a spatial resolution of 0.3 m (1 foot) and are adequate to interpret well pads and well locations (Figure 8). By comparing the two images, the new wells in well pads were interpreted and counted as the active producing wells between 2006 and 2014. More particularly, these interpreted wells can represent the active unconventional wells between 2011 and 2014 since there was no unconventional well been permitted and drilled before 2011 in Carroll County (Liu 2021). These interpreted wells were used as ground truth. The number of dated active unconventional wells in the Section 2.2 data set was compared with the number of interpreted wells to assess the temporal accuracy of the wells. On the other hand, the locations of dated active unconventional wells in the Section 2.2 data set were compared with the interpreted well locations to access the spatial accuracy of the wells. 97 well pads were interpreted and 352 corresponding unconventional wells were verified in Carroll County from the OSIP 2006/2014 aerial image pair. The number of dated active unconventional wells in the Section 2.2 data set is 385 in this county from 2006 to 2014 (actually from 2011 to 2014). The temporal accuracy of the used data set is 92.15%. Moreover, the completion year was considered as the start time node to define the active producing wells, which may lead to limited bias. For instance, the environmental impacts such as the land cover disturbance may already start in several months before the completion of the well. The spatiotemporal waves were formed with a 1-year temporal resolution in the past one century and a half. In this case, the limited bias of several months has tolerable impacts on the results and conclusions. On the other hand, the average spatial displacement between the interpreted well locations and the well locations in the Section 2.2 data set is 1.51 m, which is close to the georeferencing accuracy of the aerial images (1.5 m, 5 feet). The spatial displacement is far less than the spatial resolution of the density surfaces (1 km). The small spatial bias has an ignorable impact on the spatial accuracy of the wave. In one word, the spatial and temporal accuracy of the used data set is adequate for the density estimation of the active shale oil and gas extraction activities. The high accuracy of the used data set and the density surfaces lay down the solid foundation for the reliable formation of the spatiotemporal waves.

Driving force of the propagation of extraction activities
Fossil fuel consumption in the U.S. has a long history to rely on imports and the U.S. is still on the way to reach a balanced import and export through the development of unconventional shale oil and gas extraction activities (Neff and Coleman 2014). This consumption pattern indicates that the U.S. energy security is easily vulnerable to the dynamic international fossil fuel market, which turns back to trigger the domestic development of shale oil and gas extraction activities. The international crude oil price was collected from Macrotrends (https://www.macrotrends.net/1369/ crude-oil-price-history-chart) as the indicator to represent the international fossil fuel market. The numbers of newly drilled wells each year in the Marcellus Shale were plotted and overlapped with the international oil price curve (Figure 9). An apparent coincidental trend is observed from the plot. Two booms of the U.S. domestic shale oil and gas extraction activities are coupled with two skyrocketing oil price trends. The first one was triggered by the first and second oil crises (Venn 2016) and returned to inactive after OPEC's alter on oil production policies (Hamilton 1996;Griffin and Neilson 1994). The second one was triggered by the massive rise in the oil price since 2002 and dramatically declined after the 2008 global economic crisis (Thompson 2017).
The global dynamic fossil fuel market also has apparent impacts on the spatiotemporal propagation of these shale oil and gas extraction activities in the Marcellus Shale (Table 6). The booming of wells would like to form larger and longer waves. The average wave extent sizes and wave durations declined after the boom of the extraction activities. The latest period (2008 ~ 2018) showed a further declining trend of drilling wells, but it also contributed to the largest and longest wave in the Marcellus Shale. The recession of drilling activities mainly has two reasons. The first is the revolutionary development of unconventional shale oil and gas extraction activities. The new complicated technologies dramatically increased the efficiency in shale oil and gas production, and it encouraged a limited number of wells while yielding sufficient fossil fuel products. On the other hand, the U.S. government recently keeps showing a negative attitude on shale oil and gas extraction activities, e.g., Obama's permanent ban on the Arctic and Atlantic offshore oil and gas drilling activities (2016) and Biden's moratorium on oil and gas leasing on federal lands and waters for 60 days (2021). Due to these reasons, an extended period of inactive shale oil and gas extraction activities is expectable in a short term in the Marcellus Shale.
The primary objective of this basic research is to uncover the spatiotemporal evolution patterns of the chaotic dense shale oil and gas extraction activities. The spatiotemporal wave model is successfully applied to visualize these activities as wave propagation. The results help people understand the story of shale oil and gas development straightforwardly within a large spatial scope and a long temporal scope. The outcomes indeed help the shale oil and gas industry examine the regions of lack exploration in different technical generations. On the other side, the clear illustration of the shale oil and gas extraction activities can also help environmentalists to bond the environmental degradation with these activities within a reliable spatial domain and period. Moreover, these waves have not only well demonstrated the free transmission of these activities in space-time, but also indicate apparent correspondence to exterior driving force (Table 6). There is then a great potential to predict the shale oil and gas extraction activities using the wave model. Though beyond the scope of this study, it would be greatly encouraged for any further contribution and application on the wave representation to predict the new shale oil and gas extraction activities and other geographic phenomena.

Environmental and health impacts along with the propagation of extraction activities
The shale oil and gas drilling activities have negative impacts on local environments and they elevate the risks to human health, especially these unconventional shale oil and gas extraction activities that incorporate well pad, horizontal drilling, and hydraulic fracturing (Liu 2021). The leakage of gaseous materials while drilling may lead to air pollution issues, e.g., volatile organic compounds (VOCs) and nitric oxide (NO); hydraulic fracturing may break up the confining beds of aquafers and it may result in groundwater contamination; improper treatment and disposal of produced wastewater release toxicants and total dissolved solids (TDS) into the surface water, and they induce strong stress on surface water security (Annevelink, Meesters, and Hendriks 2016). The drilling activities and the construction of associated facilities produce odors, noise, and light pollution. These negative factors pose mental stress on residents in the neighborhood (Ferrar et al. 2013). These impacts may directly cause health symptoms in residents who live near the unconventional shale oil and gas extraction activities, e.g., sleep disruption, headache, and throat irritation (Weinberger et al. 2017). Moreover, unconventional shale gas extraction activities have also caused inevitable landscape changes. The landscape changes  are not limited to deforestation and the land cover conversion into the impervious surface through the construction of well pads and pipeline right-of-ways (ROWs). These land cover conversions are not only a disturbance to local biological communities (Brittingham et al. 2014) but also an alteration of local hydrological functions (Brantley et al. 2014). Without the interception of the tree canopies and the infiltration into the ground, more precipitation becomes surface runoff. The extra surface runoff increases drainage discharge and elevates the risks to downstream riverine buildings. The latest wave in the Marcellus Shale ( Figure 6, wave ID: 461) has witnessed the recession of drilling activities along with the booming of unconventional shale oil and gas extraction activities. This long and large wave shows the general retreat of active producing wells from the northern and western edges of the Marcellus Shale to the center of the shale. The retreat is from the highly populated regions to the mountainous rural areas of the deep Appalachian regions. This retreat may indicate the mitigation of the health impacts on human communities from these shale oil and gas extraction activities. Instead, stronger stress on local biological communities is expectable, especially on those forest-raised wild species, e.g., songbirds and reptilians. These issues call for more attentions and investigations.

Conclusions
This study has investigated the spatiotemporal development of shale oil and gas extraction activities in the Marcellus Shale as wave propagation. Density surfaces of active producing wells were estimated in each year from 1868 to 2018. The density curve at each location (surface cell) was plotted and a 15-year focal statistic filter was Figure 8. An illustration of the wells and well pad that are interpreted from the Ohio Statewide Imagery Program (OSIP) aerial image. The high spatial resolution (1 foot) is adequate to locate the wells on the image. The well pad boundary is determined by tracing the fence that encloses the well pad facilities. employed to identify temporal objects. Temporal objects were then used to form waves by expanding spatiotemporal adjacent crest year cells of the temporal objects while maximally allowing 1 gap year. Comprehensive attributes were derived for the wave. The general patterns of waves were analyzed to summarize and inspect the development of the shale oil and gas extraction activities in the Marcellus Shale. 527 waves were formed in the Marcellus Shale. Only a few of them are large and long waves that show the general patterns of these activities. The west and north edges of the Marcellus Shale have chaotic and fragmentary wave propagation, while the center in the Marcellus Shale shows several simple and large waves. The large waves of shale oil and gas extraction activities have a general trend in southward propagation. The trend follows the Appalachian region's mountainous topography. The largest and longest wave was from 2005 to 2018 and it partially covered all the four states in the study area. This wave witnessed the boom of the unconventional wells and the retreat of shale oil and gas extraction activities from Marcellus Shale's edges to its core. Diversified propagation patterns are found in different types of shale oil and gas extraction activities and the conventional gas extraction shows the most chaotic propagation.
The wave representation and corresponding analysis method developed in this study have been successfully applied to the spatiotemporal evolution of shale oil and gas extraction activities. The developed method has two critical parameters: the period P of a temporal window to identify temporal objects and the gap year G to form the spatiotemporal waves. These two critical parameters make it possible to control the formation of wave propagation on multiple scopes to meet diversified demands. The method developed in this study is also transferable and applicable to other dynamic geographic phenomena.

Disclosure statement
No potential conflict of interest was reported by the author(s).