Mapping and understanding the vulnerability of northern peatlands to permafrost thaw at scales relevant to community adaptation planning

To update permafrost peatland vulnerability maps at local scales within the discontinuous permafrost zone of the NWT, the 372 220 km² study area was mapped using Sentinel 2 imagery (Sentinel-2A, B04 (Red), B03 (Green), B02 (Blue), 10 m, 2016 and 2017, July and August). Peatland complex area was visually mapped using a 3.75 × 3.75 km grid cells using the percent cover of permafrost peatland complex area per grid cell. The grid cell classes include none (0% coverage), negligible (<2%), low (3%–25%), moderate (26%–50%), high (51%–75%), and very high (76%–100%). A subset of study area (∼5% of the mapped area) was assessed by two mappers in order to test the guidelines on interpretation and the accuracy of the approach. With respect to percent area estimates, overall mapper accuracy is 89%. The data was compiled and all GIS analysis was completed using ArcGIS (ESRI, 2014, version 10.2.2, Redlands, CA, USA). For complete methods and data see Gibson et al (2020). Given the cyclical lifecycle, and high-ice content properties of peatland permafrost (Zoltai 1993), we consider any intact permafrost peatland to be predisposed to thermokarst formation.

To assess how the degree of thermokarst formation within permafrost peatland varies across the climatic gradient, we visually quantified the proportion of peatland complex area that had undergone thermokarst formation using the ESRI World Imagery map downloaded from ArcGIS.com (World Imagery (arcgis.com)). Mosaiced satellite imagery used in this interpretation were acquired during the growing season (May to September) from 2009 to 2019. All data were provided to the World Imagery archive by Maxar Inc. Data were acquired at varying pixel resolutions between 0.31 and 0.6 m collected using GeoEye-1, WorldView2, WorldView3 and Quickbird-2. Visual estimations are possible due to the distinct vegetation differences between thermokarst bogs and intact permafrost peat plateaus that are clearly discernable on RGB (optical) high resolution satellite imagery (figure 1(b)).

To determine how the proportion of peatland complex that has thawed varied across the study region, 3.75 × 3.75 km grid cells with high and very high estimates of permafrost peatlands were first identified. This ensured that we were comparing similar peatlands areas across a latitudinal gradient and that differences in the amount of thaw are not being driven by other factors that are occurring at fine scales in small peatlands. Furthermore, from a community perspective, it is likely that thawing of large permafrost peatlands will have the most impact on their land use regarding travel and changes in hydrology (Quinton et al 2011a).

A subset of these identified cells where then randomly selected (n = 70, or ∼5%, figure 2). Selected grid cells spanned the entirety of the study area and had high resolution ArcGIS DigitalGlobe, Geo Eye basemap (ArcGIS version 10.3) imagery available. Selected grid cells where then overlain with a 10 × 10 grid (creating 100 sub-grid cells with size equal to 375 × 375 m). Of the 100 sub grid cells, those contained within the peatland complex area were identified visually and ten random cells of those were selected to determine the degree of thermokarst formation (figure 2). In the selected sub grid cells, we visually estimated the percent of the peatland complex with thermokarst formation. Estimates were made in intervals of 10 (i.e. 0%, 10%, 20%, 30% coverage etc). The mean percent of peatland complex area that was thawed and standard deviation for each 3.75 × 3.75 km grid cell was calculated.

Zoom In Zoom Out Reset image size

To assess how the proportion of peatland complex that has thawed varied across a climatic gradient, mean annual air temperature was estimated from WorldClim 2.1 climate model (Fick and Hijmans 2017). This dataset is a grid (resolution = 1 km²) of average monthly temperature interpolated from weather station data (1970–2000). Long time series of historical observations of climate and hydrology are scarce in the NWT, therefore gridded datasets have been used as alternatives to instrumental observations for climate analysis (Persaud et al 2020, Segal et al 2016). The mean annual air temperature from the WorldClim 2.1 climate model was assigned to each grid cell using the zonal statistics tool in ArcGIS. Given the high collinearity between mean annual air temperature and latitude (R² = 0.95) and the greater certainty in latitude compared to mean annual air temperature, latitude was used in all subsequent analyses.

