Remote Sensing of Multitemporal Functional Lake‐To‐Channel Connectivity and Implications for Water Movement Through the Mackenzie River Delta, Canada

The Mackenzie River Delta in Canada is a mediator of hydrological transport between the expansive Mackenzie River watershed and the Beaufort Sea. Within the delta, lakes frequently act as water and sediment traps, limiting or delaying the movement of material to the coastal ocean. The degree to which this filtering takes place depends on the ease with which sediment‐laden water is transported from distributary channels into deltaic lakes, referred to as functional lake‐to‐channel connectivity, which varies both spatially and temporally. Tracking of connectivity has previously been limited to either small regions of the delta or has focused on a snapshot of connectivity at a single instance in time. Here we describe an algorithm that uses Landsat imagery to track summertime functional lake‐to‐channel connectivity of 10,362 lakes between 1984 and 2022 on an image‐by‐image basis. We calculate a total average connected lake area of 1400.7 km2 during the 2 weeks after peak discharge, 763.6 km2 higher than previous estimates, suggesting a larger influence of connected lakes on water movement through the delta than previously estimated. We also identify water level thresholds that lead to the initiation of high sediment river water movement into 5,989 lakes (908 lakes with uncertainty ≤±0.5 m), and identify an additional 2899 lakes whose connectivity does not vary at all. As the Arctic hydrological cycle responds to climate change, this work lays a foundation for tracking the movement of water, and the matter it carries, from the Mackenzie River watershed to the Beaufort Sea.


Introduction
Arctic delta lakes and channels are mediators of sediment and water transport between terrestrial landscapes and the Arctic Ocean.As river water flows into lakes, velocity slows, material falls out of suspension, and the lakes behaving functionally as sediment and water traps (Emmerton et al., 2008;Marsh et al., 1999;Piliouras & Rowland, 2020).Additionally, lake-to-channel connectivity modifies lacustrine biogeochemical properties and ecological function (Cunada et al., 2018;Gareis & Lesack, 2018;Lattaud et al., 2021;Michelutti et al., 2001;Scott et al., 2020;Spears & Lesack, 2006;Squires & Lesack, 2002;Tank et al., 2011;Vonk et al., 2015).Lake connectedness can be defined as "structural" or "functional."Structural connectivity refers to physical connectedness of the lake, often based on connected channel presence/absence or the water elevation needed to induce flow of water into a lake (Lesack & Marsh, 2010;Marsh & Hey, 1988, 1989;Passalacqua, 2017;Piliouras & Rowland, 2020).Functional connectivity specifically refers to the movement of water and sediment from a river channel into a lake (Dolan et al., 2021;Long & Pavelsky, 2013).For example, a lake with a consistent but indirect structural connection to the channel network (e.g., a lake that is structurally connected to the distributary channel network via a chain of other structurally connected lakes) could be high functional connectivity when freshet-related elevated discharge carries high sediment water into the lake but low functional connectivity during other periods.
Extensive research on structural connectivity and its implications has focused on the Mackenzie Delta region in the Northwest Territories of Canada (Lattaud et al., 2021;Lesack & Marsh, 2010;Marsh & Hey, 1988, 1989;Michelutti et al., 2001;Scott et al., 2020;Tank et al., 2011;Vonk et al., 2015).This body of literature defines connectivity based on lake sill elevations, the lowest point along the thalweg of the connection channel above which water needs to rise to flow into the lake.Sill elevation has historically been measured using high spatiotemporal resolution aerial surveying to detect the temporal onset of lake flooding paired with temporally coincident water level data from nearby Water Survey of Canada (WSC) gages (Lesack & Marsh, 2010;Marsh & Hey, 1988, 1989).Using this method, researchers have been able to track connection durations of 141 lakes near Inuvik over the past 40 years, defining lakes as high-closure (>4 m asl), low-closure (2.5-4 m asl), and no-closure (sill elevation <2.5 m asl).
It is through this framework that we understand the role of connectivity in lake biogeochemistry in the Mackenzie Delta.High closure lakes with short and variable connection times are extremely individualized, in part because they are impacted by legacy effects of river water inflow sporadically through time (Lesack & Marsh, 2010;Scott et al., 2020).Typically, these lakes experience elevated macrophyte (plant) productivity as well as long water residence times (Lesack & Marsh, 2010;Michelutti et al., 2001;Tank et al., 2011), leading to high sedimentary organic carbon levels made up of young plant-derived carbon (Vonk et al., 2015).These high-closure lakes bury very little petrogenic (rock-derived) hydrocarbon and sequester carbon less efficiently than their low-and noclosure counterparts (Lattaud et al., 2021).Shorter connection times associated with high-closure lakes additionally lead to higher rates of methane production (Cunada et al., 2018;Michelutti et al., 2001).While Dissolved Organic Carbon (DOC) tends to be higher in high-closure poorly connected lakes (Wolfe et al., 2007), the source of DOC, which varies across the connectivity spectrum, impacts how that DOC is used (Tank et al., 2011).For example, macrophytic (plant) DOC found predominately in high-closure lakes is quickly processed by bacteria, so it is less likely to play a major role in physiochemical processes that impact water chemistry and light penetration.Contrastingly, allochthonous DOC from river water more commonly associated with low-to noclosure lakes is more readily useable for physiochemical processes (Tank et al., 2011).
Connectivity in the Mackenzie Delta is linked to phytoplankton and bacterioplankton productivity.Phytoplankton productivity can be limited either by nutrient availability or light availability (e.g., Squires & Lesack, 2002).Water in highly connected lakes is elevated in suspended sediments; therefore, light was previously hypothesized to be the factor limiting productivity (Squires & Lesack, 2002).However, while suspended sediment is the dominant attenuator of light in Mackenzie Delta lakes during the spring freshet, chromophoric dissolved organic matter becomes the dominant control on light attenuation during the rest of the open water season (Squires & Lesack, 2002).Additionally, light limitation does not appear to be a major limiting factor in phytoplankton productivity in highly connected lakes, likely due to regular water mixing during the summer and ample nutrients in riverine water inputs (Emmerton et al., 2008;Squires & Lesack, 2002).However, nutrient limitation becomes a factor as connectivity decreases.Conversely, bacterioplankton productivity increases with decreasing flood frequency (Spears & Lesack, 2006).
While extensive prior research indicates the importance of connectivity in the Mackenzie Delta, there is limited knowledge regarding spatiotemporal patterns of connectivity throughout the delta as a whole.The sill elevation method generally requires aerial or field-based observations and has only been applied within small areas of the delta near the town of Inuvik. Piliouras and Rowland (2020) developed a channel detection method using 30 m Landsat imagery and found a total connected lake area of 637.1 km 2 .While this method is widely applicable across deltas without training or in situ data, the technique is limited to lakes whose connection channels are wide enough to be visible in Landsat imagery and only has been applied to the late summertime period in a single year.
In this analysis, we use Landsat imagery to investigate patterns, trends, and drivers of lake-to-channel functional connectivity in the Mackenzie River Delta between 1984 and 2022.Long and Pavelsky (2013) and Dolan et al. (2021) demonstrated the efficacy of tracking high sediment river water recharge to detect functional connectivity, and it is on this idea that we build our methods.Our method does not require channel detection, nor does it require frequent aerial monitoring of lakes, and it is spatially applicable across all 10,362 Landsat-observable Mackenzie Delta lakes.We additionally quantify the ability of Mackenzie Delta lakes to attenuate flood water compared to other large Arctic deltas (following Piliouras & Rowland, 2020) and calculate water level thresholds that must be reached in order to initiate transport of high sediment river water into each lake.

