Detecting and Mapping Gas Emission Craters on the Yamal and Gydan Peninsulas, Western Siberia

Rapid climate warming at northern high latitudes is driving geomorphic changes across the permafrost zone. In the Yamal and Gydan peninsulas in western Siberia, subterranean accumulation of methane beneath or within ice-rich permafrost can create mounds at the land surface. Once overpressurized by methane, these mounds can explode and eject frozen ground, forming a gas emission crater (GEC). While GECs pose a hazard to human populations and infrastructure, only a small number have been identified in the Yamal and Gydan peninsulas, where the regional distribution and frequency of GECs and other types of land surface change are relatively unconstrained. To understand the distribution of landscape change within 327,000 km2 of the Yamal-Gydan region, we developed a semi-automated multivariate change detection algorithm using satellite-derived surface reflectance, elevation, and water extent in the Google Earth Engine cloud computing platform. We found that 5% of the landscape changed from 1984 to 2017. The algorithm detected all seven GECs reported in the scientific literature and three new GEC-like features, and further revealed that retrogressive thaw slumps were more abundant than GECs. Our methodology can be refined to detect and better understand diverse types of land surface change and potentially mitigate risks across the northern permafrost zone.


Introduction
Rapid climate warming [1] and intensifying hydrologic cycles [2] at northern high latitudes are driving a suite of landscape changes, including changes in vegetation productivity [3], thawing of perennially frozen soils (permafrost) [4], and lake expansion and drainage [5]. These changes may pose a hazard to humans and infrastructure and alter ecosystems [6], particularly where permafrost is relatively ice-rich [7] and, therefore, may induce greater terrain consolidation following thaw (i.e., thermokarst). Detecting and mapping northern land surface change is thus emerging as a top research priority for understanding northern ecosystem change and mitigating potential risks associated with intensifying geomorphic activity (e.g., [8]).
For decades, remote sensing has enabled broad-scale assessment of environmental change at northern high latitudes. For instance, data from Landsat missions 5-8 (1984present) have been used to classify northern land cover and vegetation change at 30 m resolution across regional to continental scales [9][10][11]. Landsat time-series data were used by the European Commission Joint Research Centre (JRC) to generate the first global estimates of surface water change at 30 m resolution [12]. In the northern permafrost zone, where thermokarst causes vertical land surface displacement, multi-temporal digital surface models and digital elevation models stand to enhance the detection of thermokarst at its onset by tracking vertical changes at the land surface. Such data are typically acquired by drone or aircraft, and can support detailed investigations from local to regional scales (e.g., [8]). For instance, 4-D data derived using high-resolution (2 m) multi-temporal strip data from the ArcticDEM digital surface model [13] enabled estimates of the volume of recent lava flows over~1000 km 2 of active volcanic terrain in Alaska [14]. The pan-Arctic ArcticDEM data are freely available on the cloud computing platform Google Earth Engine [15], enabling rapid computations on terabytes of data for investigating geomorphic change at northern high latitudes.
Gas emission craters (GECs), which are only known to occur in the Yamal and Gydan peninsulas in western Siberia [16], are among the most abruptly forming geomorphic features in the northern permafrost zone. GECs are thought to form as a result of the migration of subterranean biogenic or thermogenic gases, from unfrozen saline deposits (i.e., cryopegs) or deeper geogenic sources via faults, and their accumulation beneath relatively impermeable layers of tabular ground ice and overlying clay-rich permafrost [17,18]. GECs show common geomorphic traits during three general stages of evolution [19] that can be readily observed with moderate-to high-resolution satellite imagery [16,20]. First, GECs are preceded by a pingo-like feature up to ca. 6 m in height and 20-55 m in width, resulting from the mechanical deformation of near-surface permafrost caused by the build-up of excessive gas [16]. This initial stage of terrain disturbance may develop over several to ten years. Second, during the eruption stage, up to hundreds of cubic meters of ground material is exploded [16]. Blocks of ground material several meters in diameter can be ejected radially tens of meters, scattering permafrost and unfrozen soils amongst peripheral vegetation and forming small thaw ponds when ejecta are sufficiently ice-rich [20]. Calderas of newly formed GECs (stage 2) are commonly ca. 35 m in diameter, ranging from 10-90 m [16,21] and reaching depths up to 70 m [17]. Internally, funnel-like upper caldera walls and nearly vertical lower walls of the evacuated GEC expose ice-and sediment-rich permafrost. In the final stage, the crater structure degrades into a sub-or non-circular turbid lake over the course of one to several years. Inputs of precipitation, surface water, and meltwater and solids from exposed ice-and sediment-rich permafrost, contribute to the infilling and shallowing of GECs, resulting in a turbid lake which can bear visual, morphologic, and biogeochemical resemblance to thermokarst lakes that are ubiquitous in the region [17]. Warming permafrost temperatures [4], which are associated with GEC formation [19], may also facilitate the development of thermokarst features such as retrogressive thaw slumps (RTS) [22,23]. Despite the abrupt geomorphic change caused by GECs, little is known about the potential risks to infrastructure and human populations at regional scales, or the implications for carbon cycling and feedbacks to global climate. This is because research on GEC abundance and distribution started only recently [19] and few have been studied remotely or on the ground. Only seven GECs have been reported in the scientific literature thus far [16,17,24,25] and there are no existing algorithms that can accurately detect GECs.
In this study, we use freely available moderate-to high-resolution satellite remote sensing data to detect changes in vegetation (LandTrendr), elevation (ArcticDEM), and surface water (JRC) on a cloud computing platform (Google Earth Engine) to detect and map land surface change across 327,000 km 2 of the Yamal and Gydan peninsulas. Objectives of this research were to develop a methodological framework to detect land surface change; to quantify the extent of land surface change and the ecosystems in which it occurred; and to identify and characterize new potential GECs within the context of land surface change across the study region.

