Climate change and permafrost thaw-induced boreal forest loss in northwestern Canada

Permafrost distribution throughout the western Canadian subarctic is not well understood due to the remoteness and size of the region, its spatial and temporal heterogeneity, limited data availability, and sparse monitoring networks. These factors severely challenge investigations of how climate warming might affect the distribution of permafrost and provide strong justification for new methods of evaluating permafrost extent using remote sensing platforms. This study quantifies forest loss at ten subarctic boreal sites in the southern Northwest Territories and northeastern British Columbia between 1970 and 2010. Historical air photos and optical remote sensing images were assessed using a change detection approach over the ten sites, each 10 km2 spanning a north/south transect of 200 km. This study is the first to apply change detection methods to a large-scale gradient and spans the southern margin of discontinuous permafrost where results demonstrate variable patterns of net forest loss at each site ranged from 6.9% to 11.6% over the 40 year study period. Here we show that these differential rates of landcover change can be explained in part through climatic and environmental factors that vary latitudinally across the selected sites. Change statistics—net change, forest gain and forest loss were significantly correlated with an assortment of factors that varied across the ten-site transect.


Introduction
Up to 80% of the world's boreal forests are located in environments underlain by permafrost (Helbig et al 2016). Rapid climate warming throughout the vast, highly-sensitive boreal region, has led to changes in both geomorphology (Olefeldt et al 2016) and ecology (Baltzer et al 2014) across the landscape, which have proceeded unpredictably due to complex and poorly understood feedback processes (Chasmer and Hopkinson 2017). Permafrost thaw is one of the most important and dramatic manifestations of climate warming in subarctic boreal ecoregions. In northwestern Canada, the most substantial permafrost thaw-induced landscape transformations are observed between 59°N and 62°N (Johannessen et al 2004, Quinton et al 2011, where the boreal forest overlaps with sporadic-discontinuous permafrost . In the sporadic-discontinuous zone (10%-50% areal permafrost coverage), permafrost is on the order of 10 m thick and typically isothermal at, or very close to, melting point temperature throughout its thickness (Kwong and Gan 1994), signifying it is in disequilibrium with current climatic conditions (Helbig et al 2016).
In the sporadic-discontinuous permafrost zone, permafrost occurs preferentially below forested peatlands such as peat plateaus, while adjacent peatlands occupying lower topographic positions including collapse scars and channel fens are permafrost-free (Zoltai and Tarnocai 1975, Hayashi et al 2004. As such, thaw of this permafrost results in its disappearance and the transformation of the overlying forested plateaus into permafrost-free, treeless wetlands Moore 2000, Quinton et al 2011). Boreal ecosystems with (e.g. peat plateaus) Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.
Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. and without (e.g. collapse scar wetlands) permafrost have contrasting ground surface energy balances (Ketteridge et al 2013, Chasmer and Hopkinson 2017) and hydrological regimes (Connon et al 2014). Consequently, a widespread transition from landcovers underlain by permafrost to those that are not, has the potential to fundamentally alter the direction, pattern, and rate of surface-atmosphere exchanges of mass and energy at regional scales, with feedbacks to the atmosphere (Hinzman et al 2005) and greenhouse gas fluxes (Prowse andFurgal 2009, Schuur et al 2015). The rates and patterns of this transformation have not been well characterised over the sporadic-discontinuous permafrost zone, and as such, the impacts of this landcover transition on surface-atmospheric fluxes remain poorly understood despite the demonstrated pan-arctic occurrence of this process (Camill and Clark 2000).
An average annual air temperature 0°C is an approximate indicator of the presence of permafrost (Halsey et al 1995, Beilman andRobinson 2003). However, dry peat is an extremely effective thermal insulator (Beilman et al 2001, Turetsky et al 2007, and as such, permafrost can persist in areas where the average annual air temperatures exceeds 2°C (Camill 1999). However, local variations of soil moisture, ground surface albedo and tree canopy cover result in large spatial variations in ground thermal regimes including the presence and characteristics of permafrost (Hayashi et al 2007).
Peatlands occupy approximately 12% of Canada's land area and almost exclusively (97%) occur in the subarctic boreal region (Tarnocai 2006, Tarnocai 2009 where cycles of permafrost aggradation and degradation over periods of decades to centuries are discernable from peat profiles (Tarnocai 2006, Turetsky et al 2007. However, warming-induced landcover change appears to have accelerated over recent decades (Halsey et al 1995, Quinton et al 2011 and disrupted the balance between permafrost aggradation and degradation such that the latter is favoured (Beilman et al 2001, Schuur andAbbott 2011). The most dramatic landcover change impacts appear in areas with rising air temperatures, particularly those with annual average temperatures nearing 0°C, indicating a climatic effect (Halsey et al 1995, Beilman and Robinson 2003, Prowse and Furgal 2009. Substantial permafrost thaw and degradation has led to forest loss at an unprecedented rate and much faster than new peat plateaus can form (Hayashi et al 2004. Two principal thawinduced mechanisms are understood to be responsible for this transition from forest to wetland; plateau inundation from the surrounding wetlands (Hayashi et al 2004) causing waterlogging and loss of the tree cover (Islam andMacdonald 2004, Iwata et al 2012) and/or peat over-maturation, as observed through fissures or cracks forming in plateaus (Zoltai andTarnocai 1975, Zoltai 1993). The latter further predisposes forested plateaus to collapse due to subsequent inundation (Beilman and Robinson 2003). This transition from forest to wetland is hastened by positive feedbacks that increase the rate of permafrost thaw (Connon et al 2018). As the ground surface transitions from a relatively dry forest floor to wetland, its thermal conductivity increases (Wright et al 2009) as does its receipt of shortwave radiation .
Landcover change is a striking and increasingly apparent consequence of climate change in the sporadic-discontinuous permafrost zone, and because of the direct connection between permafrost and landcover, remotely sensed landcover change observations are an excellent proxy for underlying permafrost thaw. This study examines ten areas of interest (AOIs) over a 200 km north-to-south transect extending over the sporadic-discontinuous permafrost zone and aims to (1) characterise the rate and pattern of forest to wetland conversion; and (2) quantify the variable changes to landcover across the AOIs with respect to climatic and environmental factors. These variables included annual and seasonal components of: average temperature, change in temperature, temperature range, average precipitation, change in precipitation, as well as changes to environmental elements such as fragmentation and changes to growing season length. These climatic and environmental factors can be used to assess a site's vulnerability to not only transition from forest to wetland but also has implications for the permafrost loss that is occurring over the same time period.

Study area
The sporadic-discontinuous permafrost zone across the southern Northwest Territories (NWT) and northeastern British Columbia (BC) is a peatland-dominated landscape and represented by a patchwork of forested peat plateaus underlain by permafrost interspersed by permafrost-free and treeless wetlands (Wright et al 2009). The study area is located near the isothermal boundary between the extensive-discontinuous and sporadic-discontinuous zones where the presence of peat-rich organic soils is a key control in determining the distribution of permafrost (Vitt et al 1994, Halsey et al 1995 and is almost exclusively responsible for any permafrost found in this zone (Wright et al 2009, Helbig et al 2016. Despite the thermal insulating properties of peat, the area underlain by permafrost has significantly decreased across the study area, and the amount of surface water on the landscape has increased Moore 2000, Wright et al 2009). This has been observed through a transition from permafrost plateaus and boreal forest dominating the landscape to an increase in overall wetland cover (Quinton et al 2011, Baltzer et al 2014. The ten AOIs (figure 1) were selected by targeting peatland complexes (i.e. subarctic lowlands) containing forests, bogs and fens with sufficient remote sensing coverage (i.e. approximately a 40 year historical record). The 10 mapped areas of interest cover a latitudinal gradient of 59.96°N to 61.30°N and are characterised by a dry continental climate with short summers and long, cold winters with mean annual air temperature (MAAT) ranging from −3.4°C to −2.5°C across the ten AOIs (Natural Resources Canada 2017). Over the 1970-2010 study period, average annual air temperatures increased approximately 0.9°C, which is part of a longer-term trend across these subarctic AOIs where up to 1.4°C in warming has occurred over the past century (Natural Resources Canada 2017). Precipitation also varies across the ten sites, generally increasing with decreasing latitude and ranging from 346 to 417 mm, annually (Natural Resources Canada 2017). Annual precipitation has increased by approximately 50 mm across all sites over the study period and while the majority of precipitation in this region falls as snow, effectively all of the observed increase was recorded as rainfall (Natural Resources Canada 2017).

Methods
Remote sensing methods Panchromatic Worldview 1 and 2 imagery captured in 2010 and 2011 were combined with 1970 and 1971 historical aerial photographs over each of the ten AOIs. The historical imagery was georeferenced using multi-temporal stationary marks as tie points. Given the remoteness of the region and limited permanent anthropogenic features, seismic line intersections and distinct lake shore boundary forms were selected in each image (Chasmer et al 2010). The georeferencing process used no fewer than 30 tie points at each site and the historical image was adjusted using a first order polynomial. The adjusted historical aerial photographs (∼1-1.5 m resolution) and satellite imagery (0.3 m resolution) were then subset to a 2.5 km by 4 km area, which was delineated using NAD 1983 Zone 10N over each study area.
Landcover classifications were then completed by manually digitising the 10-square kilometre subsets in ArcGIS (ESRI, Redlands, California) into two classes: 'forest' and 'wetland', representing the landcovers that are thought to be indicative of permafrost and permafrost-free ground, respectively. Lakes, streams, seismic lines and drill pads were also digitised at AOIs when applicable and all landcover types were classified and distinguished based on unique spectral signatures. To determine the amount of landcover change that occurred over the approximately 40 year study period, the total area of forested landcover as a proportion of the total 10 km 2 area was calculated for both the historical and present-day imagery at each site and areas of forest gain or loss were then identified and enumerated. To account for differences in proportional forest area between AOIs, spatial changes to the distribution of treed landcovers were reported as percentages of only the forested land. These change detection analyses included three measurements: net-forest change between the historical and satellite imagery (reported as a percentage), and its two components-forest loss and forest gain. The net change to forested area (C F ) was computed from: where A 1 is the area within an AOI occupied by forest in 1970, and A 2 is the area within the same AOI occupied by forest in 2010. In order to focus on changes to landcover that could be attributed to permafrost thaw, areas that experienced forest loss or gain that could be confidently explained by forest-fires (as per the Canadian National Fire Database 2015) or industrial development were removed from these change detection calculations (Burton et al 2008). These seemingly thaw-induced changes to the areal distribution of treed landcovers were then reported as net-forest loss between the historical and satellite imagery and were quantified as percentages for each site.
Climate data methods Modelled MAAT as well as monthly air temperature and precipitation grids were sourced from Natural Resources Canada (NRCAN) (2017). All temperature and precipitation grids were produced using the thin plate smoothing splines in the ANUSPLIN suite of programs ( 2006,2011). Despite this, the model appeared to avoid bias as both positive and negative residuals were observed, meaning these errors are more likely due to temperature variability induced by the northern latitude and proximity to complex topography (McKenney et al 2006). Additionally, despite the challenges associated with measuring winter precipitation, higher errors relative to summer precipitation were not observed given the relatively larger amounts of precipitation collected during summer months (McKenney et al 2006). The obtained climate grid was then re-projected to centre on the ten areas of interest included in this study. The data were then reanalysed as a mean annual air temperature dataset covering a variety of time periods within the entire 1901-2010 span of the data, while more seasonal factors-both temperature and precipitation-were confined to the period of time covered by this study . The centre point of each boreal peatland study area was then used to extract the climate variables for each year from 1901 to 2010. From this re-analysis product, changes in temperature and precipitation across the same time periods were also calculated. Despite the errors recognised in northwestern Canada, these climate data products have been previously applied to climate change and forest-related research throughout Canada (McKenney and Pedlar 2003, Price and Scott 2006, McKenney et al 2007a, 2007b.

Statistical methods
The landcover change measures of net forest loss, forest loss and forest gain are the residual values after the areas directly impacted by forest fire or industrial clearing were removed. The remaining landscape changes (net change, loss and gain) were significant at all ten AOIs across the 40 year study period (p<0.05).
To explore the potential mechanisms associated with these landscape change statistics, a variety of variables were explored. These factors included: a measure of fragmentation in the form of perimeter-area ratio, two average annual temperature measurements over different time periods, the average annual temperature range during the study period, six change in temperature measurements with varying time periods and lags, average annual temperature range over the study period, seasonal temperature and precipitation components-recorded as both averages and changes over the study period, and finally, length of the growing season. Three landscape-change metrics (net loss, loss and gain) from all ten AOIs were measured for strength and direction of association with 20 identified monotonic factors using a Spearman rank-order correlation coefficient (two-tailed, significant at 0.05 level).
Given the possibility of the multiple comparisons or multiplicity problem emerging (i.e. a greater likelihood of false positives or Type I errors occurring when a greater number of inferences or comparisons are made), the Benjamini-Hochberg (1995) procedure was applied. The false discovery rate method described by Benjamini and Hochberg (1995) controls the proportion of rejected null hypotheses by employing more stringent significance thresholds for individual correlations. This procedure is widely used as it is less conservative than the original Bonferroni correction, which has the opposite affliction where the number of false negatives or Type II errors can be falsely inflated (Benjamini and Hochberg 1995).

Landcover change
Changes to forested landcovers were assessed for all areas except for those directly affected by forest fire or by linear disturbances such as seismic lines. All AOIs experienced a net forest loss across the 40 year study period (figure 2). Net forest loss decreased with latitude from a maximum value of 11.6% at AOI A (Scotty Creek, NWT) to a minimum value of 6.9% at AOI I (Hossitl Creek, BC) (figure 3). Aside from the net change in forest cover, the gross values of forest loss and gain between 1970 and 2010 were also quantified for each AOI. At seven of the 10 sites, forest loss over the 40 year study period follows a similar trend to net forest loss as forest gain is occurring at negligible rates ranging from 0.2% to 1%. However, the two southernmost AOIs (I and J) experienced the maximum values of both gross forest loss (16.3%) and gain (9.4%). This pattern limited the net change in forest cover at these more southerly sites.

Predictors of landscape change
The three landscape change components: net forest change, gross forest loss, and gross forest gain were assessed for strength of correlation with 20 factors that were identified as potential predictors. Of all the factors and the 60 (3×20) potential relationships, 27 emerged as significant controls on one or more of the three landscape change statistics after the false discovery rate correction as determined by a Spearman rank-order correlation coefficient (ρ) (table 1).
Average temperatures, both over the length of the record (1901-2010) and the study period, are  . Relationships between net forest loss, forest loss, and forest gain with latitude over the 40 year study period across the ten selected study areas of interest (AOIs) (p<0.05). Table 1. Non-parametric correlation associations between net forest loss, forest loss and forest gain with (a) original (1970/71) perimeter-area ratio, and all annual climatic variables: average temperatures (1901-2010 and 1970-2010), change in temperature (1901-2010, 1920-1970, 1930-1970, 1940-1970, 1950-1970, 1970-2010), average annual temperature range, and average annual precipitation; and (b) all seasonal factors over the study period : average summer and winter temperature, average summer and winter precipitation , change in summer and winter temperature . Significant (p<0.05) correlations (ρ) have been highlighted by bold text.  positively correlated with both forest loss and gain between 1970 and 2010. Change in temperature over the length of the modelled record was also significantly correlated with forest loss, while change in temperature over the study period, as well as the four selected time periods before the first image, were found to be significantly correlated with both net loss and gain. The four time periods encompassing 50, 40, 30 and 20 years before the start of the study period attempt to account for any lags in landscape change associated with changes to temperature (Camill andClark 2000, Chasmer andHopkinson 2017). The magnitude of change in temperature over each of these time periods was significantly negatively correlated with net loss and positively correlated with any forest gain. Changes in temperature confined to the study period showed an opposite relationship where increases in temperature corresponded with increases in net forest loss while any decreases in MAAT over the study period corresponded with decreases in forest gain. A broader annual temperature range corresponded with further conversion to wetland (i.e. a higher net forest loss rate) and less forest gain on similar magnitudes. Average winter temperature over the study period was negatively correlated with net forest loss and positively correlated with forest gain between 1970 and 2010, while average summer temperature did not emerge as a significant contributor. The change in both summer and winter temperatures over the 40 year study period were both positively correlated with net forest loss indicating a greater positive change in temperature over the two seasons would lead to greater net forest loss over the same time period. Average annual precipitation was significantly negatively correlated with net forest loss and positively correlated with forest gain. Seasonally, average summer precipitation emerged as a negative correlation with net landcover change and positive with forest gain, while average winter precipitation did not emerge as a significant factor. Changes to summer precipitation did not appear as significant with any landcover change measurement. However, changes to winter precipitation emerged as significant in both directions, where a negative relationship emerged with net forest change and a positive correlation with forest gain.

Discussion
Previous studies have demonstrated significant permafrost degradation has manifested as forest loss and plateau collapse as well as wetland expansion (Wright et al 2009, Quinton et al 2011, Baltzer et al 2014 and these changes appear across all ten AOIs along the latitudinal transect of this study. Forest to wetland conversion is the dominant landscape succession trend observed across peatlands within the sporadicdiscontinuous permafrost zone, as higher average air temperatures corresponded with greater forest loss rates. This thaw-induced forest loss, whether through plateau inundation (Hayashi et al 2004 or peat over-maturation (Zoltai 1993, Beilman andRobinson 2003), has been well established as the dominant mechanism responsible for landcover change in these environments (Baltzer et al 2014, Quinton et al 2011. However, this study also introduces an additional positive correlation with forest gain. The relationship between increased air temperatures and greater afforestation is less well established but forest expansion and increased productivity have been proposed as potential responses to climate warming (Soja et al 2007, Baltzer et al 2014. However, it is important to note that no permafrost has returned to these areas despite the afforestation that has occurred, given the generally warmer climatic conditions (Robinson and Moore 2000, Turetsky et al 2007, Jorgenson et al 2010. This study also suggests that warming air temperatures may be responsible for forest expansion in addition to forest loss via thaw-induced mechanisms ( Northwestern Canada's sporadic-discontinuous permafrost zone has not only experienced increased climatic warming, but climate change has also impacted regional precipitation patterns. In addition to the increased rainfall observed over the study period (Natural Resources Canada 2017), there is also evidence that climate warming has impacted seasonal precipitation and water budget patterns via increased rain-onsnow events (Putkonen et al 2009), earlier spring melt (Pachauri and Reisinger 2007), and overall decreased snowpack (Hinzman et al 2005)-all of which appear to increase permafrost degradation. These hydroclimatic changes are representative of how seasonallydependent (and thus temperature-dependent) the impacts of precipitation are throughout this region.
The relationships between landcover change and summer precipitation found in this study can be explained through seasonally variable characteristics unique to peatland hydrology throughout the sporadic-discontinuous permafrost zone. Climate change has led to increases in warm-season temperature and precipitation (Natural Resources Canada 2017) as well as a thickening of the suprapermafrost layer, as observed by a deeper permafrost table (Price et al 2005, Connon et al 2018. While shallower frost tables leave the active layer in the highly porous portion of the profile, deeper frost tables are associated with lower porosity and hydraulic conductivity (Quinton and Gray 2003). This new and more dominant condition results in decreased subsurface runoff and an increase in overflow to the surrounding wetlands (Quinton and Gray 2003). Increases to overland flow may result in rainfall having a lessened impact on forest losses as the elevated plateau features are predominantly saturated from surrounding wetlands rather than the combined effects of both direct rainfall infiltration and wetland inundation. Despite this, both processes ultimately lead to forested plateau collapse as the associated tree species begin to decompose in their newly inundated environment due to their relative intolerance to waterlogging (Islam andMcDonald 2004, Iwata et al 2012).
The relationships found between landcover change and winter climate, both temperature and precipitation, can be attributed to the additional insulation snow accumulation provides the ground and underlying permafrost (Osterkamp 2007, Chasmer andHopkinson 2017). During mid-winter, when air temperatures reach their annual minimum, greater snow depth shields the ground from the extreme cold keeping the subsurface relatively warmer and further degrading the contained permafrost (Hinkel and Hurd 2006, Osterkamp 2007, Chasmer and Hopkinson 2017. However, greater snow cover also insulates the ground and the underlying permafrost from warming air temperatures, particularly through high albedo decreasing shortwave radiation absorption and delaying meltwater movement and the associated thermal conduction via latent heat transfer (Wright et al 2008(Wright et al , 2009. Similar northern regions have experienced decreases to snowpack depth (Hinzman et al 2005), earlier end-of-winter melt (Pachauri and Reisinger 2007) and increases in the frequency of midwinter melt events through warm periods and rainon-snow events (Putkonen et al 2009). As such, decreases in winter precipitation appear to lead to greater forest-to-wetland conversion, through permafrost thaw, while greater snow accumulation may have the potential to increase forest gain, though reduced plateau saturation following snowmelt.
Permafrost thaw-induced landcover change is occurring across a broad latitudinal span of the sporadic-discontinuous permafrost zone. Prior studies in the sporadic-discontinuous permafrost zone have identified fragmentation as a significant factor influencing landcover change (Payette et al 2004, Baltzer et al 2014, and specifically forest loss, however this relationship was not revealed in this study. Anthropogenic developments such as seismic lines, winter roads and drill pads for the oil and gas industry, while observed throughout many of the AOIs, were not directly included in this analysis, but may have impacts on these landscapes beyond forest clearing and fragmentation. Most immediately, seismic lines have been shown to result in surface subsidence that promotes the accumulation of water and the expansion of hydrologic features such as wetlands where permafrost (and their resulting black spruce canopies) cannot regenerate (Williams et al 2013). Ultimately the increased hydrological connectivity arising from permafrost thaw has the potential to reduce water storage on the landscape to the extent that black spruce forest may re-establish (Camill 1999, Murphy et al 2009, Ketteridge et al 2013, Waddington et al 2015. These patterns of reforestation appear to be particularly prevalent in areas that have already experienced significant permafrost thaw (Chasmer and Hopkinson 2017). These environments, including AOIs I and J, may be further along in the anticipated landscape succession and may be moving beyond the initial phase of landcover change where wetland expansion is the dominant process (Murphy et al 2009, Ketteridge et al 2013, Waddington et al 2015.

Conclusion
The results presented in this study established that widespread changes to forest cover and underlying permafrost have occurred over a 40 year period across northwestern Canada's vulnerable sporadic-discontinuous permafrost zone. A variety of mechanisms were explored to help explain how these relationships change across this large-scale gradient and 27 latitudinally-varying climatic and environmental factors emerged as significant across the ten AOIs. These results suggest that landcover in the sporadic-discontinuous permafrost zone is not only closely related to the climatic changes that are occurring but also that the more southerly AOIs appear to be governed by more complex and climatically and environmentally interdependent relationships as compared to the more linear changes occurring further north. This paper provides important insights into the trajectory of change in the sporadic-discontinuous permafrost zone as this study indicates that the established wetland-dominated conversion prevailing to the north may not necessarily represent a climax community, as formerly predicted, as there may be opportunity for reforestation to occur, even without permafrost returning.
(Consortium for Permafrost Ecosystems in Transition -CPET). We also acknowledge Geoscience BC for their support of CPET.