Study Area
In The Mackenzie River Delta is in the Northwest Territories of Canada at the coastal intersection between the Mackenzie River and the Beaufort Sea (Figure 1).The expansive (∼13,000 km 2 ) delta is fed primarily by the Mackenzie River (90% of flow), with additional inflow from the Peel River and other smaller creeks draining the Richardson Mountains to the west (Emmerton et al., 2007).The delta is lake-rich, with more than 45,000 mostly small and shallow water bodies larger than 0.0014 km 2 (Emmerton et al., 2007).An average discharge volume of 316 km 3 /yr ± 10% (Stubbins et al., 2015) moves from the delta into the Beaufort Sea each year.On the way to the coast, the water and its suspended materials move through a complex network of lakes and channels, leading to net deposition of materials within the delta.An estimated 128 Mt of sediment enter the delta each year and 85 Mt of sediment reach the Beaufort Sea, with an additional 43 Mt internally deposited on levees and lake beds (Carson et al., 1999).However, the actual flux of sediment is likely much higher due to bank erosion and deposition taking place within the delta, particularly within the Middle Channel (Carson et al., 1999).Approximately 60% of suspended sediment enters the delta in late May through early June due to river ice breakup and the spring freshet (Carson et al., 1998).Sedimentation rates in lakes vary based on lake connectivity and range from 1 to 10 mm annually (Marsh et al., 1999).While the delta has not experienced statistically significant trends in total freshet volume or annual discharge volume, the freshet has trended earlier in concert with local spring warming during the ice melt period and earlier snow melt (Lesack et al., 2014;Rood et al., 2017;H. Yang et al., 2015).
The Mackenzie Delta is the traditional lands and the current settlement area of the Inuvialuit and Gwichya Gwich'in First Nations.Traditional Gwich'in land use areas are south of the Mackenzie River, as far south and east as the Thunder River, throughout the Arctic Red River, Cranswik River, and Snake River watersheds, as well as seasonally within the delta ("Eedyee" in Gwich'in) itself (Kritsch & Andre, 1994).Historically, Eedyee was a place for spring muskrat trapping and summer fishing; permanent settlements were relatively uncommon until recently.Similarly, Inuvialuit settlements were historically transitory.Inuvialuit traditional land use areas are near the Beaufort Sea Coast, centered around hunting for Bowhead and Beluga whales (Friesen, 2015).However, Inuvialuit settlements began to move inland as early as the fourteenth and fifteenth centuries CE.The two largest present-day permanent settlements in the delta are Aklavik, founded in 1912 as a Hudson Bay trading post (Gwich'in Social & Cultural Institute, 2016), and Inuvik founded in the late 1950s as a less flood-prone alternative to Aklavik ("Inuvik," 2020), both located on the boundaries between Inuvialuit and Gwich'in settlement areas (Krause, 2022).

Methods
In this research, we used optical Landsat imagery in conjunction with a Random Forest model to classify functional connectivity of Mackenzie Delta lakes during the open water season (late May through September) between 1984 and 2022 on an image-by-image basis.