Study Area
The northern and southern boundaries of our 327,000 km 2 study area were defined using the Circumpolar Arctic Vegetation Map (CAVM) bioclimatic zones and included Sub-zones C, D, and E [26,27]. These subzones are characterized by summer air temperatures of 9-35 • C and vegetation ranging from 5 to 100% vascular plant cover in open patches to closed canopy dwarf shrubs with moss layers 3-10 cm in thickness. The east and west boundaries were defined by the geographic extent of the Yamal and Gydan peninsulas themselves. Situated within the continuous permafrost zone of western Siberia [28], the relatively low-relief and low-elevation Yamal and Gydan peninsulas contain Quaternary alluvial-lacustrine-marine plains and terraces hosting abundant lakes, drained lake basins, rivers, and thermokarst features (e.g., retrogressive thaw slumps, thermo-erosion gullies) [29]. Permafrost in the Yamal-Gydan region ranges in thickness from 100 to >400 m and in temperature from −2 to −8 • C [29,30]. Tabular ground ice is widespread in the region, exceeding 20 m in thickness and extending for hundreds of meters ( [29] and references therein) and is commonly overlain by fluvial late Pleistocene-Holocene-aged sand, silt, and clay interbeds, and underlain by marine late Pleistocene sand deposits containing saline cryopeg lenses [31]. The co-occurrence of near-surface tabular ground ice, cryopegs [17], and relatively high gas content (including methane) in permafrost and deep deposits [32] are thought to be unique to the Yamal-Gydan region and the requisite conditions for GEC formation [16]. Consequently, GECs in the Yamal-Gydan region are increasingly sought, studied, and observed by the scientific community [17,19,24,29,[32][33][34][35]. While not necessarily considered the primary source of GEC methane [32,36], deep hydrocarbon reservoirs in the Yamal-Gydan region are among the largest in the world [37]. Infrastructure for oil and gas exploration and extraction, including roads, pipelines, processing plants, and workers barracks, and activities of the Yamal Nenets indigenous people are among the primary human uses of the region [38,39]. At the city of Salekhard (66.31 • N, 66.67 • E) in the south of the study area, mean annual and summertime (June-September) air temperatures from 1984 to 2018 were −5.3 • C and 10.5 • C, respectively [40].

Multivariate Change Detection Algorithm for Detecting and Mapping Land Surface Change
To determine the potential distribution of GECs across our study area, we developed a multivariate land surface change detection algorithm (CDA) using geospatial data within the high-performance cloud computing platform Google Earth Engine [15]. The CDA identified hotspot areas, where abrupt changes in vegetation, elevation, and surface water reflected terrain disturbance that is characteristic of GEC activity. The CDA used freely available, pan-Arctic geospatial data that were derived from moderate-to high-resolution satellite imagery (≤30 m pixels) and represented changes in (i) land surface reflectance [41], (ii) elevation [13], and (iii) surface water extent and distribution [12]. Because the CDA did not ascribe the type of land surface change, expert verification was required to confirm the presence or absence of GECs using high-resolution (sub-meter) satellite imagery. This approach allowed us to understand the distribution of GECs and also general landscape change across our study area. The workflow to develop the CDA, identify and map land surface change, and verify GECs is summarized in Figure 1. As described below, we first compiled the three geospatial datasets for land surface reflectance, elevation, and surface water.

Change in Surface Reflectance
The first dataset was based on the normalized difference vegetation index (NDVI), which can be derived from Landsat (30 m) imagery and is commonly used as a proxy for vegetation productivity in boreal and tundra environments [42]. Higher NDVI values reflect more productive vegetation; thus, decreased NDVI can indicate the decline or loss of vegetation. Because GEC formation causes the rapid replacement of surface vegetation with water and occurs across areas that are typically larger than a Landsat pixel, we hypothesized that changes in NDVI would be sensitive to GEC formation. Collecting imagery approximately every two weeks since 1984, the Landsat 5-8 campaigns provide a multi-decadal record of surface reflectance that is well-suited to identifying abrupt land surface change associated with GEC formation. To test this, we used the LandTrendr algorithm in Google Earth Engine, which ingests a series of user-defined algorithm parameters and determines trends and identifies break-points in surface reflectance within Landsat time-series stacks [41,43]. LandTrendr considers the entire Landsat time-series when identifying breakpoint(s). Landsat surface reflectance data made available in Google Earth Engine are atmospherically corrected using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) (Landsat 5-7) or the Land Surface Reflectance Code (LaSRC) (Landsat 8). We iteratively tested LandTrendr input parameters to optimize the detection of rapidly forming (GEC-like) permafrost thaw features. The full list of input parameters is presented in Table A1. Because any given Landsat pixel's timeseries may contain more than one change event, the user must specify in LandTrendr the type of change to query. Briefly, our final set of input parameters queried the greatest decrease in NDVI during the summer thaw season (May 1-September 30), when GEC formation is likely to be more common [16], from 1984 to 2018. Clouds, shadows, and snow were masked from the Landsat imagery based on CFMask from the United States Geological Survey [44].

