Functional Lake‐to‐Channel Connectivity Impacts Lake Ice in the Colville Delta, Alaska

Within Arctic deltas, surficial hydrologic connectivity of lakes to nearby river channels influences physical processes like sediment transport and ice phenology as well as biogeochemical processes such as photochemistry. As the Arctic hydrologic cycle is impacted by climate change, it is important to quantify temporal variability in connectivity. However, current connectivity detection methods are either spatially limited due to data availability constraints or have been applied at only a single time. Additionally, the relationship between connectivity and lake ice is still poorly quantified. In this study, we present a multitemporal classification of functional lake connectivity in the Colville River Delta, Alaska. We introduce a connectivity detection algorithm based on remote sensing of high sediment river water recharge that is expandable to other high sediment Arctic deltas. We compare results to three existing data sets, producing 64.4%, 75.8%, and 85.2% lake classification accuracy. Mismatches between the data sets are often due to inconsistencies in methodology or definition of connectivity. Connectivity varies temporally in about 10% of studied lakes and correlates strongly with discharge and lake elevation, supporting the idea that future changes in discharge will be a driver of changes in connectivity. Highly connected lakes start and end ice break up an average of 26 and 17 days earlier, respectively, compared to lakes that are poorly connected. Because spring and summer ice conditions drive Arctic lake photochemistry processes, our research suggests that surface connectivity is an important parameter to consider when studying biogeochemistry of Arctic delta lakes.

Functional connectivity also impacts lake ice. While research on the relationship between lake ice and connectivity is in its early stages, several observations suggest that connectivity may be a driving factor in intraregional lake ice timing variability. For example, as observed in Canada's Mackenzie River Delta during the fall freeze-up period, lakes with channel connections are more likely to be impacted by water movement from lakes to channels. This water movement results in lower water levels, increased likelihood of bottom-fast ice, and decreased lake bottom temperatures during the winter (Ensom et al., 2012). In a study of Alaskan lakes, Arp et al. (2013) found that while air temperature and lake area explained 80% of the variance in ice-out timing, they hypothesize that connectivity may be a driving force behind additional variability in ice-out timing within regions but emphasize the need for further studies focusing directly on investigating the topic. Lake ice is particularly important from a biogeochemical standpoint because it is an effective blocker of light (Vione & Scozzaro, 2019). Earlier ice-out timing is very effective at increasing light availability to the lacustrine water column because it usually occurs during the late spring/early summer maximum in solar radiation (Vione & Scozzaro, 2019). A relationship between ice timing and connectivity would thus drive differences in aquatic photochemistry based on lake connectivity. The relationship between connectivity and lake ice is also potentially important from the perspective of methane production in lake sediments. Within ice-covered lakes, greenhouse gasses, particularly methane, accumulate under and within ice and are released to the atmosphere as a large pulse during breakup (Denfeld et al., 2018;Karlsson et al., 2013;Phelps et al., 1998;Vachon et al., 2017). Connectivity-impacted changes in ice timing could influence the timing of this pulse. Analysis of four lakes in Sweden (Boereboom et al., 2012) found that a lake connected via channel to the river network had almost no methane bubbles trapped in the ice. This contrasted with the other three lakes that were rich in methane bubbles and that were either completely isolated or only connected via channel to another lake. Thus, understanding a delta lake's functional connectivity is essential to understanding its physical characteristics, and, by extension, its biogeochemistry.
Despite the importance of connectivity in deltaic ecosystems, previous studies have focused largely on either relatively small scales, a single instance in time, or both. Prior research has used in situ data (e.g., Remmer et al., 2020), remote sensing (e.g., Piliouras & Rowland, 2020), or a combination of the two (e.g., Lesack & Marsh, 2010;Long & Pavelsky, 2013;Marsh & Hey, 1989) to identify lake connectivity. One common method of connectivity identification involves pairing discharge records with lake sill elevations calculated via aerial imagery (e.g., Lesack & Marsh, 2010;Marsh & Hey, 1989). While this method effectively combines current lake sill elevations and historical discharge records to study past connectivity, it assumes that sill elevations are stable over time, which may not hold true given increased erosion, deposition, and subsidence due to permafrost loss (e.g., Lauzon et al., 2019;Nitzbon et al., 2019). Additionally, because historically this method has required high resolution aerial surveys or fieldwork to determine lake sill elevations, it has been challenging to apply on large scales across multiple deltas. However, in the future, ArcticDEM (Porter et al., 2018) may be used to spatially expand the viability of this method. In contrast, Piliouras and Rowland (2020) used Landsat (30 m resolution) 10.1029/2021JF006362 3 of 17 imagery to detect the presence of channel connections between lakes and the main channel network. While this works well for lakes connected via large channels and is scalable across deltas and time periods, it is unable to detect connectivity driven by channels narrower than 30 m without future incorporation of higher resolution imagery such as Sentinel-2 or Planet.
In this paper, we seek to understand the variability of lake-to-channel connectivity over decadal scales, including lakes with sub-Landsat-pixel channels, and analyze the relationship of that connectivity to lake ice breakup and freeze-up phenology. We develop and test our algorithm within the Colville Delta, Alaska because there is ample comparison data available there.