Lake Selection
We focused on lakes in the delta larger than 8100 m 2 (3 × 3 Landsat pixels) as defined by the UCLA Lake Database (J.Wang et al., 2023) that also contained at least 100 pixels that were classified as water in at least 80%  et al., 2021), study lakes in white (UCLA Lake Database; J. Wang et al., 2023), Water Survey of Canada gauges used for water level (black triangle) and discharge (black square) analyses, and (inset) the location of the delta relative to North America. of Landsat observations since 1984 according to the 30 m JRC Global Surface Water Mapping Layers v1.3 (Pekel et al., 2016).This analysis focuses on the 10,362 lakes that meet these criteria.

Landsat Pre-Processing
We used Landsat 5 (Thematic Mapper), 7 (Enhanced Thematic Mapper Plus), and 8 (Operational Land Imager) Collection 2 Tier 1 Surface Reflectance imagery during the open water season (late May through September) that were less than 75% cloudy (n = 3,581 images) to investigate long-term spatiotemporal patterns and trends in Mackenzie Delta connectivity.Detailed discussion of Landsat pre-processing steps including atmospheric correction, cloud filtering, harmonization of imagery between Landsat missions, and both the identification of water pixels and removal of snow/ice pixels can be found in Text S1 in Supporting Information S1.
For each cloud-free and ice-free preprocessed Landsat image, we calculated the mean reflectance of each calibrated band within each lake polygon as well as a number of other optical parameters (Table S1 in Supporting Information S1).If fewer than 50% of the lake's open water pixels were cloud-free in a given Landsat observation, we removed that observation from analysis.While Dolan et al. (2021) used ratios of lake reflectance to channel reflectance to determine connectivity, we found that this method overestimated the movement of high sediment water into lakes in the Mackenzie Delta, specifically during periods of low sediment transport.Therefore, this analysis instead uses the optical properties of each lake as input to a random forest model.

Model Development
This analysis relies on Random Forest classification, a supervised classification method that requires training data.A suitable training data set was not previously available and therefore needed to be developed.We randomly chose ∼20% of Landsat-observable lakes (n = 1815) and manually classified their functional connectivity during two distinct periods representing a high sediment period in the delta (June 10th-30th 2020) and a low sediment period in the delta (August 22nd-10th September 2020), as defined by visual observation of imagery.Given that high connectivity is commonly only spatially pervasive in the delta during this early summer period as part of the spring freshet, a random sampling scheme throughout the summer would have resulted in an undersampling of high connectivity observations.We first built a Google Earth Engine User Interface displaying side-by-side highresolution Sentinel-2 (10 m) and Planet imagery (5 m) from the high sediment period and the low sediment period.Through this user interface, we visualized one lake at a time and classified each lake (Figure S1 in Supporting Information S1) based both on structural (channel presence/absence) and functional connectivity during each period, defined by expert assessment of the visible movement of high sediment river water into a lake, referred to as high sediment river water recharge.These 14 classes were then combined into three final connectivity classes specifically focused on functional connectivity: 0-low functional connectivity, 1-some functional connectivity, and 2-high functional connectivity (Figure S2 in Supporting Information S1).Of these 3,630 observations of 1,815 lakes, 982 observations were cloudy, ice-covered, visually influenced by coastal processes (lakes that visually appeared to receive more water from the coastal ocean than river channels), or the connectivity was uncertain.An additional 887 lacked coincident Landsat observations, leaving a remaining 1,759 observations of 1,257 distinct lakes for training and testing purposes (n 0 = 778, n 1 = 494, n 2 = 487).
While we obtained Landsat observations of each lake on an image-by-image basis, the training data were collected using imagery averaged and mosaicked during two distinct time periods (June 10th-30th 2020 and August 22nd-10th September 2020).Therefore, we averaged together Landsat observations during each of the two time periods and used those optical values in training and testing our model.We used the Boruta method (Kursa & Rudnicki, 2010) to identify which optical properties in Table S1 in Supporting Information S1 were important for distinguishing between connectivity classes.As all properties were found to be important, we kept the most important property, Band Combination 1 (NIR-SWIR1), and removed other properties that were either highly correlated (ρ > 0.8) to NIR-SWIR1 or to each other.The remaining optical properties were NIR-SWIR1, NSMI, Blue/NIR, and Green/Blue.We split the resulting observations into 75% training (n = 1,318, n 0 = 583, n 1 = 370, n 2 = 365) and 25% testing (n = 441, n 0 = 195, n 1 = 124, n 2 = 122), stratified by connectivity classification, and tuned a 500-tree random forest model on the training data set via 5-fold cross-validation.We selected the model with the highest Receiver Operating Characteristics Area Under the Curve metric, which had the following parameters: mtry = 1 (number of randomly selected variables sampled at each split) and min_n = 28 (number of data points in a node required for ).We assessed model performance by applying the final model to the testing data set and analyzed both the resulting confusion matrix and Gwet's AC 2 metric weighted linearly (Gwet, 2015(Gwet, , 2021).Gwet's AC 2 accounts both for classification agreement due to random chance and for the fact that certain misclassifications are better or worse than others.For example, misclassifying a class 0 observation (low functional connectivity) as a class 1 observation (some functional connectivity) is less problematic than misclassifying a class 0 as a class 2 (high functional connectivity).We then applied the classification algorithm to all open water Landsat observations between 1984 and 2022 (n = 2,831,027 total observations of 10,362 lakes).

Data Availability
Cloud-free and ice-free data availability peaked in June, followed, in decreasing order of availability, by July, August, September, and May (Figure 2).May imagery was limited by ice cover and September imagery was limited by both ice cover and clouds.For the remainder of the analysis, the month of May was excluded due to a lack of ice-free imagery.