Change in Elevation
To detect abrupt elevation change-a hallmark of GEC formation (e.g., [16])-we used the ArcticDEM, a high-resolution (2 m) pan-Arctic digital surface model generated by the Polar Geospatial Center and freely available in Google Earth Engine ( Figure 1) [13]. We used ArcticDEM strip data, comprising a time series of stereo-imagery collected from

Change in Surface Reflectance
The first dataset was based on the normalized difference vegetation index (NDVI), which can be derived from Landsat (30 m) imagery and is commonly used as a proxy for vegetation productivity in boreal and tundra environments [42]. Higher NDVI values reflect more productive vegetation; thus, decreased NDVI can indicate the decline or loss of vegetation. Because GEC formation causes the rapid replacement of surface vegetation with water and occurs across areas that are typically larger than a Landsat pixel, we hypothesized that changes in NDVI would be sensitive to GEC formation. Collecting imagery approximately every two weeks since 1984, the Landsat 5-8 campaigns provide a multidecadal record of surface reflectance that is well-suited to identifying abrupt land surface change associated with GEC formation. To test this, we used the LandTrendr algorithm in Google Earth Engine, which ingests a series of user-defined algorithm parameters and determines trends and identifies break-points in surface reflectance within Landsat timeseries stacks [41,43]. LandTrendr considers the entire Landsat time-series when identifying breakpoint(s). Landsat surface reflectance data made available in Google Earth Engine are atmospherically corrected using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) (Landsat 5-7) or the Land Surface Reflectance Code (LaSRC) (Landsat 8). We iteratively tested LandTrendr input parameters to optimize the detection of rapidly forming (GEC-like) permafrost thaw features. The full list of input parameters is presented in Table A1. Because any given Landsat pixel's time-series may contain more than one change event, the user must specify in LandTrendr the type of change to query. Briefly, our final set of input parameters queried the greatest decrease in NDVI during the summer thaw season (1 May-September 30), when GEC formation is likely to be more common [16], from 1984 to 2018. Clouds, shadows, and snow were masked from the Landsat imagery based on CFMask from the United States Geological Survey [44].

Change in Elevation
To detect abrupt elevation change-a hallmark of GEC formation (e.g., [16])-we used the ArcticDEM, a high-resolution (2 m) pan-Arctic digital surface model generated by the Polar Geospatial Center and freely available in Google Earth Engine ( Figure 1) [13]. We used ArcticDEM strip data, comprising a time series of stereo-imagery collected from 2008 to 2017, to quantify elevation change over time. We recognize that elevation change could have occurred prior to 2008. Therefore, because our change detection product accounts for changes in surface reflectance and water from 1984 to 2017, but not elevation over this entire time span, our change product is likely a conservative estimate of change. Sensor model errors can result in geolocation offsets between individual ArcticDEM strips on the order of several meters [14], which must first be corrected for assessment of change across broad scales (>1000 km 2 ). We first corrected for elevation artifacts in the ArcticDEM strip data, which can result in overlapping strips being anchored to different elevation values and potentially result in erroneous change values ( Figure A1). Artifacts were corrected using offset values from the Ice, Cloud and land Elevation satellite (ICESat) LiDAR altimetry data available in the ArcticDEM strip metadata. Of the 5783 ArcticDEM strips in our study area, 4523 had metadata information that allowed us to correct the values to ICESat data (Appendix A). We filtered elevation change to within ±20 m, the typical range of geomorphic activity associated with GEC formation [16], to omit change associated with variation in fog and/or cloud cover over time. Using the corrected and filtered ArcticDEM strips, we calculated change as the maximum difference in elevation change for each 30 m pixel in our study area.

Change in Surface Water
Precipitation, surface water, and thawed ground ice rapidly fill newly formed GECs, creating a lake within~1 to 2 years of formation [16]. To limit change in surface water to within our study period, surface water that was present both in 2018 and before the study period was omitted from the analysis. This may have resulted in the omission of change associated with GECs which formed within drained lake beds (e.g., [36]) and GECs which formed prior to 1984. Surface water presence during the study period was obtained from a dataset of global surface water (30 m resolution) derived from Landsat data [12].

