Long-term climate-influenced land cover change in discontinuous permafrost peatland complexes

The discontinuous permafrost zone is undergoing rapid transformation as a result of unprecedented permafrost thaw brought on by circumpolar climate warming. Rapid warming over recent decades has significantly decreased the area underlain by permafrost in peatland complexes. It has catalyzed extensive landscape transitions in the Taiga Plains of northwestern Canada, transforming forest-dominated landscapes to those that are wetland- 15 dominated. However, the advanced stages of this landscape transition, and the hydrological and thermal mechanisms and feedbacks governing these environments, are unclear. This study explores the current trajectory of land cover change across a 300,000 km 2 region of northwestern Canada’s discontinuous permafrost zone by presenting a north-south space-for-time substitution that capitalizes on the region’s 600 km latitudinal span. We combine extensive geomatics data across the Taiga Plains with ground-based hydrometeorological measurements collected in the 20 Scotty Creek basin, Northwest Territories, Canada, which is located in the medial latitudes of the Taiga Plains and is undergoing rapid landscape change. This data is used to inform a new conceptual framework of landscape evolution that accounts for the observed patterns of permafrost thaw-induced land cover change, and provides a basis for predicting future changes. Permafrost thaw-induced changes in hydrology promote partial drainage and drying of collapse scar wetlands, leading to areas of afforestation forming treed wetlands without underlying permafrost. 25 Across the north-south latitudinal gradient spanning the Taiga Plains, relatively undisturbed forested plateau-wetland complexes dominate the region’s higher latitudes, forest-wetland patchworks are most prevalent at the medial latitudes, and forested peatlands are increasingly present across lower latitudes. This trend reflects the progression of wetland transition occurring locally in the plateau-wetland complexes of the Scotty Creek basin and informs our understanding of the anticipated trajectory of change in the discontinuous permafrost zone.