Trend Analysis
We calculated trends in connectivity during the month after high freshet discharge begins at the WSC station at Arctic Red River (Station Number: 10LC014;data available 1984-2019).We defined annual freshet initiation as the first day when discharge increases exceeded 3% (Lesack et al., 2013) and defined the start of the high discharge freshet as the first date after that initiation in which discharge stopped increasing and was also higher than 10,000 m 3 /s.We then selected all connectivity observations within 4 weeks of this date and calculated the average connectivity class in this period for each lake in each year.We then used a non-parametric Mann-Kendall test (Yue et al., 2002) to analyze trends in connectivity within three periods: 1984-2001, 2002-2019, and 1984-2019.The trend analysis only included lakes that had at least 10 years of observation in both shorter temporal periods (n = 6,457 lakes).

Connected Lake Area and Storage
To calculate the mean area of connected lakes across study years, we selected all lakes that experienced a mean connectivity class of greater than or equal to one during the 2 weeks after the initiation of the annual spring freshet and summed their areas.Following Piliouras and Rowland (2020), we multiplied an assumed average flood water depth of 1 m, a common freshet water level increase in the Mackenzie Delta (Lesack & Marsh, 2007;Marsh & Lesack, 1996), by the connected lake area to capture an approximate estimate of flood storage in lakes that are commonly connected during the spring freshet.We then compared mean annual flood volume (Piliouras & Rowland, 2020) to the estimated connected lake storage volume to better understand the capabilities of connected lakes to attenuate freshet water and sediment transport.We also calculated the fraction of total lake area greater than varying connectivity class thresholds, ranging from zero to two, to assess the sensitivity of that threshold, both during the 2 weeks after the start of the high discharge freshet, and within 2 week increments throughout the entire summer.

Functional Connectivity Elevation Threshold Calculations
We combined connectivity classifications with WSC water level data at six stations throughout the delta (Table S2 in Supporting Information S1; accessed on 3/23/2023) to estimate the water level needed for high sediment river water to be moved into a lake, which we refer to as the functional connectivity elevation threshold.First, we corrected water levels at each station to CGG05, which is thought to be the most accurate vertical datum for regions close to the Beaufort Sea (Véronneau et al., 2006;Table S3 in Supporting Information S1).As water level is measured at 6 discrete locations, to better understand water level variability in the delta, we first calculated both the distance between all station pairs as well as the difference in average water level between all station pairs.Next, we used a linear model to define the relationship between both distance and water level parameters (Figure S5 in Supporting Information S1).Under approximately 50 km, the difference in average water level between stations was less than 0.5 m.
To calculate the functional connectivity elevation threshold for a lake, we paired same-day water level with sameday connectivity classifications from any WSC station within 50 km of the lake.We then selected connectivity observations of end-member classes (0 and 2) and assessed the water level necessary for the lake to transition from class 0 to class 2. We first completed a t-test (Kalpić et al., 2011) to investigate whether or not there were significant (p < 0.05) differences in water level when a lake was at class 0 compared to when it was at class 2. If there was a significant difference, then we proceeded to complete one of two calculations depending on the data.For lakes with any overlap in water level between class 0 and class 2 observations, we defined a possible range of functional connectivity elevation thresholds bounded by the minimum water level that we observed at class 2 on the lower end and the maximum water level that we observe as class 0 on the higher end (Figure 3b).If there were no overlaps in water level, then the range of possible functional connectivity elevation thresholds was bounded on the lower end by the maximum water level observed at class 0 and by the minimum water level observed at class 2 on the upper end (Figure 3a).We added additional uncertainty onto these ranges based on how far the lake was from the WSC station, using the distance from the station multiplied by the slope of the linear relationship between water level difference and distances described above and shown in Figure S5 in Supporting Information S1 (slope = 1.15 × 10 5 m/m).We did not calculate functional connectivity elevation thresholds for lakes with connectivity class of 0 (i.e., lakes nearly always disconnected) or of 2 (lakes nearly always connected) in >95% of all observations.Given that some lakes were within 50 km of two or three gauge stations, we calculated functional connectivity elevation threshold ranges based on all relevant stations independently and assessed overlap between the ranges.If ranges did overlap, we then narrowed the range of possible elevation thresholds to that region of overlap.If the elevation threshold ranges did not overlap, we did not calculate a final functional connectivity elevation threshold range for that lake.If, for example, a lake's connectivity was not related to water levels at one station but was at another station, we discarded the former measurement and calculated an elevation threshold using the latter.The final connectivity elevation threshold for each lake was defined as the midpoint of the possible range, and the uncertainty of that measurement is described as the total range divided by 2.
While field-based observations of functional connectivity elevation threshold do not exist for validation purposes, structural sill elevations were measured in a small area of the delta in the mid-1980s using methods similar to those used in this paper (Marsh & Hey, 1988).During fieldwork described by Marsh and Hey (1988), researchers paired daily WSC water level data at Inuvik with observations of floodwater intrusion onset in 141 nearby lakes.Similarly, flood onset was determined by visually identifying the intrusion of high sediment floodwater into each lake.These methods specifically focused on the water levels that were required to overtop structural sills and described the state of the delta during a single period in the 1980s, rather than an average delta state over more than 35 years, as is done in this paper.However, these data provide a useful comparison tool with which to evaluate our functional connectivity elevation thresholds.

