Vulnerability assessment of peatland complexes in the Hudson Plains (Ontario, Canada) to permafrost‐thaw‐induced landcover and hydrological change using a multiscale approach

The Hudson Plains, Canada, is one of the largest, undisturbed peatland regions (370,000 km2) in the world. Air temperature in the Hudson Plains is increasing rapidly leading to unprecedented permafrost thaw. The region's remoteness has hindered our knowledge of how permafrost thaw alters peatland land cover and hydrological response at multiple scales. To assess which landscapes in the Hudson Plains are vulnerable to such disturbances, we analysed latitudinal distributions of land cover over a 300‐km transect spanning the sporadic (<30% areal) to continuous (>80% areal) permafrost zone in northern Ontario and quantified land cover changes over 40 years using multiple remote sensing datasets (lidar, air photographs, and high‐resolution satellite imagery). We then evaluated these landscapes at a fundamental hydrological unit, the peatland complex, identified five peatland complex types, and conceptualized their potential hydrological response using circuitry analogues. Over four decades, we found that permafrost peatlands declined by 4%, 8.5%, and 2% areal in the sporadic, discontinuous, and continuous permafrost zones, respectively. Circuitry analogues partitioned peatland complexes into their component peatland forms (e.g., permafrost peatland, bog, and fen) and represented each component's hydrological function using an electrical equivalent (e.g., generator, switch, and conductor). When interpreted at the landscape scale, circuitry analogues demonstrated latitudinal patterns in landscape structure (i.e., circuitry wiring) and indicated where permafrost thaw will have the greatest impact on landscape structure (i.e., rewiring) and therefore hydrological response. Based on these analyses, we suggest a 60‐km latitudinal segment (54.5°N to 54.9°N) where peatland complexes are most vulnerable to permafrost‐thaw‐induced land cover and hydrological change and should therefore be the focus of future research and monitoring efforts.


