Changes in surface water dynamics across northwestern Canada are influenced by wildfire and permafrost thaw

The abundance and distribution of surface water at high latitudes is shifting rapidly in response to both climate change and permafrost thaw. In particular, the expansion and drainage of lakes and ponds is widespread but spatially variable, and more research is needed to understand factors driving these processes. In this study we used medium resolution (30 m) remote sensing data to analyse changes in lake area in permafrost-rich lowland regions across northwestern Canada. First, we used the Global Surface Water Dataset developed by the GLAD research group to map the absolute area of different land–water transitions across a 1.4 million km2 study domain. Next, we selected six regional study areas representing a range of climatic, geologic and hydrologic conditions. Within these regional study areas, we used the Landsat satellite archive to map annual trends in the area of 27 755 lakes between 1985 and 2020. We trained a random forests model to classify lakes exhibiting significant increasing or decreasing trends in area, and assessed the relative importance of climate, disturbance and environmental variables in determining the direction of change. Our analysis shows that significant increases in lake area were 5.6 times more frequent than decreases during the study period. Wildfire and ground ice abundance were the most important predictors of the direction of change. Greater ground ice content was associated with regions that experienced increases in lake area, while wildfire was associated with regions that experienced decreases in lake area. The effects of climate, including trends in mean annual temperature and total annual precipitation were smaller than disturbance and environmental factors, indicating that climate has likely had indirect effects on lake area changes over our period of study.