Average Annual Functional Connection Duration
We compared the range of possible functional connectivity elevation thresholds to water level data at the six WSC stations to calculate average annual functional connection times for lakes with low uncertainty elevation thresholds (<±0.5 m; n = 908).Given that each station observed water level data over a different window of time, we selected years of data with temporal overlap between stations.More specifically, for each station, we included only years that were missing fewer than 20 days of water level data between May and September.We also only kept years of water level data where at least five out of the six stations had a relatively complete record (Table S2 in Supporting Information S1).This filtering process resulted in 15 years of water level data (2001-2012 and 2014-2017).Using these data, we created an average daily water level timeseries between May 1st-September 30th for each station.Then, for each lake, we identified the stations used to calculate its elevation threshold-if there was more than one station within 50 km, we averaged together those stations' water level timeseries.We then compared that timeseries to the range of possible functional connectivity elevation thresholds for that lake, observing how many days the water level at the gage was above the minimum, median, and maximum possible functional connectivity elevation threshold.

Resuspension
High sediment lake water is derivable both from an influx of riverine water and from resuspension of existing lacustrine sediment.To assess the frequency of resuspension, we compared lake reflectance to same-day nearby channel reflectance in the near-infrared band, looking for observations when lake reflectance was much higher than channel reflectance.Nearby channel reflectance was calculated by pairing lake polygons with delta island polygons (Piliouras & Rowland, 2020) that they are contained within, and then analyzing channel pixels surrounding that island (isolated from other waterbodies using a channel mask raster from Piliouras & Rowland, 2020) For each Landsat observation of channels surrounding each island, we calculated the mean reflectance in the near-infrared band.
Next, we paired same-day near-infrared lake reflectance with same-day nearby channel near-infrared reflectance, resulting in a total of 1,694,790 matchups (60% of all observations) of 10,024 lakes (97% of all lakes), and divided lake reflectance by channel reflectance, which we refer to as the lake-to-channel ratio (Dolan et al., 2021).We flagged any observations in which both the lake-to-channel ratio was greater than 1.3 and the connectivity classification was a class 2, indicating high lake reflectance relative to nearby channel reflectance.If there was more than one flagged observation of a lake, we randomly chose one observation per lake, and then imported all flagged data into a Google Earth Engine API for further investigation.We manually inspected each Landsat image associated with each observation to derive the cause of the high lake-to-channel ratio.While it is challenging to visually determine whether or not we are observing resuspension, it is often more possible to pick out when something is definitely not resuspension; for example, a small delta-adjacent channel that was not in the channel water mask might bring sediment into lakes but would not be captured as connected in our primarily algorithms.Therefore, we classified each observation as either "definitely not resuspension," "possible resuspension," "cloud" (issues caused by thin cirrus clouds or odd atmospheric corrections), or "ice" (ice that wasn't properly filtered from the image).We note that this method likely results in an overestimation of resuspension in our observations.

Accuracy Assessment
Classification accuracy for the testing data set is high (Gwet's AC 2 = 0.79).Misclassifications typically occurred along class boundaries, which is to be expected given that functional connectivity is a continuous spectrum that we divided into three discrete classes (Table 1).Classification accuracy is not meaningfully different between the two periods that were used in training/testing data set creation.
We also investigated the occurrence of misclassifications due to sediment resuspension by flagging observations where lake reflectance was much higher than channel reflectance.We initially flagged 0.1% of all lake-to-channel ratios, 1955 observations of 581 distinct lakes, and manually investigated one flagged observation per lake.Of the 581 observations, 230 (40%) were classified as possible resuspension, 210 (36%) were definitely not resuspension, 134 (23%) were caused by issues with cloud-filtering or atmospheric correction, and 7 (1%) resulted from ice filtering errors.Assuming that the resuspension analysis of the individual observation of each lake matches all other observations of the lake that were flagged (total flagged observations = 1,955), a total of 535 observations could be at risk for possible resuspension, which is 0.03% of all lake-to-channel ratio observations.Given this result, we conclude that resuspension is not a major issue for functional connectivity classification of Mackenzie Delta lakes.

Spatiotemporal Patterns of Functional Connectivity
Lakes in the delta experienced the highest functional connectivity in the month of June, corresponding to the spring freshet, with connectivity increasing down-delta toward the coast (Figures 4a-4d; summarized by island in Figures S6a-S6d in Supporting Information S1).Connectivity then decreased across much of the delta during July through September.Functional connectivity variability through time followed a similar pattern, with June experiencing the highest variability, likely due to variations in the timing of freshet arrival (Figures 4e-4h; summarized by island in Figures S6e-S6h in Supporting Information S1).Functional connectivity and its variability remained relatively high near the coast throughout the open water season, possibly due to the influence of storm surge-driven connectivity events (Kokelj et al., 2012;Marsh & Schmidt, 1993;Vermaire et al., 2013).Statistically significant (p < 0.05) trends in freshet connectivity depend on the time period but were generally uncommon (Table 2).The total lake area observed during the 2 weeks after peak discharge in this study was 2664.4 km 2 .An area of 1400.7 km 2 was associated with lakes of an average connectivity of class 1 or higher during that same 2-week window.However, connected lake area was highly dependent on the class threshold used to define "connected" (Figure 5).Using the class 1 connectivity threshold, we found that approximately 53% of observed deltaic lake area was available to receive high sediment river water during the spring freshet.Given a freshet water level increase of 1 m throughout the delta, we estimated an approximate freshet volume storage of 1.38 km 3 in connected lakes.Based on a mean annual flood volume of 19.9 km 3 (Piliouras & Rowland, 2020), connected lakes can, therefore, hold approximately 7% of mean annual flood volume.Additionally, it is notable that the relationship between connectivity threshold and lake area is linear during the 2 weeks after peak discharge, indicating the continuous spectrum of connectivity in the Mackenzie Delta, different from more bimodal connectivity systems like the Colville Delta in Alaska (Dolan et al., 2021).However, as the summer progresses, lakes skew toward lower functional connectivity (class 0; Figure 5).

Functional Connectivity Elevation Thresholds
The majority of lakes within 50 km of a WSC gauge had temporally variable connectivity (n = 6,325, 69%).Of those lakes, only 305 did not experience statistically significant differences in connectivity based on water level.Out of the remaining 6,020 lakes, we were able to calculate ranges of functional connectivity elevation thresholds for 5,989 lakes, and we acquired ranges narrower than 1 m (uncertainty of ≤±0.5 m) for 908 lakes, narrower than 1.5 m (±0.75 m) for 2,377 lakes, and narrower than 2 m (±1.0 m) for 3,224 lakes.For these lakes, elevation thresholds peaked in upstream regions of the delta and decreased downstream (Figure 6).
We assessed the accuracy of our functional connectivity elevation thresholds by comparing them to sill elevations from Marsh and Hey (1988).Unfortunately, many of the lakes studied by Marsh and Hey (1988) were too small to be observed using remote sensing or were difficult to locate on a map (n = 87).Of the remaining 54 lakes contained in both data sets, we had calculated elevation threshold ranges for 35 lakes, only seven of which had  Connectivity Following Spring Freshet in Three Distinct Time Periods: 1984-2001, 2002-2019, and 1984- ranges narrower than 1 m (±0.5 m; Figure 7), which is similar to the ±0.5 m estimated error of Marsh and Hey (1988).Out of these 35 lakes, the possible elevation ranges of 26 lakes overlap between the two data sets; and out of the 7 lakes with lower uncertainty, five ranges overlap.

Table 2 A Summary of Trends in
The remaining 19 lakes were either always low functional connectivity (n = 15), always high functional connectivity (n = 3), or had no significant relationship between functional connectivity and discharge (n = 1).We compared sill elevations from Marsh and Hey (1988) between the former two consistent functional connectivity groups.While there are limited observations, we note that low functional connectivity lakes have a much higher sill elevation compared to those lakes that are always high functional connectivity, which is to be expected given that lower sills lead to more frequent inputs of high sediment river water (Marsh & Hey, 1988; Figure 8).Water Resources Research 10.1029/2023WR036614

Average Functional Connectivity Duration Based on Water Level Climatology
Given the low topographic relief of the delta (Hopkinson et al., 2011), average functional connectivity duration was highly dependent on which part of the elevation threshold range was used in the calculation.Specifically focusing on lakes with narrow elevation threshold ranges (<1 m; n = 908 lakes), we observed that using the upper end of the range to define connectivity duration resulted in much shorter connectivity durations compared the low end of the elevation threshold range.While many lakes are connected between 15 and 60 days per summer regardless of assumption, those lakes connected from 0 to 14 days decreased from 40 lakes using the high elevation threshold value to only 6 lakes using the lowest possible elevation threshold value.Additionally, 52 lakes were connected for more than 60 days using the highest possible elevation threshold value, which increased to 644 lakes using the lowest possible elevation threshold value (Table 3).

Discussion
This study represents a meaningful step toward understanding both water movement through and sediment storage in the Mackenzie River Delta on large spatiotemporal scales.Our work builds on a long and in-depth lineage of structural connectivity work in the delta based on sill elevations (e.g., Lesack & Marsh, 2010;Marsh & Hey, 1988, 1989) and channel absence/ presence detection via Landsat (Piliouras & Rowland, 2020).Structural connectivity is a particularly useful metric for understanding the movement of both water (Marsh & Hey, 1988, 1989) and aquatic animals (Amoros & Bornette, 2002) and correlates to the input of sediment and nutrients into deltaic lakes (Marsh et al., 1999;Spears & Lesack, 2006;Squires et al., 2002;Squires & Lesack, 2002, 2003).However, given that sediment-driven light limitation along with concurrent nutrient inputs are particularly important because they determine primary production sources in Mackenzie Delta lakes (Squires et al., 2002(Squires et al., , 2009;;Squires & Lesack, 2002, 2003), functional connectivity is a more direct metric with which to track the movement of this high sediment water.However, given the lack of prior functional connectivity data in the delta, we compared large-scale observations of functional connectivity to these previously developed structural connectivity data sets.
From a structural standpoint focused on lakes with channel pixels detectable in Landsat imagery, Piliouras and Rowland (2020) found that the Mackenzie Delta had an average connected lake area of 637.1 km 2 , compared to a total lake area of 3201.9 km 2 .Using an estimated water level increase of 1 m and a calculated average annual flood volume of 19.9 km 3 , they found that connected lakes could store 3.2% of the Mackenzie River's flood volume, higher than many other major Arctic deltas except the Colville Delta.While our method observed slightly fewer lakes overall and therefore a slightly smaller total lake area (2,664 km 2 ), our functionally connected lake area, 1,400.7 km 2 , is 763.6 km 2 higher than Piliouras and Rowland (2020)'s estimate.This difference is likely explained by contrasting methodologies, in which Piliouras and Rowland (2020) define connectivity structurally based upon detection of channels in 30 m resolution Landsat imagery, which may under detect connectivity that occurs as overbank flooding or within narrow sub-Landsat pixel channels.Methods described in this paper are based on lake reflectance as opposed to channel detection and are therefore more adept at handling these more challenging scenarios.Additionally, using the same lake water volume estimation technique as Piliouras and Rowland (2020), we found that connected lakes stored approximately 7% of mean annual flood volume, making the delta's storage percentage higher than even the Colville (Piliouras & Figure 7. Median functional connectivity elevation thresholds compared to sill elevations from Marsh and Hey (1988).Gray error bars reflect expertestimated uncertainty for sill elevation measurements and the full possible range for functional connectivity elevation threshold values.Black dots reflect those lakes with elevation threshold ranges narrower than or equal to 1 m, whereas red dots are those lakes with ranges wider than 1 m.Marsh and Hey (1988) summertime sill elevations between those lakes that do not have functional connectivity elevation thresholds because they are always high functional connectivity (n = 3) and those that are always low functional connectivity (n = 15).

Water Resources Research
10.1029/2023WR036614 Rowland, 2020).This finding highlights the importance of connected lake storage in sediment and water movement through the delta toward the coast and indicates that lakes with either no structural connection or a sub-Landsat pixel width connection may be important attenuators of sediment during the month preceding peak freshet.While methods used to calculate connected lake volume are comparatively simplistic, given available data, they provide a consistent framework through which to compare floodwater attenuation across multiple deltas.
Additionally, structural sill elevation measurements were made in 141 lakes during frequent aerial surveys of the region in the late 1980s with an accuracy of ±0.5 m (Marsh & Hey, 1988, 1989).The functional connectivity elevation thresholds in this paper represent a 767-lake expansion of sill elevation-like observations, also within an uncertainty of ±0.5 m or less.Adding to that number are those lakes that are always either high or low functional connectivity, and therefore do not have meaningful elevation thresholds (n = 2,899), expanding information useful to high sediment water movement through the delta to a total of 3,807 lakes.For example, it would be possible to predict which lakes would receive high sediment water based on water level in surrounding channels, which would be relevant for light availability (Squires & Lesack, 2002), phytoplankton productivity (Emmerton et al., 2008;Squires & Lesack, 2002) and bacterioplankton productivity (Spears & Lesack, 2006).
This work highlights the sensitivity of average functional connectivity duration both to water level and to the accuracy of sill elevation or elevation threshold measurements.Increased elevation threshold accuracy could come from two avenues: (a) High temporal frequency water level measurements at each individual lake and (b) More frequent remotely sensed functional connectivity observations.The recently launched Surface Water and Ocean Topography (SWOT) satellite uses Ka-band radar altimetry to provide lake and channel water surface elevations four to six times every 21 days in the Mackenzie Delta, regardless of cloud cover (Altenau et al., 2021;Biancamaria et al., 2016).SWOT measurements could be used to better calibrate daily WSC station data to the water level at or near each individual lake.While the recent launch of the Landsat 9 mission increases temporal observation frequency at 30 m resolution (Lulla et al., 2021), Sentinel-2 now also provides 10 m 5-day temporal frequency observations of lakes, not accounting for the impact of clouds.While the Sentinel-2 record is still relatively short, going forward, its data will allow for more frequent observations of lakes, as well as of smaller lakes.
We detect few statistically significant trends in connectivity during the period following annual freshet, a finding that is well-aligned with nonexistent trends in discharge volume over a similar period (Lesack et al., 2014;Rood et al., 2017).Freshet timing upstream of the delta at Arctic Red River has shifted from June into May due to earlier snowmelt (H.Yang et al., 2015), but detection of connectivity in May remains challenging due to the presence of melting lake ice.The sensitivity of connectivity to water level indicates that predicted future changes to the northern Mackenzie River hydrologic cycle, such as increasing winter and spring streamflow at the expense of summer streamflow (Scheepers et al., 2018) and increasingly frequent and intense future precipitation events (Kuo et al., 2020), are cause for future monitoring.Additionally, permafrost thaw (Lantz & Kokelj, 2008) and resulting increased erosion rates (Rowland et al., 2023), and subsidence (Forbes et al., 2022) are predicted to influence the Mackenzie River Delta region going forward.Each of these projected changes may influence functional and structural lake-to-channel connectivity in different ways; with changes to streamflow resulting in higher and earlier springtime connectivity but decreased summertime connectivity, increased erosion leading to Note."Highest elevation threshold" refers to the connectivity duration calculated using the upper end of each lake's elevation threshold range."Median elevation threshold" and "Lowest elevation threshold" refer to the connectivity duration calculated using the mid-point and lower end of the elevation threshold range, respectively.
higher TSS and therefore higher functional connectivity, and increased subsidence leading to shifting lake sill elevations and functional connectivity elevation thresholds.
Outside of the Mackenzie Delta, future discharge and suspended sediment concentrations are predicted to vary regionally in response to climate change.Controls on these processes include shifts in precipitation, temperature, and permafrost degradation (Feng et al., 2021;Li et al., 2020;Pavelsky & Smith, 2006), along with deforestation (increases discharge and suspended sediment) and dam construction (decreases suspended sediment, modulates discharge) (Döll et al., 2009;Hou et al., 2023;Li et al., 2020).Across the arctic, particularly in Northern Eurasia, discharge has increased (Durocher et al., 2019;Feng et al., 2021;Li et al., 2020) and is expected to continue increasing over the next century (Moragoda & Cohen, 2020;Stadnyk et al., 2021), with suspended sediment concentrations following suit in many places (Moragoda & Cohen, 2020).Whereas in places like Central Asia, dam construction and decreasing precipitation have had the opposite effect (Li et al., 2020;Moragoda & Cohen, 2020).Given that we find that floodplain lake connectivity duration is highly sensitive to water levels, these predicted shifts in discharge and suspended sediment concentrations will likely impact connectivity processes.However, because we do not currently have functional connectivity elevation thresholds or sill elevations across geographically diverse landscapes, it is challenging to quantify the influence of climate change on lake-tochannel connectivity more globally.
Methods described in this paper may be applicable to other lake-rich deltaic or floodplain regions depending on whether suspended sediment concentrations are high enough to distinguish between riverine water and isolated lake water.Pinpointing a specific sediment concentration threshold required for this method is challenging given the lack of in situ lacustrine TSS observations in the Mackenzie Delta.However, given the Mackenzie River's average TSS during the months of June and July taken at the water surface (222 mg/L; ArcticGRO-McClelland et al., 2024), we can infer that the Yukon River (305 mg/L; ArcticGRO-McClelland et al., 2024) may also be a good candidate for application of these methods.Additionally, prior research in the Colville Delta (Dolan et al., 2021) and the Peace-Athabasca Delta (Long & Pavelsky, 2013) used lacustrine water reflectance for functional connectivity detection, suggesting that these locations may also be ideal sites for use of these methods.Application of the methods described in this paper may be more challenging in deltas associated with rivers with lower TSS concentrations over the same period such as the Lena (56 mg/L), the Kolyma (81 mg/L), the Ob' (47 mg/L), and the Yenesey (9 mg/L) (ArcticGRO-McClelland et al., 2024).Other regions globally with floodplain lake-river interaction may also be viable for these methods depending upon the river's sediment levels.Such viable regions may include the Mompox Depression of the Magdalena River in Colombia (Restrepo & Kjerfve, 2000;Restrepo & Syvitski, 2006) or the Poyang Lake-floodplain system in China (Tan et al., 2021).
While not applicable to all lake-rich deltas and floodplain regions, the functional connectivity elevation methods described in this paper provide an opportunity to track water and sediment movement through systems that are often too large or remote to study in their entirety in the field.Specifically, there is an opportunity to pair functional connectivity elevation thresholds with future predicted water levels to better understand how connectivity and therefore sediment and water fluxes may shift in response to climate change.

Conclusion
In this paper we shared a novel data set describing functional lake-to-channel connectivity over a 38-year period across 10,362 lakes in the Mackenzie River Delta.Using this data set, we demonstrated the importance of lake-tochannel connectivity for attenuating floodwater (7% of average annual floodwater) and calculated functional connectivity elevation thresholds that can be used to model sediment movement through the delta (n = 908 lakes with low uncertainty elevation threshold and n = 2,899 lakes with no meaningful elevation threshold because their connectivity is unchanging through time).Going forward, water surface elevations from SWOT can be used to constrain uncertainty bounds and increase accuracy of estimates of functional connectivity elevation thresholds across the delta.Additionally, incorporating high spatiotemporal resolution Sentinel-2 optical imagery into the connectivity detection algorithm described above could allow us to better study both smaller lakes and how individual high frequency events, such as coastal storm surge or upstream rainstorms, impact functional connectivity at shorter timescales.

Figure 1 .
Figure 1.Map of major Mackenzie Delta channels in blue (SWORD, Altenauet al., 2021), study lakes in white (UCLA Lake Database; J.Wang et al., 2023), Water Survey of Canada gauges used for water level (black triangle) and discharge (black square) analyses, and (inset) the location of the delta relative to North America.
be split further).Overall model accuracy and Receiver Operating Characteristics Area Under the Curve were relatively stable regardless of model parameters.The optical properties used in the final model, in order of decreasing importance, were NIR-SWIR1, NSMI, Blue/NIR, and Green/Blue (FigureS4in Supporting Information S1

Figure 2 .
Figure 2. (a-e) Total number of Landsat-derived connectivity observations of each lake from May through September.

Figure 3 .
Figure 3.A conceptual figure demonstrating two methods of functional sill elevation range calculation; (a) shows a scenario where water levels associated with class 0 and class 2 connectivity observations do not overlap while (b) shows a scenario where water levels do overlap.

Figure 4 .
Figure 4.A summary of (a-d) mean monthly functional connectivity classifications and (e-h) the variability of that connectivity over the entire study period of 1984-2022.

Figure 5 .
Figure 5.Total area of lakes greater than a certain mean connectivity class threshold during consecutive 2-week periods after peak discharge.

Figure 6 .
Figure 6.(a) Median functional connectivity elevation threshold defined as the midpoint of the range of possible elevation thresholds.Inset-location of lakes with uncertainty ≤±0.5 m.(b) Uncertainty in elevation threshold as defined as the range of possible elevation threshold values divided by two.(c) Connectivity classifications for lakes with no calculable elevation threshold.

Table 1
Confusion Matrix for Testing Data Classification DOLAN ET AL. 2019 DOLAN ET AL.

Table 3 A
Summary of Average Summertime Functional Connectivity Durations Calculated Using Varying Functional Connectivity Elevation Thresholds