Introduction
Northwestern Canada is one of the most rapidly warming regions on Earth (Vincent et al., 2015;Box et al., 2019) and it is transitioning to a warmer state at a rate that appears to have no analogue in the historical record (Porter et al., 2019). This transition includes region-wide thaw and disappearance of permafrost at unprecedented rates (Rowland et al., 2010). The Taiga 45 Plains ecoregion of northwestern Canada extends from 55º to 68º N and as such, encompasses the spectrum of permafrost cover, from continuous to sporadic. Permafrost thaw in the Taiga Plains ecoregion is especially pronounced in its lower latitudes where the permafrost is relatively thin and warm, often already at the thaw-point temperature (Biskaborn et al., 2019), indicating a state of disequilibrium with the current climate (Helbig et al. 2016a). For example, Kwong & 50 Gan (1994) repeated the permafrost surveys of Brown (1964) in northern Alberta and the southern Northwest Territories (NWT) and found that the southern limit of permafrost occurrence had migrated northward by about 120 km over a period of 26 years. Beilman & Robinson (2003) estimated that 30-65% of the permafrost has disappeared from the southern Taiga Plains in the preceding 150 years, most of which disappeared in the latter 50 years. The 55 accelerated rates of permafrost warming and thaw observed in recent decades throughout the circumpolar region (Biskaborn et al., 2019), including all of northwestern Canada Holloway & Lewkowiz 2019), have dramatically transformed land covers in the southern Taiga Plains (Chasmer & Hopkinson, 2017).
Much of the southern Taiga Plains is occupied by peatland-dominated lowlands, a 60 landscape of raised, black spruce (Picea mariana) tree-covered peat plateaus overlying thin (<10 m), ice-rich permafrost interspersed with permafrost-free, treeless wetlands. These permafrostfree wetlands are predominantly classified as channel fens and collapse scar wetlands, the latter of which is developed from thermokarst erosion of the plateaus (Robinson & Moore, 2000). Peat plateaus and collapse scar wetlands are typically arranged into distinct "plateau-wetland 65 complexes," which are separated by channel fens. Each of these major land cover types in the lowlands of the southern Taiga Plains, have contrasting hydrological functions (Hayshi et al., 2004) and therefore changes to their relative proportions on the landscape can affect water flux and storage at the basin scale . Permafrost thaw underlying plateaus is driven by horizontal conduction and advection from adjacent wetlands, and vertical heat flows 70 from the ground surface (Walvoord & Kurylyk, 2016). As this permafrost thaws, the overlying plateau ground surface subsides and is engulfed by the surrounding wetlands (Beilman et al., 2001;Quinton et al., 2011;Helbig et al., 2016a). As such, permafrost thaw in this environment transforms forests to treeless, permafrost-free wetlands (Robinson & Moore, 2000). In the process, this also changes the hydrological function of the transformed land cover, in part due to 75 a change in surface water-groundwater interactions (McKenzie & Voss 2013). Such a transformation can profoundly affect local drainage processes and pathways (Connon et al. 2014; with implications to regional hydrology (St. Jacques & Sauchyn, 2009;Korosi et al. 2017;Connon et al., 2018), ecology (Beilman, 2001) biogeochemical processes (Gordon et al., 2016) and carbon cycling (Vonk et al., 2019;Helbig et al. 2016a). 80 Zoltai (1993) described a perpetual cycle of permafrost development and thaw in which permafrost evolves from perennial ice bulbs that form below Sphagnum hummocks in permafrost-free treeless wetlands (i.e. collapse scars). Such hummocks expand and coalesce eventually forming tree-covered plateaus. However, over time plateaus experience a disturbance (e.g. fire, disease) that initiates the development of collapse scars and as a result, the plateau or 85 portions of it revert to a permafrost-free wetland. In a stable climate, the permafrost and permafrost-free fractions of a landscape are assumed to remain relatively consistent. Zoltai (1993) estimated that the time required to complete this cycle is approximately 600 years. Treat and Jones (2018) indicated time scales for forest recovery following permafrost thaw in the range of 450 to 1500 years. However, there is growing evidence throughout the southern Taiga Plains 90 that the climate warming of recent decades has disrupted the cycle of permafrost thaw and redevelopment such that the rates of permafrost loss greatly exceed those of permafrost development (e.g. Halsey et al., 1995;Robinson & Moore, 2002;Quinton et al., 2011).
The accelerated rates of permafrost thaw and resulting land cover change described above call into question the utility of existing concepts (e.g. Zoltai, 1993) as a means to estimate the 95 current trajectory of land cover change since such concepts were developed from analyses of geological sediments (e.g. peat cores) which generally lack the resolution needed to identify land cover change sequences over relatively short (i.e. decadal) periods. Moreover, it is uncertain whether the current rates of climate warming are represented in the sediment record. As a result, there remains considerable uncertainty on the trajectory of permafrost thaw-induced land cover 100 change in this region, including possible end-members and intermediate stages. Because of the close connection between land cover type and hydrological function in this region, the uncertainty related to possible land cover change trajectories also raises new uncertainties in regards to the region's water resources.
In addition to unprecedented climate warming in the North, accelerated permafrost thaw 105 is also driven by positive feedbacks including increased fragmentation of forested peat plateaus with increasing thaw , a process which increases the length of interface between permafrost and permafrost-free terrain, and therefore also increases the overall flux of energy into the remaining permafrost bodies . Connon et al. (2018) demonstrated that talik layers situated between the active layer and underlying permafrost are 110 widespread in thawing peatland-dominated terrains and their occurrence increases with increasing permafrost thaw. Devoie et al. (2019) demonstrated that once a talik forms, the rate of permafrost thaw can increase 10-fold.
Since the permafrost table beneath peat plateaus rises above the water surface of the adjacent wetlands, plateaus function as "permafrost dams" that prevent wetlands from draining. 115 Permafrost thaw therefore removes this effect and enables previously impounded wetlands to partially drain until the hydraulic gradient driving their partial drainage reaches an equilibrium state (Haynes et al., 2020). The slow release of water from the long-term storage of wetlands no longer impounded by permafrost changes the physical and ecological characteristics and hydrological function of these wetlands (Haynes et al., 2020). Such drainage transforms the 120 uniformly wet Sphagnum lawns that characterise impounded wetlands, into hummocky surfaces that provide a wider range of near surface moisture conditions including those sufficiently dry to support the re-growth of trees (Haynes et al., 2020). There is also evidence that when black spruce forest is lost due to permafrost thaw and plateau inundation, forest regeneration does not depend on the regeneration of permafrost (Haynes et al. 2020;Chasmer & Hopkinson 2017). For 125 example, treeless collapse scars have transformed into black spruce forest within two to three decades after the permafrost dams disappear (Haynes et al. 2018).
In addition to the transient drainage process described above that may occur following the removal of the impounding permafrost (Haynes et al., 2018), such removal also increases the hydrological connectivity of basins through the incorporation of wetlands that were previously 130 impounded and therefore hydrologically isolated from the basin drainage network (Connon et al., 2015). This process of "wetland capture" expands the runoff contributing areas of basins, a process that increases their runoff potential. Connon et al. (2018) attributed the trends of increasing runoff ratio (i.e. fraction of basin runoff per unit input of precipitation) in basins throughout the Taiga Plains to this permafrost thaw-induced process of runoff contributing area 135 expansion.
The transition of one type of ground cover to another as a result of permafrost thaw or a subsequent process such as partial wetland drainage and re-establishment of forest also results in a change of surface energy balance Devoie et al. 2019). Insight into the nature of such changes can be obtained through comparing the energy regimes of the existing 140 suite of land covers including the end-members of land cover change. For example, the incoming solar radiation measured at a height of 2 m above the ground surface is highest in the treeless wetlands and lowest in areas of peat plateaus with dense forest . The average shortwave radiation flux density of treeless wetlands is approximately twice of that measured below dense forest . Plateau areas with moderate or sparse tree 145 canopies have incoming solar radiation values intermediate between these two end members . The ground surface albedo varies over the narrow range of 0.15 to 0.19 (Hayashi et al., 2007) among the ground surface types discussed here, the exception being the late snowmelt period while plateau ground surfaces are still snow covered and the treeless wetlands are snow-free (Disher et al. 2021;Connon et al., Submitted). 150 The nature of changes to a land cover's surface energy balance is governed by the properties of its subsurface, ground surface, and the overlying tree canopy, all of which change as one land cover type transitions to another (Helbig et al. 2016b). The reduction in the areal cover of forested plateaus and concomitant increase in the coverage of treeless wetlands indicates that in the first instance, permafrost thaw increases the incoming shortwave flux to the 155 transformed land cover Devoie et al. 2019). Chasmer et al. (2011) found that this thaw-induced transition and associated increase of incoming shortwave radiation occurs over several years as tree mortality decreases the density of tree canopies. However, the forests that subsequently re-establish in partially drained wetlands may have an energy balance that shares some characteristics of the forested peat plateaus, where insolation is relatively low, and 160 the low albedo (and therefore high energy adsorption) of trunks, branches and stems result in relatively high long-wave and sensible heat compared to the treeless wetland surfaces (Helbig et al. 2016b).
Unprecedented climate warming and the feedbacks to thaw and land cover change are new factors not accounted for in current theories on permafrost degradation-aggregation cycles 165 based on the analysis of peat cores. As a result, the time scales for land cover transformations derived from such theories cannot account for the current rates and patterns of all thaw-induced land cover change. This study examines peat plateau-wetland complexes along a latitudinal gradient through the Taiga Plains to improve the understanding of permafrost thaw-driven land cover change in this region as well as to advance the ability to predict land cover changes over 170 the coming decades. This overall objective will be accomplished by: (1) delineating the current extent of peatlands and forest distribution along the latitudinal span of discontinuous permafrost; (2) characterising the end-members and intervening stages of land cover transition; (3) providing an interpretation of the hydrological and ground surface energy balance regimes for each stage of land cover transition based on twenty years of field studies at the Scotty Creek Research Station; 175 and (4) presenting a conceptual framework of peatland transition during and following permafrost thaw.