The proportion of thawed peatland complex area was linearly regressed against latitude. The scatterplot was visually assessed for trends in thermokarst formation to make inferences about how thaw may progress in a warming climate. This data was also binned into three groups (59.9–62° N, 62–64.1° N, and 64.1–66.1° N) and statistically tested for differences in the proportion of peatland complex thawed using an ANOVA. The binned data was also statistically tested for differences in variance using a Fligner–Killeen test of homogeneity of variances for non-normally distributed data (Williams et al 1981).

To assess the potential for elevational controls on thermokarst formation, elevation data was derived from the 0.75 arcsec (20 m) Canadian Digital Elevation Model (CDEM—Natural Resources Canada 2015). Individual CDEM tiles were mosaiced to the extent of the study area and mean elevation was assigned to each grid cell using the zonal statistic tool in ArcGIS. Visual assessments and interpretations were made to determine how elevation influenced the proportion of thermokarst formation within an individual bin and between bins. This was assessed using regression analysis with a model of the proportion of thermokarst with the main effects of latitude and elevation.

The mapping of a 372 220 km² area of northwestern Canada confirmed widespread coverage of permafrost peatlands. In total, 53% of the grid cells contained permafrost peatlands (figure 3) with 16% classified as negligible, 21% low, 9% medium, 5% high, and 2% very high cover. The proportion of peatland complex already containing thaw ranged from 3 ± 3% to 77 ± 12% within the study area. The proportion of peatland complex thawed also varied along a latitudinal gradient with greater proportions of thermokarst formation in peatland complexes in the south compared to the north (figure 4(a)). The proportion of peatland complex thawed varied between latitudinal bins (figure 4(a) inset; ANOVA, p < 0.001, F = 73.89). Additionally, there was greater variability in the proportion of peatland complex thawed in southern permafrost peatlands compared to northern ones (Fligner–Killeen test of homogeneity of variances, p < 0.001, figure 4(a) inset). This was also apparent from increased scatter around the regression line.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Elevation appears to govern the proportion of thermokarst formation within peatland complexes at lower latitudes (figure 4(b)). Our results suggest that latitude and elevation together improved prediction of thermokarst formation compared to latitude or mean annual air temperature alone (p < 0.001). Furthermore, the proportion of thermokarst formation was significantly correlated to elevation for the lower (linear model, p = 0.004) and mid latitude (linear model, p = 0.006) bins, but not for the high latitude bin (linear model, p = 0.8). We infer that the increased variance in the proportion of thermokarst formation at lower latitudes is being driven by this elevation effect, whereby higher elevation peatlands remain protected from increasing temperatures.

One of the driving forces for this study was to create permafrost thaw vulnerability maps at a scale appropriate for community-use, as coarse circumpolar scale assessments can lead to feelings of eco-anxiety (Cunsolo and Ellis 2018). If thaw probability maps are coarse in scale and large areas of community territories are indicated to be 'at high risk' of abrupt thaw, it may contribute to feelings of hopelessness in the face of climate change (Cunsolo and Ellis 2018). The spatially explicit nature of our permafrost thaw dataset provides a more manageable perspective of risk and allow more effective identification of 'hot spots', and conversely also 'cold spot' areas that are deemed less vulnerable to landscape change in the face of warming and permafrost thaw. Additionally, the permafrost peatland and thermokarst dataset described here helps to address the subjective concept of permafrost risk (Aven and Renn 2010) in which community members are likely to carry their own risk narratives (including past experiences with permafrost thaw) and apply it to any mapping product (Flynn et al 2019).

