Surface water area in a changing climate: Differential responses of Alaska’s subarctic lakes

Lake surface area in arctic and sub-arctic Alaska is changing in response to permafrost deterioration, changes in precipitation, and shifts in landscape hydrology. In interior Alaska, the National Park Service’s Central Alaska Network Shallow Lakes program studies lakes and ponds in a wide range of geomorphological settings ranging from alpine lakes to low lying lakes on fluvial plains. The purpose of this study was to determine if and how lake area was changing across this diverse environment. Using the USGS Dynamic Surface Water Extent product, we tested landscape-scale trends in surface water area from 2000–2019 in 32 distinct ecological areas, or ecological subsections, within the three parks. Surface water area declined in 9 subsections, largely in glaciated landscapes with coarse substrates and areas underlain by ice-rich permafrost. Surface water increase was seen in one subsection in the floodplain of the Copper River in Wrangell-St. Elias National Park. No net change was observed in many subsections, but individual lake analysis showed that within several ecological subsections some lakes were increasing in area while others decreased in area, masking changes in lake surface area within the subsection. Over the course of the study period, surface water area in all parks experienced similar fluctuations, likely due to oscillations in regional climate. Periods of high surface water area coincided with relatively warm, wet periods. Climate change models project increases in both temperature and precipitation in Alaska; our results suggest periods of regional wetting may mask longer-term declines in surface water area in some geomorphological settings. Overall, lake surface area declined over the study period; declines were greatest in the Glaciated Lowlands in Denali National Park and Preserve.


Introduction
The hydrologic environment in subarctic Alaska is undergoing rapid change [1], due in part to changes in precipitation patterns [2], groundwater flowpaths [3] and catchment characteristics [4]. A growing number of records indicate that lake surface area has declined over the past 50-100 years in the Subarctic [5][6][7] and low latitude regions of the Arctic [8]. Changes in lake surface area are attributed to a host of climate-related drivers including permafrost thaw [9], thermokarst [10], and shifts in precipitation regimes [11] as well as terrestrialization, or the infilling of lakes to peatland [7]. Lake surface area at high latitudes is highly dynamic and has historically been regulated by overwinter snowpack [12]. Overwintering conditions are changing in Alaska with warmer winter temperatures [13], a shorter snow season [14], increased seasonal thaw depth [15,16], and in some areas increased precipitation [2]. However, not all lakes appear to be declining in surface area; for example, lakes in the Arctic found in moraine settings currently appear to be stable, with lake levels oscillating with precipitation [8]. In some regions, where glaciers are actively retreating, lakes are expanding [17,18]. These differential responses to climate change, combined with the large acreage of subarctic shallow lake habitat (550,000 km 2 and >57,000 lakes), make understanding lake succession in interior Alaska an important but challenging task.
Interior Alaska lakes are concentrated in three major geomorphological settings: glacial moraines, thermokarst plains and floodplains. Source water inputs (precipitation, groundwater, glacial meltwater, etc.) vary among these three settings, and catchment characteristics (eg. vegetation, topography and permafrost) influence the timing and duration of surface water influx. Research in the low Arctic [8] shows that in glacial deposits lake surface area is correlated to precipitation, whereas lakes in thermokarst depressions are influenced by developmental stage: initial inception, expansion, and reduction by either catastrophic drainage or by becoming connected to groundwater [19][20][21][22]. Surface area extent in riverine systems is governed by groundwater or flooding and is dependent upon stage, sediment porosity and hydraulic conductivity [23]. These lakes tend to have discontinuous permafrost and once isolated can transition to thermokarst lakes or infill with vegetation [24]. These geomorphological differences create a complex mosaic of lakes-often within close proximity to each other-that will respond differently to change.
Shallow lakes provide critical habitat for wetland obligate species and understanding longterm patterns in distribution and abundance of lakes is critical to land management agencies.
To better understand the current status of lakes and to identify important drivers of lake change we have developed a study to track changes in lake surface area across three national parks and preserves in subarctic Alaska. We identified four primary objectives of this study: 1) detect trends in lake surface area annually from 2000-2019; 2) identify geomorphological patterns of change; 3) relate surface area change patterns to environmental climate variables; and 4) evaluate fine scale changes in lake surface area. Results from this study will be used to enhance our understanding of lake surface water dynamics in the parks and in subarctic Alaska.