The Taiga Plains Ecozone
Much of northwestern Canada's boreal region is located within the discontinuous 180 permafrost zone, which ranges latitudinally from extensive-discontinuous (50-90% areal permafrost coverage) in the north to sporadic-discontinuous (10-50%) in the south. Within this region, the Taiga Plains ecozone contains a patchwork of mineral and organic terrain. This study examines the peat plateau-collapse scar wetland complexes that dominate the lowlands of this ecoregion (Wright et al. 2009;Helbig et al. 2016a). While air temperature is the predominant 185 control on permafrost, relatively dry peat at the ground surface can allow permafrost to exist where mean annual air temperatures (MAATs) are at or even above 0°C due to thermal insulation (Vitt et al. 1994;Camill & Clark 1998). Permafrost is therefore largely restricted to below peat plateaus since only these features contain unsaturated layers sufficiently developed to insulate permafrost (Zoltai & Tarnocai 1975;Hayashi et al. 2004;Quinton et al. 2009). The areal 190 coverage of permafrost in the discontinuous zone has significantly decreased in recent decades due to increasing MAATs and has resulted in a shift towards more wetland-dominated landscapes (Thie, 1974;Robinson & Moore, 2000;Wright et al. 2009;Quinton et al., 2011;Olefeldt et al. 2016).
The discontinuous permafrost zone of the Taiga Plains ecozone covers 312,000 km 2 and, 195 for the purposes of this study, is divided into the areas of extensive-discontinuous permafrost (151,000 km 2 ) and sporadic-discontinuous permafrost (161,000 km 2 ) (Brown et al., 2002; Figure   1). The Taiga Plains, bounded by the Taiga Cordillera to the west and Taiga Shield to the east, has a dry continental climate with short summers and long, cold winters with MAATs ranging from -5.5°C to -1.5°C (Vincent et al. 2012). MAATs have increased across the Taiga Plains over 200 the past 50 years (1970 -2019) (Vincent et al. 2012) in a manner consistent with panarctic warming (Overland et al., 2019). This is largely due to increases in average winter and spring temperatures of approximately 3°C over this period (Vincent et al. 2012). However, there has been no consistent trend in mean annual precipitation over this period in the Taiga Plains (Mekis & Vincent, 2011). 205 Figure 1: The Taiga Plains ecoregion with the discontinuous permafrost zones (coloured) defining the study region (Brown et al. 2002). The location of Scotty Creek Research Station (SCRS) is also indicated. Contains information licensed under the Open Government Licence -Canada.

Scotty Creek, Northwest Territories
Scotty Creek (61.3ºN, 121.3ºW) has been the focus of field studies and monitoring since the mid-1990s and as such, the long-term and detailed data archive at Scotty Creek  provide a unique opportunity to evaluate land cover changes over a period that 215 coincides with rapid climate warming. Scotty Creek therefore also provides a reference to interpret land cover changes for terrains that are also present throughout the region. Scotty Creek lowlands of the Taiga Plains, and as such its landscape is dominated by complexes containing tree-covered peat plateaus overlying permafrost alongside treeless and permafrost free collapse scar wetlands. Such plateau-wetland complexes are separated by channel fens that collectively function as the basin drainage network (Hayashi et al. 2004;Quinton et al. 2009). This type of land cover not only dominates the lowlands of the Taiga Plains but is also found extensively 230 throughout northwestern Canada and across the circumpolar subarctic (Olefeldt et al. 2016).