Mapping Hotspots of Land Surface Change
Geospatial data for changes in NDVI (1984-2017), elevation (2008-2017), and surface water (1984-2017) were combined by scaling each metric from 0 to 1, assuming these three metrics of land surface change carried equal weight for detecting GECs, and then adding them together for every 30 m pixel. Maximum possible change, therefore, for any given pixel would have a value of 3. While elevation change was computed in the native resolution of the ArcticDEM (2 m), we calculated the mean elevation change over 30 m to combine the three datasets at the same spatial resolution. Hereafter, land surface change refers to change within the context of the CDA.
Additional filters were applied to the change product to mask out mountainous terrain and human infrastructure. Topography-driven landcover masking was used to omit parts of the change product that was in mountainous terrain, where avalanches or other forms of mass-wasting may induce large scale elevational or landcover changes unrelated to GEC activity. This mask was defined by thresholds of elevation (>150 m) and slope (≥2 • ), and CAVM category (mountain complexes, both carbonate and non-carbonate). Human development in the study area, which includes mines, roads, pipelines, oil and gas platforms and machinery, and urban infrastructure, was expert-identified and digitized by hand using ArcGIS Pro, and masked out of the change product. Open Street Map data were used to enhance and validate the human infrastructure mask. To limit potential effects of human infrastructure on natural change in surface reflectance, we applied a buffer of 100 m to digitized infrastructure.

Validation of the Land Surface Change Map
In order to validate CDA identified GECs versus other forms of change (e.g., thermokarst such as RTSs), we developed a three-stage validation protocol based on expert knowledge. First, initial validation efforts explored six cells around GEC-1, near the Bovanenkovo gas field. Additional validation cells were developed by dividing the study domain into 604 cells (~500 km 2 /cell). Based on the majority CAVM landcover type within each cell, we used random stratified sampling to select an additional 31 cells for validation. Next, a team of GEC 'crater hunters', trained by permafrost experts, independently manually interpreted high-resolution (sub-meter) satellite imagery available for the entire study region (Maxar, www.maxar.com). Manual interpretation included labeling features which appeared to be GECs within the subset of validation cells. RTSs were also noted and differentiated from other cryogenic landslide features (e.g., active layer detachments or ground flows; [22]) on the basis of a near-vertical headwall of exposed permafrost, readily identified by a shadow cast onto the RTS floor. Landcover data were obtained from the CAVM [27]. Crater hunters were assigned the randomly selected cells, which contained a grid of smaller cells (~5 × 5 km 2 ) for visual interpretation. All crater hunters independently interpreted the same two larger cells, containing known GECs, to ensure rigor and accuracy of the validation process. Second, features of interest identified by the crater hunters were verified by the permafrost experts. Third, for features resembling GECs, we interpreted high-resolution, multi-temporal imagery to confirm the feature type. GECs were positively identified when features showed: (i) a pingo-like mound prior to eruption; (ii) a debrisencircled crater containing turbid water; and (iii) growth of a turbid lake one or more years post-eruption [20]. Additionally, spectral and morphological characteristics of potential GECs were extracted from remote sensing data and compared with GECs previously identified in the literature (Section 3.3), to support our hypotheses about the origin of the GEC-like features. Due to the low number of GECs reported in the literature (n = 7), GEC validation was not evaluated using accuracy statistics [45] or omission/commission errors.

Statistics
To assess variability in landscape change by ecosystem and land cover type, we report change summary statistics for the categories in CAVM [27]. The summary statistics were derived using Google Earth Engine. To test whether change pixel values associated with GECs were significantly different from the median of non-GEC change pixels, we used a one-sample Wilcoxon Signed Rank Test. The Wilcoxon test was performed using the stats package in R software v.3.6.2 [46]. We assessed long-term trends  in mean annual and mean summertime air temperature for the city of Salekhard, to place land surface change within the context of regional climate change. Long-term climate trends were evaluated using the non-parametric Mann-Kendall test from the R software package zyp [47], following the trend pre-whitening approach to account for serial autocorrelation [48].

Land Cover Types and Change
Across the 327,000 km 2 study area, land cover was primarily classified as upland/nonwetland tundra (81%), with lower amounts of wetlands (11%), water (7%), and mountains/barren terrain (1%) (Table 1, Figure A2). Approximately 5% (16,600 km 2 ) of the land surface showed changes in vegetation, elevation, and/or water extent (Figures 1 and 2). Land surface change was generally greater along the shorelines of water bodies, along coasts and in northern regions, and in the Ural Mountains in the southwest of the study area ( Figure 2). Land surface change occurred among 14 categories of vegetation of the CAVM, which could be generalized into four land cover categories ( Table 1). The total area of change was greatest in upland/non-wetland tundra (11,600 km 2 ), but this category showed the lowest amount of change relative to total area (4%). Compared to upland/nonwetland tundra, proportional change was roughly two times higher in wetland complexes, areas of surface water, and mountains/barren terrain (Table 1). We recognize that masking some mountainous terrain (Section 2.2.2) likely resulted in a conservative estimate of proportional change. Change was often relatively more intense (higher change pixel values) for upland thermokarst (Figure 3a,b), thermo-erosional features along lake shorelines (Figure 3c,d), exposed littoral zones (Figure 3e,f), and, in some regions, for headwater stream networks (Figure 3g,h).     Pixel change values ranged from 0.10 to 1.62 and were generally low (median = 0.31, IQR = 0.29-0.33) relative to the maximum potential change (3.0). The distribution of values was right-skewed, indicating a smaller proportion of high-magnitude change values ( Figure 4). Generally, previously identified GECs were associated with moderate values of land surface change (median = 0.78, n = 7). Values of land surface change associated with GECs were significantly greater than non-GEC change (p < 0.01, one-sample Wilcoxon Signed Rank Test). All GECs were located within upland/non-wetland tundra except for SeYkhGEC, which was located within a river in a wetland complex. Three additional GEClike change detection features (CDF-1, CDF-2, CDF-3), which have not been previously reported, were detected during validation of the CDA ( Figure 5). These three features were located in the central Yamal Peninsula, two of which were located approximately 40 km south of the Bovanenkovo gas field, in relatively close proximity to GEC-1, GEC-2, and GEC-3. Across the study area, RTSs were also prevalent. Within the 28,300 km 2 of validated terrain, we found 150 terrain disturbance features that were likely to be RTSs. RTSs were present within 500 m of GEC-1, GEC-2, CDF-2 and not as close to the other GECs. AntGEc, SeYkhGEC, YeniGEC, and ErkutaGEC were located within~400 m of terrain that appeared to be polygonal tundra.

