Investigating lake-area dynamics across a permafrost-thaw spectrum using airborne electromagnetic surveys and remote sensing time-series data in Yukon Flats, Alaska

Lakes in boreal lowlands cycle carbon and supply an important source of freshwater for wildlife and migratory waterfowl. The abundance and distribution of these lakes are supported, in part, by permafrost distribution, which is subject to change. Relationships between permafrost thaw and lake dynamics remain poorly known in most boreal regions. Here, new airborne electromagnetic (AEM) data collected during June 2010 and February 2016 were used to constrain deep permafrost distribution. AEM data were coupled with Landsat-derived lake surface-area data from 1979 through 2011 to inform temporal lake behavior changes in the 35 500- km2 Yukon Flats ecoregion of Alaska. Together, over 1500 km of AEM data, and roughly 30 years of Landsat data were used to explore processes that drive lake dynamics across a variety of permafrost thaw states not possible in studies conducted with satellite imagery or field measurements alone. Clustered time-series data identified lakes with similar temporal dynamics. Clusters possessed similarities in lake permanence (i.e. ephemeral versus perennial), subsurface permafrost distribution, and proximity to rivers and streams. Of the clustered lakes, ∼66% are inferred to have at least intermittent connectivity with other surface-water features, ∼19% are inferred to have shallow subsurface connectivity to other surface water features that served as a low-pass filter for hydroclimatic fluctuations, and ∼15% appear to be isolated by surrounding permafrost (i.e. no connectivity). Integrated analysis of AEM and Landsat data reveals a progression from relatively synchronous lake dynamics among disconnected lakes in the most spatially continuous, thick permafrost to quite high spatiotemporal heterogeneity in lake behavior among variably-connected lakes in regions with notably less continuous permafrost. Variability can be explained by the preferential development of thawed permeable gravel pathways for lateral water redistribution in this area. The general spatial progression in permafrost thaw state and lake area behavior may be extended to the temporal dimension. However, extensive permafrost thaw, beyond what is currently observed, is expected to promote ubiquitous subsurface connectivity, eventually evolving to a state of increased lake synchronicity.