| INTRODUCTION
The 370,000-km 2 Hudson Plains ecozone in north-central Canada is the world's second-largest permafrost-affected peatland region ( Figure 1). This ecozone extends from continental boreal environments in permafrost-free terrain to low-arctic tundra underlain by continuous permafrost (>80% areal) along the south-western Hudson Bay coast (Harris et al., 1988). Between these southern and northern limits lies an expansive peatland-dominated region of sporadic (<30% areal) and discontinuous (>30% to <80% areal) permafrost, with relatively minor occurrences of mineral wetlands, deeper waterbodies (>2 m deep), upland forests, and barren lands (Riley, 2011). Peatlands consist of two predominant wetland classes, bogs and fens (National Wetlands Working Group [NWWG], 1988). The primary source of water to bogs is precipitation, whereas, for fens, precipitation is augmented by surface and or subsurface inputs (NWWG, 1988). The large assemblages of bogs and fens in the Hudson Plains are typically arranged into peatland complexes (Cajander, 1913), with each multisquare kilometre complex containing its own internal topographical, hydrological, and ecological gradients.
With increasing latitude, peatland complexes in the Hudson Plains are increasingly composed of nonwetland peatland forms underlain by permafrost (Riley, 2011). Such 'permafrost peatlands' include polygonal peat plateaus, peat plateaus, palsas, and inter-ridge bogs/fens (Riley, 2011;Sjörs, 1959;Zoltai & Tarnocai, 1974), which typically rise 1-2 m above adjacent permafrost-free terrain. Detailed descriptions of permafrost-free and permafrost peatlands are found in the Supporting Information (Table S1). As permafrost peatlands thaw, these features transform into permafrost-free peatlands because of the simultaneous processes of ground surface subsidence and inundation.
This transformation changes the flux and storage of water within and from the transformed peatland as preferential flow is no longer impeded by permafrost barriers. To address the question of the impact of climate warming on basin run-off in the Hudson Plains, a reasonable first step is to investigate the region's distribution of peatland complexes and how permafrost thaw is changing their component peatland forms and therefore hydrological function.
Little is known of the present distribution of permafrost peatlands in the Hudson Plains, how the distribution of these features changes with increasing latitude (e.g., McLaughlin & Webster, 2014;Riley, 2011), and how their presence and distributions have changed over recent decades (e.g., Kirkwood et al., 2019;Ou et al., 2016;Pironkova, 2017). Even less is known of how permafrost thaw affects the form and hydrological function of the peatland complexes of which they are a part. This is a significant dearth of knowledge considering that the Hudson Plains occupies the southernmost extent of nonalpine permafrost (Brown et al., 2002) and is therefore highly vulnerable to hydrological change resulting from permafrost thaw (Spence et al., 2020).
Several decades of field (Connon et al., 2014(Connon et al., , 2015Hayashi et al., 2004;Quinton et al., 2003), remote sensing (Carpino et al., 2018(Carpino et al., , 2021Chasmer et al., 2014;Chasmer & Hopkinson, 2017), and numerical modelling (Devoie et al., 2019;Devoie & Craig, 2020) studies in the southern Taiga Plains ecozone of north-western Canada have characterized specific permafrost thaw processes, defined their rate and pattern, and characterized their effect on the form and hydrological function of the peatland complexes in permafrost terrain. Permafrost thaw reduces the ability of F I G U R E 1 Map depicting the Hudson Plains study region, permafrost zonation, and transect locations. The lidar transect is shown in green and the gridded transect is shown in grey. Shown in red are the six 36-km2 areas of interest (AOIs). From south to north, AOI labels refer to Koper Creek (KPR), Muketei River (MUK), Ekwan River (EKW), North Washagami River (NWG), Sutton Lake (SLK), and Kinusheo River (KNS). The base map is a 'true colour' composite Sentinel-2A mosaicked image. Permafrost zonation is based on Brown et al. (2002) Circum-Arctic Map of Permafrost and Ground-Ice Conditions, Version 2. The subset map depicts the entire Hudson Plains domain (dark grey) within Canada and the nearest climate stations (blue dots) to our study area with 1981-2010 averages. From bottom left to bottom right in a clockwise direction, those climate stations are Pickle Lake, Gillam, Kuujjuarapik, and Moosonee. the landscape to obstruct drainage and impound water (Connon et al., 2014;Hayashi et al., 2004;Quinton et al., 2003) and introduces new surface (Connon et al., 2015) and subsurface Devoie et al., 2019) flow paths. As such, permafrost thaw increases the hydrological connectivity of peatland complexes (Connon et al., 2015), increases the hydrological response of streams (Connon et al., 2014;St. Jacques & Sauchyn, 2009), and partially drains areas of the landscape previously impounded by 'permafrost dams' (Haynes et al., 2018). However, a study by Mack et al. (2021) demonstrated that between 1970 and 2016, increased hydrological response was largely not seen in a cohort of peatland-dominated basins found in the sporadic and discontinuous permafrost zones outside the Taiga Plains, including four basins in the Hudson Plains.
Regionally slower permafrost thaw rates and more diverse land cover assemblages were possible explanations for the absence of any change. The knowledge gained from these studies provides valuable and detailed information directly applicable to the interpretation of changes to landscapes in the Hudson Plains observed through remote sensing.
The few remote sensing observations (e.g., Kirkwood et al., 2019;Ou et al., 2016;Pironkova, 2017) suggest that permafrost thaw is changing peatland forms throughout the Hudson Plains, and as a result, the region is potentially in a state of ecohydrological transition without precedent in the historical record. Aerial and/or satellite images provide snapshots of the spatial arrangement of peatlands and peatland complexes, and the comparison of images representing different times can therefore show permafrost thaw-induced changes to peatland forms and hydrological function.
As permafrost thaw changes the relative proportion of permafrost peatland forms and therefore the composition of peatland complexes, the flux and storage of water within and from such complexes must also change. In this study, we examine a latitudinal segment through the Hudson Plains using circuitry analogues and our knowledge of peatlands in the Taiga Plains to characterize peatland complexes and how their form and hydrological function are changing due to permafrost thaw. This will be accomplished through the following specific objectives: (1) identify permafrost peatland forms, their position within peatland complexes, and contemporary distribution; (2) compare present-day and historical imagery over four decades, analyse land cover changes, and determine a latitudinal distribution of peatland complexes; (3) perform a systematic landform survey to characterize present-day peatland types, forms, and complexes over our entire latitudinal transect; (4) using the results of objectives 1-3, interpret the hydrological function of peatland complexes using circuitry analogues to generate an initial map of peatland complexes and landscapes vulnerable to change due to thaw.

| CONCEPTUAL APPROACH
Several authors have characterized the hydrological functioning of peatlands based on their physical form (e.g., Galkina, 1964;Galkina & Kiryushkin, 1969;Ivanov, 1975;Sjörs, 1959). Ivanov (1975) extended this causal relation between form and hydrological function to multiple spatial scales including what he referred to as the relatively loworder 'microforms' within individual peatlands (e.g., hummocks and hollows), the middle-order 'mires' (i.e., individual peatlands), and highorder 'mire massifs' (i.e., peatland complexes). Accordingly, the hydrological functioning of an individual peatland is affected by the form and therefore function of its component microforms, and the hydrological functioning of a peatland complex (>10 0 km 2 to <10 1 km 2 ) by the form and function of its component peatlands. Likewise, the hydrological functioning of a peatland-dominated landscape (<10 2 km 2 ) by that of its component complexes, and the hydrological functioning of a region (>10 3 3 km 2 ) by the landscapes it contains.
By assigning specific hydrological functions of individual peatland types based on their form, we can infer hydrological flux and storage processes within and between peatlands and at the scales of peatland complexes and of landscapes containing multiple peatland complexes (i.e., 'compound mire massifs'; Ivanov, 1975). It has long been recognized that peatland complexes are deterministic features (Cajander, 1913;Ivanov, 1975;Sjörs, 1959;Vitt et al., 1994;Zoltai & Tarnocai, 1974) to the extent that their forms are shaped by the environmental conditions in which they developed (Waddington et al., 2015;Zoltai & Vitt, 1995). Therefore, as such conditions change, their forms also change (e.g., Kuhry, 2008), and presumably so too would their hydrological function. However, beyond consideration of how hydrological processes may change at specific locations is the need to understand how permafrost thaw affects the routing of water through a landscape of contrasting peatland forms and hydrological functions. Such scaling problems remain central questions in hydrology more broadly (Blöschl et al., 2019).
Hydrological functions are often conceptualized using electrical circuitry analogues (e.g., Dingman, 2002;Freeze & Cherry, 1979;Kirkham, 2005), because the flow of water and electricity are both governed by linear transport equations. Water flux is analogous to amperage in electrical terminology, storage to capacitance, and hydrological systems to circuits. Circuitry analogues offer insights into how water is routed through peatland landscapes, and how such routing will change as the landscape transitions in response to permafrost thaw ( Figure 2). The electrical analogue has little practical value in most applications since it is always preferable to characterize and predict water flow and storage directly from knowledge of hydrological processes and pathways. However, defining the broad circuitry over landscapes is a good first step to understanding hydrological functioning at larger scales and will also help identify knowledge gaps and needs. The electrical analogue also provides a framework to evaluate changes over time from historical aerial and satellite images. However, it must be stressed that while the analogue can ascribe general hydrological functions based on peatland form over large space and time scales, it is not intended to substitute information on specific processes and pathways arising from hydrological field studies. Rather it is used to identify general hydrological functions based on peatland form and position relative to other peatland forms. As such, it is F I G U R E 2 Components of peatland complexes in the Taiga Plains and Hudson Plains ecozones. For each landform, we show their hydrological function using a circuitry component symbol. For the Taiga Plains, all the components of the peatland complex (location: 61.36 N, 121.36 W) and circuitry analogue are known. For the Hudson Plains, there are additional peatland and circuitry components that occur. Hence, the peatland complexes and circuitry analogues are unknown. Note that for the Taiga Plains, both the peatland complex and circuitry analogue are present-day representations. Presumably, as the proportions of the peatland components change so too will the circuitry analogue.
intended to complement detailed information on peatland hydrological functioning derived from field studies.
By occupying the highest topographic position in peatland complexes and owing to their impermeable frost table, permafrost peatlands function as run-off generators with an active layer (depth) controlled resistor (Quinton & Baltzer, 2013;Wright et al., 2008).
While ephemeral connections between down-gradient collapse scar bogs, pools, and ponds are governed by water table thresholds, hence these forms function as water level switches (Kirpotin et al., 2009;Price & Maloney, 1994). In largely permafrost-free peatland complexes, bogs are often situated at the highest topographic position and are therefore run-off generators. However, the efficiency by which a bog generates run-off can be differentiated by whether it is patterned or nonpatterned. Patterned bogs have raised ridges and sunken hollows or pools arranged in repeating mosaics (Ivanov, 1975;NWWG, 1988;Riley, 2011) and therefore function as ridge-controlled run-off resistors. While nonpatterned bogs, which have microtopographic surfaces that are more similar to peat plateaus, function as water table (depth) controlled resistors.
The lower topographic position of fens relative to bogs ensures that run-off from peatland complexes flows through fens before entering the open-channel network (Price, 1987;Quinton et al., 2003). Like bogs, fens may be patterned or nonpatterned (Ivanov, 1975;NWWG, 1988;Riley, 2011). With more substantive microtopography, run-off from patterned fens encounters more barriers as it flows down-gradient, hence their hydrological function as run-off resistors (McCarter & Price, 2017;Quinton & Roulet, 1998).
By contrast, with fewer barriers to down-gradient flow, nonpatterned fens function as run-off conductors (Quinton et al., 2003;Roulet & Woo, 1988). As a result, assuming similar flow path lengths (Beven & Kirkby, 1979), drainage densities (Dingman, 2002), and meteorological conditions, fen form governs the efficiency of water fluxes from permafrost peatlands to the open-channel network and onward to the basin outlet. Together, the occurrence and proportion of each peatland form (e.g., permafrost peatland, bog, fen) and their specific types (e.g., peat plateau, collapse scar bog, and patterned [net] fen) ultimately determines the hydrological function (i.e., wiring) of a peatland complex.

| STUDY REGION
Our primary Hudson Plains study region consists of a 300 km (8550 km 2 ), from 52.8 N to 55.1 N latitudinal segment, which samples the sporadic, discontinuous, and continuous permafrost zones  Figure 3a. Briefly, we aimed to delineate elevated surfaces (≥ 1 m above surroundings), which are indicative of permafrost peatlands in an otherwise low-relief terrain.
We refined our initial classified image by manually removing elevated surfaces underlain by mineral substrates based on the Quaternary geology of Ontario shapefile (Pala et al., 1991). We then summarized our final classified image by 0.2 N latitudinal increments to facilitate our interpretation of the effects of latitude. We also used the classified lidar segment to cross-reference the visual appearance of permafrost peatlands on the other spectral imagery used in this study.
Accuracy was evaluated by counting intersections between our classified image and 'ground truths'. These consisted of 100 elevation profiles distributed along the entire transect where an elevated surface occurred concurrently with a landform that upon visual inspection was indicative of a permafrost peatland.

| Imagery acquisition and classification: six areas of interest (AOIs)
We selected six AOIs, evenly spaced along a 300-km north-to-south transect from Kinusheo River in continuous permafrost to Koper Creek in sporadic permafrost to focus our analyses (Figure 1). Sutton Lake, North Washagami River, Ekwan River, and Muketei River com- Canada. Each AOI included distinct watercourses and/or waterbodies to aid with georeferencing photos in landscapes without static manmade structures (see Chasmer et al., 2010). The metadata for each satellite image and air photo is found in Table 1.
Initially, we manually classify land cover forms in each air photo and satellite image ( Figure 3b). Where manual classification was deemed impractical (e.g., numerous waterbodies) or unsatisfactory (producer defined), we proceeded with algorithm-based classifications. First, we applied an iso-cluster unsupervised classification using 10-20 initial classes. If results were still unsatisfactory, a maximum likelihood supervised classification using between 10 and 50 training areas was implemented. For Kinusheo River, we used a segment mean shift object-based classification with spectral detail and spatial detail parameters both set to 20. We then manually edited all manual and algorithm-based classifications until we attained a satisfactory result for each photo and image. A complete accuracy assessment for each AOI classification is found in the Supporting Information document and includes an estimate of overall accuracy, the Kappa coefficient, confusion matrix, ground truths, commission, omission, user's, and producer's accuracies. To assess land cover changes, rather than directly subtracting our two classified images (i.e., exact spatial difference), we quantified proportional changes to each land cover form, and therefore did not calculate statistical significance. This approach was chosen because the incongruent alignment of air photos and satellite imagery where no static man-made structures exist results in greater georeferencing error (Chasmer et al., 2010). Potential land cover forms for each AOI were water, forest, nonpatterned fen, pat-

| Imagery acquisition and survey: gridded transect
To better understand the land cover patterns we observed among our AOIs, we also estimated land cover types across an overlapping transect using 10-m resolution Sentinel-2A imagery (Table S3) acquired using Google Earth Engine (Gorelick et al., 2017) (Figure 3c). Separately, in consultation with the Thermokarst Mapping Collective (see Gibson et al., 2021), we extended their 7.5 km Â 7.5 km mapping grids covering the Northwest Territories (NWT), Canada to the Hudson Plains. The collective aims to create a systematic and broadscale database describing the variation of thermokarst and periglacial features across the NWT using a grid-based framework and Sentinel-2A imagery. We concluded that applying a similar grid-based approach was favourable to performing a large-scale algorithm-based classification because our study was specifically interested in understanding F I G U R E 3 Simplified workflows depicting the major processing steps for analysing the (a) lidar, (b) six areas of interest, and (c) gridded transect.
the Hudson Plains landscape at a higher-order (i.e., peatland complexes) and structural perspective (i.e., patterned or nonpatterned) for which such classification techniques are not well suited.
Our gridded transect consisted of 152 grid cells centred on and inclusive of our six 36-km 2 AOIs (Figure 1). For each grid cell, we systematically moved around a 1% quadrat (750 m 2 ) as a reference, then using the Sentinel-2A imagery we estimated the proportion of each land cover type (Figure 3c). Land cover types were determined using available literature (Sjörs, 1959;Ivanov, 1975;NWWG, 1988;Radforth, 1969;Riley, 2011;Zoltai & Tarnocai, 1974;Zoltai & Vitt, 1995) and, where available, georeferenced oblique air photos (N = 522). Furthermore, our classified lidar image was used to crossreference the visual appearance of permafrost peatlands relative to uplands on Sentinel-2A imagery. The minimum proportion assigned to any land cover type was 1%. Land cover types with greater proportions were estimated to the nearest 5% increment (i.e., 5% and 10%).
Descriptions of each peatland form and type are provided in Table S1.

| Mapping peatland complex and landscape vulnerability
We defined vulnerability based on Chapin et al. (2004), as systems more sensitive to change following subjection to a disturbance. Therefore, in each grid cell, we evaluated the vulnerability of the most prevalent peatland complex to increased hydrological response (i.e., more sensitive to change) as permafrost thaws (i.e., a disturbance) over the coming decades. The vulnerability was ranked into four categories, most, more, less, and least using a decision tree (Figure 4). Grid cells were ranked as most vulnerable where peatland complexes were composed of ≥35% permafrost peatland areal and where nonpatterned fens were the most prevalent fen form (i.e., Taiga Plains like). Grid cells were given the rank of least vulnerable where permafrost peatlands account for <35% areal. We chose 35% as a conservative first approximation for minimum permafrost peatland area for landscapes to be vulnerable to hydrological change in the Hudson Plains because the percentage reflected estimates from the southern sporadic permafrost zone of the Taiga Plains where such changes occurred (Carpino et al., 2021). We assigned intermediary vulnerability ranks based on whether the most prevalent fen forms in each grid cell were patterned (terraced) fens or patterned (net) fens, which corresponded to the ranks of more and less vulnerable, respectively. Assigning a higher vulnerability to patterned (terraced) fens acknowledges their steeper slopes and more efficient flow paths relative to patterned (net) fens (Table S1).
T A B L E 1 Description of each 36-km 2 Hudson Plains area of interest (AOI) along the 300-km interrupted transect. Names are based on the nearest known waterbody or watercourse.

| Latitudinal trends of permafrost peatland area estimates
Our lidar classification demonstrated that permafrost peatlands were present throughout the sporadic ( Figure 5a) and discontinuous permafrost zones (Figure 5b) but occupied only a minor area within the overall landscape. We observed that permafrost peatlands were positioned either as 'islands' within larger permafrost-free peatland complexes or on the periphery of watercourses and waterbodies. For

| Latitudinal and temporal variation of peatland forms
Classifications of recent (2014, 2016, and 2017) imagery demonstrated several distinct latitudinal patterns at the peatland form hierarchical level. Bogs were the most prevalent peatland forms in four out of six AOIs ( Lake and North Washagami River, but without an ability to confirm F I G U R E 4 Decision tree used for mapping peatland complex and landscape vulnerability to permafrost-thaw-induced hydrological change.
using ground verification or lidar classification, those landforms were classified as a nonpatterned bog.
Since the mid-1970s, forest cover (a proxy for permafrost peatlands) decreased and nonpatterned bogs increased by less than 5% of the total areas, in the sporadic zone and southern part of the discontinuous zones ( Figure 6). However, in Koper Creek, our most southern AOI, several treed permafrost peatlands present in 1975 were completely degraded by 2014 ( Figure S1f). In contrast, forest cover decreased by 9% to 16% in the northern discontinuous permafrost zone. In Kinusheo River, our comparison showed losses to water (À1%) and nonpatterned fen ( Figure S2).
In addition to permafrost-thaw-induced disturbances, we also observed the impacts of wildfires in Sutton Lake and North Washagami River. Using Landsat 5 TM imagery to highlight recently burned vegetation and the most up-to-date National Fire Database (NFDB, 2022) we estimated the total area burned and dates for each wildfire ( Figure S3) and amended Table 2 and Figure 6 to include burned area (black). In Sutton Lake, a 1.9-km 2 wildfire occurred in Based on our burned area analyses, the remaining 78% (2.5 km 2 , 7% areal) and 60% (6.4 km 2 , 10% areal) of all forest losses in Sutton Lake and North Washagami River, respectively, resulted from permafrostthaw-induced disturbances. For each classified AOI image (N = 12), the overall accuracies and Kappa coefficients were between 79% and 86% (Table 3).

| Identification of peatland complex types
Using the most prevalent run-off generating and run-off conducting (see Figure 2) peatland forms (see Table S1) classified in each AOI, we identified the following peatland complexes (Figure 7a

| Latitudinal land cover trends and comparison to classified AOIs
Within our 8550-km 2 gridded transect open, nonpatterned fens (see Table S1), rapidly transitioned from more than 50% of the total area at 55.0 N to less than 5%at 54.7 N. (Figure 7). Open, nonpatterned fens were surrounded by open, inter-ridge bogs/fens and peat plateau fields. These peatland types, along with the numerous pools and ponds form widely dispersed often multisquare kilometre heterogeneous peatland complexes. Open, inter-ridge bogs/fens and peat plateau fields contributed nearly 50% to the total land cover area and were most common near the continuous-discontinuous transition zone (54.7 N), and steadily declined southward, stretching into the sporadic permafrost zone (53.8 N). Over the same latitudinal range, we observed large areas of forest and peatland, which showed signs of past wildfires, accounting for nearly 50% of the total land cover at 54.3 N. Treed peat plateaus contributed nearly 20% of the total land cover area at 54.3 N and were the only permafrost peatland type to occupy areas >1% south of 53.8 N. Thaw disturbances/collapse scars (Harris et al., 1988) occurred throughout the transect adjacent to many permafrost peatlands and occupied a maximum of 10% of the total land cover at 53.5 N. Overall, we found that open, nonpatterned bogs and treed, patterned (terraced) fens (see Table S1) were the most prevalent permafrost-free peatland types and, therefore, were important components of our identified peatland complexes. Open, nonpatterned bog area contributed nearly 45% to the total land cover area at 53.6 N, while further south, treed, patterned (terraced) fen area contributed nearly 20% to the total land cover area at 53.3 N. Comparing all surveyed and classified land cover forms (i.e., one hierarchal level above land cover type) and within each AOI, correlations were all greater than 0.76 (Table 4).

| Apparent hydrological function
The hydrological function of peatland complexes has previously been conceptualized by applying the fill-and-spill concept to permafrost (Connon et al., 2015) and permafrost-free (McCarter & Price, 2017) settings. Fill-and-spill described run-off controlled by a stepwise expansion of contributing areas as storage thresholds in individual F I G U R E 7 A 'true colour' (red, green, blue) composite satellite image for each Hudson Plains area of interest (see Figure 1) using the multispectral bands from their respective satellite platform (see Table 1). Each (6 km Â 6 km) AOI is labelled by the primary peatland complex. Note that Sutton Lake, North Washagami River, and Ekwan River images were each rotated 2 , 1.5 , and 1.5 counterclockwise, respectively, to align with the other AOIs.

T A B L E 4 Correlation summary between classified and surveyed
Hudson Plains land cover forms and among areas of interest. Land cover types from survey results were combined into their land cover forms prior to calculating correlations. Correlations were statistically significant at p values <0.01 (*), <0.001 (**), and <0.0001 (***).  (Spence & Woo, 2003. In addition to the abovementioned fill-and-spill limitations, fieldbased hydrological processes studies are spatially limited, hence, we suggest applying circuitry analogues to understand the potential hydrological vulnerability of peatland complexes as permafrost thaws at the landscape and regional scale. Although circuitry analogues cannot be used to provide insights into specific hydrological processes or interannual hydrological variability within peatland complexes and their component forms, they do demonstrate latitudinal patterns in landscape structure (i.e., circuitry wiring) and how such structure may respond to a disturbance (i.e., permafrost thaw).
The high correlations between our classified and surveyed results (Table 4)

| Vulnerability map
Using the decision tree outlined in Figure 4, our initial map of landscape vulnerability to permafrost-thaw-induced hydrological change in the Hudson Plains is shown in Figure 10a and suggests a 60-km latitudinal segment (54.5 N to 54.9 N) where peatland complexes are most vulnerable. We assert that each 56.25-km 2 grid cell represents the landscape scale because they contain multiple peatland complexes ( Figure 1). Landscapes where complexes are characterized as permafrost peatland-fen-pond or permafrost peatland-bog-fen corresponds to the most, more, and less vulnerable ( Figure 10b). All other landscapes are ranked as least vulnerable and occur mainly in the sporadic permafrost zone and parallel the lidar transect ( Figure 5).
Because fen form controls run-off efficiency (Quinton et al., 2003;Quinton & Roulet, 1998), among the permafrost peatland complexes the most prominent fen form further discriminates between landscapes most vulnerable (non-patterned), more vulnerable (patterned and terraced), and less vulnerable (patterned and net) (see Table S1).
An exception is found in the continuous permafrost zone where nonpatterned fens are the most prevalent of all peatland forms (Figure 8; F I G U R E 8 Area plot of Hudson Plains land cover types. Land cover is summarized by grid row (four grid cells per row). The locations of the six AOIs are shown with solid black line boxes and the approximate boundaries between permafrost zones are shown with dashed lines. An estimated 5%-10% uncertainty exists across all grid cells (white area). Figure 10b). In such landscapes, we substantially reduced vulnerability from most to less (Figure 10a), because the development of new flow paths as permafrost peatlands thaw would be spatially limited (e.g., Figure 7a; Figure S2).

| Latitudinal land cover trends and temporal variation
Concerning the distribution of permafrost peatlands in the Hudson Plains and quantifying the rate and pattern at which these landforms are thawing, previous studies and this study suggest that in the ecoregion's sporadic permafrost zone, permafrost peatlands occupy <5%-15% areal ( Figure 5), consist primarily of treed permafrost peatlands ( Figure S1), and have thawed between 0.08%-0.79% decade À1 ( Figure 6). One study by Pironkova (2017) (Figure 6), which fall within the sporadic permafrost zone.
We posit that the range of permafrost area estimates reflects the uneven distribution of permafrost peatlands across the sporadic permafrost zone as seen in our lidar analysis ( Figure 5). However, given that the longitudinal width of our lidar surface ranged from 600 to 900 m, sampling error was certainly a contributing factor to our estimates. The larger losses in forest cover in this study may reflect a sampling error as well because we selected AOIs near distinct water bodies for georeferencing purposes where thaw rates were highest in the Pironkova (2017)  For example, in the sporadic and discontinuous zones of Québec's Taiga Shield and Southern Arctic ecozones, Payette et al. (2004) found that 80% of the permafrost in an inter-ridge bog/fen had thawed over 47 years , and further north, Vallée and Payette (2007) found that 23% of palsa areas had thawed over 44 years . Similar latitudinal patterns and even greater permafrost thaw rates have been described in Lapland (Luoto & Seppälä, 2002. In the Taiga Plains, Carpino et al. (2018) found that forest area losses between 1970 and 2010 along a 200-km latitudinal transect in the sporadic permafrost zone ranged from 6.9% to 11.6%, with greater losses observed further north.  (Figure 10a).
F I G U R E 1 0 (a) Vulnerability map of peatland complexes in the Hudson Plains to permafrost thaw-induced hydrological change annotated with the locations of the lidar transect, six 36-km 2 AOIs, and permafrost zones. (b) Distribution of the most prevalent peatland complexes, fen types, and peatland types in each grid cell, which together with the decision tree (see Figure 4) were used to determine vulnerability in each grid cell.

| The Hudson Plains and Taiga Plains compared
Studies in Taiga Plains have shown that to increase hydrological response, peatland complexes must be initially composed of a substantial area of permafrost peatlands, permafrost thaw must form enough new hydrological pathways to increase the contributing area of the landscape, and fens must not inhibit new pathways from reaching the basin outlet. For example, Connon et al. (2014) showed that a 33% increase in contributing area (from 31.1% to 41.8%) and a 137 mm increase in run-off occurred over two decades in a landscape initially composed of 54.6% permafrost peatlands rapidly decreases to 47.5% (À7.1%) over 34 years . Areal permafrost losses in Sutton Lake (À10%) and North Washagami River (À7%) were like those measured by Connon et al. (2014); however, neither landscape was initially composed of >40% permafrost peatlands and each contained more patterned than nonpatterned fens ( Figure 6). These structural differences highlight their reduced vulnerable rankings (Figure 10a). Despite only a small proportion of peatland complexes in our study area resembling those in Taiga Plains, that is, permafrost peatland-bog-fen (concentrated amperage) complexes, we assert that this does not preclude permafrost peatland-fen-pond (dispersed amperage) complexes from undergoing increased hydrological response as permafrost thaws, so long as they meet the conditions described above. However, for either peatland complex, what remains uncertain and what circuitry analogues cannot determine is which complexes will connect to the larger drainage network and which to isolated waterbodies as permafrost thaws (see Kirpotin et al., 2009).

| CONCLUSION
We characterized peatland complexes along a 300-km latitudinal segment in the Hudson Plains using multiple remote sensing approaches, conceptualized their hydrological function using circuitry analogues, and assessed their vulnerability to permafrostthaw-induced hydrological change based on those conditions and estimated land cover changes. Out of the five peatland complexes we identified, two had sufficient permafrost peatland proportions (≥35%, areal) to suggest some vulnerability to permafrost-thawinduced hydrological change, permafrost peatland-fen-pond (dispersed amperage) and permafrost-bog-fen (concentrated amperage). We suggested the most vulnerable landscapes occur where nonpatterned fens were the most prevalent fen form, not the most prevalent peatland form. Based on these parameters, we identified a 60-km latitudinal segment (54.5 N to 54.9 N) where landscapes are most vulnerable and should therefore be the focus of future research and monitoring efforts in the Hudson Plains of northern Ontario. Conceptualizing the hydrological function of peatland complexes using circuitry analogues indicated which complexes and landscapes in the Hudson Plains were most vulnerable to increased hydrological response as permafrost thaws.
However, because each circuitry analogue represents a nonspecific peatland complex, the analogues cannot be used to determine if individual peatland complexes were vulnerable or used as routing networks in hydrological models. This study has highlighted an ongoing challenge, which is the juxtaposition between the large scale (10 2 km 2 ) of the Hudson Plains which necessitates remote sensing approaches and the heterogeneity of peatland components over small scales (10 À3 km 2 ) which impedes classification of individual permafrost peatlands without the use of high-resolution satellite imagery. Furthermore, because patterned peatlands are composed of several land cover types, algorithm-based classifications are ineffective at delineating their bounds. We know that as permafrost thaws, the hydrological function of individual permafrost peatlands changes. At some critical threshold, hydrological changes aggregate so that the hydrological response of entire peatland complexes, landscapes, and regions changes. Therefore, future work should aim to decipher peatland complexes that may connect to the drainage network as permafrost thaws, and which will not.