Study area
The study was conducted by the Central Alaska Network (CAKN) Inventory and Monitoring program which consists of three national park units located in interior Alaska: Yukon-Charley Rivers National Preserve (YUCH), Denali National Park and Preserve (DENA), and Wrangell-St. Elias National Park and Preserve (WRST, Fig 1). These subarctic parklands contain numerous lakes distributed across a large latitudinal and longitudinal gradient that extends from southeast Alaska to eastern interior Alaska, crossing a diverse range of permafrost conditions (continuous to sporadic), geologic settings, and climate regions. Surface water area trends were tested in pre-defined ecological subsections [25][26][27][28] rich in surface water, which were previously assigned based on geology, landforms, soils and vegetation.
Lakes in these areas are concentrated in glacial moraines, thermokarst plains and floodplains. The Yukon River Valley and its associated floodplains contain most of the lakes found in YUCH, making up a large portion of riverine-influenced lakes within this study. Numerous lakes are also found in the Copper and Bremner River basins in WRST, and scattered along rivers throughout DENA. Glacial processes formed most lakes within DENA and WRST including large glacial scour lakes such as Tanada, Copper and Wonder Lakes, as well as small lakes distributed throughout the glacial moraines. Permafrost is present in all three parks, and thermokarst lakes are common in YUCH and DENA.
The parks span a range of localized interior climates including the Northeast Interior, Central Interior, Southeast Interior, Cook Inlet, and Northeast Gulf climate regions (Fig 1) [29]. These climates vary from extreme annual temperature swings and low precipitation in the Interior to milder temperatures and high precipitation in the southern climate zones near the coast. However, the climate in these parks is changing. Mean annual temperature (MAT) increased significantly at long-term National Weather Service stations near YUCH and DENA (p<0.01, Mann-Kendall trend test) between 1980-2019, and the growing season (June-September) temperature has increased in WRST (p = 0.04) although MAT has not changed significantly (S1 Fig). However, 2014-2019 was significantly warmer than the long-term average [16]. No significant trends in precipitation have been observed between 1980-2019; however, potential evapotranspiration rates increased in all parks (p�0.01, Mann-Kendall trend test, S2

Mapping surface water extent
We used the USGS Landsat Dynamic Surface Water Extent or DSWE product [31] to study surface water change in YUCH, DENA, and WRST from 2000-2019. The DSWE was developed as a tool to identify surface water using Landsat 4-8 data. Briefly, Landsat -detected surface reflectance and elevation data are used in a decision-rule style model to classify pixels. Water was mapped at 30 m resolution and the layer is classified by pixel to five levels of inundation, coded by categorical raster values: Not Water (raster pixel value 0), Water-High Confidence (raster pixel value 1), Water-Moderate Confidence (raster pixel value 2), Potential Wetland (raster pixel value 3), and Water or Wetland-Low Confidence (raster pixel value 4) [31]. DSWE tiles representing the growing season months (June-September) with less than 40% cloud cover and cloud shadow were downloaded from the US Geological Survey Earth Explorer (https://earthexplorer.usgs.gov/) for all three parks 1980-2019. Tiles were visually compared to satellite and aerial imagery (including Alaska SDMI SPOT 5 imagery, Alaska Worldview imagery, and State of Alaska Open Data Geoportal services, https://geoportal. alaska.gov/arcgis/rest/services/ahri_2020_rgb_cache/MapServer) to assess the accuracy of DSWE pixel classifications. "Water-High Confidence" and "Water-Moderate Confidence" accurately identified surface water in the study area. However, the other classified categories ("Potential Wetland", "Water or Wetland-Low Confidence") proved unreliable in our study ecosystems. This was likely influenced by vegetation type-specifically conifers-which often classified as false positives in these categories as described in Jones 2019. We excluded these classifications from our statistical analyses. The tiles were reclassified into three categories by pixel value: high confidence water (used in analyses), wetland, and low confidence water (Table 1), using the Reclassify tool (Spatial Analyst) in ArcMap 10.7.
Reclassified rasters were summarized by year to provide annual estimates of surface area extent for each category. Each pixel was summarized as the mode (most common value) in a three-year moving window to create a composite image that was used to calculate annual surface water extent (Fig 2).
We used this method because the annual surface water maps rendered from the DSWE periodically contain terrain shadows, cloud shadows (etc.), and Landsat mechanical failures [32] that affect the quality and accuracy of surface water estimates. This method was adopted from Swanson (2019), who demonstrated the effectiveness of this technique in arctic parks. Very little satellite coverage exists for the 1980s and 1990s (not included in our analysis but used for comparison purposes only in Table 2), so images were summarized by the decade (mode taken from all tiles within the 1980s, all within the 1990s). Summarized tiles were merged to create a complete surface map for each park per year. Evaluation of water surface area change over time by ecological subsection was accomplished by extracting water pixel value counts with the Zonal Histograms tool in ArcMap 10.6.1. Surface water area trends were tested by subsection, using Mann-Kendall trend tests with Theil-Sen slope in R open source software, using Kendall and TTAinterfaceTrendAnalysis packages [33][34][35]. Environmental variables were summarized by a 3-year sliding average to match the averaged surface water area data. Stepwise linear model selection by AIC evaluation was used to inform which environmental variables related most to surface water area, using the MASS package in R. The best model selected was compared to a mixed model with the same variables with subsection added as a random intercept using a likelihood ratio test. After determining that the mixed model better fit the data, we elaborated on individual environmental variable effects using a mixed effects hypothesis testing approach. Ecological subsection was used as a random effect (random intercept) to determine study area-wide trends with environmental variables as fixed effects. The lmer function within the lme4 package was used to test linear mixed effects fit by Restricted Maximum Likelihood (REML). Fixed effects were evaluated via type II analysis of variance, reporting F-tests using Kenward-Rogers approximation for denominator degrees of freedom using the pbkrtest package in R (https://CRAN.R-project.org/package=pbkrtest). P values were reported for ease of interpretation, derived from the F-statistic, with significance reported at α = 0.05.General linear models were used to explore correlations between surface water area and environmental variables per each subsection. Packages used include lme4 [36] and corrplot [37].

Mapping ecological subsections
The change in DSWE-classified water area was summarized by ecological subsection within the parks. Previous work delineated ecological subsections within each park based on soils, vegetation, geology, and other factors [26]. Briefly, ecological areas were assigned in ArcView 3.2a using land cover maps, satellite imagery, aerial photography, geologic and soil maps, and US Geological Survey topographic maps and delineated using principles presented by Bailey (1996) [38] and Wertz and Arnold (1972) [39]. Where ecological subsections were found to be too coarse, a finer ecological division, or "detailed subsection" was used. For example, the entire Yukon River Valley was classified as one ecological subsection but contained six ecologically distinct areas (Yukon River Active Floodplain; Wet Terraces with Oxbows; Wet Terraces with Few Ponds; Wet Terraces with Thermokarst Lakes; High Terraces, Undulating; Nation/ Kandik/Bonanza Valleys). In this instance, we used the detailed subsection instead of the larger subsection delineation. Subsections with >0.1% surface water were included in the analysis (Fig 1, S5-S7 Figs). Prior to analysis, the ecological subsections maps were manually edited to exclude rivers, streams, and glaciers, resulting in the removal of approximately 4.8% of the area to be analyzed within the three parks.

Individual lake surface area trends
In addition to ecosystem-wide surface water area change, we analyzed surface water trends in 365 individual lakes that are monitored as part of the CAKN long term monitoring program. Lakes were selected for monitoring in the early 2000s using a general randomized tessellation stratified (GRTS) sampling approach [40] to select a spatially balanced sample of lakes [41]. Shoreline edges were delineated using ArcMap editing tools for each of the lakes in the mid-2000s using the USGS National Hydrography Dataset, IKONOS, and high-resolution aerial photography. Using ArcMap 10.7, each lake shape was buffered by 100 m to account for lake area increase. The Zonal Histogram tool was used on 2000-2019 DSWE layers to extract the number of pixel values classified as water for each year per lake. Surface water area trends were tested by individual study lake, using Mann-Kendall trend tests with Theil-Sen slope in R open source software using Kendall and TTAinterfaceTrendAnalysis packages [33][34][35]. Each study lake was visually inspected using the DSWE layers in GIS. Lakes that were too small to be recognized by the DSWE were identified and removed, resulting in the omission of 21 study lakes from further analysis totaling 344 study lake results reported.

Surface water trends in ecological subsections
Over the course of the study period (2000-2019) total surface water area remained stable in 22 of 32 ecological subsections and detailed subsections, although individual study lakes within some of the subsections did experience increases or decreases in surface water area (Table 2). Lake surface area declined in 9 of 32 ecological subsections and increased in one ecological subsection. On a park-wide basis, DENA lost approximately 33 km 2 of surface water between 2000 and 2019, WRST lost 4 km 2 , and YUCH lost 0.05 km 2 . Together, the subsections significantly losing surface water account for approximately 3,905 km 2 of parkland, and together contain roughly 36 km 2 of surface water in lakes and ponds. Lake surface area declines occurred most frequently in glaciated regions of the parks where coarse substrate is prevalent; five of the nine subsections undergoing declines were in glacierinfluenced regions of DENA and WRST (Table 3, Figs 3 and 4). Of these, the Glaciated Lowlands of DENA contains the greatest number of lakes (approximately 1,251) and occupies 861 km 2 in land area. Surface water decreases in this subsection could have the greatest impact on habitat loss in this study, along with the Northern Chugach Foothills in WRST, which contains 131 lakes and occupies a land area of 1,350 km 2 .
Declines in lake surface area also occurred in a subset of permafrost-affected landscapes; three of the 9 subsections where declines occurred are in areas that contain ice-rich permafrost or thermokarst influence in YUCH and WRST. In YUCH, these were Wet Terraces with Thermokarst Lakes and High Terraces, Undulating-the latter of which contains areas of yedoma, an ice-rich form of permafrost that is laden with carbon ( Fig 5). Carden Lakes in WRST is one of the few subsections in this study that contains lakes with volcanic ice-rich basins [41]; these lakes are significantly losing water including Carden Lake, the largest lake in the subsection. The final subsection losing water, Little Black River Hills in YUCH, burned in wildfires within the past 40 years.
Lake surface area increased in one ecological subsection in WRST-the Middle Copper River Floodplain and Terraces, accounting for 89 km 2 in land area and approximately 3.6 km 2 in surface water area. This subsection is partially within the active floodplain of the highly dynamic Copper River on a terrace characterized by finer substrates with thicker organic surfaces.

Surface water trends in long-term study lakes
Overall, 96 of the 344 individual study lakes showed a significant trend over the study period; 77 decreased in surface area while 19 increased in area. Of the 77, approximately 8 lakes fell within subsections associated with riverine processes, 52 fell within subsections associated with permafrost, and 28 fell within a glacially-influenced area. Although no trend was detected in many ecological subsections, individual study lakes within some of the subsections did experience significant change in lake surface area. Many lakes in the Eolian Lowlands and Glaciated Uplands in DENA, for example, significantly lost surface water area over the study period even though no significant change was detected at the subsection level (Table 2). Numerous lakes within the Wet Terraces with Oxbows in the Yukon River Valley, and the Chitina Valley Moraines and Hills subsection in WRST, experienced change-some increasing and some decreasing in surface area. Subsections which contain both increasing and decreasing lakes, leading to a net "no change" in surface water area within a subsection, experienced what we hereafter refer to as "masking" of change.

Environmental variables
All three parks experienced similar patterns in lake surface area throughout the study: lake surface area peaked in the early 2000's and declined during the mid-2000's, followed by an increase again in 2015. When compared to regional climate variables, surface water highs broadly appear to occur during warm, wet periods with gradual declines occurring between (Fig 6). The three best fit linear models (lowest AIC values, ΔAIC<10) included mean annual temperature, mean annual precipitation, potential evaporation, water balance, and maximum annual snowpack as predictive environmental variables (adj R 2 = 0.26). Individual environmental variable effects are reported here, with subsection as a mixed effect to account for geographic differences in climate effects. (Table 4). Although lake surface area across the parks fluctuated in a similar pattern, the subsections experiencing significant declines departed from the pattern around 2015. When the second warm, wet period began around 2014, the surface water in these subsections continued to decline. Explanatory environmental variables varied among parks and ecological subsections. Mean annual temperature, average summer temperature, and PET correlated with surface water area in DENA subsections; precipitation, water balance (precipitation-PET), PET, and average summer temperature correlated with surface water area in WRST subsections; and PET, Yukon River discharge, precipitation, snowpack, and average minimum winter temperature correlated with surface water area in YUCH subsections (S8-S10 Figs).

Discussion
For central Alaska parklands, long-term regional climate trends may be more important than short term variations in precipitation and temperature in determining lake surface area extent. Between 2000 and 2019, we observed large scale oscillations in lake surface area throughout the entire study region. Surface water area increased during warm and wet periods, peaking at the end of the warm period, and decreased during cool periods, regardless of precipitation patterns (Fig 6). Surface water dynamics at the subsection-level, however, appear more dependent  Table 4. Results of mixed effects modeling investigating individual environmental variable effects on surface water area in Central Alaska Network parks. Climate variables are averaged over a 3-year moving window to match the DSWE data. Subsection was used as a random effect.

PLOS CLIMATE
upon complex interactions among climate variables, local groundwater, glacial meltwater supplies, river discharge, etc. Several subsections within the parks undergoing surface water declines departed from this pattern in 2014 and continued to decline, despite the warming and wetting trend occurring on the landscape. In other areas, some individual lakes lost surface water while others gained surface water, likely resulting in a masking of change in some subsections. Overall, decreases in surface water area exceeded gains in both the ecological subsections and individual study lakes, corroborating recent findings of overall lake area loss in the Subarctic and Arctic over a similar time period [42]. Surface water area decreases occurred most often in glacier and thermokarst-affected landscapes.

Glacial lake succession
Significant changes in lake surface area occurred in 10 of the 32 ecological subsections. Changes were most apparent in glaciated landscapes of DENA and WRST. In agreement, Riordan et al. (2006) [43] found surface area declines in interior Alaska including DENA and WRST between the 1950s and 2002. Arctic and subarctic landscapes, including DENA, are greening, with a documented increase in shrub and tree growth in glaciated regions [44][45][46][47][48]. Encroachment of shrubs and trees is a natural part of the succession of glaciated landscapes and is accelerated by warming temperatures. Trees and shrubs can intercept 15-30% of water moving across the landscape, increasing transpiration and biological water demand within watersheds and decreasing the volume ultimately received by water bodies [49]. Lake surface area declines in glaciated regions may be influenced by greening that is occurring here. Additionally, numerous historic lake records demonstrate that warming coincides with an intensification of lake terrestrialization and infilling rates in glacially-formed lakes [50][51][52][53]. This could be happening in the glaciated lowlands in DENA where the last glacial maximum was 26,500-19,000 years ago [54], as well as in subsections not influenced by glacial activity.

Thermokarst lake succession
Significant declines in surface water area also occurred in three subsections containing thermokarst lakes or ice-rich permafrost. Wetting conditions in other regions have been shown to degrade permafrost in the long-term [15], and warm periods are directly responsible for the degradation of permafrost and glaciers [16]. More individual study lakes in YUCH increased in surface water area than in both WRST and DENA. This may be due to the more extensive ice-rich permafrost found here; lake expansion during this period may be an indication of additional permafrost degradation. At least four other subsections with discontinuous permafrost exhibited masking of change, including a thaw lake basin along the Yukon River and the Eolian Lowlands in DENA (Fig 7). Because thermokarst activity can expand lakes and create new ponds and wet areas, thermokarst-affected areas may be the most likely to mask overall change trends as increases and decreases in surface water area create no net surface water change over a geographic area. Masking has likely caused a significant underestimate of the changes occurring in permafrost affected landscapes in this study.

Riverine lake succession
Fourteen subsections within this study, containing at least 1,770 lakes, fell in areas of riverine influence. Of these, five subsections decreased and one subsection increased in surface water area over the study period. Subarctic riverine lakes are generally surrounded by discontinuous permafrost and typically in contact with groundwater for at least part of the year, which can drive lake surface area change [7]. The Yukon River Valley in YUCH contains approximately 567 lakes, over half of which are surrounded by permafrost. As permafrost thaws, groundwater connectivity to rivers including the Yukon is increasing [55][56][57]. This increased connectivity could potentially 1) increase underground recharge to floodplain lakes or 2) drain the landscape, exporting water to the river and carrying it downstream [3,58]. These differential effects may be why masking of lake area change was seen in some riverine subsections in this study.

Other factors
Lake surface areas were stable in much of the landscape; we detected no significant long-term trend in surface water area in 22 ecological subsections. The continued widespread presence of permafrost and glaciers still support a wet landscape via the prevention of subsurface drainage and addition of meltwater, respectively. Montane snow and glaciers hydrologically feed the landscape throughout the growing season and likely have an important impact on the parks' surface water dynamics. However, glaciers throughout Alaska are receding and glacial extent has declined significantly over at least the last 60 years [59]. Lack of understanding of these flow paths impair our ability to determine if glacial meltwater is compensating for possible long-term landscape drying trends. As climate models project Alaska to become warmer and wetter, it is possible that we could see a landscape wetting before we cross a permafrost/glacial deterioration threshold and landscape-wide drying. Additionally, we did not consider the effects of wildfire history on lake surface area in this study, which has been shown to affect rates of lake surface area change in permafrost areas [60].

PLOS CLIMATE
The accuracy of the DSWE product could influence surface water estimates; for example, we noticed that dense floating vegetation (e.g. pond lilies) obscured the detection of surface water in a few instances. However, the DSWE is available for the entire state which is difficult to attain in such a large landmass, and provides a promising resource for surface water change. The influence of any of these factors on lake area change in the Subarctic should continue to be pursued in future study.

Conclusions
In Alaska's interior parks, more ecological areas and individual lakes are losing surface water area than are gaining. Lakes found to be most susceptible to surface water area decreases appear to be related to the evolution of glacial and permafrost landscapes, although masking of overall change likely occurred in thermokarst-affected areas. Since 2014, mean annual air and ground temperature trends have dramatically increased, which will lead to significant permafrost degradation in Alaska's national parks [16,61]. Precipitation rates are also expected to increase in most of Alaska [2], which may accelerate permafrost thaw [15]. In this study, it appears that the interaction of these two climate variables (warming, wetting) could lead in the short-term to a wetter landscape overall as lakes endure a lag effect before expected surface water loss. As the climate continues to warm, enhanced groundwater connectivity could lead to both the stabilization or growth of some lakes and the drainage or infilling of others. The thaw of areas with coarse substrates (gravel, sand) are more likely to lose surface water area or perhaps transition to peatland.
Whether situated in glacial, thermokarst, or river-affected landscapes, shallow lake ecosystems provide critical habitat for wildlife, including many species of migratory birds. The Park Service also values these habitats as subsistence, cultural, recreational, and aesthetic resources within its jurisdiction. The greatest surface water losses appear to be occurring in DENA, which is the most well-known and visited parkland in Alaska. Loss of aquatic habitat could affect important breeding populations of the many species that use shallow lakes and parkland users alike.
As landscape infrastructures of the Subarctic-permafrost and glaciers-degrade with a changing climate, differential responses of surface water area to warming temperatures and increased precipitation will make surface water area predictions difficult. Modelling efforts to predict landscape change tend to generalize ecological shifts, and it is important to continue monitoring change in variable subarctic settings in a quantifiable manner. Landscape and water area change is complex and updating models and scenarios with observational data is essential for proactive management and adaptation. increased significantly between 1980-2019 near all three parks. Precipitation data obtained from the National Weather Service Cooperative Observing Network. Potential evapotranspiration was calculated using the Hargreaves-Samani method [62]. Data are presented with a loess curve for better visualization of trends (degree = 2, span = 0.5).