Recent warming reverses forty-year decline in catastrophic lake drainage and hastens gradual lake drainage across northern Alaska

Lakes represent as much as ∼25% of the total land surface area in lowland permafrost regions. Though decreasing lake area has become a widespread phenomenon in permafrost regions, our ability to forecast future patterns of lake drainage spanning gradients of space and time remain limited. Here, we modeled the drivers of gradual (steady declining lake area) and catastrophic (temporally abrupt decrease in lake area) lake drainage using 45 years of Landsat observations (i.e. 1975–2019) across 32 690 lakes spanning climate and environmental gradients across northern Alaska. We mapped lake area using supervised support vector machine classifiers and object based image analyses using five-year Landsat image composites spanning 388 968 km2. Drivers of lake drainage were determined with boosted regression tree models, using both static (e.g. lake morphology, proximity to drainage gradient) and dynamic predictor variables (e.g. temperature, precipitation, wildfire). Over the past 45 years, gradual drainage decreased lake area between 10% and 16%, but rates varied over time as the 1990s recorded the highest rates of gradual lake area losses associated with warm periods. Interestingly, the number of catastrophically drained lakes progressively decreased at a rate of ∼37% decade−1 from 1975–1979 (102–273 lakes draining year−1) to 2010–2014 (3–8 lakes draining year−1). However this 40 year negative trend was reversed during the most recent time-period (2015–2019), with observations of catastrophic drainage among the highest on record (i.e. 100–250 lakes draining year−1), the majority of which occurred in northwestern Alaska. Gradual drainage processes were driven by lake morphology, summer air and lake temperature, snow cover, active layer depth, and the thermokarst lake settlement index (R 2 adj = 0.42, CV = 0.35, p < 0.0001), whereas, catastrophic drainage was driven by the thawing season length, total precipitation, permafrost thickness, and lake temperature (R 2 adj = 0.75, CV = 0.67, p < 0.0001). Models forecast a continued decline in lake area across northern Alaska by 15%–21% by 2050. However these estimates are conservative, as the anticipated amplitude of future climate change were well-beyond historical variability and thus insufficient to forecast abrupt ‘catastrophic’ drainage processes. Results highlight the urgency to understand the potential ecological responses and feedbacks linked with ongoing Arctic landscape reorganization.