Introduction
The abundance and distribution of lakes and ponds are changing rapidly across the circumpolar north (Smith et al 2005, Watts et al 2012, Nitze et al 2017, Pastick et al 2019. Lake drainage and expansion are widespread and one of the most well-documented hydrologic changes in regions underlain by permafrost (Grosse et al 2013, Kokelj and Jorgenson 2013). Several remote sensing studies have attributed decreases in the area of lakes and ponds to changes in temperature and precipitation (Smol and Douglas 2007, Plug et al 2008, Lantz and Turner 2015, Campbell et al 2018, Swanson 2019). In the high-Arctic, rising air temperatures and increasing evaporation have driven drying within small ponds (Smol andDouglas 2007, Campbell et al 2018). In other regions, rapid thermokarst lake drainage has been attributed to intense precipitation events, as high water levels can lead to bank overtopping and the formation of outlet channels via ground-ice thaw (Plug et al 2008, Lantz and Turner 2015, Swanson 2019. Drivers of lake expansion are generally less studied, but likely include increased precipitation and changes in water balance (Olthof et al 2015, Korosi et al 2017. Additional investigations are needed to understand the heterogeneous response of lakes to ongoing climate change.
The direct effects of climate on lake area (including increasing evapotranspiration and precipitation), are likely mediated by a range of lake and environmental factors including lake size, lake depth, catchment characteristics Bigras 1988, Campbell et al 2018), underlying ground-ice conditions (Yoshikawa and Hinzman 2003, Swanson 2019, Jones et al 2020, surficial geology (Nitze et al 2017, Carroll and Loboda 2018, Wang et al 2018 and vegetation cover (Connon et al 2014, Turner et al 2014. These environmental factors likely influence surface water dynamics by controlling the potential for thermokarst activity , Jones et al 2020, the permeability of surficial materials (Wang et al 2018) and runoff and snow melt intensity (Connon et al 2014, Turner et al 2014. However, the relative importance and interaction of different climate and environmental factors at broad spatial scales is unclear. To address this uncertainty, trends in the number and size of lakes have been investigated using medium (30 m) and high (<1 m) resolution remote sensing imagery (Jones et al 2011, Bouchard et al 2013, Nitze et al 2017, Cooley et al 2019. Recent studies showing high interannual and seasonal variation in lake area as well as complex non-linear changes highlight the need to conduct analyses at multiple spatial and temporal scales.
In regions with ice-rich permafrost, thermokarst processes can also have a significant impact on the area of lakes and ponds (Grosse et al 2013, Kokelj and Jorgenson 2013, Lantz and Turner 2015, Fraser et al 2018, Nitze et al 2020. Rising air temperatures, changes in the amount and intensity of precipitation, as well as more frequent wildfires, are increasing the frequency of thermokarst disturbances across the circumpolar (Kasischke and Turetsky 2006, Hu et al 2015, Segal et al 2016, Fraser et al 2018. Wildfire can impact permafrost stability by increasing thaw depth and ground temperature, which in turn can alter flow paths and hydrological fluxes (Yoshikawa et al 2002, Jones et al 2015, Narita et al 2015, Zipper et al 2018, Holloway et al 2020. In northern boreal regions, wildfire has been shown to initiate a range of thermokarst processes including the formation of thermokarst bogs, fens, and the degradation of icewedges, which all impact the hydrological connectivity of the landscape (Jones et al 2015, Gibson et al 2018, Zipper et al 2018. Previous studies have also observed decreases in lake area following wildfire (Roach et al 2013, Frost et al 2019. In this study, we mapped changes in surface water across northwestern Canada and examined the effects of climate, disturbance and environmental factors on increases and decreases in surface water area. To accomplish this, we used the Global Surface Water Change dataset developed by Pickens et al (2020) to summarize surface water change across a 1.4 million km 2 study domain encompassing the western Canadian Arctic. Next, we selected six study areas representing ecological and geologic gradients across lake-rich lowlands likely to be sensitive to changes in permafrost stability. In these study areas, we mapped interannual changes in the area of 27 755 lakes and used a Random Forests classification model to investigate the relative influence of disturbance, environmental, and climatic factors on significant increases and decreases in lake area. Due to the increasing frequency of wildfire in these areas (Kasischke and Turetsky 2006), we also investigated the impact of fire on the timing and magnitude of change in lake area. To our knowledge, this represents one of the first studies documenting interannual variation in lake area across climatic and environmental gradients at broad spatial scales in northwestern Canada.

Study domain
Our study domain covers approximately 1.4 million km 2 of northwestern Canada, including the Northwest Territories, Yukon and Nunavut (figure 1). Climate, vegetation and surficial geology vary across the study domain, which is divided into four primary ecoregions: the Southern Arctic, the Taiga Plains, the Taiga Shield and the Taiga Cordillera.
The Southern Arctic extends northward from the limit of continuous forest to the coast of the Beaufort Sea and Amundsen Gulf. It is characterized by dwarf and low-shrub tundra and is primarily underlain by a combination of glacial till, glaciofluvial and lacustrine sediments (Ecosystem Classification Group 2012). The permafrost is continuous and typically ice-rich due to the abundance of wedge, relict and segregated ground ice (O'Neill et al 2020b). The Taiga-Plains covers the central portion of the study area and is characterized by spruce and pine forests; forest fire is an important driver of vegetation change in this ecoregion (Ecosystem Classification Group 2009). Extensive wetland complexes containing peatlands and thaw lakes are also common (Ecosystem Classification Group 2009). There is an abundance of fine-grained tills across this region, which contains ice-rich continuous and discontinuous permafrost. The Taiga Shield is located in the northeastern part of the study domain and is characterized by a mosaic of Precambrian bedrock, glacial till and glaciolacustrine deposits. Much of this region was once covered by glacial Lake McConnell, which drained and deposited fine-grained lacustrine sediments throughout the region (Smith 1994). Areas with these deposits have the potential to host ice-rich permafrost, whereas bedrock regions are not thaw sensitive (O'Neill et al 2020b). Vegetation in this area varies from closed spruce, pine and aspen forests in the south, to isolated spruce stands in the north  Within the larger study domain we selected six smaller study areas (figure 1). Regional study areas were selected across a range of climatic, geologic and permafrost conditions within the ecoregions described previously (table 1). We restricted study areas to regions that were lake dense and close to a meteorological station. Shadowing effects in steep mountainous regions of the Cordillera make it difficult to detect surface water, especially when lakes are sparse, thus we focused our analysis on low lying, lake-rich regions. Detailed descriptions of each study area are presented in the supplemental materials.

Climate data
We acquired daily meteorological data from 1985 to 2020 from Environment and Climate Change Canada (2021) using the R package weathercan   across the study domain in the western Canadian Arctic and Sub-Arctic. Change in mean annual temperature was assessed for the full study period , but change in total rainfall was limited to available data from 1985-2015 at stations marked with an asterisk ( * ). Temperature and precipitation data from Tuktoyaktuk also begins in 1995. All historic climate data is from Environment and Climate Change Canada (2021).

Station
Temperature (LaZerte and Albers 2018). Average temperature and total rainfall were aggregated annually for each station, and years were excluded if they were missing more than one month of data. Trends in average annual temperature were calculated using ordinary least squares (OLS) regression (table 2). Scatterplots showing OLS regression are presented in supplemental figure 1. Climate normals for mean annual temperature and total annual rainfall based on the 1981-2010 period were also acquired from Environment and Climate Change Canada (2021).

Changes in surface water area 2.3.1. Broad-scale patterns
We used the Global Surface Water Dynamics (1999-2020) data product developed by the Global Land Analysis and Discovery group (Pickens et al 2020) to summarize patterns in surface water dynamics across the study domain. This dataset was developed using the Landsat satellite archive at 30 m resolution and classifies pixels into the following categories: stable water, stable land, seasonal water and one of five dynamic land-water transitions (Pickens et al 2020).
The dynamic classes, represent pixels that change between water and land cover based on interannual changes in water percent. Pickens et al (2020) define water percent as the number of clear-sky observations classified as water in a given year. Annual water percent was normalized by season to ensure that seasons with more observations did not bias the results.
In this analysis, we focus on the classes representing land-water transitions, which include (a) water loss (decreasing water percent), (b) gain (increasing water percent), (c) dry period (a bi-directional change representing decreasing then increasing water percent), (d) wet period, (a bi-directional representing increasing then decreasing water percent) and (e) high frequency, (representing more than three transitions between land and water). To map spatial patterns in these water transition classes, we calculated the absolute area covered by each class within 10 × 10 km grid cells covering the study domain. Because water bodies within the same grid cell may exhibit divergent types of change (i.e. drying in one part and wetting in another), a single grid cell may be characterized by multiple water transition types.

Regional patterns
We chose six regional study areas in lake-rich lowlands underlain by continuous and discontinuous permafrost that are likely to be sensitive to climate driven changes in permafrost stability (figure 1).
Within each study areas we used the Landsat satellite archive to map interannual changes in the area of 27 755 water bodies larger than 0.03 km 2 .

Landsat processing and sub-pixel water fraction
In each study area we acquired all available July and August top of atmosphere Landsat scenes from 1985 to 2020 with less than 30% cloud cover (Dwyer 2019). We restricted our acquisitions to Tier 1 Collection 1 images that are most suitable for time series analysis (Dywer 2019). These images have minimal geometric error and have been radiometrically calibrated for consistency across different Landsat missions (Dwyer 2019). We used the quality assessment band to mask clouds and cloud shadows and each scene was visually inspected and discarded if smoke or clouds remained in the study area. We used the histogram breakpoint method (Olthof et al 2015) to map annual coverage of surface water. This method uses thresholds in the shortwave infrared reflectance band (SWIR1) to quantify the proportion of a pixel covered by water (sub-pixel water fraction). Details about this method can be found in the supplementary materials. We created a raster of maximal water coverage in each study area using pixels where sub-pixel water fraction was greater than 50% in at least two scenes. We created a mask of all lakes in the study area by converting  (1995) Fire presence/absence as intersection of lake boundary with fire perimeter Fire perimeters below treeline in NWT and Yukon Geomatics Yukon (2019), NWT Centre for Geomatics (2020) Lake size between 1985 and 1990 Initial lake area the water coverage raster to polygon features. We removed ponds smaller than 0.03 km 2 and water bodies that were not completely inside the study area from analysis. The lake mask was buffered by 30 m and lake area was calculated as the sum of sub-pixel water fraction (ranging from 0 to 1) within each lake polygon multiplied by the pixel resolution (900 m 2 ). This method of mapping surface water has been previously validated within the Lower Mackzenzie Plain study area (Travers-Smith et al 2021) and in two tundra sites within the Tuktoyaktuk Coastal Plain and on Banks Island (Olthof et al 2015, Campbell et al 2018. Lakes were also classified into the following size categories (a) large lakes (>0.5 km 2 ), medium lakes (0.05-0.5 km 2 ), and small lakes (0.03-0.05 km 2 ) based on the average area between 1985 and 2020 (Supplemental table 1).
To quantify interannual change in total lake area we used a subset of lakes where Landsat imagery allowed estimates of area in more than half of the available Landsat acquisitions between 1985 and 2020. To minimize the impact of missing data in the time-series due to clouds or smoke, we calculated a five-point rolling mean of the area of each lake. We used this data to estimate annual change in total lake area, and change by size class. Next, we classified each lake as increasing, decreasing or non-trended in area over time using a Mann-Kendall Trend test and a significance threshold of 0.01.

Random forests
We used a random forests classification model (Breiman 2001) to classify lakes exhibiting significant increasing and decreasing trends in annual area within the six regional study areas. The explanatory variables used in this model are described in table 3. Random forests models were created in R using the random forests package (Liaw and Wiener 2002). Before training the model, we used the caret package in R to upsample our training dataset to ensure equal representation among the six study areas Upsampling randomly replicated the data so that each study area had an equal number of lakes. Each decision tree in the random forests model was then trained on an equal number of lakes showing increasing and decreasing trends. After training the model, we estimated the relative importance of climate, disturbance and environmental variables using the unscaled mean decrease in accuracy, which is calculated as the mean decrease in model accuracy when a variable is excluded. We also generated partial dependence plots showing the marginal probability of a lake belonging to the increasing or decreasing class across the range of the explanatory variables using the pdp package in R (Greenwell 2017).

Impact of fire on lake area
To explore the timing and magnitude of any changes in lake area following fire, we selected one fire that burned between 1985 and 2020 that affected an area at least 10 000 m 2 from four study areas below the treeline (Lower Mackenzie Plain, Central Mackenzie Valley, Bulmer Plain, Great Slave Upland regions); (figure 2). Most large fires above the treeline were not included in the fire database (Geomatics Yukon 2019, NWT Centre for Geomatics 2020), so we were unable to assess the possible impacts of fire in the Tuktoyaktuk Coastlands and Eagle Plain study areas. For each study area, we used boxplots to compare annual changes in lake area for lakes inside the burned area Fire history within the four study areas located below the treeline. Fire history data is from the NWT Centre for Geomatics (2019). Dates show the timing of fires within each study area. Fires outlined in red were used to determine the timing and magnitude of changes in lake area. and lakes outside the burned area, but within 10 km of fire boundaries.
It is possible that increases in SWIR reflectance driven by post-fire vegetation recovery contributed to decreases in sub-pixel water fraction along lake margins (Frazier et al 2018). To test the effect of fire on lake area trends independent of mixed pixels (pixels containing both water and land), we calculated annual lake area only using pixels with sub-pixel water fraction greater than 0.99. We identified lakes showing significant trends in area (p < 0.01) using Mann-Kendall trend detection, and then used a chi-square test to compare the frequency of lake area trends inside and outside of fire boundaries. We conducted this test using lake trends and fire data from the Lower Mackenzie Plain.

Broad scale patterns
Over the entire study domain, pixels representing water losses covered 1121 km 2 (0.4% of the area of permanent water) and water gains covered 3002 km 2 (1.2%). The area representing wet periods (increasing water percent followed by decreasing water percent) was 18 557 km 2 (7.4%) and the area representing dry periods (decreasing water percent followed by increasing percent) was 1536 km 2 (0.6%). Maps showing the spatial distribution of these classes show that water losses and dry periods were clustered in the Yellowknife area, parts of the Mackenzie Delta and within the Mackenzie River Valley (figure 3). Water gains were generally distributed throughout the eastern Taiga Shield, the Mackenzie Delta and within the northern and central Mackenzie River Valley (figure 3). The class representing wet periods was the most common and was distributed among most of the study domain, except the Cordillera ecoregions where lake density and absolute change in water area was low among all classes.

Regional change in lake area
Between the beginning and end of the study period, four of the six regional study areas showed an increase in the total area covered by lakes, with an average change of +1.3% (table 4). We observed the greatest percent increase in lake area in the Tuktoyaktuk Coastal Plain (2.7%), while the Lower Mackenzie Plain was the only study area that showed a decrease in lake area (table 4). Across all six study areas 65% of lakes were classified as non-trended in area, 29% were classified as increasing and 6% were classified as decreasing (table 4). All study areas Table 4. Absolute and relative change in total lake area and the proportion of lakes showing significant increasing, decreasing and non-trended patterns in each study area. Absolute change in area was calculated as the difference between total lake area in the first and last five years of the time series and relative change was calculated as the absolute change divided by total area in the first time period. The direction and magnitude of interannual change also varied across different study areas and by lake size (figure 4). Small lakes generally showed the largest absolute change in area, but in most study areas, the overall trend in small lakes was mirrored in larger lakes. In the Eagle Plains, total lake area increased by 2.4%, while the area of small lakes increased by approximately 10% (figure 4). In the Lower Mackenzie Plain and Great Slave Upland study areas, the area of small lakes increased by 5% and 7%, respectively, while the area of large lakes decreased by 1% and 0.4% ( figure 4). Overall, the patterns in interannual lake area reflected in the analysis of broad-scale patterns showing that most regions gained surface water area. While the broad-scale patterns suggested that wet periods were more common than dry periods, our regional analysis suggests that both types of change have occurred, and that decreases in lake area frequently follow periods of increases.

Drivers of change
The random forests model classifying significant increases and decreases in lake area within the six study areas had an out-of-bag error rate of 17.16% and a kappa score of 0.60. Unscaled mean decrease in accuracy shows that the presence or absence of fire was the most important predictor across all study areas (figure 5). The next most important predictors were ground ice abundance in the surrounding terrain, and trends in mean annual temperature (figure 5). Partial dependence plots of fire history show that lakes in the increasing class were more likely to be outside of fire scars (figure 6), while lakes in the decreasing class were more likely to be within fire scars (figure 6). Partial dependence plots of ground ice abundance suggest that lakes in regions with ice-rich terrain and in regions that have experienced the most warming are also more likely to be in the increasing class.

Impact of fire on lake area
In all four study areas south of the treeline, lake area changes demonstrate a consistent trajectory following fire ( figure 7). Post-fire lake area was consistently lower for lakes within fire scars compared to nearby lakes outside of fire scars. In all study areas, differences between these groups were apparent for at least 15 years after the fire ( figure 7). These patterns reflect real changes in surface water area and not shifts in vegetation composition or soil moisture in mixed pixels along lake shorelines. When only considering pixels with water fraction >0.99, we found that lakes with a significant decreasing trend were also more likely to be located within fire scars (p < 0.01, supplemental table 2).

Discussion
Our analysis shows that the direction of change in lake area across northwestern Canada is more strongly associated with wildfire and ice-rich permafrost terrain than with the direct effects of climate (i.e. evaporation and precipitation). We observed that lakes in regions with higher ground ice abundance were more likely to increase in area, while lakes impacted by wildfire were more likely to decrease. At both broad and fine scales of analysis, we found that surface water gains and increasing trends in lake area were more common than water losses and decreasing lake area. This finding contrasts with studies in Alaska, where gradual and catastrophic lake drainage are the primary drivers of changes in lake area (Nitze et al 2020, Lara et al 2021. In the western Canadian Arctic, decreases in lake area associated with catastrophic drainage have been observed in the Old Crow Flats (Labrecque et al 2009) and the Lower Mackenzie Plain (Travers-Smith et al 2021). While drainage is an important process, our analyses suggest that localized instances of lake area losses are outweighed . Changes in lake area in the six study areas relative to the average area in the first time period (t1: [1985][1986][1987][1988][1989][1990] grouped by lake size: large lakes (>0.5 km 2 ), medium lakes (0.05-0.5 km 2 ), and small lakes (0.03-0.05 km 2 ). The lines show a smoothed fit using a generalized additive model. by widespread increases in lake area at larger spatial scales.
We observed increases in lake area associated with higher ground ice content, which may reflect ongoing changes in permafrost stability (Segal et al 2016, Kokelj et al 2017. In regions that are rich in groundice, thermokarst can lead to lateral lake expansion and the formation of new ponds (Jorgenson and Shur 2007, Kokelj et al 2017, Roy-Leveillee and Burn 2017. Permafrost degradation around existing lakes is likely being accelerated by increases in air temperature and precipitation that have been reported across the circumpolar region (Biskaborn et al 2019, Farquharson  et al 2019). Increases in lake area in northwestern Canada have been documented in the Tuktoyaktuk Coastlands (Olthof et al 2015) and Nunavut (Carroll and Loboda 2017). Among our regional study areas, the Tuktoyaktuk Coastlands showed the greatest number of lakes with significant increasing trends in area (57%) and the greatest relative change over the study period (+2.7%). Olthof et al (2015) found that increases in lake area in the Tuktoyaktuk Coastlands were driven by sub-pixel changes along lake margins, characteristic of thermokarst shoreline expansion.
Over the time-scale of analysis, our results also indicate that fire is an important determinant of decreases in lake area. In all four of the boreal  study areas (Lower Mackenzie Plains, Central Mackenzie Valley, Bulmer Plain and Great Slave Upland), lakes impacted by recent wildfire showed persistent decreases in area compared to nearby unimpacted lakes. These results are consistent with previous studies showing that lakes affected by wildfires were more likely to decline in surface area (Roach et al 2013, Frost et al 2019. Decreases in lake area following fire are also likely driven by permafrost thaw. However, compared to undisturbed areas, thaw driven by fires may be deeper and more rapid, resulting in greater active layer thicknesses and talik development, altering flowpaths and increasing hydrological connectivity (Woo and Guan 2006, Rocha and Shaver 2011, Nossov et al 2013. The removal of vegetation and surface organic cover following fire also increases ground temperatures and annual thaw depths by changing albedo and increasing ground heat flux (Nossov et al 2013). Changes in ground temperature surrounding lake basins may in turn cause taliks expansion and the loss of permafrost (O'Neill et al 2020a), which can promote surface drainage through coarse-textured soils (Nossov et al 2013). Although fire-induced ground warming could cause shoreline instability and expansion in areas of icerich permafrost, a reduction in lake levels following fire would decrease the potential for thermokarst lake expansion (Kokelj et al 2009, Roy-Leveillee and. The response of surface water to fire may also depend on the extent of underlying permafrost. Within continuous permafrost, Rocha and Shaver (2011) observed that a thinner organic layer following fire contributed to surface ponding, however in organic-rich terrain with discontinuous permafrost, Connon et al (2014) showed that thaw-driven landcover changes after fire increased surface runoff and streamflow. We were unable to fully assess the impact of fire in tundra study areas (Tuktoyaktuk Coastlands and Eagle Plain) as few large fires were recorded above the treeline. Although fire is still likely to influence surface permafrost conditions in tundra areas, the influence of fire in regions with cold, thick permafrost is likely different than regions where permafrost is warm and thin (Jorgenson et al 2010, Kokelj et al 2017. Larger and more frequent tundra fires as a result of a warming climate also demonstrates the need for additional research on lake area changes following fire across a range of permafrost conditions (Kasischke andTuretsky 2006, Hu et al 2015).
Our observation that increases in lake area were more likely in regions that have warmed more rapidly, suggests that rising temperatures are not driving widespread evaporative losses across our study area, as has been reported at some high Arctic sites (Smith et al 2005, Smol and Douglas 2007, Campbell et al 2018. It is likely however, that rising air temperatures have indirectly impacted lakes by thawing ice-rich permafrost, leading to lake expansion and ponding (Kokelj et al 2009, Roy-Leveillee and Burn 2017, Fraser et al 2018, Farquharson et al 2019. This explanation is consistent with our observation that greater increases in lake area were observed in regions with higher ground ice abundance, despite increased air temperatures across the entire study domain. In our statistical modelling, the rate of change in mean annual temperature was the third most important predictor of lake area, and trends in total annual rainfall was one of the least important variables. Between 1985 and 2020, trends in mean annual temperature and total annual precipitation across the study domain were generally increasing, but the statistical significance of these trends were variable (table 3). Snow cover, which was not included in this analysis due to lack of data at relevant spatial resolution, may also play a role in lake area change through spring snowmelt events and by regulating permafrost temperatures in the winter (Burn 2002, Arp et al 2016. Decadal trends suggest that snow cover is declining at high latitudes and more research is needed to understand interactions between topography and snowmelt at fine spatial scales (Callaghan et al 2011, Mudryk et al 2018). Previous research indicates that lake area responses to snowmelt and precipitation are controlled by fine-scale vegetation cover within lake catchments (Connon et al 2014, Turner et al 2014. We also observed that smaller lakes showed greater relative changes in area than larger lakes, suggesting that smaller water bodies may be more responsive to changes in precipitation and evaporative loss given their smaller surface area to volume ratios Bigras 1988, Arp et al 2011). Taken together, this suggests that, over our period of study, lake responses to climate have been indirect and are mediated strongly by lake size and catchment characteristics.
Our work highlights the importance of conducting hydrological studies at different spatial and temporal scales of analysis. Ecologically significant changes evident in our regional scale analyses, including declines in the area of large lakes in the Lower Mackenzie Plain and widespread increases in lake area in the Tuktoyaktuk Coastlands, were not clearly identified in our subcontinental-scale analysis. Analyses at both spatial scales show that non-linear changes in surface water over decadal periods are common. Wet periods, representing bi-directional increasing then decreasing trends in annual water fraction, were the most frequent and widely distributed water transition class across the study domain. Dry periods representing decreasing and then increasing water fraction were also more common than unidirectional water losses. Non-linear change in lake area may be linked to short periods of high-intensity rainfall and drought events, rather than longer term variation in climate , Lantz and Turner 2015, Nitze et al 2020. Anomalously high air temperatures may also decouple permafrost stability from longer term trends (Farquharson et al 2019). This highlights the sensitivity of surface water change to the temporal resolution of the study. Corroborating previous work using different sources of data and at different time scales is vital to predicting future lake area responses to climate change. Field-based studies will also be crucial for testing these hypotheses and understanding the underlying processes.

Conclusions
In this study, we found that lakes showing significant increasing trends in area were 5.6 times more frequent than lakes showing decreasing trends in area. This finding was corroborated by our broad-scale analysis, which showed that the area representing gains in water coverage was 3002 km 2 compared to 1121 km 2 of water losses. We found that the presence of wildfire, modelled ground ice abundance and the change in mean annual temperature were the most important predictors of change. Lakes located in regions with higher ground ice content were more likely to have experienced increases in area, while lakes impacted by wildfire were more likely to have experienced decreases in area. Trends in precipitation were less important in our model, suggesting that over the period of study changes in rainfall and summer temperature had indirect effects on the direction of change in lake area. Taken together this indicates that the trajectory of change in lake area and its regional drivers are complex and spatially variable. Remote sensing studies are useful in identifying broad-scale patterns and generating hypotheses, however, fieldbased studies are needed to test these hypotheses. Future studies should examine how terrain, permafrost conditions and fire severity interact to influence changes in lake area. Understanding the variable sensitivity of northern lakes will be important in predicting future impacts of lake expansion and drainage on wildlife habitat, permafrost stability and freshwater resources.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.