Geomatics Methods
To place Scotty Creek into a regional context, geomatics methods were applied to both zones of discontinuous permafrost within the Taiga Plains to quantify the areas occupied by each 235 of the major land covers of all areas identified as peatland-dominated lowland. Multispectral Landsat 8 imagery (30 m resolution; Figure 2a) was acquired across an area of over 300,000 km 2 totalling 70 Landsat scenes. Of these, 59 scenes were used to construct the base of the mosaic and 11 were used as secondary data to patch and minimize cloud cover. The 59 primary scenes were acquired in 2017 and 2018 while the 11 secondary scenes were acquired between 2013 and 240 2016 as data of suitable quality was unavailable during the preferred time period. Acquiring imagery during the snow-free season was prioritized and as such, all 70 Landsat tiles were acquired in June, July, or August, rendering the coniferous forest cover seasonally comparable and allowing for a more streamlined mosaicking process. A colour infrared mosaic (Landsat 8 bands 5, 4, 3 displayed as R, G, B; Figure 2b) was created across the study region in ArcGIS 245 (ESRI, Redlands, California) using a Lambert Conformal Conic projection. The mosaic dataset was colour balanced and the boundary was amended to the Taiga Plains ecozone including the delineations dividing the sporadic and extensive discontinuous zones (Brown et al. 2002). The saturated soils dataset is part of a larger digital cartographical project of Natural Resources Canada, CanVec. The CanVec dataset is a vector format dataset, which can be 270 downloaded by province/territory or Canada-wide and includes over 60 features organized into 8 themes, including land features. Land features in this dataset, such as the distribution of saturated soils, were originally digitized at a scale of 1:50000 (Natural Resources Canada 2017). The NCSCD is also a polygon database developed by the Bolin Centre for Climate Research through synthesizing data from numerous regional and national soil maps alongside field-data collected 275 across Canada, USA, Russia, and the European Union. The NCSCD includes data on the fractional coverage of different soil types and stored soil organic carbon (Hugelius et al. 2013a;Hugelius et al. 2013b). In the present study, the layer containing information on the fractional coverage of soil types was used. While the original format of the NCSCD is a vector of delineated zones, gridded data is also available at resolutions varying from 0.012° to 1° 280 (Hugelius et al. 2013b). The NCSCD is comprised of a circumarctic dataset as well as countrywide and regional datasets, including one of Canada (Hugelius et al. 2013b).
The NCSCD is a widely used dataset (Olefeldt et al. 2014;Gibson et al. 2018;Stofferahn et al., 2019;etc.) but the zones do not map specific locations of peatland-dominated terrain ( Figure 2d). The locations of peatlands is helpful for work in regions such as the Taiga Plains,  285 where the landscape is a patchwork of both organic and mineral terrain. The saturated soils dataset and the NCSCD were then both masked to the Taiga Plains boundaries in ArcGIS, where over 26,000 saturated soil polygons and 572 NCSCD zones were contained within the study region. The saturated soils dataset was mapped to display probable peatland terrain across the study region (Figure 2c). The areas of each saturated soil polygon were calculated alongside the 290 areas for each NCSCD zone using the boundaries in the dataset. As the fractional coverage product from the NCSCD was used in this study, the fractional area of probable peatland terrain within the same NCSCD zone was calculated. The fractional areas of organic soils reported in the NCSCD were then compared to the fractional areas of probable peatland terrain from the saturated soils dataset within the same NCSCD zone boundary (Figure 2e). 295 The Landsat mosaic dataset (Figure 2b) was then combined with the resultant product displaying peatland terrain (Figure 2e). An unsupervised land cover classification was subsequently completed on the Landsat mosaic across the areas identified by the saturated soils and NCSCD datasets to identify and classify the land covers within these peat plateau-wetland complexes ( Figure 2f). The first iteration of the unsupervised classification (Iso Cluster 300 classification approach) targeted 50-75 classes (72 created). The original 72 classes were then aggregated into 12 final classes within the peatland terrain outlined across the Taiga Plains study region. The final 12 aggregated classes include: coniferous (dense and sparse), mixed (dense and sparse), and broad leaf forests stands (dense and sparse), collapse scar, fen, open water, bare ground, cloud, and cloud shadow. 305 Forested peatlands are particularly indicative of landscape change in this region (Quinton et al. 2010;Baltzer et al. 2014;Chasmer & Hopkinson 2017) and as such, identifying the forested areas within the already identified peatland-dominated terrain was the focus of the Landsat classification. Specifically, the proportion of coniferous forested area within the total peatland area was quantified across the region's latitudinal span (Figure 2g). Fractional 310 coniferous forested area was selected rather than total forested area to account for the observed spatial differences in peatland distribution across the Taiga Plains. For each degree of latitude, a bin was created for fractional forested area and the median was calculated alongside upper (i.e. 75 th percentile) and lower (i.e. 25 th percentile) quartiles. This data was plotted as a function of latitude across the Taiga Plains ecozone. This generated a dataset of forest cover across the 315 peatland-dominated regions of interest that was subsequently complemented by field data collected in the Scotty Creek basin to guide the proposed conceptual framework.

Scotty Creek Imagery
To

Hydrological Data
A comprehensive archive of hydrometeorological measurements was used in this study to examine the temporal variation in hydrological characteristics as land cover transition from one 335 stage to another. The form and hydrological function of the major land cover types of permafrost plateau, collapse scar, and channel fen are well understood from numerous studies at Scotty Creek since the 1990s . Field studies and monitoring at Scotty Creek over this period have also provided firsthand accounts of how permafrost thaw changes land covers ). In the present study, we examined how runoff, evapotranspiration, and 340 water storage are affected as land cover changes. In addition, we examined the precipitation data collected from 2008 to 2019 (Geonor, Model T200B) in relation to the three hydrological components listed above to gain insights into how changes in land cover affect the water balance for each stage in the land cover transition. These stages will be presented in detail in section 4.2.
The Geonor precipitation data include both rain and snow measurements logged at 30 minute 345 intervals ( Table 1). Monitoring of discharge from Scotty Creek by the Water Survey of Canada began in 1996. For this study, annual basin runoff (mm year -1 ) between 1996 and 2015 was calculated and used in the basin runoff component of the conceptual framework (Table 1) (Connon et al. 2014;Haynes et al. 2018). Given that this period of discharge monitoring coincided with a period of considerable climate warming and documented land cover change at 350 Scotty Creek, the trend in calculated runoff over the period of record reflects a shift from a permafrost plateau-dominated landscape to one increasingly influenced by hydrologicallyconnected wetlands. Therefore, the temporal trend of runoff from the Scotty Creek basin is driven by permafrost thaw-induced land cover change (Connon et al. 2014;Haynes et al. 2018).  Warren et al. (2018) were converted to annual ET