Introduction
Arctic and boreal regions are warming at a faster rate than the global average, enhancing vulnerability to permafrost thaw (Lawrence and Slater 2005, Jorgenson et al 2006, Serreze et al 2009, Screen and Simmonds, 2010. In turn, permafrost thaw has the potential to alter the routing and storage of water above and below ground (Walvoord and Kurylyk 2016). Of particular concern is the impact of permafrost thaw on surface-water dominated regions such as boreal lowlands, which comprise ∼30% of earth's forested area and contain more surface-water area than any other terrestrial biome (Burton et al 2010, Brandt et al 2013. Northern lowland lakes provide critical habitat for migratory arctic shorebirds and waterfowl (Alerstam 2001), fish, and wildlife (Prowse and Brown, 2010), thereby supporting critical ecosystem services to indigenous communities. Lakes also generate feedbacks to the climate system by partitioning sensible and latent heat fluxes and regulating carbon cycling at the landscape level (Tranvik et al 2009, Abnizova et al 2012.
Decadal-scale changes in the abundance and distribution of lakes have been observed in arctic and boreal lowlands. Previous studies identified regional patterns that reflect a general shift towards increasing net lake coverage in continuous permafrost (>90% spatial coverage), such as northern Siberia (Smith et al 2005, Sannikov 2012); and decreasing net lake area in discontinuous permafrost (50%-90% spatial coverage), as in Canada (Labrecque et al 2009, Carroll et al 2011, interior Alaska (Rover et al 2012, Chen et al 2014 and southern Siberia (Smith et al 2005, Vonk et al 2015. Yet recent studies observe patterns in surface-water distribution that do not strictly adhere to these broad generalizations due to complexities in local conditions and processes (Pastick et al 2018, Karlsson et al 2012. This variability highlights the need for the identification of processes that drive surface-water dynamics in heterogeneous permafrost regions. Though all of the above studies acknowledge permafrost thaw to be a driving mechanism for lakesurface area changes, hydroclimatic factors and other influential processes, such as ice-jam flooding and catchment vegetation (Chen et al 2014, Turner et al 2014, Jepsen et al 2016, also require consideration. Quantifying surface-water redistribution and understanding the underlying mechanisms of change in arctic-boreal regions are required for accurate climate predictions, ecosystem assessments, and judicious land management decisions in northern latitudes (e.g. Abnizova et al 2012, Gauthier et al 2015.
The Yukon Flats in Alaska is a lake-rich lowland that encompasses the Yukon Flats National Wildlife Refuge and has been identified as an ecologically-sensitive region poised for major hydrologic change with permafrost degradation (Walvoord et al 2012). Despite low precipitation (∼207 mm with ∼45% falling as snow) associated with a continental climate, a high abundance of lakes persist (∼5.5% by surface area), partially due to the presence of permafrost, which limits subsurface connectivity between lakes and river networks (Arp andJones 2008, Chen et al 2013). The region spans a broadly-mapped transition from continuous to discontinuous permafrost (Jorgenson et al 2008), and hosts a shallow, continuous, fluvial gravel layer (Minsley et al 2012). Thawing of the gravel layer could dramatically enhance subsurface connectivity among inland waters in the Yukon Flats, altering lake distribution and chemistry, although this process is not expected to be spatially uniform, and few studies exist on the topic. A regional remote sensing analysis by Rover et al (2012) identifies lakes in the Yukon Flats that are shrinking, expanding, and stable over the past several decades, revealing a net decline in lake-area extent. Jepsen et al (2013a) combined those observations with airborne electromagnetic (AEM) data from Minsley et al (2012) to evaluate connections between below-lake permafrost conditions and long-term (∼30 yr) lake-area trends. They found no statistically relevant relationships between open vertical taliks (unfrozen conduit connecting sub-permafrost groundwater and the overlying lake) and lake dynamics. Rather, the findings highlighted the importance of shallow permafrost thaw state and suggested that enhanced connectivity in the shallow subsurface (upper tens of meters) may promote lake/groundwater exchange, thereby contributing to lake shrinkage and expansion in thawed areas. In addition to permafrost distribution, inter/ intra-annual variability in ice-jam flooding, precipitation, and evapotranspiration likely influence the nonuniform spatial and temporal behavior of adjacent lakes within the Yukon Flats (Rover et al 2012, Chen et al 2013, 2014, Jepsen et al 2013a. However, the contribution of each of the above processes' spatial and temporal influence on lake distribution remains largely unknown (Chen et al 2014, Jepsen et al 2016, limiting our ability to project regional lake response to ongoing climate change and permafrost thaw.
The goals of this study are to (1) elucidate potential mechanisms that drive lake-area dynamics in regions of variable permafrost and lake talik conditions, and (2) infer likely trajectories of lake-area change expected with permafrost degradation in the Yukon Flats. We address these goals by performing a detailed cluster analysis of Landsat-derived lake surface-area timeseries data, using previous and newly attained airborne geophysical data to constrain spatially variable permafrost conditions. The spatial time-series clusters were co-located with AEM surveys in places, allowing for a subsurface characterization of the permafrost distribution under individual lakes, an in-depth comparison not previously possible with satellite imagery or field measurements alone. This approach provides improved understanding of the physical mechanisms and characteristics governing lake-area behavior. New AEM data acquired in 2016 (Minsley et al 2017) provide an ice-rich, continuous permafrost complement to results from a nearby AEM survey conducted in 2010 that covered discontinuous permafrost in the Yukon Flats (Ball et al 2010, Minsley et al 2012. Comparing mechanisms that drive lake surface-area dynamics in areas containing a spectrum of local permafrost conditions yields insight into the future of permafrost-supported lowland lakes in transition towards an increasingly thawed subsurface

Study area
The Yukon Flats ecoregion is roughly 35 500 km 2 , 225 km northeast of Fairbanks in the Alaskan interior, and is physically bounded by the Yukon-Tanana Plateau to the south, the southern foothills of the Brooks Range to the north, the Hodzana Highland to the west, and the Porcupine Plateau to the east (figure 1; Williams 1962). The Yukon Flats basin is comprised of a broad Cenozoic sedimentary lowland (∼23 300 km 2 ) with a marginal upland (∼12 200 km 2 ) that provides a transition into the surrounding highlands comprised of Paleozoic bedrock with metamorphic and intrusive igneous rocks (Williams 1962). Mid-Tertiary to early-Quaternary clay, silt, and sand underlie surficial deposits within the Yukon Flats. Surficial deposits are comprised of fluvial gravel in the lowlands, and ice-rich loess in the uplands. Fluvial gravel deposits range in thickness from ∼50 to 75 m east to west (Williams 1962, Clark et al 2009 and form alluvial fans, terraces from Yukon River channel migration, and floodplains of Pleistocene to present age. Little topographic relief exists within the lowlands, with elevations that range from approximately 100 m in the west to just over 150 m in the east and an average slope of 0.15 m per km (Gesch 2007). (Gesch 2007) Many rivers and streams are present in the area, including the main channel of the Yukon and several tributaries-most notably, the Porcupine River, Chandalar River, and Birch Creek. Lakes are abundant within the area, ranging in size from tens of m 2 to over 5 km 2 , with vegetation consisting primarily of hardwood and low-growing spruce forests (Viereck and Little 2007).
The Yukon Flats has a continental boreal climate. Seasonal air temperatures reach ∼22°C in the summer and ∼−33°C in the winter, with a mean annual air temperature of −6°C (Nakanishi and Dorava 1994). Mean annual precipitation totals 170 mm throughout the region, with the uplands receiving more precipitation than the lowlands during the summer months (Nakanishi and Dorava 1994). The Yukon Flats encompasses the gradational permafrost boundary from continuous to discontinuous permafrost, resulting in a spectrum of thaw states across the region (Jorgenson et al 2008). Regional active layer thickness varies from 0.2 to 1.3 m according to Pastick et al (2014), however, localized field studies in the Yukon Flats have noted measured depths to permafrost of several meters (Jepsen et al 2013b), emphasizing a high degree of variability in the potential for lateral flow through the shallow subsurface. Permafrost thickness is also highly variable throughout the Yukon Flats with a maximum permafrost thickness of up to ∼100 m in areas devoid of existing or previously existing surface-water bodies (Williams 1962, Minsley et al 2012, Pastick et al 2014.

AEM surveys
Geophysical surveys such as AEM measure physical properties that can be used to help identify the presence or absence of permafrost efficiently over much larger areas than can be surveyed using conventional methods (i.e. push-probe, auger), especially in remote and inaccessible areas. AEM systems transmit an electromagnetic signal from an airborne platform, inducing secondary currents in the Earth that are subsequently measured by receiver coils in the sensor housing (Paine and Minty 2005). The measured secondary magnetic fields are sensitive to subsurface electrical conductivity structure, which is controlled primarily by lithology, liquid water content, salinity, and other factors such as metallic mineralization in the material (Siemon et al 2009). Contrasts in subsurface conductivity in artic-boreal environments stem primarily from frozen material having very low conductivity values and unfrozen materials having relatively high conductivity values. Consequently, AEM surveys have been shown to be a productive tool for mapping permafrost distribution at regional scales (Minsley et al 2012, Mikucki et al 2015 and determining thaw conditions beneath northern lakes (Minsley et al 2015).
Two surveys that crossed more than 6000 lakes were conducted within the Yukon Flats: (1) a survey flown in February 2016 totaling 300 line km in the western Yukon Flats near Stevens Village, and (2) a survey flown in June 2010 covering a densely spaced block of roughly 300 km 2 near Fort Yukon and widely spaced reconnaissance lines totaling an additional 900 line km in length (figure 1). Data were collected at a nominal height of 30 m above ground surface with the RESOLVE frequency-domain system (https://cgg. com, previously https://fugroairborne.com/), using six coil pairs that measure signals at different frequencies (400 Hz, 1800 Hz, 3300 Hz, 8200 Hz, 40 kHz, and 140 kHz), with raw data output (∼ 3 m spacing) averaged into soundings output every 15 m. The 3300 Hz coil pair is in a horizontal magnetic dipole oriented in the flight direction with ∼9 m spacing; other coil pairs are all vertical magnetic dipoles with ∼7 m spacing. The AEM average depth of investigation is roughly 100 m for both the 2010 and 2016 surveys. Raw in-phase (IP) and quadrature (Q) data expressed in parts per million (i.e. ppm=10 6 * secondary magnetic field/primary magnetic field) from the airborne surveys were imported into the Aarhus Workbench software (Auken et al 2014). Upon import, all negative values were removed from the raw data (in ppm), and each frequency was averaged with a depth-dependent running mean used to populate each 15 m sounding: (1) 100 m for 400 Hz; (2) 80 m for 1800 Hz; (3) 75 m for 3300 Hz; (4) 50 m for 8200 Hz ; (5) 30 m for 40 kHz; and (6) 15 m for 140 kHz. Absolute error tolerance for the raw data was set to 5 ppm, and relative error tolerance for data was set to 5%. Error tolerance for measured elevation data was set to 2 m.
A laterally constrained inversion of the processed data was conducted in the Aarhus Workbench software using the 1D stitched nonlinear least-squares AarhusInv algorithm (Auken et al 2014), providing values of electrical resistivity with depth beneath each sounding. A smooth inversion was run with 25 layers of fixed, but variable, thicknesses that extend 150 m below the ground surface. General inversion structure was not strongly sensitive to the number of layers, layer thicknesses, or regularization parameters that were allowed to vary in the inversion. Layer thickness increased with depth to accommodate the loss of resolution with depth for AEM measurements. Preliminary inversions were conducted with varying homogeneous prior-model values; final inversions used a homogenous half-space prior model with resistivity values of 500 ohm m. The final model parameters described above were selected because they best represented our physical understanding of the system and minimized data misfit. A more detailed description of inversion parameters, final inverted resistivity values for each layer, layer thicknesses, and the uncertainty associated with these values are provided in Minsley et al (2017).
In order to quantitatively analyze similarities of medians, distributions and variances in the raw 2010 and 2016 survey data, Wilcoxon ranksum tests, Kolmogorov-Smirnov tests, Levene's tests, and Kruskal-Wallis tests were conducted. Tests were run with raw data and inverted models; both yielded the same results.

Landsat data and water body classification
Multi-date measurements spanning 1979-2011 from Landsat Multispectral Scanner, Thematic Mapper, and Enhanced Thematic Mapper Plus (ETM+) were used to determine water body extents, procedures and decision tree model accuracies are outlined in Rover et al (2012). Landsat data for analysis were selected based on the absence of snow, cloud contamination, and aquatic vegetation sometimes present later in the growing season. All pixels were classified as water or not water using a decision-tree model developed from all Landsat spectral data including indices (i.e. normalized difference vegetation index (Rouse et al 1974)), normalized difference water index (McFeeters 1996), and principle components. Training data were selected for water and non-water pixels in each Landsat scene. After applying each model, results were visually inspected and when needed, additional training data were added. Classified data were then converted to polygons, and the areas containing clouds, cloud shadows, and snow were identified as missing observations. Lake polygons were created to eliminate bias from shrinking and coalescing lakes as described in Rover et al (2012). The area of each water body (i.e. polygon) for remaining observations was calculated and normalized as a percentage of the maximum extent the lake reached. This allowed for the direct comparison of a range of lake sizes found to occur across the landscape. The water body extents were considered a proxy for water budget and lake volume changes over the observational period (Jepsen et al 2013a). The final, normalized, lake observations were represented as a vector time series for clustering. In total, 22 Landsat scenes that included ∼6830 lakes and spanned a duration of 32 years were processed and used to estimate lake areas, calculate the distance to nearest rivers or streams (http://nhd.usgs.gov), and create the time-series of lake areas for clustering.
3.3. Clustering of lake-area time series data Clustering was preformed to identify lakes with similar lake-surface area behavior over the study period. Individual lakes were clustered based on their surfacearea time-series using a hierarchical correlation-based clustering algorithm. Only lakes with valid data (i.e. clear observations) for at least 80% of the 22 available Landsat scenes were included in the clustering analysis. This restriction eliminated ∼2600 lakes, resulting in 4230 of the original 6838 lake time series being used for analysis. Lake-area time series were represented as numeric vectors containing the normalized lake extent over time and were used to create a matrix of Pearson correlation coefficients comparing each lake-surfacearea time series to the others. The lake correlated with the largest number of other lakes above a correlation threshold of 0.8 was identified Preliminary algorithm iterations were run to gauge sensitivity to different correlation values. A correlation value of 0.8 was determined to best optimize the balance between not fitting environmental noise and capturing time-series that exhibit similar trends. All lakes correlated with the identified lake above the correlation threshold were removed from the algorithm and grouped into a cluster. This was repeated iteratively until lake cluster sizes diminished to 1% of the total available lakes. Clusters that included less than 1% of the sample population were not considered as spatially representative and thus not incorporated in further analysis. To simplify time series comparison between clusters, time-series data within each cluster were averaged to attain a mean lake surface-area time-series. All inverted AEM profiles overlapping the clustered lakes were extracted, allowing for comparison of permafrost distribution between the clusters.

Permafrost distribution in the yukon flats
Results from the 2016 AEM survey (data available in Minsley et al 2017) confirm that the western area of the Yukon Flats represents a system with a thicker gravel layer, and more spatially extensive permafrost than areas to the east and south covered by the 2010 survey (data available in Ball et al 2010). This assertion is qualitatively supported when visually comparing geophysical results between the two surveys (figure 2), as well as quantitatively by previous remote sensing analysis (Pastick et al 2014) and when the raw data or inverted models of the AEM surveys were compared statistically. Abundant permafrost in the form of frozen gravel is visible in the 2016 survey as laterally and vertically extensive resistive (i.e. dark blue) segments of the inverted cross-sections (figure 2). In the 2016 survey, these highly resistive areas appear to be 25-50 m thicker on average than in the 2010 survey. All statistical hypothesis tests conducted provided p-values below 0.01, so at the 99% confidence level, we rejected the null hypotheses of equivalent sample medians, distributions, and variances for each survey. These results support the perception that the 2016 survey is different than the 2010 survey. In this case, the 2016 dataset is more electrically resistive, and thus more ice-rich and less thawed than its 2010 analog. Comparison of crosssections from the 2010 (i.e. A-A' and C-C') and 2016 (B-B' and D-D') surveys demonstrates that the eastern portion of the Yukon Flats has thinner, less resistive, and less pervasive permafrost distribution than the western portion of the Yukon Flats, offering a spectrum of permafrost conditions for further examination.

Landsat lake area clustering results
Of the 4230 lakes included in the clustering analysis, 2200 lakes were grouped into one of 14 unique clusters. Clusters ranged in size from ∼550 to ∼40 lakes. Between some clusters, the frequency and magnitude of lake-area changes were similar, but the timing of peak lake-area changes differed. Cluster time-series similarities included: (1) flashy time series characterized by episodic spikes in lake area; (2) stable time series that do not experience abrupt lake area changes; (3) perennial lakes characterized by time series that have no dry observations (i.e. time series that never reach 0% lake area); and (4) ephemeral lakes that have dry observations in their time series. Lake clusters with surface-area dynamics governed by similar physical characteristics (e.g. subsurface connectivity, surface-water connectivity, proximity to rivers and streams, maximum surface area-extent) manifest similar lake-area time-series patterns (e.g. flashy, stable), with differences in maximum and minimum lake-area timing likely caused by different clusters' spatial distribution. For example, two small ephemeral lakes within the Yukon Flats may both have a water budget primarily driven by periodic flooding, resulting in flashy time series with multiple observations where there the lake was dry. However, if the two lakes are located near different streams, local conditions may cause floods to occur at different times, creating similar time-series behavior (i.e. flashy with some dry observations) but different peak lake-area timing.
Each of the 14 clusters were placed into 1 of 5 groups (table 1, figures 3, 4, 5) to capture similarities in lake permanence (i.e. ephemeral versus perennial), subsurface permafrost distribution, and inferred connectivity to adjacent hydrologic systems (i.e. surface, subsurface, or without connectivity to either). The resulting groups represent a continuum between hydrologically connected versus non-connected and ephemeral versus perennial lakes. The following data were used to identify lake characteristics and mechanisms responsible for lake-area time-series patterns: (1) inverted AEM data over the clustered lakes (figures 2 and 4); (2) the median size of lakes in each cluster as estimated from Landsat ( figure 4); (3) the median distance of a lake cluster from the nearest stream as measured from Landsat and NHD (figure 4); and (4) lake surface-area time-series behavior as calculated from Landsat (figure 5).

Cluster grouping
Each of the five groups predominantly clustered in unique landscape positions (figure 3). Clusters with inferred surface-water connectivity were located near rivers and streams. Specifically, perennial lakes with surface-water connectivity were predominantly located near the Yukon River (i.e. SurC-Surface connectivity), while ephemeral lakes with surfacewater connectivity clustered closer to smaller rivers and streams and farther from the Yukon River (i.e. SFT-Shallow, flashy, thawed; and SFF-Shallow, flashy, frozen). Ephemeral lakes with inferred shallow subsurface connectivity (i.e. SubC-Subsurface connectivity) clustered farther from rivers and streams and south of the Yukon River, while ephemeral lakes inferred to have no connectivity (i.e. NoC-No connectivity) are located exclusively in the 2016 survey extent characterized by more continuous permafrost.
Of the 14 original clusters identified, clusters with similar time-series behavior but different peak surface-area timing tended to plot together when maximum lake-area extent and the distance of each lake polygon from the closest river or stream (i.e. a proxy for surface connectivity) were compared (figure 4). Lakes with flashy time series (i.e. SurC, SFT, and SFF groups) showed closest proximity to rivers and streams. Lakes with time series indicating perennial water (i.e. NoC, SubC, SurC groups) had the largest maximum observed lake-area size when compared with lakes whose time series indicated ephemeral behavior (i.e. SFT, SFF). The SubC group displayed the greatest likelihood of subsurface connectivity as inferred from low resistivity values, while the NoC group, hypothesized to have no surface or subsurface connectivity, showed the most consistently high resistivity values with depth, consistent with their location in continuous permafrost areas as identified via AEM. Both the SFT and SurC groups had AEM-associated resistivity values higher than the SubC group (i.e. least resistive end-member), but lower than the SFF and NoC groups (i.e. highly resistive end-members).
Ephemeral clusters' time series were identified by an abundance of dry observations, indicating no detectable water in the lakes at the time of the observation, with occasional large spikes in area. Clusters identified as ephemeral were predominantly located adjacent to rivers or streams. Clusters with ephemeral lakes were placed into two mechanistic groups: (1) small, flashy, and thawed lakes (SFT); and (2) small, flashy, and frozen lakes (SFF). The main difference between the SFT and SFF groups is subsurface permafrost distribution as inferred by AEM data, where lakes within the SFT appeared highly thawed, and the SFF group appeared predominantly frozen. Given small lake sizes and long periods without thermal energy inputs from perennial water features, small ephemeral lakes were expected to have permafrost below them. Deviation from this expectation (i.e. permafrost absence beneath small ephemeral lakes) most likely represents a thermal legacy effect from drained lakes (previously larger in size than present) or river channel migrations (Minsley et al 2012).
Clusters of permanent lakes adjacent to rivers and streams exhibit flashy time-series behavior but have no dry observations. These lakes show larger median maximum lake surface-area observations when compared to the ephemeral SFT and SFF mechanistic groups ( figure 4). AEM data co-located with these clusters are characterized by lower resistivity values than ephemeral lakes with similar time-series behavior Table 1. Characteristics of the 5 mechanistic lake groups derived from the 14 original lake clusters. The maximum lake extent and river proximity rank metrics were defined by the median value for all lake polygons contained within the cluster.

Original clusters included
Maximum lake extent (Rank 1-5, largesmall) Distance from nearest river (Rank Clusters with the above characteristics comprise a mechanistic group (SurC) of permanent lakes, which likely have active hydrologic connection to adjacent rivers or streams. In addition to the SurC group, two other types of perennial lake clusters were identified (i.e. SubC and NoC groups). Clusters placed into the SubC and NoC groups are characterized by large median maximum lake surface area observations, greater distances from rivers or streams than other clusters, and relatively stable time series. Time series for clusters in both groups demonstrate stability or gradual change with occasional peaks that are muted compared with time series associated with lakes connected to rivers and streams. AEM data differentiates the SubC and NoC groups, showing a stark contrast in their local permafrost distribution as inferred from inverted resistivity values. AEM data show the more frozen NoC group located entirely in areas of continuous permafrost (i.e. 2016 survey). The NoC cluster likely has no surface or subsurface connectivity deduced from the cluster's distance to streams, spatial location exclusively in thick, continuous permafrost, stable time-series behavior, and co-located AEM data, suggesting each of these lakes to be a closed talik. Clusters placed in the SubC group are characterized as having a high potential for shallow subsurface connectivity and were located exclusively south of the Yukon River. This region of the Yukon Flats coincides with the topographic transition from upland plateau to lowland that supplies a hydraulic gradient to drive shallow subsurface lateral flow through thawed pathways and drive upward flow through open lake taliks (Walvoord et al 2012). The SubC group exhibits a long-term decreasing trend in area, while the NoC group shows a relatively stable time series, oscillating about a relative baseline likely in response to variability in precipitation and evaporation.

Discussion
Mechanistic groups identified by this analysis represent a continuum from hydrologically connected to not connected, thawed to frozen, and perennial versus ephemeral. Two major attributes divided the mechanistic groups, (1) lake permanence and (2) hydrologic connectivity. Clustered lakes were ∼37% ephemeral Figure 4. Mechanistic groups identified by color (SurC, SubC, NoC, SFT, and SFF; table 1) in the cross-plot with distance from the nearest river or stream versus the maximum extent of the lake polygon in log-log space. 1D-resistivity profiles from individual lakes within each mechanistic group, and a normalized histograms of resistivity values for the group with associated cumulative distribution function (CDF) are plotted adjacent to the cross-plot. The histograms are normalized so that their areas sum to one, and the CDF shows the cumulative probability of a particular group having a specific resistivity value. The dashed line represents a range of values outside of which the material is likely thawed or frozen. Values to the left of the vertical red line (resistivity of 300 ohm m) are interpreted to most likely be thawed, while values to the right of the line were interpreted to most likely be frozen. and ∼63% perennial surface-water features. Among the clustered lakes ∼66% were assigned to a mechanistic group driven primarily by surface-water connectivity (i.e. SurC, SFF, SFT), ∼19% fall into a group driven primarily by shallow subsurface connectivity (i.e. SubC), and ∼15% have no connectivity (i.e. NoC). Groups with inferred surface-water connectivity clustered closest to rivers and streams supporting the hypothesis that flashy time-series behavior is likely the result of periodic input from ice-jam and open-water flooding. In this region, ice-jam flooding during spring break up of large rivers contributes to occurrences of more intense (higher stage) flooding than open-water flooding (USACE 1992, Jepsen et al 2016. Conversely, the SubC group's muted response to hydrologic inputs could be a result of enhanced shallow subsurface connectivity supporting regular inflow and outflow, as determined from the close spatial grouping of the clustered lakes. A highly thawed subsurface and the presence of open taliks inferred from AEM, stable time-series behavior, and large topographic gradient in areas where these clusters occurred further support the SubC group's potential for subsurface connectivity. We have developed a mechanistic description of the differences in these lake systems (figure 5), considering two different hypothetical times in the life of a lake. During time 1, ephemeral lakes have no water present, while time 2 represents post-flood conditions. For the perennial SurC group, time 1 represents a post-flood period, while time 2 illustrates a period characterized by recent or active flooding. Shallow connectivity exists during both time 1 and time 2 for the perennial SubC group, resulting in more stable fluxes into and out of the lake, reflected in the group's more stable time series. The time series for the perennial NoC group shows lake-size stability relative to other mechanistic groups; it is likely that fluctuations in size are primarily due to precipitation and evapotranspiration in the absence of surface and subsurface hydrologic connectivity.
Changes in stream discharge and ice-jam flooding have the potential to influence the greatest number of lakes (∼66%) within the Yukon Flats because surfacewater features are likely their primary hydrologic input other than precipitation (i.e. SurC, SFT, and SFF groups). Climate warming will affect the timing of river-ice breakup, potentially trending towards early breakup timing and less severe break-up events (e.g. Smith 2000, Shiklomanov andLammers 2014); however, other factors also contribute to regional ice-jam flooding propensity (Cooley and Pavelsky 2016). The effects of increasing temperatures on changes in winter snow accumulation, solar radiation-induced weakening of ice strength, Figure 5. A conceptual diagram of mechanistic lake groups. Two times from each of the time series highlight processes governing the hydrologic mass balance of each group. Mean time series for each of the 14 clusters placed into the mechanistic groups above are depicted as a separate time series, elucidating similar responses to fluxes in and out of similar clusters. Solid and dashed lines represent the mean time series for each of the unique clusters put into a particular mechanistic group. changes in early season freeze-up regimes, and increasing frequency of mid-winter break-ups complicate future ice-jam flooding predictions (Prowse and Beltaos 2002). Due to the multitude of variables controlling ice-jam flooding, in the trajectories of ice-jam flooding frequency and intensity remain uncertain. An increase in flooding events may cause perennial lakes with surface connectivity (i.e. SurC) to grow in size, and ephemeral lakes with surface connectivity (i.e. SFF, and SFT) to transition towards perennial surface-water features (i.e. SurC), unless offset by increasing lake evaporation. Consequently, an increase in ice-jam flooding propensity may cause a net increase in surface-water area in lakes with surface-water connectivity. If ice-jam flooding propensity decreases, net lake-surface area will likely decrease in lakes connected to rivers and streams, as evaporation increases with warmer temperatures and the lakes' primary hydrologic influx decreases. Constraining the future propensity for ice-jam flooding will be required to predict surface-water distribution with greater certainty.
Lakes with minimal or no surface-water connectivity (∼34%) may be the most sensitive to climateinduced permafrost degradation due to the greater influence of subsurface hydrologic connectivity on the water budgets of these lakes. The NoC group exists exclusively in areas of more continuous permafrost relative to other areas observed via AEM in the Yukon Flats. In contrast, the SubC group is located in discontinuous permafrost areas that appear to be the most thawed compared to other areas observed via AEM within the Yukon Flats. Assuming the NoC to SubC mechanistic lake group spatial pattern reflects a permafrost thaw continuum that is transferrable to a temporal thaw analog, we can infer that a shift toward increasing thawed conditions will lead to a reduction in the percentage of lakes being driven primarily by direct hydrologic input and outputs (i.e., precipitation-evaporation) and a concomitant increase in the percentage of lakes influenced by lateral hydrologic inputs and outputs through shallow subsurface connective pathways. Based on the broadly observed decline in lake area of SubC lakes, we may further infer that continued permafrost thaw in the Yukon Flats will result in a net decrease in lake surface area for lakes not connected to rivers or streams.
The AEM data indicate thick, continuous permafrost (i.e. frozen gravel) in the western portion of the Yukon Flats, and the clustering results for lake-area time series in this area showed the greatest degree of homogenous lake behavior (synchronous time series and same mechanistic groups over large area), with less diversity in processes governing lake surface area. This homogeneity may be a result of extensive permafrost. Consequently, as thaw progresses, driving hydrologic processes may change in the western portion of the Yukon Flats such that lake-area dynamics may exhibit more heterogenous patterns, similar to the eastern portion of the Yukon Flats presently, where hydrologic pathways through the shallow gravel layer connecting inland water bodies are non-uniformly distributed. As the shallow gravel layer thaws in the central and eastern Yukon Flats, we can expect shallow subsurface hydrologic connectivity to become increasingly more uniform, thus promoting a shift back toward homogenous lake-area dynamics across the Yukon Flats as lake levels reflect subtle fluctuations in the regional water table. Homogeneity in lake behavior may reflect homogeneity in the ground, whether the permafrost is continuous, or mostly absent.
Due to the logistical challenge of collecting field data in such a remote region, all the measurements used in this study were either a satellite remote sensing product or collected via helicopter (i.e. AEM). Some Landsat scenes were acquired at different times during the same year. Image acquisition schedules (i.e. every 8 to 16 d)and observation frequency due to clouds, cloud shadows, and snow limited the opportunity for attaining a more complete inventory of the lakes. A selection of high-quality images were used to provide the clustering algorithm with a multi-decade record, although intra-annual lake surface area variations could influence results. Additionally, observed lakes varied in size, meaning lake responses to hydrologic mechanisms could be influenced by lake volume (lake surface-area was used as a proxy) and bathymetry undetectable with AEM or Landsat images. Collected AEM data lack the near-surface resolution to image shallow flow paths directly, despite being the best option for mapping regional permafrost distribution. Consequently, synoptic ground-based, high resolution shallow monitoring of thawed regions within the SubC clusters could provide useful information to further evaluate our hypothesis on evolving shallow subsurface flow paths. Landsat images were spatially but not temporally-extended beyond the processed data set presented in Rover et al (2012). As such, this study provides a basis for future testing of proposed Yukon Flats lake-area trajectories using temporally-extended Landsat data, repeat AEM surveys, and additional information on lake bathymetry as remote sensing technologies continue to improve.

Conclusions
This study presents a synthesis of new data products and analyses from the Yukon Flats of Alaska including: (1) permafrost mapping via AEM; (2) clustering of Landsat derived lake-area time series and lake landscape position; and (3) exploration of hydrologic mechanisms and characteristics that influence lake-area dynamics. This is the first study to simultaneously constrain the effects of deep permafrost distribution via AEM and the heterogeneity of lake dynamics based on lake surfacearea time series and landscape position. Comparing regions with permafrost conditions along a continuous-discontinuous permafrost distribution spectrum provides a space-for-time proxy in which the mechanisms locally controlling lake-surface area dynamics may be transitioning towards those observed in areas of more pervasive thaw. Mechanisms and characteristics driving lake surface areas were elucidated by constraining permafrost distribution via AEM data, observing lake time-series behavior with Landsat, and clustering lakes via time-series correlation and landscape position.
Based on the findings of this study, permafrost distribution has an influential effect on regulating lake surface-area dynamics. Specifically, areas with continuous permafrost tend to exhibit relatively synchronous lake-area dynamics, whereas areas of discontinuous permafrost tend to reflect more asynchronous lake-area behavior among lakes in close proximity. This regional pattern may be explained by the evolution of subsurface connectivity along a thaw spectrum. In continuous permafrost, general lack of subsurface connectivity among lakes results in lake water budgets that are controlled primarily by hydroclimatic factors (and river discharge conditions for those lakes connected to active flow systems). In discontinuous permafrost, differential thawing leads to variability in subsurface connectivity and thus variability among lakes in the influence of groundwater fluxes through thawed pathways of the shallow gravel layer. With continued permafrost degradation in the Yukon Flats resulting in a ubiquitous (versus differential) enhancement in subsurface connectivity between water bodies, we expect lake-area fluctuations to become increasingly synchronous. We also expect lakearea fluctuations to become more muted as the expanding thawed subsurface acts as a low-pass filter for hydroclimatic variability. Further refinement in time and space of these general projections require improved shallow subsurface characterization throughout the Yukon Flats, more certainty in future predictions of icejam flooding events, and integrated (climate-surface water-subsurface) hydrologic modeling.