Study Area
We acknowledge that, first and foremost, the Colville River Delta is home to the village of Nuiqsut, an Iñupiat, or, more specifically, a Kuukpikmiut (people of the Colville River) community (population 492 as of 2020; Department of Commerce, Community, and Economic Development, 2021), which was founded in 1974 when 27 Utquiavik families moved to the area (http://www.north-slope.org/our-communities/nuiqsut). Subsistence activities are one of the primary industries of the village, including fishing, hunting, and whaling (Brubaker et al., 2014). Both community knowledge and quantitative climate data paint a picture of a changing delta, with temperatures rapidly increasing (Applied Climate Information System, 2021;Brubaker et al., 2014). During the study period for this analysis (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019), mean annual air temperature in Nuiqsut increased by an average 0.15°C per year (same results using both a linear model and Sen's slope, p < 0.05 for both), with larger temperature changes in the late winter months such as February (Linear Slope: 0.37°C/yr, Sen's Slope: 0.42°C/yr) (Applied Climate Information System, 2021). Residents share that thaw of permafrost below lakes has led to lake draining and drying (Brubaker et al., 2014). This change in permafrost, in combination with later ice freeze-up and earlier ice breakup have reportedly led to dangerous travel conditions both on land and over ice roads (Brubaker et al., 2014). Residents also express concern with loss of habitat due to lake draining as well as loss of accessibility to certain fishing locations due to channels drying up and reducing connectivity (Brubaker et al., 2014). However, more frequent extreme events, such as elevated flood water levels (Brubaker et al., 2014), could potentially result in increased short-duration connectivity.
The Colville River originates in the Brooks Range, flows west to east across the North Slope, and then about 20 km downstream of Umiat, Alaska, takes a northward turn and proceeds ∼140 km to drain into the Beaufort Sea ( Figure 1). The entire Colville River Basin is underlain by continuous permafrost (Mikhailova, 2009;Walker, 1998), except for a few taliks under deep lakes and deeper parts of distributary channels within the delta (Walker et al., 1987). The delta itself is lake-rich, and estimates of the delta's size range from 550 to 665 km 2 (Mikhailova, 2009). While recent estimates of sediment concentration are limited, Walker (1998) indicates that the Colville River Delta has one of the highest suspended sediment concentrations compared to other large Arctic rivers, with a 1998 estimate of 362 g m −3 , compared, for example, to the Yukon (882 g m −3 ) and the Mackenzie (149 g m −3 ) river deltas. Estimates of suspended load grain size range from 3.1 to 12.1 μm, with a median size of 7.6 μm (Arnborg et al., 1967), or very fine silt to fine silt according to the Wentworth Grain Size Classification (Wentworth, 1922). The delta channels remain frozen from approximately late September/early October to late May/mid-June, when the spring discharge maximum drives river ice breakup (Mikhailova, 2009). While the delta itself only receives 300-400 mm of precipitation each year (Mikhailova, 2009), historical records show that river water level can increase by 3-5 m during breakup, inundating at least 65% of the delta's total area (Walker, 1975).

Development of Landsat Water Masks
To create water masks based on Landsat imagery, we use summertime (June to September) 30 m Landsat 5/7/8 Tier 1 Surface Reflectance imagery from 2000 to 2019. To reduce noise from clouds and snow/ice, we use only those Landsat scenes that are less than 75% cloudy overall. We further remove pixel-level cloud and cloud shadow influence using bits 3 (cloud shadow) and 5 (cloud) of the quality band ("pixel_qa"), and we use bit 4 to remove snow/ice. Additionally, we use the Dynamic Surface Water Extent (DSWE) algorithm to select only high confidence water pixels in each image (Jones, 2019) and to remove pixels with unreasonable reflectance values (>10,000 or <0). Lastly, because the delta is sometimes observed by multiple Landsat tiles per day, we create a daily mean composite image. We refer to images that have gone through this process as "filtered Landsat images." Additionally, to create a water mask that only contains channel pixels, we begin with the JRC Global Surface Water Mapping Layers (v1.3, Pekel et al., 2016), 30 m resolution, to select pixels that have a water occurrence value of 60% or greater over the Landsat record. We select this relatively low initial threshold to keep channel connectivity intact through the following steps. Next, to remove lakes from the mask, we buffer lake polygons (Alaskan Lake Database, J. Wang et al., 2011) by 60 m (two Landsat pixels) and remove these regions from the channel mask. However, this leaves behind waterbodies too small to be included in the Alaskan Lake Database. To remove these lakes, we perform a connected components analysis to remove any disconnected features of the mask (e.g., lakes) that are smaller than 300 pixels. Lastly, we apply an 80% water occurrence threshold mask to keep only those pixels that are frequently classified as water over the Landsat record. However, certain lakes are still included in this mask due to their small size and proximity to the channel network, so we manually remove these lakes from the final channel mask.

Lake Boundaries
To identify lakes for connectivity analysis, we select lakes from the Alaskan Lakes Database (J. Wang et al., 2011) that contain at least 100 pixels classified as water 90% of the time and are within the bounds of the Colville Delta (v1.3, Pekel et al., 2016, n = 120 lakes). Note that we have split up several large lake complexes that were originally grouped into one polygon into single lake polygons.

Comparison Data Sets
We compare lake connectivity classification results against three different data sets to assess the strengths and weaknesses of our algorithm.
First, we use the high resolution (<10 m) Google Earth Composite Image (GECI) (a composite of 2013-2016 imagery) to manually identify those lakes that have a visible channel connection. Lakes whose channel connection is unclear (e.g., lakes that look like they have a dry channel that perhaps could fill up during floods) are marked as "uncertain channel connection." This data set is used for both algorithm training and testing. The other two data sets define connectivity slightly differently, hence we use the terminology of comparison, rather than validation.
We compare connectivity classifications to a data set from Piliouras and Rowland (2020) that is based on summer Landsat images from 2014 that defines connectivity based on Landsat-visible channel presence. This data set is originally stored in an image in which waterbodies connected to the channel network have values of 1 and waterbodies that are not connected are masked. To convert this image to a connectivity classification for each lake, we erode the study lake polygons by one Landsat pixel (30 m), and then count connected water body pixels inside each polygon. Polygons with 10 or more connected pixels are considered connected in the comparison data set. If there are fewer than 10 connected pixels, the lake is not considered to be connected. A threshold of 10 is used to further account for any differences in lake shorelines between our study lake polygons and the comparison data set. However, there are no differences in the number of lakes classified as connected using thresholds ranging from zero to 18 pixels.
Lastly, to compare our results to a ground-based data set, we use connectivity classifications from a 1992 survey of the Colville Delta (Jorgenson et al., 1997). Using the print map contained in the manuscript, we manually match connectivity classifications to our study lakes. Lakes in this data set that correspond with our study lakes are classified into the following categories (full descriptions can be found in Table 4 of Jorgenson et al., 1997): Medium-high connectivity 1. Deep connected lake: Lakes ≥1.5 m in depth with long, narrow, and persistent connections to the delta channel network. 2. Deep tapped lake with low-water connection: Same as prior, but connection channels are filled with water, even at low discharge. 3. Deep tapped lake with high-water connection: Lakes ≥1.5 m in depth that have been "tapped" by erosion of adjacent river banks. Connecting channels are dry unless discharge is high. 4. Shallow tapped lake with low-water connection: Same as prior, but connection channels are filled with water even at low discharge. 5. Shallow tapped lake with high-water connection: Lakes <1.5 m in depth that have been "tapped" by erosion of adjacent river banks. Connecting channels are dry unless discharge is high.
No connectivity 1. Deep isolated lake: Lakes ≥1.5 m in depth that have no distinct outlets and are not connected to the channel network. 2. Shallow isolated pond: Lakes <1.5 m in depth with no distinct outlets or channel connections. 3. Shallow isolated pond, riverine: Lakes <1.5 m in depth with no distinct outlets or channel connections that are associated with old river channels.
No connectivity data available 1. Brackish ponds: Shallow ponds near the shoreline that may be impacted by very high tides or storm surges. May or may not have connections to the channel network.

ArcticDEM
We investigate patterns between connectivity and lake surface elevation using 2 m-resolution ArcticDEM strip data (Porter et al., 2018) in 2015 and 2017. ArcticDEM data is only available when there is sunlight. Additionally, availability of data in the summer months is limited because the processing method is not effective over open water. As a result, there is a small window of time each year in which lake elevation data from ArcticDEM is available in this region-the month of April. Therefore, the lake surface elevation in this study represents the elevation of the top of the snow/ice on the lake as opposed to the water surface. For each observed year, we create an April mosaic of the strip data. Next, to remove noise from non-water lake edge pixels, we erode each lake polygon by two Landsat pixels (60 m). This removes six lakes from analysis that are less than two Landsat pixels wide (remaining lakes used for elevation analysis, n = 114). Lastly, for each lake, we calculate the mean 2015 and 2017 elevation within each eroded lake polygon. The vertical reference for the elevation data is the height above the WGS84 ellipsoid.

River Discharge Data
To analyze the impact of discharge on connectivity, we use daily discharge data from the USGS Gage No. 15875000 at Umiat, Alaska between 2003, when the station began collecting data, and summer 2019. We order the years based on maximum summer discharge, defined as the highest daily summer (1 June to 30 September) discharge, and then group years by discharge (Table 1). While Umiat is ∼125 km upstream from the delta, the Colville River has already drained approximately 67% of its catchment area (U.S. Geological Survey, 2021; Walker & Hadden, 1998). We assume that extreme discharge years at Umiat, Alaska coincide with high discharge years within the Colville Delta.

Detection of Connectivity
In this study, we differentiate functionally connected and functionally disconnected lakes based on the relationship between their optical properties and the optical properties of nearby channels in the delta. We extract Landsat reflectance measurements for lakes and river channels, calculate lake to channel optical ratios for each observation, and then use the median and standard deviation of the ratio between lake and channel color/band combinations within four-to-five-year periods (Section 3.2.1.). These ratio values are then used to classify lakes as connected or disconnected within each time period (Section 3.2.2.). We then compare the method to three data sets (Section 3.3) and investigate the relationship between connectivity, discharge, lake surface elevation, and lake ice phenology (Sections 3.4 and 3.5).

Extraction of Paired Lake and Channel Water Color Observations
For each Landsat observation of a lake on a given day, we calculate mean optical property values within the lake polygon and filter out those days when there are fewer than 100 cloud-free pixels observing the lake (see Table 2 for optical values calculated). Next, we search for same-day channel pixels within a one-kilometer buffer of the lake using the channel water mask. If at least 200 pixels are not present, we look within a 2-km buffer, then 5 km, then 10 km, until the 200-pixel requirement is satisfied. Once the requirement is satisfied, we calculate the mean optical property values for those channel pixels within the buffer. While there is some spatial reflectance variability within the delta, the sensitivity of the classification to that variation is small. For each optical property, we then take the ratio between the same-day lake observation and the corresponding channel observation, here on referred to as lake-to-channel ratios.
In Figures 2a and 2b, which cover the same geographic extent on two different days, we show examples of two lakes, one with an obvious channel connecting the lake and the river (A) and one with no connective channel (B). The lake with obvious channel presence routinely has daily Green/Blue ratios and NSMI ratios of close to  one during the period between 2000 and 2004, indicating the similarity of lake-to-channel ratios through time (Figure 2c). The lake with no visible channel presence has a much wider distribution of NSMI ratios and a slightly wider distribution of Green/Blue ratios (Figure 2d). When the river is transporting little sediment, the low connectivity lake and the channel look quite similar (Figure 2a), but when the river water is high in sediment, the low sediment lake looks very different from the nearby channel ( Figure 2b). As a result, we find that it is more effective to use lake-to-channel ratios to classify functional connectivity within a time period as opposed to at a single instance in time. Therefore, within a given time period, we calculate the median, standard deviation, and kurtosis of single-day lake-to-channel ratios within that period and use those values for classification. The different time periods used in this study are described in Sections 3.3. and 3.4.

Variable Selection and Decision Tree-Based Classification
We use a supervised classification method in conjunction with the GECI validation data set, excluding lakes that have uncertain channel connectivity. For each lake, we calculate the median, standard deviation, and kurtosis of each type of lake-to-channel ratio during the period corresponding to the GECI data set (2013-2016).
Next, we use the Boruta method in R (Kursa & Rudnicki, 2010) to select the lake-to-channel ratio variables with the highest ability to differentiate between the high connectivity and low connectivity lakes. Those properties are: NSMI ratio (median/standard deviation), Green/Blue ratio (median/standard deviation), dominant wavelength ratio (median/standard deviation/kurtosis), Red band ratio (median/standard deviation), Green band ratio (median/ standard deviation), and Blue band ratio (median/standard deviation).
Next, we develop a decision tree model using the tidymodels package in R (Kuhn & Wickham, 2020) using the variables selected in the Boruta process. The lakes are split into 25% (testing, n = 27) and 75% (training, n = 77) groups, and the training data are further split into five groups for five-fold cross validation. Using five-fold cross validation, we tune the hyperparameters of the model and select hyperparameters from the model producing the highest overall accuracy, with the following results: cost complexity (0.1), the maximum depth the tree can reach before it is forced to stop growing (8), and the number of data points in a node that are required for the node to split further (2). Using these hyperparameters, we fit the decision tree to the training data set, resulting in the decision tree shown in Figure 3.

Comparison to Existing Connectivity Data Sets
To evaluate the accuracy of our classification algorithm, we apply decision tree to the testing portion of the GECI validation data set and compare the classification from time periods temporally closest to the collection of the Jorgenson et al. (1997) data set and the Piliouras and Rowland (2020) data set (Table 3). We then calculate the agreement of each classification and use visual inspection of GECI and Landsat 5/7/8 imagery to investigate disagreements.

Examination of Temporal, Discharge, and Elevation-Related Patterns in Connectivity
To better understand the variability of functional lake connectivity in the Colville River Delta, we group single day lake-to-channel ratios for each lake into four time periods: 2000-2004, 2005-2009, 2010-2014, and 2015-2019. We also group single day ratios for each lake into four different groups corresponding to increasing levels of maximum summer discharge (Table 1). For each of these periods, we calculate median, standard deviation, and kurtosis values for each of the ratios and classify the connectivity of each lake using the decision tree ( Figure 3). We only analyze lakes with at least ten single-day observations within each of the periods (n = 120 lakes). We also investigate patterns between connectivity and lake surface elevation by using a Mann-Whitney U-test to compare lake surface elevations from 2015 and 2017 (Section 3.1.4) to classified connectivity in the 2015-2019 time period.

Relationship Between Connectivity and Lake Ice
To better understand the influence of connectivity on the deltaic environment, we compare functional connectivity to Landsat-derived lake ice climatology. The lake ice data set was developed using a logistic regression based on a series of 2000-2019 Landsat 5/7/8 lake ice fraction observations via the Yang et al. (2021) lake ice detection algorithm. For each lake, the method pools together all single-day ice fraction data between 2000-2019 and divides the data into two subsets-one for freeze-up period (day of the year 230-365 and 1-58 of the following year) and one for the breakup period (day of the year 59-229). Then two separate logistic regressions were fitted to the two subsets of data to model freeze-up and breakup ice fraction dynamics. These models together classify each day of the calendar year (DOY) as likely ice covered (typical ice fraction >80%), transitional (typical ice fraction 20%-80%), or likely ice free (typical ice fraction <20%). We compare our functional connectivity results to ice phenology events including DOY of breakup start, DOY of breakup end, breakup transition duration (number of days in the spring between 80% ice cover and 20% ice cover), DOY of freeze-up start, DOY of freeze-up end, freeze-up transition duration (number of days in the fall between 20% ice cover and 80% ice cover), ice cover duration (number of days classified as likely ice covered), and ice free duration (number of days classified as likely ice free). Based on the temporal lake connectivity classifications, lakes are divided into three connectivity groups for analysis: those that are always classified as high functional connectivity, those that are always classified as low functional connectivity, and those that have variable functional connectivity through time. We use a Mann-Whitney U-test to compare ice climatology between high and low connectivity groups.

Connectivity Classification
Functional connectivity detection via optical properties of water has certain limitations and potential sources of error. First, this method exclusively works when there is a visible difference in color between river water and lake water during at least some portion of the open water period. While this is true in the Colville Delta, deltas that export higher sediment overall, but have lower sediment concentrations, are difficult to analyze using this method because of indistinct optical differences between river water and lake water.  Lakes may also experience non-SSC-driven changes in color, for example, via high chlorophyll concentrations (Bukata et al., 1983;Jerome et al., 1994). While the DSWE algorithm used in this analysis to identify water pixels also filters out vegetated land pixels as well as submerged vegetation, it may not filter out other high chlorophyll water pixels. However, based on visible inspection of imagery, this seems to be a limited issue in the Colville Delta, with only 2-3 lakes potentially impacted. Lakes near the delta shoreline present challenges in validation because they are likely both low-elevation and brackish, which may impact water color. These lakes commonly do not have clearly visible channels, so are not classified as connected in the GECI validation, and are simply labeled as 'brackish', with no additional connectivity information, in the land survey validation data set (Jorgenson et al., 1997). They are also close enough to the shoreline that backwater may impact the color of these lakes.

Lake Ice Climatology
While individual Landsat-derived ice fraction observations that go into creating the ice climatology model are dense during the breakup period, they are less frequent during the freeze-up period due to low light and cloud conditions ( Figure 4). Therefore, there is higher uncertainty/inaccuracy in climatological freeze-up dates compared to breakup dates. More specifically, for one week prior to the breakup start date through one week after the breakup end date, there are an average of 23 ice fraction observations used in the model, compared to an average of eight for the freeze-up encompassing period. Also, because we are using a climatological timeseries based on 20 years of data, distinguishing between a longer freeze-up/breakup period and a more variable freeze-up/breakup period is challenging.

Comparison to Existing Data Sets
Compared against channel absence/presence classifications from the testing portion of the GECI data set, functional connectivity classification matches with the channel presence classification in 85.2% of lakes (n = 23/27 lakes). Of the incorrect classifications, only one lake classified in GECI as connected is classified as low functional connectivity in our analysis, whereas three lakes without channels are classified as high functional connectivity ( Figure 5a). However, in the case of the lake that we classify as low functional connectivity, the GECI validation data set was actually in error, and the lake did not have a visible channel. In the case of the three lakes that we classified as connected that do not have a GECI visible channel, two may have dry channel beds upon closer inspection, and one is at low elevation near the delta shoreline.
Next, we compared our 2013-2016 connectivity results to the 2014 classification of the Colville River Delta by Piliouras and Rowland (2020) (Figure 5b). This comparison data set is similar to the GECI-derived validation data set because it defines connectivity via channel presence, in this case determined using coarser resolution 30 m Landsat imagery. Our classification agrees with this data set for 75.8% of lakes (n = 29/120). Of the mismatched lakes that we label low functional connectivity but are connected in the Piliouras and Rowland (2020) data set (n = 12 lakes), seven are near the delta shoreline, which is mostly classified as water by Piliouras and Rowland. Upon visual inspection, one of the five remaining lakes appears to be connected via a channel inundated only at high water, a second is the third lake in a series of lakes leading off of the channel so it is hydrologically connected but does not regularly receive high sediment water, and the final three lakes fall within one to two Landsat pixels of larger, connected lakes but do not have their own channels. Of the lakes that our algorithm classifies as high functional connectivity but are not connected in the comparison data set (n = 17 lakes), the majority (n = 15 lakes) have channel connections narrower than one Landsat pixel in GECI imagery.  the lakes with mismatching classifications, 35 are predicted to be low functional connectivity, but are of a high connectedness class in the comparison data, whereas two are differently classified in the opposite direction. However, 23 of these differently classified lakes are deep tapped lakes with high water connections, which accounts for all but one of the lakes in this class in this validation data set. Similarly, nine of the 35 differently classified lakes are deep connected lakes or deep tapped lakes with a low-water connection that either have very long and narrow connections or are several lakes deep within a lake complex, so it takes higher energy to move high sediment water from the channel into these lakes. While we lack discharge records prior to 2003, it is possible that these lakes were rarely connected during our period of analysis. Alternatively, these results may indicate that our algorithm is most suited for detecting connectivity in those lakes that are connected to the channel network even at low water periods. Additionally, the 1992 survey was eight years prior to our first observation, and both lake tapping and lake or channel infilling are possible within this period.

Functional Connectivity Variation With Time, Discharge, and Elevation
Within the Colville Delta, functional connectivity remains consistent in the majority of lakes throughout the four temporal periods between 2000 and 2019 ( Figure 6). Of the observable lakes in the delta that had stable connectivity through time, more lakes were classified as always low functional connectivity (n = 79) than always high connectivity (n = 28). Of those that change connectivity through time (n = 13), three transition from high to low connectivity over time, one moves from low to high connectivity over time, and nine lakes transition between states multiple times. Nanuk Lake, highlighted in Figure 6, switches from high connectivity in the first three periods to low connectivity in the last period. This pattern may be related to the lake's complex infilling processes in conjunction with the timing of Landsat observations during the latter period. An image of period-by-period lake classifications can be found in Figure S1.
Overall, 14 lakes vary their functional connectivity within different discharge periods (Figure 7). We identify eight lakes that always shift from low functional connectivity to high functional connectivity as discharge increases. Of these, seven have narrow channels, visible channels that appear dry in the GECI data set, or connect to the channel networks via a series of lakes, hindering the routine movement of high sediment river water. The remaining lake was tapped during the study period. Additionally, there are five lakes that experienced multiple shifts in connectivity amongst the discharge periods as well as one lake that experienced decreasing connectivity with discharge. This likely indicates the importance of non-discharge related drivers of functional connectivity, such as wind, infilling, and backwater (e.g., Geleynse et al., 2015). Of these lakes that experienced any change in connectivity with discharge, 11 overlap with those lakes that change connectivity over time, six of which show increasing connectivity with increasing discharge.
As expected, high functional connectivity is also associated with significantly (p < 0.01) lower wintertime lake surface elevation in both 2015 and 2017, compared to low functional connectivity lakes (Figure 8, 2.2 m lower and 1.7 m lower, on average, in 2015 and 2017, respectively). This makes sense from a definitional standpoint, as low-elevation lakes are more likely to be flooded during periods of elevated river water level. Additionally, because the measured lake surface elevation reflects wintertime snow or ice cover elevation, the lower elevations of connected lakes may also reflect riverine drawdown of water from during the fall freeze-up process.

The Relationship Between Connectivity and Lake Ice
Within the Colville Delta, functional connectivity strongly correlates to lake ice climatology (Figure 9). While there are statistically significant differences (p < 0.01) between many of the studied ice parameters, the largest differences are found in ice breakup start date and ice breakup end date. Ice breakup in connected lakes begins a mean of 25.7 days earlier and ends 17.3 days earlier than their disconnected counterparts, leading to a difference in total ice duration of 22.0 days, typically during the months of May and June. Connected lakes often start freezing up in mid-to-late September compared to late September to early October for low functional connectivity lakes, though the mean difference in freezeup start and end times are only 8.4 and 3.6 days, respectively.

Discussion
The primary takeaways of this paper are threefold. First, we have demonstrated the efficacy of a widely applicable algorithm for detecting temporal variations in functional connectivity in a high sediment Arctic delta. Second, we find that over a 20-year period, connectivity in the Colville is variable through time in 10% of the studied lakes and is related to discharge as well as lake elevation. Third, we find that functional connectivity may control lake ice timing, especially ice breakup, which has implications for photochemistry and methane production.

Efficacy of the Functional Connectivity Algorithm and Comparison to Existing Data Sets
This functional connectivity algorithm provides a new alternative to other published deltaic lake-to-channel detection methods. Each method is useful in specific situations and has unique pros and cons. First, the "Sill Elevation Method" used by Marsh and Hey (1989) and Lesack and Marsh (2010) uses lake sill elevations in combination with historical discharge observations to assess functional connectivity (by channel or overland flow) duration through time. This method is excellent in situations in which lake sill elevation data, collected via aerial or field survey, is available, and in situations when lake sill elevations are not assumed to change over the study period. While the Figure 7. Lake connectivity classification with increasing discharge. X-axis numbers correspond to discharge groups in Table 1. reliance on aerial imagery or field data currently limits the wider spatial and temporal applicability, this method allows for the most temporally fine-grained analysis of functional connectivity duration and allows calculation of connectivity prior to Landsat observations, if discharge data is available. Second, Piliouras and Rowland (2020) developed a method that uses Landsat imagery to detect direct connections between lakes and channels. The challenge with this method is that it only detects connectivity via channels that are visible within Landsat imagery, which excludes lakes connected by sub-pixel channels or by overland flow. The lake-to-channel ratio method described in our analysis builds on work by Long and Pavelsky (2013) that suggested the usefulness of suspended sediment as a tracer for river water recharge within six floodplain lakes in the Peace-Athabasca Delta. In our laketo-channel ratio method, similar to the Sill Elevation Method, we are able to detect functional connectivity via channels that are sub-Landsat pixel in width, but we are less likely to capture short-duration connectivity due to overland flow, particularly due to the training data set that was used. Our method does not require aerial or field surveys to define sill elevation, and therefore is easier to expand to multiple deltas. Additionally, our method does not make assumptions about the consistency of sill elevations through time. However, unlike the Sill Elevation Method, the lake-to-channel ratio method does not allow calculation of yearly connectivity durations. Compared to the Piliouras and Rowland (2020) method, our method allows for sub-pixel channel connectivity; however, the method relies on river water with high sediment concentrations, limiting application to high sediment deltas. Our method currently produces connectivity classifications within coarser five-year periods. However, going forward, high temporal and spatial resolution optical satellites such as Sentinel-2 and Planet may be used to increase the temporal resolution of these connectivity classifications.

Discharge and Connectivity
Discharge is a major control on changes in functional connectivity, with high discharge corresponding to increased functional connectivity, and low discharge corresponding to reduced connectivity. These discharge-related shifts in connectivity are concentrated within lakes that are connected to the main channel network indirectly (d) Total ice duration (days with >80% ice cover) and total ice-free duration (days with <20% ice fraction). *significant difference at p < 0.01 between always high functional connectivity lakes and always low functional connectivity lakes. Dots represent outliers.
via large lake complexes, long and narrow channels, or occasionally channels that are only water-filled at high discharge. While the North Slope region overall has not experienced a statistically significant trend in discharge over the past 41 years (1975( ) (Durocher et al., 2019, future changes in precipitation type, timing, and amount, along with warmer air temperatures, are predicted to increase discharge on the North Slope over the next hundred years (Bring et al., 2017). Based on the pattern we have observed, projected increases in discharge could result in increased functional connectivity either via increased overland flow, higher sediment river water being transported further into lake complexes via channels, or lake tapping (Walker & McGraw, 2015). This increase in connectivity, particularly near the shoreline, will take place in conjunction with sea level rise, which will likely increase intrusion of saltwater into low-elevation lakes, impacting the ecology and optical properties of those systems (Jorgenson et al., 1997). Alternatively, increased SSC associated with elevated discharge will also result in faster infilling of highly connected lakes (Walker & McGraw, 2015), leading to decreased connectivity. Any future changes in discharge will also be concurrent with degradation of permafrost . Permafrost degradation impacts Arctic lakes in several ways, including increased erosion rates (Piliouras et al., 2021), increased subsurface connectivity due to active layer thickening (Connon et al., 2014;Walvoord & Kurylyk, 2016), and resulting lake area loss and drainage (Brubaker et al., 2014;Marsh et al., 2009;Nitze et al., 2017;Smith et al., 2005) or expansion (Nitze et al., 2017;Smith et al., 2005). Residents of Nuiqsut are already reporting drying tributaries and loss of connectivity due to increased infilling and drainage as well as resulting concerns about fish migration and travel hazards (Brubaker et al., 2014).

Connectivity and Lake Ice
Air temperature is one of the dominant controls of lake ice phenology in conjunction with, to a lesser extent, lake morphometry, latitude, and elevation (Arp et al., 2013;Sharma et al., 2019;Warne et al., 2020;Williams et al., 2004;Šmejkalová et al., 2016). However, this work demonstrates the importance of functional connectivity as a driving factor of lake ice timing. Additionally, our results fit well with initial findings by Arp et al. (2013) that suggests hydrologic connectivity may be responsible for intraregional lake ice variability in Alaska not explained by temperature or morphometric controls. As lake ice continues to be modeled and interpreted both as a proxy for historical climate change (Assel & Robertson, 1995;Magnuson et al., 2000;Robertson et al., 1992;Šmejkalová et al., 2016;Warne et al., 2020) and as a signal of and response to current and future climate change (Arp et al., 2018;Sharma et al., 2019), and as we seek to quantify the impact of climate change on Arctic transportation (Prowse et al., 2009, Prowse, Alfredsen, Beltaos, Bonsal, Bowden, et al., 2011Stephenson et al., 2011), it is important to consider functional connectivity as an input to lake ice models. Additionally, because ice is an effective blocker of light (Vione & Scozzaro, 2019), connectivity-driven differences in ice timing during the breakup period, which experiences very high solar radiation (Figure 10), may also be an important control on photochemistry in these systems.

Conclusions and Future Work
This research demonstrates the efficacy of using remotely sensed water color to detect spatial and temporal variations in functional lake connectivity in the Colville River Delta, Alaska. Results suggest that connectivity is variable in about 10% of delta lakes. Discharge is an important control of connectivity, particularly within those lakes that have indirect channel connections. Additionally, our results suggest a strong relationship between functional connectivity and lake ice phenology, especially observed earlier breakup timing in connected lakes. More fieldwork is needed to quantify the impact of this relationship on photochemistry and carbon processing in these lakes. In the future, we can use methods described in this paper to track connectivity in the Colville and other Arctic deltas and to understand how patterns of connectivity interact with the changing physical and biogeochemical environment of the Arctic.