355
(mm year -1 ) for the purpose of the conceptual framework water balance in the present work 365 (Table 1). The different land covers monitored by Warren et al. (2018) are representative of the land cover end-members identified in our conceptual framework. Therefore, the annual ET values based on the data collected by Warren et al. (2018) for the black spruce forests, open collapse scar wetlands and the integrated landscape were associated with the appropriate stage along our proposed trajectory of change. Stages of the trajectory for which representative 370 measurements were not collected are interpolated between the land cover end-members for which ET was measured.
Given the insignificant changes in annual precipitation over the period of measurement (Connon et al. 2014;Haynes et al. 2018), annual storage was calculated as the residual of annual precipitation inputs, and annual evapotranspiration and runoff outputs for the conceptual 375 framework water balance (Table 1).

Radiation Fluxes
Four meteorological stations at Scotty Creek were selected for use in this study ( Figure   3). This included a station installed in a collapse scar wetland in 2004 (hereafter "wetland station") followed by a second station on a densely forested peat plateau in 2007 (hereafter 380 "dense forest station"). The radiation flux data from these two stations are representative of the collapse scar wetland and permafrost plateau land cover types, respectively. Two additional stations located on forested plateaus were also used to represent tree canopy densities different from that of the dense forest station. These stations were installed on a sparsely forested peat plateau in 2015 (hereafter "sparse forest station") and a forested plateau with a canopy of 385 intermediate density between that of the dense and sparse stations in 2014 (hereafter "intermediate forest station"). All radiation measurements were made below the tree canopy at a height of 2 m above the ground surface. Four component radiation data were collected at the dense forest, sparse forest, and wetland meteorological stations, while only shortwave radiation was collected at the intermediate forest station. The reader is directed to Haynes et al. (2019) for 390 full descriptions of the radiation instrumentation within the Scotty Creek basin. Radiation was measured every minute, and averaged and recorded every 30 minutes from which daily (24 hourly) averages were computed. The daily averages were then used to compute annual average radiation for each station. While these computations defined some of the variability of radiation fluxes among land cover types, they do not account for flux variations over short temporal and 395 spatial scales (Webster et al. 2016). To address these, the daily average four component radiation data from each station were compared on a monthly time step. The monthly averages were calculated and compared across the land covers represented by each of the four meteorological stations using a one-way analysis of variance (ANOVA) with Tukey post-hoc test (a = 0.05), thereby testing the effect of land cover on monthly shortwave and longwave incoming and 400 outgoing radiation.

Peatland and Forest Occurrence
The type of peatland-dominated terrain composed of peat plateau-wetland complexes separated by channel fens as described for Scotty Creek, occupy approximately 35% of the discontinuous permafrost zones of the Taiga Plains (Figure 4a). Large peatland clusters are located in lowland areas with high histel or histosol soil percentages. In the extensive-415 discontinuous permafrost zone, peatlands are clustered to the west near to the Mackenzie River, and are largely absent from the eastern portion of the study area in the region bounded by Great Bear Lake to the north and the Taiga Shield to the east. In the sporadic-discontinuous zone however, the peatland clusters are more longitudinally dispersed.