Validation of Land Surface Change
We validated 37 fishnet cells (6.1% of total), accounting for 28,300 km 2 within the study area ( Figure A3). The randomly selected fishnet cells represented 14 of the CAVM land cover categories (Table A2). Validation efforts revealed that all seven of the GECs previously reported in the literature ( Table 2) were detected using the CDA ( Figure A4). In addition to detecting land surface change associated with GEC formation, validation efforts further revealed that the CDA detected change that appeared to be due to declines

Validation of Land Surface Change
We validated 37 fishnet cells (6.1% of total), accounting for 28,300 km 2 within the study area ( Figure A3). The randomly selected fishnet cells represented 14 of the CAVM land cover categories (Table A2). Validation efforts revealed that all seven of the GECs previously reported in the literature ( Table 2) were detected using the CDA ( Figure A4). In addition to detecting land surface change associated with GEC formation, validation efforts further revealed that the CDA detected change that appeared to be due to declines

Validation of Land Surface Change
We validated 37 fishnet cells (6.1% of total), accounting for 28,300 km 2 within the study area ( Figure A3). The randomly selected fishnet cells represented 14 of the CAVM land cover categories (Table A2). Validation efforts revealed that all seven of the GECs previously reported in the literature ( Table 2) were detected using the CDA ( Figure A4). In addition to detecting land surface change associated with GEC formation, validation efforts further revealed that the CDA detected change that appeared to be due to declines in NDVI, loss or gain of surface water extent, and non-GEC geomorphic processes (e.g., lake shoreline erosion, RTS activity). Additionally, validation efforts identified 23 lakes having a visual appearance (e.g., morphometry, turbidity) similar to that of previously identified GECs. While 19 of these lakes were associated with change pixels, inspection of high-resolution satellite imagery archives found that none of the lakes were likely to be GECs. Table 2. Characteristics of gas emission craters (GECs) on the Yamal and Gydan peninsulas. CDF = change detection feature. Latitude (Lat) and longitude (Long) are reported in decimal degrees. ∆ = surface change value derived from the change detection algorithm (CDA). Approximate formation dates reported as month, when known, or a range of months/years. * No change pixels in area due to omission by water mask. 1 [24], 2 [25], 3 [17], 4 [16], 5 This study.

Characteristics of GEC Formation from the Change Detection Algorithm and High-Resolution Satellite Imagery
High-resolution satellite imagery together with elevation data from the ArcticDEM strips showed the geomorphic development of GEC-2 prior to (Figure 6a-h) and after its detection by the CDA (Figure 6i,j). CDA detection of GEC-2 was approximately coincident with its formation in late 2012. The CDA demarcated a cluster of pixels with a high magnitude of change centered on GEC-2 and the peripheral ejecta. Subsequent degradation of the crater parapet and apparent infilling of the GEC lake by precipitation, surface water, ground ice, and soil sediments resulted in the formation of a turbid, sub-circular lakẽ 0.1 km 2 in area by August 2020. At least one dozen small ponds (2-10 m diameter) around the GEC lake periphery were visible in satellite imagery from August 2020. NDVI generally increased from 1985 until 2012. An abrupt decrease in NDVI showed as a distinct breakpoint in the time-series, signaling the formation of GEC-2 ( Figure 6b). After the formation of GEC-2, NDVI values decreased until~2019. Similarly, elevation values from the ArcticDEM decrease by~5 m after GEC-2 formed, and decrease thereafter until 2018.  Image on the bottom right shows a cluster of change pixels superimposed on the lake formed by GEC-2 (i). Chart at the bottom (j) shows the original trends and fitted trends (by LandTrendr) in Landsat-derived NDVI (scaled from −1000 to 1000), the breakpoint identified by LandTrendr, and corresponding trends in elevation extracted from the corrected ArcticDEM strips. Imagery © 2020, DigitalGlobe.
The landscape change features CDF-1, CDF-2, and CDF-3 have not been previously reported in the literature. CDF-1 was located 53 km south of SeYkhGEC. Satellite imagery of CDF-1 showed an elongate lake adorned by what appeared to be a sandy parapet encircling a lake with relatively low turbidity. There was no high-resolution, multi-year satellite imagery of CDF-1 with which to characterize its geomorphic evolution. CDF-2 was located 10.8 km due north of GEC-1 and was situated within the edge of a partially drained lake basin (Lake Khalevto), within 300 m of RTSs and polygonal tundra. Satellite imagery of CDF-2 showed a pingo-like feature in 2011, which had begun to degrade by August 2015 (Figure 7a,b). ArcticDEM elevation data indicate that the collapse of the feature apex occurred from ca. 2015 to 2017, resulting in a GEC-like crater by July 2018 with a diameter of 12-15 m (Figure 7c,d; Table 2). There was no evidence from the high-resolution satellite imagery of ejected debris around the periphery of the feature, nor signs of newly formed ponds associated with the thawing of ice-rich ejecta. Degradation of the feature parapet and infilling of the crater formed a~600 m 2 circular lake by July 2020 and the development of secondary thermokarst adjacent to the lake. Unlike GEC-2, NDVI was relatively constant from 1984 to 2002, increased from 2002 to 2015, and then decreased in conjunction with the collapse of the feature and formation of a lake (Figure 7e).