Previous research used circumpolar data to quantify how much lowland organic-rich permafrost (mainly peatlands) is predisposed to thermokarst, including peatlands that may have already experienced thermokarst (Olefeldt et al 2016). Because of its coarse scale, Olefeldt et al (2016) caution against regional or local applications of the product. Rather than relying on proxies known to be important in triggering thermokarst, such as topography or ground ice information (Olefeldt et al 2016), here we directly quantified the area of peatland complexes and, in a subset of grid cells, the proportion of thermokarst. Because peatlands form in particular topographic and ground ice settings, our mapping product has strong correlations with the input data used in Olefeldt et al (2016) coarser product, despite the differences between the two studies. The Olefeldt et al (2016) classifies 'High' vulnerability as 30%–60% coverage (equivalent to this study's 'moderate' category) and assumes 50% of their 'very high' coverage to be permafrost-free given it is located in the discontinuous permafrost zone. Though these two mapping products have different approaches and vary in how the thermokarst predisposition classes were binned and classified, we were interested in comparing the two products in terms of how the area of land deemed most vulnerable to thermokarst varied (i.e. the difference in area of 'high' and 'very high' categories between the two products).

Using the circumpolar mapping approach outlined in Olefeldt et al (2016) applied to our study region, the proportion of total land associated with 'high' or 'very high' predisposition to lowland/peatland thermokarst was 61%. Our results indicate that only 6% of this same region has 'high' or 'very high' predisposition to peatland thermokarst based on our mapping of permafrost peatland complexes (figure 5). This reduction in the area classified as being highly predisposed to peatland thermokarst possibly is attributed to the coarse scale of input data utilized by Olefeldt et al (2016). The scale of circumpolar data could be leading to an overestimation of peatland area or an overestimation of the potential for peatland thermokarst. Here we assume that all permafrost peatlands are vulnerable to thermokarst. The large differences between our results and the Olefeldt et al (2016) approach is likely due more to differences in classifying lowland organic-dominated terrain (Olefeldt et al approach) versus actual peatland area (our approach). We conclude that the Olefeldt et al (2016) approach is suitable for general applications that need non-spatially explicit information on thermokarst potential, while this study's approach is appropriate for regional or fine-scale applications that require a spatially explicit understanding of where thermokarst potential areas truly are.

Zoom In Zoom Out Reset image size

For purposes of community use and planning, the mapping products described in this study provide spatially explicit empirically-based estimates of permafrost peatland area and thus a more direct assessment of thermokarst potential. Given that proactive planning of infrastructure and land use based on the potential for permafrost thaw is typically more cost-effective than retrofitting (Melvin et al 2017), appropriately-scaled thaw vulnerability maps are critical for supporting landscape planning, cumulative effects assessments within environmental assessments, and for environmental management through range planning. Our assumption that all permafrost peatland area is predisposed to future thermokarst is reasonable for permafrost thaw in boreal and subarctic peatlands but may not be appropriate for other types of thermokarst occurring in Arctic tundra. This means that other approaches may be required for creating community-relevant thaw vulnerability maps in regions where hillslope or lake thermokarst is likely to be a dominant type of landscape change.

We note that there are challenges in using the information from this study portrayed in figure 3 to identify future 'hazard potential' from thawing permafrost. A key limitation of this current dataset is that it assesses permafrost peatland complex area and does not identify the amount of thermokarst, or 'thaw to date' within the complex. 'High' classified grid cells in the north versus south can contain similar amounts of permafrost peatlands; however, peatlands located further south in our region have already experienced a significant amount of thaw (figure 4). Understanding future 'thaw hazard' will require a consideration of the period in which 'hazard' is being assessed. Southern peatlands, though they have less permafrost left to thaw (figure 4), are going to lose the remaining permafrost in much shorter timeframes than northern peatlands. Therefore, hazard assessments with shorter time frames (years to decades) would suggest higher hazard in the south, where thawing is expected to be most pronounced over this short time period. On the other hand, northern peatlands have significantly larger areas of intact permafrost remaining to thaw (figure 4), suggesting that the total magnitude of the hazard will be greater over the long term. Thus, for hazard assessments with a longer time frame (decades to next century), the potential for future thaw or thaw-hazards may be more pronounced in northern peatlands. For a more detailed description of the dataset's limitations see Gibson et al (2020).

The grid-based mapping approach used in this study allowed for improved spatial resolution and continuous coverage while balancing the time required to analyse and interpret the satellite imagery. Rather than 'mapping' with points, lines and polygons, grid-based mapping allowed us to effectively record the locations of permafrost peatlands and identify high density area of peatlands. Grid-based mapping provides an efficient solution to the problems of mapping small landforms over large areas, by providing a consistent and standardized approach to spatial data collection (e.g. Segal et al 2016). The simplicity of the grid-based mapping approach makes it extremely scalable and workable for group efforts, requiring minimal user experience and producing consistent and repeatable results (Ramsdale et al 2017). Although the grid-based assessment cannot identify the specific locations that thaw has or is likely to occur within the permafrost peatlands (i.e. does not identify specific thermokarst locations), what it does provide is a higher order assessment of where vulnerable areas are located. These vulnerable areas can then be assessed against important and traditional areas of communities to help direct and inform where finer scale studies and efforts should be applied (Andrews et al 2016). This approach is feasible because thermokarst in peatlands leads to surface changes that are easily detected. In other situations (such as active layer thickening), permafrost thaw may not be easily detected from surficial changes.

The latitudinal effect on the proportion of thaw in permafrost peatland complexes illustrates the potential for continued widespread thawing across the discontinuous permafrost zone of the NWT. A near 70% difference in the proportion of thermokarst area in our study region occurs across a ∼3°C–4°C difference in air surface temperatures (Fick and Hijmans 2017). In response to these air temperatures, mean annual ground temperatures range from >0 °C in the southern portion to −2°C in the northern portion of our study region (Smith et al 2010). Mean annual air temperatures in our study area are expected to rise by 3°C by 2100 (IPCC 2018). Altogether, this suggests that the ground temperatures and thaw-extent in the northern extent of the study area in the future will be similar to those currently observed in the southern extent of the study area today. Jorgenson et al (2020) concluded that permafrost thaw in interior Alaska, with mean air temperatures of −2.4 °C, has already reached a tipping point with irreversible thaw. If this is true in the NWT, we speculate that by 2100 as much as 70% of northern permafrost peatlands will have thawed permanently within the study area.

We were surprised by the linear relationship between latitude and proportion of thermokarst. We had predicted there would be nonlinear evidence of abrupt ecological change, defined as a substantial change in ecosystem states over relatively short periods of time when compared to typical rates of change (Ratajczak et al 2018). Abrupt ecological changes are increasing being reported in nature due a warming climate and spans diverse ecosystems and scales (Bestelmeyer et al 2011, Cloern et al 2015, Rocha et al 2015, Thomson et al 2015, Westerling 2016). In this study, it was expected there would be a nonlinear response in the proportion of thermokarst area across the latitudinal gradient as Baltzer et al (2014) showed an exponential increase in the rate of thaw (i.e. percent plateau loss per year) with climate warming using a time series analysis. Although this study shows no evidence of nonlinearity in the proportion of thermokarst area within permafrost peatlands, it could be occurring within discrete ecological areas as opposed to the larger ecological gradient we used in our space-for-time substitution. If so, additional work is needed using time series analysis to test for non-linearities within permafrost peatlands.

Our latitudinal gradient and opportunity to think about a space-for-time substitution allowed us to speculate about how northern climate-protected permafrost peatlands may be impacted by future warming. We do, however, acknowledge that our approach assumes that northern permafrost peatlands will respond to climate warming in the same way southern permafrost peatlands have. This will be complicated by complex interconnected controls on thermokarst formation including but not limited to permafrost thickness, subsurface condition, drainage and more (e.g. Quinton et al 2011b, Quinton and Baltzer 2013, Baltzer et al 2014). However, despite these cautions and caveats, space-for-time approaches are commonly used as one approach for providing insights into potential ecosystem changes associated with climate change (Dieleman et al 2020, Pokrovsky et al 2020).

One of the key challenges for predicting how and when permafrost peatlands will respond to warming is that projections often depend on accurate mean annual air temperature data. The NWT is data sparse. Mean annual air temperature data for the study region is only based on two weather reporting stations (Environment Canada 2016). While the resulting interpolation may be sufficient for large scale climate and permafrost modeling, predicting finer scale patterns and processes is difficult. If we are to support community adaptation and planning to changing permafrost conditions, we will require a better understanding of regional differences in mean annual air temperature (i.e. more long-term climate monitoring).

This study also demonstrates the importance of considering fine-scale regional differences in elevation when assessing trends in thermokarst formation. Our results show substantial variance in thaw-extent at lower latitudes. We attribute this to differences in elevation, in which higher elevation peatlands are more protected from increasing temperatures than lower elevational peatlands and are generally far better drained (figure 5). We assumed that permafrost in northern peatlands would be more climate protected than southern peatlands in our study region, but these results suggest that peatland permafrost also can be resistant to change due to high elevation. This finding, coupled with the known mean annual air temperature sparsity within the study area, introduces the need more fine-scale data. We recommend that fine-scale data collection of mean annual air temperatures and mean annual ground temperatures is prioritized in order to make more valid predictions of future permafrost thaw in and around communities.