425
Comparing the fractional areas of probable peatland terrain from the saturated soils dataset to the NCSCD showed the saturated soils dataset was more likely to overstate the distribution of probable peatland terrain compared to the NCSCD maps. Approximately 20% of the fractional areas were exact matches between the two datasets, 20% were lower in the 430 saturated soils dataset, and 60% were higher in the saturated soils dataset. However, despite these disagreements, 79% of the fractional areas determined using the saturated soils dataset were within 15% of the fractional areas in the NCSCD. This suggests that using the Natural Resources Canada saturated soils dataset may be an appropriate method of mapping probable peatland terrain in the Taiga Plains on a finer scale (Figure 4a) compared to the broad zones 435 presented by the NCSCD (Figure 4b). Only 11 of the 572 NCSCD zones (~2%) had disagreements over 25% when comparing the fractional areas between both datasets. The majority of these zones of disagreement were located along the Slave River, in the far southeast of the Taiga Plains study region.
A latitudinal trend in land cover percentage was found for the mapped peatland-440 dominated terrain ( Figure 5). Along the boundary between the extensive-discontinuous and sporadic-discontinuous permafrost zones near the centre of the study region, collapse scar wetland features are most prevalent. Median fractional forest cover in peatlands (i.e. peat plateaus or treed wetlands) reaches its minimum value of 33% within the 61º N bin, near the latitude of Scotty Creek. The proportion of forested peatlands remains relatively low throughout 445 the transitional zone between sporadic and extensive discontinuous permafrost (approximately 61º to 62º N) where the median forest cover does not exceed 34%. The widespread occurrence of collapse scars suggest that permafrost thaw and the resulting processes of ground surface subsidence and inundation are particularly active in this zone compared to the extensive discontinuous permafrost zone to the north and the sporadic discontinuous permafrost zone to the 450 south, where the fractional forested areas are higher. The median proportion of forested peatlands in the extensive-discontinuous zone (63º to 460 66º N) ranges from approximately 35 to 45%, indicating that permafrost thaw is less prevalent over the landscape than in the transition zone immediately to the south. However, the median fractional forested area is at its greatest (52%) south of the transition in the sporadicdiscontinuous zone (59º to 60º N), where about half of the peatland area is forest covered.
Expansion of forest cover in this zone, especially in the areas of north-eastern British Columbia 465 and north-western Alberta, has been reported by others (e.g. Zoltai 1993, Carpino et al., 2018. A pattern of forest expansion over permafrost-free terrain is consistent with the observation of a northward-moving southern limit of permafrost reported by Kwong and Gan (1994) and with the process of tree re-establishment following permafrost thaw-induced partial drainage of wetlands described by Haynes et al. (2020). 470

Conceptual Framework of Land Cover Change
From the remote sensing and field-based hydrological studies at Scotty Creek since the mid-1990s, key insights into incremental land cover changes initiated by permafrost thaw have emerged. Using this knowledge as a foundation, the present study examines both the hydrological and radiation regimes of each incremental stage and the land cover changes over the 475 larger region in which permafrost thaw is known to be widely occurring and within which the southern edge of permafrost is migrating northward (Kwong & Gan, 1994). From this approach, a new conceptual framework is presented (Figure 6 (VI) Hummock growth with forest establishment; and (VII) Forested peatlands ( Figure 6). In the following sections, the biophysical, hydrological and radiation regimes of each of these stages are presented and discussed, drawing on several investigations in the region.

Biophysical Characteristics
The present study found that the early land cover stages presented in Figure 6 are more 505 prevalent at the higher latitudes of the study region and the later stages at the lower latitudes.
Considering that permafrost thaw is more advanced in the lower latitudes and that the southern limit of permafrost is advancing northward (Kwong & Gan, 1994), it is reasonable to expect that the more advanced stages presently characterising the lower latitudes will, in the future, characterise the higher latitudes, assuming a continuation of climate warming induced permafrost 510 thaw. Approaching a change in latitude through the study region as analogous to a change in land cover stage, or more specifically, to a change in time, is supported by studies that examined land cover change over the last half-century at Scotty Creek (Chasmer & Hopkinson, 2017;Quinton et al., 2019) and along a north to south transect extending from Scotty Creek to northeastern British Columbia (Carpino et al., 2018). 515 The Scotty Creek basin, located near the northern limit of sporadic discontinuous permafrost ( Figure 1) is characterised mainly by stages III and IV. However, examples of all seven stages can be found in local areas at Scotty Creek. For this reason, the long term monitoring and research programs at Scotty Creek involving each of these land cover types contributes detailed information on their form and functioning, and on their transition from one 520 to another. For example, permafrost thaw changes a landscape dominated by forested plateaus ( Figure 6I) to one with small, suprapermafrost taliks ) and isolated collapse scars ( Figure 6II; Quinton et al. 2011). Continued thaw expands isolated collapse scars ( Figure   6III; Devoie et al. 2019) enabling them to coalesce to form interconnected wetlands (Connon et al. 2015), a process that then leads to a landscape of expansive wetlands dotted with isolated 525 plateau "islands" (Figure 6IV; Baltzer et al. 2014;Chasmer & Hopkinson 2017). Since the process of wetland expansion removes peat plateaus, a land cover type that impounds wetlands and obstructs drainage (Connon et al., 2014), this process enables the landscape drain more efficiently (Haynes et al., 2018). As wetlands drain, hummock micro-topography develops in their relatively drier interiors ( Figure 6V; Haynes et al. 2020), which allows black spruce to 530 colonise the wetlands on the relatively dry hummock surfaces ( Figure 6VI; Iversen et al. 2018;Dymond et al. 2019). Continued drainage and drying of wetlands enables the expansion of their hummocky terrain and therefore of their tree cover (Eppinga et al. 2007;Iversen et al. 2018) until the landscape returns to a more continuous forest cover ( Figure 6VII; Carpino et al. 2018).
However, the forest cover in this final stage is permafrost-free and for that reason, the conceptual 535 framework presented in Figure 6 stands in contrast to those presented by Zoltai (1993) and Camill (1999) in which the re-emergence of a forest cover relies on the re-emergence of the underlying permafrost. According to Zoltai (1993), forest re-emerges because permafrost displaces the overlying ground surface upward, resulting in the development of an unsaturated layer suitable for tree establishment. By contrast, the tree establishment described in Figure 6  540 results not from the re-emergence of permafrost, but from its continued thaw over the landscape, a process that dewaters wetlands (Connon et al., 2014;Haynes et al., 2018) to the extent suitable for tree establishment (Haynes et al. 2020).
The concept of land cover transition depicted in Figure 6 also stands in contrast to those that couple forest with permafrost thaw and re-development in terms of the time frame. 545 Permafrost-enabled re-establishment of forest requires a minimum of several centuries to complete (Zoltai, 1993;Treat & Jones, 2018) whereas the wetland drainage enabled reestablishment of forest ( Figure 6) occurs in less than half a century as suggested by the archived imagery for Scotty Creek. There, the early stages (i.e. stages I, II) represent the changes observed between 1947 and 2000 over which time the tree-covered area decreased from approximately 550 70% to approximately 50% . The fraction of forested land above permafrost-free terrain cover is unknown for this period but assumed to be negligible based on land cover descriptions for this region (e.g. NWWG, 1988;Zoltai 1993;Robinson & Moore, 2000). This period therefore generated a concomitant rise in the cover of wetlands over the landscape (i.e. 30 to 50%) since permafrost thaw transitions the tree covered plateaus to collapse 555 scars and channel fens, as confirmed by analysis of archived imagery (Chasmer et al., 2010). By 2018, tree-covered peat plateaus decreased to 40%, and treeless wetlands occupied 45% of the land cover (13% collapse scars and 32% channel fens) (Disher 2020). These studies indicate that permafrost thaw and the resulting processes have both removed forest as a result of thaw-induced subsidence and inundation of plateau surfaces, and more recently, enabled forest re-560 establishment in the form of treed wetlands (Haynes et al. 2020;Disher et al., 2021). However, the dominant land cover transition at Scotty Creek is still from forest (peat plateau) to wetland as a result of permafrost thaw, resulting in a net forest loss, a process that will continue until the later stages of Figure 6 are reached (i.e. stages VI, VII), at which point there will be a net forest gain. 565 The sequence of land cover stages following permafrost thaw observed at Scotty Creek and depicted in Figure 6 is supported by vegetation successional changes described in the literature for wetlands as they age. For example, aquatic Sphagnum species, notably S. riparium, are the first to occupy the inundated margins between thawing permafrost plateaus and developing collapse scars (Garon-Labrecque et al. 2015;Pelletier et al. 2017). Such recent areas 570 of collapse are easily identified on high-resolution RPAS imagery by the distinct bright green colour of S. riparium ( Figure 6II, III, IV; Gibson et al. 2018;Haynes et al. 2020). These wetland-plateau edges may also be identified by bare peat or moats of water (Zoltai 1993). As collapse scars expand, lawn species, such as S. angustifolium, and hummock species, such as S. fuscum, emerge, particularly in the drier interior of wetlands (Zoltai 1993;Camill 1999;Pelletier 575 et al. 2017). Hummock species, mainly S. fuscum, first emerge near the centre of collapse scars, and expand outward over time (Camill 1999;Loisel & Yu 2013). Much like S. riparium, S. fuscum is also easily identified in high-resolution imagery, where S. fuscum is distinguished by its russet colour ( Figure 6V) (Haynes et al. 2020). As the density of the S. fuscum hummocks increases, imagery and ground-based observations indicate the presence of young black spruce 580 trees (Liefers & Rothwell 1987;Haynes et al. 2020), first on isolated hummocks ( Figure 6VI) but eventually as widespread afforestation ( Figure 6VII; Camill 2000;Ketteridge et al. 2013).

Radiation Flux Characteristics
The Scotty Creek basin is a microcosm of its larger regional setting since it contains each of the land cover stages of the conceptual framework in Figure 6. As such, the micro-585 meteorological measurements made at Scotty Creek for different land cover types provide insight into how energy regimes change as one land cover stage transitions to the next. Both incoming and outgoing shortwave radiation peak at the middle stages (IV, V), where treeless collapse scars predominate. Annual incoming and outgoing shortwave radiation is lowest at the dense forest station, which represents the initial stage (I). Likewise, incoming and outgoing annual longwave 590 radiation are greatest in the early (I, II) and late stages (VI, VII) and lowest in the wetlanddominated middle stages (IV, V). Statistically significant differences were found between stations for incoming ( Figure 7a) and outgoing (Figure 7b) shortwave and incoming longwave radiation (Figure 7c), while there was no statistical differences between stations for outgoing longwave (Figure 7d). However, 595 Tukey post-hoc tests revealed variability between the two shortwave components in terms of which groups showed these significant differences. As the four meteorological stations fall along a gradient of forest density from treeless wetland to a densely forested plateau, no significant differences in incoming shortwave radiation were ever found between stations only one rank apart on that gradient. As such, measurements indicate average monthly incoming shortwave 600 radiation is significantly greater at the wetland compared to both the intermediate forest (p < 0.05) and the dense forest (p < 0.05), while no significant difference exists between wetland and the sparse forest. However, the dense forest receives significantly less incoming shortwave radiation than both the wetland (p < 0.05) and the sparse forest (p < 0.05), but this station is not significantly different from the intermediate forest. 605 No significant differences in outgoing shortwave radiation exist between any of the forested plateau stations but all are significantly different from the wetland. Specifically, outgoing shortwave radiation recorded at the wetland station is significantly greater than the sparse forest (p < 0.05), intermediate forest (p < 0.05), and dense forest stations (p < 0.05). The differences between the wetland station and the forested plateau stations are also increasingly 610 significant with increasing tree density. There was a statistically significant difference between the three stations (intermediate forest omitted due to lack of measurements) for incoming longwave radiation, while no statistically significant differences were observed in the outgoing longwave radiation component. The only significant difference in incoming longwave was observed between the wetland and dense forest (p < 0.05). No significant differences exist 615 between the sparse forest and the wetland or dense forest.
The comparison among the plateau stations of contrasting tree canopy densities provides insight into the permafrost thaw-induced progression of radiation regimes as plateaus transition to wetlands, a process involving the gradual thinning and eventual loss of the tree canopy. Wright et al. (2009) demonstrated that small-scale changes to the tree canopy density can 620 increase insolation to the ground in localized areas leading to thaw depressions in the active layer and water flows toward such depressions from their surroundings. Such areas of preferential thaw therefore develop elevated soil moisture contents, and since soil thermal conductivity increases with its moisture content, the preferential thaw process is reinforced. This is suggested as the mechanism driving the transition from stage I to II in the trajectory . 625 This feedback is present in the initial stages of the trajectory, and is often associated with talik formation and expansion into collapse scars due to localized permafrost loss (Chasmer & Hopkinson 2017;Connon et al. 2018). Such thaw can extend to the base of the active layer in which case further thaw results in permafrost loss, ground surface subsidence, waterlogging of the ground surface, local tree mortality, and therefore further thinning of the overlying tree 630 canopy and consequently more insolation at the ground surface. These processes and feedback mechanisms are critical in the generation of collapse scars in stage III.
Differing ground surface properties, particularly albedo, can amplify the differences in incoming shortwave radiation among the land covers. However, the difference in mean albedo during the snow-free season (May-September) below the plateau canopies is less than 5% and 635 displays a small increasing gradient as the canopy becomes more dense (sparse: 0.111, intermediate: 0.127, dense: 0.147). The mean wetland albedo (0.145) during the snow-free season is also similar to the plateau surfaces and most closely resembles the surface albedo of the dense plateau. The greatest contrast in albedo occurs during the period of several weeks while snow still covers the plateaus but is absent on the adjacent wetlands (Connon et al., Submitted). 640 This contrast in albedo is also evident in Figure 7b, which shows that following winter, outgoing shortwave radiation from the wetland increases before the forested stations. Helbig et al. (2016b) attributed their observed increase in landscape albedo in this late winter/spring period to the permafrost thaw-induced conversion of forest (lower albedo) to wetland (higher albedo), and suggested that this could lead to a regional cooling effect during this time of the year. However, 645 that study implicitly assumed that the wetlands were a final land cover stage rather than an incremental step toward the re-establishment of forest as depicted in the conceptual framework presented in Figure 6.

Hydrological Characteristics
As the land covers presented in the conceptual framework transition from one to the next, hydrological processes also change (Figure 6a). In the early stages (I, II), a relatively large proportion of hydrological inputs from the atmosphere are stored in collapse scars due to their impoundment by the permafrost on their margins (Connon et al. 2014). Evapotranspiration from 660 the landscape is relatively low given the high proportion of forest and relatively low transpiration by the black spruce that dominates the plateau canopies (Warren et al. 2018). In the early land cover stages (I, II) when forests predominate, understorey vegetation provide the pathway for evapotranspiration . The incremental change in land covers presented in Figure 6 involves biophysical changes that affect the partitioning of precipitation into storage or 665 runoff. By stage III to IV wetlands are interconnected and rapidly expanding, the storage of water on the landscape reaches its minimum level while runoff from the landscape is maximized ( Figure 6a). This increased runoff is enabled by the removal of permafrost barriers (Haynes et al. 2018) and areal expansion of runoff contributing areas resulting in greater hydrological connectivity and therefore drainage of the landscape (Connon et al. 2014). These observations 670 coincide with a period of steady, unchanging annual precipitation; therefore precipitation does not account for elevated basin runoff (Connon et al. 2014). A decrease in landscape drainage then follows in the subsequent stages as the transient runoff contributions from "captured" collapse scars diminishes as the importance of evapotranspiration increases as the wetlands become the predominant land cover (IV, V). The increase in evapotranspiration is due to 675 increases in evaporation from areas occupied by standing water and saturated or near-saturated wetland vegetation, including Sphagnum mosses, with losses due to transpiration driven by shrub vegetation (Warren et al. 2018). In the advanced stages (VI, VII) evapotranspiration would decrease as a result of the drier wetland surfaces as hummock microtopography replaces saturated Sphagnum lawns. The treed (afforested) wetlands (VII) have not been studied to the 680 same degree as peat plateaus or collapse scar wetlands (Haynes et al. 2020;Disher et al., 2021) and therefore ground based hydrological data specific to these features are lacking.

Conclusions
The discontinuous permafrost zone of the Taiga Plains exemplifies a landscape in transition. Coupling a broad-scale mapping initiative with the detail of site-specific data 685 collected in the Scotty Creek basin demonstrates a permafrost thaw-induced land cover transition. This transition is incremental and involves distinct land cover stages. The first and last of these is a continuous forest cover, although in the first stage the forest is underlain by permafrost while in the last stage it is not. Unlike traditional concepts of land cover change in peatland dominated regions of discontinuous permafrost in which forest re-establishment occurs 690 over centuries and is constrained by the rate of permafrost re-development, the concept presented here described forest re-establishment within decades and resulting from continued permafrost thaw, a process which allows wetlands to de-water sufficiently for tree growth. Each land cover stage has characteristic biophysical, hydrological and micro-meteorological features.
The proposed conceptual framework of landscape evolution describes the transitions 695 occurring across the Taiga Plains in peat plateau-collapse scar wetland complexes like Scotty Creek. This study also identifies the applicability of this conceptual framework across a large region of the Canadian north. We establish the likely pattern of change across these peat plateaucollapse scar wetland complexes and project their future trajectory by combining long-term field observations with analyses of contemporary and historical imagery. It is proposed that, while 700 permafrost thaw-induced land cover changes have previously been dominated by a transition from forest to wetland, this transition is not permanent and forested land covers are likely to return over time, although unlikely to be underlain by permafrost. This research improves the understanding of how peat plateau-collapse scar wetland complexes in the Taiga Plains may be impacted by ongoing permafrost thaw and these results may also be of relevance to other 705 peatland-rich permafrost environments across the circumpolar north.

Acknowledgements
We gratefully acknowledge the support of the Dehcho First Nations, in particular, the Liidlii Kue First Nation and Jean Marie River First Nation. We also thank these communities for their long-standing support of the Scotty Creek Research Station. This work was funded by 710 ArcticNet through their support of the Dehcho Collaborative on Permafrost (DCoP), and by the Natural Sciences and Engineering Research Council of Canada (NSERC). We also acknowledge the Canada Foundation for Innovation (CFI) for providing funding for infrastructure critical to this study.

Data Availability 715
The radiation flux data used in this paper are catalogued for open access in the Wilfrid Laurier University (WLU) Data Repository at https://doi.org/10.5683/SP2/JTIQDO.

Author Contributions
All authors contributed to the development of the research question and the methodological approach used in this study. OC and RC performed the analyses. OC, KH, and WQ wrote the 720 manuscript with input and editorial contributions from RC, JC, and ÉD.

Competing Interests
The authors declare that they have no conflict of interest.