Multi-Decadal Land Surface Change on the Yamal and Gydan Peninsulas
The CDA revealed widespread land surface change across the Yamal-Gydan region in recent decades. Change was primarily associated with coastal regions-which is occurring at rates of 0.2-1.7 m y −1 in the Yamal-Gydan region [49]-surface water (e.g., lakes, streams, rivers), and geomorphically active terrain (e.g., thermokarst), which is not unexpected for lake-rich terrains hosting relatively ice-rich permafrost [5]. Change along lake shorelines is taken to reflect the expansion and/or drainage of thermokarst lakes, and The development and morphometry of the feature CDF-3 appeared to be similar to that of CDF-2. CDF-3 was located 14.3 km northeast of GEC-1, within shrub-tundra. Satellite imagery showed the feature resembled a pingo in July 2010 and had degraded into a crater~20-22 m in diameter by July 2013. Similar to CDF-2, there was no evidence for the explosive ejection of ground material associated with the formation of CDF-3. Imagery from May 2020 showed no evidence of degradation of the feature or growth of the lake within the crater.

Multi-Decadal Land Surface Change on the Yamal and Gydan Peninsulas
The CDA revealed widespread land surface change across the Yamal-Gydan region in recent decades. Change was primarily associated with coastal regions-which is occurring at rates of 0.2-1.7 m y −1 in the Yamal-Gydan region [49]-surface water (e.g., lakes, streams, rivers), and geomorphically active terrain (e.g., thermokarst), which is not unexpected for lake-rich terrains hosting relatively ice-rich permafrost [5]. Change along lake shorelines is taken to reflect the expansion and/or drainage of thermokarst lakes, and change within fluvial networks may indicate changes in snow persistence or vegetation loss associated with thawing and subsiding riparian zones. These land surface changes occurred in tandem with an increase in mean summer (June-September) air temperature by 0.4 • C per decade from 1984 to 2017 (p < 0.05, Mann-Kendall test) in Salekhard, in the southern part of the study area, and warming of permafrost on the Yamal Peninsula (ca. +1 • C per decade, 2007-2016) [4]. Together, these findings indicate that geomorphic and hydrologic activity in the Yamal and Gydan peninsulas are driving widespread land surface change with largely unknown implications for regional ecosystems, biogeochemical cycles, and hazards to human populations and infrastructure.