Introduction
Lakes are widespread across northern latitudes (Smith et al 2007, Grosse et al 2013. These often shallow water bodies initially formed during warm responsible for their creation, the same processes have been implicated in their destruction (Nitze et al 2020, Lara and Chipman 2021. Although lake drainage is a natural process (Billings and Peterson 1980, Mackay 1988, Hinkel et al 2003, Yoshikawa and Hinzman 2003, Marsh et al 2009, anthropogenic climate change has recently triggered widespread patterns of lake drainage across Arctic and subarctic permafrost regions (Carroll and Loboda 2017, Veremeeva et al 2021.
Though there is a general consensus that the trajectory of northern lake area is decreasing, spatial and temporal patterns of lake area change have been highly dynamic (i.e. increasing, decreasing, or stable). For example, lake area has been relatively stable on the Arctic Coastal Plain of northern Alaska, decreasing <1% since 1985 (Hinkel et al 2007, Jones et al 2009, 2020a, similar to observations on the Tuktoyaktuk Coastlands in northwestern Canada (Plug et al 2008, Marsh et al 2009. In the Old Crow Flats of the northern Yukon Territories, lake area increased 1. 6% between 19516% between -19726% between and decreased 5% between 19726% between -20016% between (Labrecque et al 2009, with drainage further intensifying through 2010 (Lantz and Turner 2015). In addition, lake area decreased across Canada between 2000 and 2009, but rates of lake area loss were greatest at high-latitudes (Carroll et al 2011). In western Siberia, lakes decreased 11% between 1978 and 1998 (Smith et al 2005), while others reported no change (Karlsson et al 2013). Lake area declined by 7% between 1972 and 2015 in the Kotzebue Sound Lowlands of northwestern Alaska (Lindgren et al 2021), in line with the 15% decline reported between 1950 and 2007 in the northern Seward Peninsula (Jones et al 2011). However, recent observations of lake drainage in northwestern Alaska have been ten times higher than historical drainage rates (Swanson 2019, Nitze et al 2020. Understanding the underlying causes, mechanisms, and drivers of such spatiotemporal patterns and variability in lake dynamics is key to anticipating the consequences of future landscape evolution.
Temporal patterns of lake change are hypothesized to regionally vary with gradients in climate, topography, surficial geology, landscape history, permafrost conditions, and associated interactions. For example, lake drainage has been linked to extreme precipitation events, which can elevate lake water levels, increasing lateral ice-wedge degradation (Mackay 1988, Marsh et al 2009. Whereas high snow fall events can create ice dams, promote drainage channel formation (Jones and Arp 2015), and subsequent outburst floods (Arp et al 2020a). Warmer winter air temperatures have also been shown to increase the likelihood of gully formation , while thinning lake ice promotes the growth of permafrost penetrating taliks (Surdu et al 2014, Arp et al 2016. Widespread talik growth has also been implicated in the recent anomalous lake drainage activity in northwestern Alaska (Swanson 2019, Nitze et al 2020, potentially triggered by a combination of warmer winter temperatures and higher than normal snowfall. In addition, local topographic variability can influence drainage pathways as thermokarst lakes gradually expand towards these drainage gradients (Jones et al 2011). Quaternary history is also a strong determinant of lake drainage, as it is relatively common phenomenon in lakes developed in ice-rich yedoma permafrost but rare in lakes developed in ice-poor permafrost terrains (Larsen et al 2017, Swanson 2019, Lara and Chipman 2021. Due to the array of complex climategeophysical interactions that influence patterns of lake drainage across permafrost and periglacial environments, our ability to forecast future lake change and associated carbon and energy dynamics has been severely limited. Here we advance knowledge of the spatial and temporal drivers of both gradual and catastrophic lake drainage in northern and northwestern Alaska. We follow lake drainage terminology as defined by Grosse et al (2013), where gradual drainage is the progressive decrease in lake area over time as a result of evaporation, deepening of thaw bulb, gully formation, or talik formation. Catastrophic drainage is defined as the rapid lake drainage within several days to weeks as a result of external (i.e. thermal or mechanical erosion of the permafrost via coastal erosion, river tapping, lake overflow) or internal factors (i.e. thermal erosion of lake banks or talik formation). We analyze the spatiotemporal patterns of 32 690 high-latitude lakes between 1975 and 2019, determine the drivers of change, and forecast lake drainage through 2050. We use supervised support vector machine (SVM) classifiers in eight seamless Landsat image composites created every five years using Google Earth Engine (GEE) to extract lake surface water extents. An object based image analysis (OBIA) masked all waterbodies and generated the lake drainage product for driver attribution analysis. Boosted regression tree (BRT) models were used to determine the primary drivers of gradual and catastrophic lake drainage, respectively. Results improve our ability to anticipate landscape evolutionary trajectories in warmer and wetter permafrost ecosystems.

Image processing
Forty-five years (i.e. 1975-2019) of Landsat observations were used to map the spatiotemporal patterns of lake drainage in northern Alaska. All Landsat data was pre-processed by the United States Geological Survey and downloaded by GEE in a radiometrically, atmospherically, and geometrically terrain-corrected state. We used Landsat surface reflectance products acquired from the Multispectral Scanner (MSS), Terrestrial Mapper (TM), Enhanced Terrestrial Mapper Plus (ETM+), and Operational Land Imager (OLI) sensors to compute eight image mosaics for the icefree period (15 June 15 to 1 September) at five-year time-periods (figure 2). We extracted the following spectral bands from all sensors: blue (0.45-0.52 µm), green (0.52-0.60 µm), red (0.63-0.69 µm), nearinfrared (0.77-0.90 µm), and short-wave infrared 1 (1.55-1.75 µm; actual wavelengths slightly vary between sensors), with the exception of the blue band in the MSS. Clouds and shadows were masked using the Fmask algorithm and short-wave infrared thresholds (e.g. Lara et al 2019). Median pixel values were selected to generate each five-year image composite, ensuring any remaining fog, haze, and high diffused clouds were filtered from the dataset (figure 2). However, during time-periods 1990-94 and 1995-99 we replaced an average of ∼5% of the study domain with imagery from the preceeding time-periods due to poor data coverage or data quality.

Decadal lake drainage
Supervised SVM algorithms were trained on 353 reference sites in GEE to extract waterbodies from each five-year image composite. To represent the full spectral range of land cover features and improve SVM feature extraction, we selected reference sites from wet tundra (n = 33), dry tundra (n = 66), sandy barrens (n = 24), geologic outcrops (n = 30), and surface waterbodies (n = 200). Surface waterbodies included spectral reflectance from shallow ponds, turbid rivers, lakes of varying origin (e.g. thermokarst, glacial, maar), and coastal water. All reference sites were specifically selected to represent waterbodies and terrain features that did not change over the 45 year observation period (figure 2).
Surface water classifications for each five-year image composite was exported to an object based image analysis (eCognition version 9.1) for lake extraction and decadal lake drainage product creation. Following protocols described in Lara et al (2018aLara et al ( , 2018bLara et al ( , 2019, we masked ponds, streams, rivers, and coastal water using spectral and morphological (i.e. roundness) metrics. Lakes were differentiated from ponds using surface area thresholds (i.e. ponds <1 ha), resulting in 32 690 lakes considered in this assessment. However, due to the different mechanisms likely driving lake expansion versus drainage (Jones et al 2011, 2020a, Shur et al 2012, Roach et al 2013, Lara et al 2019, Veremeeva et al 2021, and the small but quantifiable contribution of lake expansion to total lake area change (Smith 2005, Nitze et al 2020, we controlled for lake expansion by restricting the outward growth of lakes beyond the initial 1975-1979 lake boundaries. We computed the lake-specific change in area and percent over time, relative to these 1975-1979 lake boundaries. Therefore, the overall patterns of lake change presented here, represent a slight overestimation of the actual lake area declines observed over time. In addition, unlike typical surface water change analysis, this approach cannot discriminate between the formation of 'new' remnant ponds or lakes within the initial lake boundary following drainage, which is reported to be responsible for the increase in number of lakes in northwestern Alaska as lake area decreases (Jones et al 2011, Lindgren et al 2021. Despite differences in the spatial resolution between Landsat MSS (∼60 m) and TM, ETM+, and OLI (30 m), lake change studies have used a variety of approaches to integrate MSS data into time-series analysis. Some of these include resampling imagery to a common spatial resolution (Smith 2005, Karlsson et al 2013, restricting interpretation to larger lakes (Smith 2005, Hinkel et al 2007, Plug et al 2008, Karlsson et al 2013, or completely bypassing this issue (Rover et al 2012, Lantz and Turner 2015, Lindgren et al 2021 as the surface water mapping uncertainties between Landsat sensors is likely negligible. Nevertheless, we developed a post-classification method to inherit lake change data from the MSS 1975-1979 scenes, but maintain the 30 m spatial resolution of the TM scenes. To accomplish this task, we selected the earliest available TM lake classification (i.e. 1985-1989) and used an object based growth function to extend the 1985-1989 TM lake boundary into the overlapping 1975-1979 MSS lake boundary. If lake area increased between TM and MSS time-periods, the growth function expanded the 1985-89 TM lake area into the 1975-1979 MSS lake extent, keeping the initial 30 m spatial resolution from the TM, but inheriting the lake drainage data from the 1975-1979 time-period. If lake area change did not increase between TM and MSS time-period, the 1985-1989 lake boundary was not updated. The newly updated 30 m spatial resolution 1975-1979 lake boundaries were used to detect patterns of lake drainage across time (figures 2 and 3).

Driver attribution analysis
Although catastrophic drainage is defined as the rapid reduction in lake area within several days to weeks (Mackay 1988, Jones andArp 2015), we identify the occurrence of catastrophic drainage in Landsat observations as the abrupt decline in lake area by 30%, 45%, or 60% of the total lake area, relative to the prior observed lake area extent (e.g. difference between the 1995-1999 and 2000-2004 image composites). Three drainage thresholds were selected to define catastrophically drained lakes as the quantitative definition of gradual versus catastrophic lake drainage is not well defined in the literature, and to eliminate any potential misinterpretation of driver attribution associated with quantitative definitions. However for simplicity, only 30% and 60% thresholds are presented throughout the manuscript, and results from the 45% threshold are presented within supplemental materials 1 and 2 (available online at stacks.iop.org/ERL/ 16/124019/mmedia). We model the drivers of catastrophic drainage using data from the time-period of drainage (table 1). If a lake catastrophically drained, it was recorded and removed from the proceeding dataset to avoid double counting of drainage events. All lakes that did not catastrophically drain were categorized as gradually draining lakes, where we modeled the drivers of long-term lake-specific change rates (i.e. including stable and gradually draining lakes).
Driver attribution analysis was implemented using a BRT model for each drainage type (i.e. gradual and catastrophic) and threshold (30%, 45%, and 60%). We used 35 static (i.e. temporally non-changing parameter) and dynamic (i.e. temporally changing parameter) morphologic, topographic, environmental, disturbance, and climate related parameters to model lake drainage (table 1). With the exception of lake-specific morphological metrics, all spatial data were extracted using 500 m lake buffers. Morphological metrics lake perimeter, area, edge-to-area ratio, roundness, and distance to drainage gradients (adjacent lake or river), were computed within the OBIA (e.g. Lara et al 2018a), while the thermokarst lake settlement index (TSI) was calculated using the Arctic digital elevation model (DEM), following Lara and Chipman (2021). Topographical metrics, topographic position index, topographic wetness index, elevation, slope, and aspect were also computed using the Arctic DEM (Danielson and Gesch 2011). Environmental metrics, day of thaw (DOT), day of freeze (DOF), length of thawing season (number of days between DOT and DOF), and snowfall equivalent (SWE) were derived from the Scenarios Network of Alaska and Arctic Planning (SNAP) (CRU 4.0; Harris et al 2014), while permafrost temperature (MAGT), ground surface temperature (MAST), active layer thickness (ALT), and minimum permafrost thickness (Pthick) were modeled by Geophysical Institute Permafrost Laboratory model (Luo et al 2014). Near surface permafrost probability (PermProb) was obtained from Pastick et al (2015) and Vegetation cover (Veg) referred to the latest National Land Cover Dataset (NLCD) Land Cover Map (Homer et al 2015). Disturbance metrics year since fire (YSF), Fire Severity (Severity), and Fire BRT algorithms were trained and validated using all lake drainage data (e.g. Elith et al 2008). We used all candidate predictors to fit the model (table 1), and removed all non-significant parameters (relative influence <5%) in a step-wise fashion to avoid overfitting following the principle of maximum parsimony and minimizing cross-dependent and crosscorrelative variables. The regularization parameters of the model were tested before being fixed at the optimal bag fraction of 0.75, learning rate of 0.0001, and tree complexity of 8, which resulted in the number of trees >5000. Model performance was evaluated using a k-fold cross-validation procedure (ten-folds), where lake drainage data was randomly divided into ten subsets (training:testing) to estimate the average accuracy of the model after each fold. All BRT model development was performed in R v3.6.1 (with package 'gbm' , e.g. Greenwell et al 2020).

Results
Forty-five years of Landsat observations were used to generate a seamless 30 m spatial resolution lake drainage map product of 32 690 lakes spanning six ecoregions in northern Alaska (figure 3). This product characterized the timing of lake area change between 1975 and 2019, and was used to differentiate gradual from catastrophic drainage patterns and processes.   '1975-79-1980-84' and '1980-84-1985-89' were recomputed by dividing the lake area change between 1975-79 and 1985-89 by two (panel (B)). The 'All Ecoregions' category represents the mean lake change by ecoregion and time-period. and the Brooks Range (−10.0%), respectively. The overall mean lake drainage rate across all ecoregions was approximately −3.6% decade −1 , but rates varied by ecoregion and over time (figure 3). Drainage rates from northern to southern ecoregions ranged from the Brooks Foothills (−3.6% decade −1 ), Brooks Range (−2.6% decade −1 ), Davidson Mountains (−3.4% decade −1 ), Kobuk Ridges and Valleys (−4.1% decade −1 ), Kotzebue Sound Lowlands (−3.8% decade −1 ), and the Seward Peninsula (−3.2% decade −1 ). The temporal variability in lake drainage rates was apparent as only about 7% of the cumulative lake drainage over time occurred between the first 15 years of observation (i.e. 1975-1989), whereas nearly half of all drainage occurred within the 1990s (figure 4).

Spatiotemporal patterns of lake drainage
Gradual lake drainage dominated the overall proportion of lake area change across time and space (figure 5), as it was responsible for between 94.0% and 97.6% of all lake area loss (range dependent on 30%-60% drainage threshold definitions). In contrast, catastrophic drainage processes were responsible for only between 2.4% and 6.0%. Interestingly, the proportion of lakes to experience gradual versus catastrophic drainage were not static over time, as the proportion of lakes to experience gradual drainage increased between 1.8% and 1.1% decade −1 as catastrophic drainage proportionately decreased (figure 4). For example, during the earliest timeperiod (1975)(1976)(1977)(1978)(1979), catastrophic drainage represented as much as 10.0%-4.5% (range associated with 30% and 60% threshold definition) of all lake area loss, but linearly decreased over time, representing only 3.0%-0.5% of all drained lake area during the most recent time-period (2015-2019).
In line with the decreasing proportion of catastrophic drainage, the number of catastrophically

Spatiotemporal drivers of lake drainage
Despite differences in our quantitative definitions of gradual and catastrophic drainage (i.e. 30%, 45%, and 60%), BRT models consistently identified the drivers of change to differ by drainage process, suggesting drainage definitions had little influence on our overall results. Gradual lake drainage models for 30 (n = 29 020 lakes), 45 (n =31 244 lakes), and 60% (n = 32 329 lakes) drainage thresholds explained 42%, 39%, and 39% of the overall variance in lake drainage in northern Alaska, respectively. All gradual drainage models identified the same six drivers: the edge-to-area ratio (E:A), mean summer air temperature (MSAT), lake skin temperature (LakeT), snow cover days (SD), ALT, and the TSI. The ten-fold cross validation procedure indicated model validation error was low (10-CV = 0.35-0.30, P < 0.001; supplemental material 1) as locally indicated by observed versus predicted gradual drainage between 1975 and 2019 (supplemental material 3).
Due to the skill of BRT lake drainage models, we forecast the potential lake drainage across northern Alaska over the next 30 years (figure 8). Models project gradual lake drainage to continue to decrease in lake area between 19% and 23% across our study domain by 2050. However, this pattern differed by ecoregion, as cumulative gradual lake  1, 3). Despite good overall model performance of catastrophic drainage (i.e. low bias and low variance), the model was unable to account for the climate variability projected by 2050. BRT models markedly over-projected the frequency of catastrophic drainage events by projecting between 68.4% and 90.0% of all lakes to catastrophically drain in the study region by 2050. Due to limited historical observations beyond 157 ThawDays, the model was highly sensitive to values beyond this range, as all northwestern ecoregions approached the 157 ThawDays threshold by 2050 (supplemental material 4), whereas the Kobuk Ridges and Valleys exceeded this threshold as early as 2020. Therefore, we did not further scrutinize these lake drainage predictions.

Spatiotemporal patterns of lake drainage
Using nearly a half-century of Landsat observations (1975-2019), we identified a 14.4% (∼3.6% decade −1 ) decrease in lake area, spanning the upland and lowland tundra region of northern and northwestern Alaska. Catastrophically drained lakes were a relatively uncommon form of lake drainage, accounting for 2.4%-6.0% of annual lake area losses. The disproportionately higher prevalence of gradual drainage was likely due to the broader influence regional climate change had on all lakes (via evaporation versus recharge). Increased gradual drainage rates across northern Alaska during the 1990s (figure 4(B)) were nearly identical to that observed on the Old Crow Flats in the northwestern Yukon, which attributed amplified drainage to regional water deficits between 1988 and 2001, associated with the Pacific Decadal and Arctic Oscillations (Labrecque et al 2009). Additionally, we observe the number of catastrophically drained lakes to decrease between 1975-1979and 2010, similar to Jones et al (2020 and Marsh et al (2009). However, this forty-year negative trend was reversed during the most recent time-period (i.e. 2015-2019) as the number of catastrophic drainage events was among highest on record (figure 6). These lake drainage observations may be even more notable as we also identify a negative trend in the proportion of gradual versus catastrophic drainage throughout all time-periods (figure 5), suggesting the most recent time period (2015-2019) was an extremely high lake drainage epoch for both gradual and catastrophic drainage. We hypothesize this was due to the progressive thawing and elimination of lakes with thin stabilizing permafrost as the climate warmed and winter snow cover increased (Nitze et al 2020). At the start of the anomalous warming periods during the 1980s and 1990s, lakes with relatively thin stabilizing ground ice are more likely to be punctured and lost to the groundwater by the progressive thaw subsidence and lake deepening. Therefore, over time the most sensitive lakes will have been lost quickly after the climate began to change, resulting in the observed decreasing trends in catastrophic drainage rates.
Interestingly, the recent high number of lake drainage events (i.e. 2015-2019) were all recorded in northwestern ecoregions in the Kotzebue Sound Lowlands, Kobuk Ridges and Valleys, and the Seward Peninsula ( figure 6). This observation confirms recent patterns of lake area losses reported within the Kotzebue Sound Lowlands (Nitze et al 2020), the Noatak National Preserve (Swanson 2019), and the Seward Peninsula (Lindgren et al 2021). In addition to the identified proximate drivers of lake drainage, the high drainage rates in this region may be ultimately influenced by the northward shifting continuous to discontinuous permafrost boundary , Lindgren et al 2021 with regional warming (Jafarov et al 2012).

Modeling drivers of lake drainage
Due to evidence suggesting that the drivers and mechanisms of gradual and catastrophic lake drainage may differ (Jones et al 2011, 2020a, Shur et al 2012, Roach et al 2013, Lara et al 2019, Veremeeva et al 2021, we evaluated the drivers of both drainage processes separately. Gradual lake drainage was controlled by lake morphology and type (i.e. E:A, TSI), temperature (i.e. MSAT, LakeT), and snow-permafrost interactions (i.e. SD and ALT). We found the more irregular the shoreline and the warmer the air/water temperatures (i.e. increased evaporation), the more vulnerable the lake to gradual drainage. In addition, the identification of the thermokarst lake settlement index confirms the predisposition that thermokarst lakes are indeed more likely to experience lake drainage than non-thermokarst lakes of glacial or volcanic origin (Swanson 2019, Lara andChipman 2021). Our identified drivers of catastrophic drainage wellrepresented the spatiotemporal patterns of drainage events (supplemental material 1-3), and correspond with previous observations (Labrecque et al 2009, Jones et al 2011, Lantz and Turner 2015, Swanson 2019, Nitze et al 2020, Lindgren et al 2021. However we identify historical patterns of catastrophic lake drainage to be regionally controlled by the length of thawing season (i.e. ThawDays), as partial dependency plots identified the magnitude of catastrophic drainage events to be closely linked to a threshold of 157 thawing days, after which the likelihood of catastrophic drainage events rapidly increased. This pattern explains the recent high catastrophic drainage events in northwestern Alaska (2014-2019), as the ThawDays threshold was often locally exceeded but regionally approached in the Kobuk Ridges and Valleys (156.6 days), Kotzebue Sound Lowlands (150.4 days), and the Seward Peninsula (152.0 days). We interpret the additional drivers of catastrophic drainage, precipitation and permafrost thickness to be indirect indicators of regional landscape heterogeneity as (a) the dominate precipitation gradient across our study domain occurs from low to high elevations associated differences in lapse rates, and (b) the modeled permafrost thickness generally decreased with latitude and increased with elevation.
Equally important are the parameters that the models did not identify as regional drivers of lake drainage, such as fire metrics, proximity to drainage gradients, and increased precipitation. Although wildfire may indirectly alter the timing and extent of ALT by modifying surface energy dynamics , linkages to local pond and lake drainage events are limited (e.g. Roach et al 2013, Frost et al 2019, Chen et al 2020, and also not supported by our results. In addition, the proximity to drainage gradient (i.e. distance to river, stream, or lake) has been reported to be a predominant drainage mechanism in low-relief coastal tundra ecosystems as thermokarst lakes laterally expand (Jones et al 2011. However, models failed to identify this well-studied mechanism of drainage reported for the northern Seward Peninsula and the Arctic Coastal Plain of Alaska, suggesting the identified drivers (i.e. warming and extended thawing season) may be robust regional indicators of lake drainage, but represent an array of local mechanisms of drainage. Increased precipitation has also been linked with lake drainage as rising lake water may facilitate bank overtopping and gully or channel formation Arp 2015, Jones et al 2020a), a phenomenon that may be exacerbated by adjacent ice-wedge degradation that create new surface hydrological flow paths for drainage (Jorgenson et al 2006, Liljedahl et al 2016. However we did not identify precipitation change as a dominant driver of lake drainage (gradual or catastrophic), which was in line with Swanson (2019) as he reported high rates of lake drainage following record high temperatures during periods of declining precipitation across northwestern Alaska.

Forecasting lake drainage
By 2050, models project gradual lake drainage to continue to decrease lake area between 19% and 23% across our study domain (figure 8). Hotspots of lake drainage are found in the intermontane boreal ecoregions (i.e. Kobuk Ridges and Valleys and the Davidson mountains; ∼34% lake area losses), as they are projected to experience among the greatest warming in MSAT and associated vertical permafrost thaw (i.e. ALT; supplemental material 4). However, it remains unclear how these climate-driven patterns in future lake loss may be further compounded by catastrophic drainage. Assuming the decreasing trend in the proportion of catastrophic versus gradual drainage will continue (figure 5), the contribution of catastrophic drainage to overall lake area losses approach zero by 2020-24 (30% drainage threshold) or 2030-34 (60% drainage threshold) across upland and lowland tundra of Alaska, suggesting gradual drainage processes will govern nearly all lake drainage by 2050. In contrast, our catastrophic drainage model artificially triggered an enormous (i.e. 65% of all lakes) wave of catastrophically drained lakes across the region. Despite our spatially and temporally robust lake drainage dataset spanning climate, environmental, topographic, and disturbance gradients, which generated models that well-represented past catastrophic drainage dynamics (supplemental material 2, 3), our observations remained insufficient to predict how abrupt thermokarst-driven lake dynamics will change over the next 30 years. This highlights the challenges of predicting abrupt permafrost degradation processes with limited knowledge of ground ice distribution and without explicit linkages with local mechanisms of drainage. To improve forecasts of thermokarstdriven lake drainage processes, future research should aim to characterize and model the various mechanisms and pathways of talik formation, sub-surface tunnel development, gully or drainage channel formation, bank-overtopping, lake expansion, and associated interactions with adjacent land cover properties (Yoshikawa and Hinzman 2003, Shur et al 2012, Jones and Arp 2015, Arp et al 2016, Liljedahl et al 2016, as they are likely to accelerate as Arctic and sub-Arctic regions become warmer and wetter (Liu et al 2012, Bintanja and Selten 2014, Arp et al 2020b, Jones et al 2022.

Conclusion
We leveraged the entire Landsat image archive (MSS, TM, ETM+, OLI) to create wall-to-wall lake drainage maps for six northern and northwestern ecoregions of Alaska. The frequency of lake drainage observations enabled the identification of the likely drivers of change, which differed between drainage processes. The drivers of gradual and catastrophic drainage corresponded with climate-driven versus thermokarstdriven parameters, respectively. We found gradual lake drainage to intensify over the 45 year Landsat time-series as increased temperatures (i.e. MSAT and LakeT) and potentially surface water deficient (e.g. Labrecque et al 2009) reduced overall lake area, which is likely to continue through 2050. By leveraging early Landsat observations we captured intense periods of catastrophic lake drainage events that (a) linearly decreased for forty-years, but was reversed by one the highest number of catastrophic drainage periods on record (i.e. 2015-2019), and (b) corresponded with periods of amplified summer warming and lengthening of the thawing season. Although the inclusion of a larger breadth of lake drainage observations and climate and environmental variability undoubtedly improved model performance, this historical variability was still insufficient for sophisticated machine learning models to forecast catastrophic drainage. Collectively, these results bridge many previous observations of high-latitude lake drainage and provide new insights into (a) how permafrost ecosystems are changing, and (b) why they are changing, however further research into the biogeophysical consequences of such transitions from waterto-land (e.g. Walter Anthony et al 2018), should be prioritized.

Data availability statement
The data that support the findings of this study are openly available at the following DOI: 10.18739/ A2BV79W8S.