Speculation on the Origins of the Three New GEC-Like Features
Of the three features, CDF-2 appeared most similar to a GEC upon first inspection of high-resolution satellite imagery acquired in July 2018 (Figure 7). Examination of historical high-resolution satellite imagery showed a mound predecessor in July 2011, within the margin of drained Lake Khalevto. Annual time-series of Landsat imagery in Google Earth Timelapse (https://earthengine.google.com/timelapse/) showed a minor decrease in water levels in Lake Khalevto in 2003, which was followed by a period of stability for five years and then lake drainage in the period 2008-2009. This aligns with NDVI trends, which started to increase around 2003 (Figure 7), perhaps reflecting the emergence and/or growth of vegetation following subaerial exposure of the littoral zone. Drained lake margins often contain the requisite geomorphic and hydrologic conditions for the formation of hydrostatic (i.e., closed system) pingos [50]. Within drained lakes, pingos typically form within residual ponds via the downward aggradation of permafrost, causing forceful expulsion of talik porewater, formation of intrusive ice, and upward heaving of the land surface [51,52]. In early stages, hydrostatic pingos are characterized by rapid early growth rates up to~1.5 m y −1 [50] and, in following decades, on the order of cm y −1 [53]. Pingo degradation may occur over decades and show dilation (i.e., tension) cracks across the dome surface [50,54], though geomorphic characteristics of pingo collapse vary across the northern permafrost zone [51]. We presume that growth of the CDF-2 mound predecessor initiated between 2003 and 2008, when Lake Khalevto water levels began to decline. Assuming the height of CDF-2 reached~3 to 4 m before collapse initiated (Figure 7), growth aligned with theoretical rates typical of early pingo formation. Thus, development of the CDF-2 mound predecessor was characteristic of both pingos and GECs. Curiously, Google Earth Timelapse reveals that CDF-2 did not originate within a residual pond. Therefore, CDF-2 may have originated as a hydraulic (open-system) pingo, provided a sufficient water source [50]. Relatively abrupt collapse of the CDF-2 dome is indicated by high-resolution satellite imagery collected in August 2015, and July 2018, resulting in the formation of a central pond and funnel-like upper crater walls (Figure 7).
We presume that the abrupt loss of soil mechanical strength and subsequent dome collapse was coupled to the loss of internal pressure within CDF-2 via non-eruptive degassing of methane or loss of subsurface water. Subsequent thermokarst and exposure of permafrost would accelerate thaw, resulting in the turbid lake. From the available data, we can neither confirm that CDF-2 was a pingo, nor reject the hypothesis that it was a GEC that did not erupt explosively.
Google Earth Timelapse revealed that, unlike CDF-2, CDF-1 and CDF-3 were not located within drained lake or pond basins. While limited elevation data and historical highresolution satellite imagery exist for CDF-1, its elongate nature, relatively low turbidity, and sand-adorned parapet were not characteristic of GECs [16,24,36,55]. We surmise it was a degraded pingo, or perhaps resulted via thermokarst lake formation. In contrast, high-resolution satellite imagery of CDF-3 showed geomorphic evolution more similar to that of CDF-2. Without additional information, we reason that CDF-3 could have formed from a collapsing pingo induced by the escape of subterranean water or thaw of ice-rich permafrost [56], or as a GEC that released its methane non-explosively. Why CDF-3 did not appear to show much, if any, degradation between 2013 and 2020 remains unknown. Relatively low change pixel values for CDF-3 and absence of breakpoints in its NDVI timeseries are likely because the feature was smaller than the resolution of a Landsat pixel (30 m). These results highlight the limits of using remote-sensing to ascribe forms of geomorphic change in permafrost-affected terrains. The CDA detected previously undocumented GEC-like features and results suggest that not all GECs may erupt explosively. Validating land surface change features in the field is a top priority to advance understanding of these phenomena.

Performance of the Change Detection Algorithm
The CDA detected all seven of the GECs previously reported in the literature [16,24,36,55] and three new GEC-like features. Therefore, despite uncertainty whether the latter were GECs, these results demonstrate that our semi-automated CDA successfully detected land surface change associated with GECs. Importantly, however, GECs had already formed prior to detection by the CDA [16]. Detection of GECs after their formation is largely due to the CDA inputs (Figure 1), which weight breakpoints associated with abrupt changes in NDVI [41]. Such changes are most likely to occur after GECs form. This was reflected by relatively large clusters of high-value change pixels centered on lakes that formed within year(s) following GEC formation, such as for GEC-1, GEC-2 ( Figure 6), and YeniGEC. Further, the CDA product simply represents net land surface change during the 1984-2017 period. The CDA does not as readily capture more subtle multi-year terrain uplift associated with mound predecessor formation [16]. For these reasons, and the laborious validation required for the CDA (Figure 1), our semi-automated approach represents a first step towards the automated detection of GECs. Further refinement of the algorithm is required to better isolate change associated with GECs from other types of land surface change, which would reduce validation time. On the other hand, although the CDA input parameters (Figure 1) were tuned for the detection of GECs, our results demonstrate that the CDA was sensitive to diverse types of land surface change that occur in the Yamal-Gydan region [19,22] and also in permafrost regions elsewhere. Thus, our CDA workflow ( Figure 1) could be further refined to detect and monitor RTSs [57], lake or coastal shoreline erosion [58], and/or lake drainage across the circumpolar north [59].

Detection, Mapping, and Monitoring Northern Environmental Change
Geomorphic, hydrologic, and vegetation change across the northern permafrost zone encompasses a suite of land surface processes, including GECs, that can be expected to intensify with future warming [60][61][62]. Given the potential hazards to humans and infrastructure posed by GECs [17], detecting mound predecessors and predicting GEC formation in the Yamal-Gydan region is of top interest. Yet, the sparsity of known GECs hinders automated and remote detection of mound predecessors, and also the develop-ment of an approach to predict GEC formation [16]. Our results indicate that elevation change, a hallmark of mound predecessor formation, can be detected across broad scales using ICESat-corrected ArcticDEM data. Remote sensing campaigns which make repeat measurements of elevation at annual or sub-annual timescales could detect and monitor future regional land surface change associated with thermokarst activity [63]. For instance, the ESA Sentinel-1 and Sentinel-2 platforms (launched in 2014/2015) collect radar and multispectral imagery up to 84 • N every 10 to 12 days at decameter spatial resolution (https://sentinel.esa.int). Sentinel-1 interferometric synthetic aperture radar (InSAR) data capture sub-decimeter vertical resolution of ground surface displacement [64] and could potentially enable detection and monitoring of GEC formation. Developing Sentinel databased methods was beyond the scope of this work, which sought to identify trends in land surface change prior to the launch of the Sentinel platforms. Identifying changes in surface water extent using multispectral satellite imagery [12] could allow for detection and monitoring of areas considered conducive to GEC formation, such as drained thermokarst lake basins [55]. Other approaches for detecting GECs may be hindered by the limited number of documented GECs [16]. For instance, machine learning approaches show promise for mapping thermokarst based on the distinct spectral texture of thermokarst features [65], but require large amounts of input data for calibration and validation. In any case, feature detection is limited by the spatial resolution of imagery. Increasing accessibility to large amounts of remote sensing data and cloud-computing (e.g., Google Earth Engine) will push forward new avenues of research on land surface change at northern high latitudes [66]. Remote sensing methods for detecting and monitoring geomorphic, hydrologic, and vegetation change, which are rapidly evolving, are of utmost importance for understanding the extent, rate, and forms of permafrost thaw, and for mitigating effects to infrastructure [7] and ecosystems [6] across the circumpolar north.

Conclusions
In this study, we developed a semi-automated workflow using satellite remote sensing imagery to characterize land surface change and detect GECs on the Yamal and Gydan peninsulas. We found that 5% of the Yamal-Gydan region experienced changes in vegetation, elevation, and/or surface water from 1984 to 2017. Our change detection algorithm successfully identified the seven GECs previously reported in the scientific literature and three previously undocumented GEC-like features. Two of the newly identified features had geomorphic similarities to-and were located in an area nearby-several known GECs in the vicinity of the Bovanenkovo gas field. While less abundant than retrogressive thaw slumps in the Yamal-Gydan region, GEC activity is challenging to predict and poses a hazard to human populations and infrastructure. Our approach is a first step to detecting and monitoring the formation of GECs, which may be further improved and perhaps automated with additional observations of GECs. Refining methods to detect geomorphic activity will help to mitigate potential hazards and better understand regional variability in terrain responses to climate change across the northern permafrost zone.

Acknowledgments:
We thank Carl Churchill for assistance with delineating human disturbance; Kim Fiske for assistance with validating the GEC map; Paul Morin and Ian Howat (University of Minnesota Polar Geospatial Center) for insightful conversations about the ArcticDEM; and four anonymous reviewers for their constructive feedback, which greatly helped to improve the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Example code for correcting ArcticDEM data to ICESat spaceborne lidar can be found here: https://code.earthengine.google.com/4e31cc160753479c7a25a435b240278a. Table A1. Input parameters used in the LandTrendr analysis. Detailed descriptions of input parameters can be found in the LandTrendr for Google Earth Engine guide (https://emapr.github.io/LT-GEE/api.html) and in the literature [41]. Appendix C  Appendix D Figure A1. Change in elevation across the study area from 2008 to 2017, derived using the Arc-ticDEM digital surface model [13]. Vertical discrepancy among ArcticDEM strips (uncorrected, top panel) was corrected using the ICESat altimetry offset (corrected, bottom panel) (see Section 2.2). Figure A1. Change in elevation across the study area from 2008 to 2017, derived using the ArcticDEM digital surface model [13]. Vertical discrepancy among ArcticDEM strips (uncorrected, top panel) was corrected using the ICESat altimetry offset (corrected, bottom panel) (see Section 2.2).

Parameter Type LandTrendr Var Value
Appendix F Figure A3. Distribution of fishnet cells across the study area (outlined in dark gray). Out of 604 fishnet cells within the study area, 37 (6.1%) were validated, accounting for 28,300 km 2 . Figure A3. Distribution of fishnet cells across the study area (outlined in dark gray). Out of 604 fishnet cells within the study area, 37 (6.1%) were validated, accounting for 28,300 km 2 . Figure A4. High-resolution satellite imagery of the seven gas emission craters (GECs) previously reported in the literature (Table 2), shown without change pixels (left panel) and with change pixels (right panel) for GEC-1 (a,b), GEC-2 (c,d), GEC-3 (e,f), AntGEC (g,h), SeYkhGEC (i,j), YeniGEC (k,l), and ErkutaGEC (m,n). Note: due to the difficulty of obtaining cloudfree, high-resolution imagery during the relatively brief period during which the characteristic crater stage of GEC evolution, most images show the GECs prior to (e.g., ErkutaGEC) or after (e.g., YeniGEC) their formation. Imagery © 2020, DigitalGlobe. Figure A4. High-resolution satellite imagery of the seven gas emission craters (GECs) previously reported in the literature (Table 2), shown without change pixels (left panel) and with change pixels (right panel) for GEC-1 (a,b), GEC-2 (c,d), GEC-3 (e,f), AntGEC (g,h), SeYkhGEC (i,j), YeniGEC (k,l), and ErkutaGEC (m,n). Note: due to the difficulty of obtaining cloud-free, high-resolution imagery during the relatively brief period during which the characteristic crater stage of GEC evolution, most images show the GECs prior to (e.g., ErkutaGEC) or after (e.g., YeniGEC) their formation. Imagery © 2020, DigitalGlobe.