Because winter is still the dominant season in the Arctic and is widely observed to be changing most rapidly, we hypothesized that winter processes driving ice growth exert the dominant control on lake temperature regimes. We suggest that lake depth relative to ice thickness should create an important threshold governing the sensitivity of MABT in the winter (ice‐covered period). In the summer (open‐water period), we expect that lakes will respond similarly to variation in air temperature. Accordingly, we partitioned our analysis of lake temperature dynamics at the water‐sediment interface by seasonal surface conditions (i.e., ice covered versus open water) as originally recommended by Brewer [ 1958 ]. We developed empirical models based on MABT observed among a range of lake gradients by depth and regional climatology in northern Alaska. We then used these models to examine the temperature sensitivity of lake beds to modes of Arctic lake change. Our overarching goal is to place projected rates of lake temperature change into the context of other Arctic systems to better understand potential responses and feedbacks to warming, specifically permafrost thaw.

Air temperatures in Arctic Alaska have warmed by 2.7°C since 1979 at Barrow [ Wendler et al ., 2014 ], partly in response to sea ice decline and associated Arctic amplification [ Serreze and Barry , 2011 ]. Warming has been most pronounced in the winter [ Wendler et al ., 2014 ], and this, along with higher snowfall, has resulted in reduced ice growth on lakes and a corresponding regime shift in shallow thermokarst lakes [ Arp et al ., 2012 ]. For example, nearly 36% of the bedfast ice lakes on the Barrow Peninsula of northern Alaska transitioned to floating ice conditions between 1992 and 2011 [ Surdu et al ., 2014 ].

Shallow, seasonally ice‐covered lakes are abundant in lowland Arctic landscapes, often covering 20–40% of the land surface [ Smith et al ., 2007 ; Grosse et al ., 2013 ]. The majority of these lakes are of thermokarst origin, forming by the degradation of ice‐rich permafrost [ Grosse et al ., 2013 ]. Thermokarst lakes progressively expand into permafrost shorelines [ Hopkins , 1949 ] and transfer heat into the underlying permafrost over time [ Brewer , 1958 ; Ling and Zhang , 2003 ]. Bed temperature plays a key role in the rate of sublake warming and determines whether sublake permafrost thaws or is maintained [ Brewer , 1958 ; Lachenbruch et al ., 1962 ]. MABTs above freezing cause sublake permafrost thaw, thereby creating taliks [ Burn , 2002 ] and potentially releasing greenhouse gases [ Walter et al ., 2006 ; Wik et al ., 2016 ]. In most shallow thermokarst lakes, however, ice grows to the lake bottom early enough in the winter to keep MABTs below freezing which preserves sublake permafrost [ Ling and Zhang , 2003 ]. Because ice‐cover duration and ice regimes (i.e., bedfast versus floating ice conditions) are responding rapidly to changing Arctic climate [ Duguay et al ., 2006 ; Arp et al ., 2012 ; Surdu et al ., 2014 ], predicting how MABTs respond during ice‐covered (winter) and open‐water (summer) periods is critical to understand how sublake permafrost may thaw below shallow Arctic lakes. However, comprehensive lake data sets and analysis of controlling processes are currently lacking [ Burn , 2002 ; Ling and Zhang , 2003 ].

The rate of terrestrial surface warming is of keen interest in understanding how ecosystems will respond to climate change, particularly in the Arctic where climate warming is amplified [ Chapin et al ., 2005 ; Serreze and Francis , 2006 ; Jorgenson et al ., 2010 ]. Thermokarst lakes are prime ecosystems to look for such responses because they are developed in and interact with permafrost, variably exchange heat and moisture with the atmosphere depending on both ice cover and stratification, and cover vast areas of many Arctic regions [ Grosse et al ., 2013 ]. Lake surface temperature trends measured from satellites and in situ using lake buoys indicate a globally averaged warming rate of 0.3°C/decade, yet with considerable variation regionally and with lake morphology, ranging from cooling trends to rates exceeding atmospheric warming [ O ' Reilly et al ., 2015 ]. Three important points have emerged from recent studies of lake temperature dynamics: (1) that rates of warming tend to be greater at higher latitudes and in lakes with seasonal ice cover [ Schneider and Hook , 2010 ; O'Reilly et al ., 2015 ], (2) that shallow lakes are most sensitive to warming, while deeper lakes often have stable to declining temperatures [ Kraemer et al ., 2015 ; Winslow et al ., 2015 ], and (3) that most analyses only consider lake surface temperature during ice‐free seasons [ Boike et al ., 2015 ]. However, for thermokarst lakes developed in permafrost landscapes, mean annual bed temperature (MABT) is the more important thermal condition [ Lachenbruch et al ., 1962 ; Burn , 2002 ] because this value predicts occurrence of perennial sublake thaw (i.e., talik formation).

The model developed from these data was used to evaluate the sensitivity of MABT to three major controlling variables: (1) summer air temperature during open‐water periods ( T a‐open ), (2) ice‐cover duration (ICD), and (3) MIT expressed as EDR in equation 1 . The sensitivity analysis was conducted using Barrow climate data (NWS 700260 27502) from 1969 to 2015 in order to provide realistic ranges of variation. Trends in MABT records were evaluated using simple linear regression and compared to mean annual air temperature (MAAT) at 3 m height and mean annual ground temperature (MAGT) at 1 m depth. MAGT records are a combination of Geophysical Institute Permafrost Laboratory Model (GIPL) [ Marchenko et al ., 2008 ] simulations (1978–2001) and station observations (2002–2015) as described in Romanovsky et al . [ 2015 ].

The measurements used for model development consisted of MABT data from 29 lakes over a 3 year, winter‐centric period from 2013 to 2015, yielding a data set of= 45. Not all lakes had continuous records during the period of interest for a variety of reasons, primarily sensor malfunction or ice drift that caused sensor loss or movement from the lake center. An additional 10% of MABT observations were randomly excluded from the calibration data set to use for model validation ( supporting information Table S1 ). The model development approach follows increasing levels of complexity in order to arrive at a sufficient level of variability explained (coefficient of variation,). We used a combination of a three‐segment piecewise linear regression model to explain the mean bed temperature during ice‐covered periods () based on EDR and a multiple linear regression model to explain the mean bed temperature during open‐water periods () based on mean air temperature during the same open‐water period () and lake depth ( supporting information Text S2 ). These models were then combined to predict MABT in the general form 2 where ice‐cover duration (ICD) is ice‐cover duration in days and OWD is open‐water duration in days and equations forandare provided in the supporting information (Text S2 ). This empirical model was then validated with the excluded set of lakes, as well as comparisons to (1) a set of lakes outside the spatial domain used for model development (Toolik Lake node), (2) a set of lakes outside the time domain used for model development (2011–2012), and (3) a lake in Barrow, AK, during the year 1954–1955 [] to evaluate the model suitability for making hindcasts well beyond the data range used in its development ( supporting information Table S1 ).

Ice thickness was measured at each study lake during the late winter (late March to middle of April). Estimation of each year's maximum ice thickness (MIT) was determined by fitting observational data to a modeled ice growth curve ( supporting information Text S1 ). The effective depth ratio (EDR) was developed as a dimensionless index of MIT relative to lake depth 1 where lake depth was determined from atmospherically corrected pressure transducer data averaged during the month of September (period prior to freezeup) and most representative of depth over the winter. MIT was estimated individually for each lake per year, except for bedfast ice lakes (i.e., MIT = lake depth), for which we applied the regional mean from adjacent floating ice lakes. Thus, EDR > 1 is floating ice and EDR ≤ 1 is bedfast ice in that year at the observation point in the lake center.

(a) The surface thermal landscape of the Alaska North Slope showing the location of lakes monitored in this study (terrestrial mean annual ground temperature (MAGT, 2000–2009) based on GIPL model output [.,] and lacustrine lake ice regimes (bedfast ice (red) and floating ice (blue) based on synthetic aperture radar analysis []). (b) Example of lake thermal regimes over the winter for a bedfast ice lake (Tes‐005) and a floating ice lake (Tes‐001) including sublake permafrost, terrestrial permafrost, and air temperatures.

Monitoring of multiple lakes in Arctic Alaska began in 2012 as part of the Circum‐Arctic Lake Observation Network (CALON, www.arcticlakes.org ) and ended in 2015. The objective of CALON was to collect data to understand patterns in lake physical processes across broad physioclimatic gradients (foothills to coastal plain) and within region scales of lake area and depth. The primary observations utilized for this study consisted of continuously recorded lake temperature (surface and bed) from 29 lakes, located at seven nodes along two ~175 km north to south transects (Figure 1 a). Additional information on the CALON project, study lakes, sensors, and data collection and processing are detailed in the supporting information (Text S1 ).

3 Results and Discussion

3.1 Winter and Summer Controls on Lake Bed Thermal Regimes Temperature at the bed (water‐sediment interface) of shallow lakes shows distinctly different patterns between open‐water and ice‐covered periods (Figure 1b). Following spring ice out, bed temperature closely tracks air temperature through the summer. Once ice cover forms in the early winter, the lake water column becomes isothermal near 0°C and then temperature at the bed rapidly rises by several degrees over days to months depending on the amount of heat stored and released from the sediment. The extent and duration of early winter warming are likely indicative of the talik depth, though this relationship has not yet been quantified [Burn, 2005; Boike et al., 2015]. Bed temperature in floating ice lakes remains above freezing during the entire ice‐covered period, whereas bedfast ice lakes often freeze solid by early winter and maintain permafrost (Figure 1b). Sublake permafrost temperatures below the example bedfast ice lake in Figure 1b averaged −4.3°C, slightly warmer than adjacent terrestrial permafrost temperatures, −6.3°C. The range of MABT measured among 29 lakes over 3 years was greater than 10°C, ranging from −5.5°C in a 0.4 m deep coastal plain lake in 2013 to 5.1°C in a 2 m deep foothills lake also in 2013. Over these same set of years and study areas, MAAT varied by only 3°C and was unrelated to MABT (Figure 2a). This high interlake variability in MABT is primarily a function of lake depth relative to ice thickness (i.e., the effective depth ratio, EDR), as has been documented by previous studies [Burn, 2002, 2005; Arp et al., 2012]. Variation in ice‐cover duration in our data set (Figure 2b) was primarily a function of ice‐out timing among study lakes, which mainly is driven by distinctly earlier ice out of bedfast ice lakes relative to floating ice lakes [Arp et al., 2015] and by interannual variation in early summer air temperature that affects both lake types [Duguay et al., 2006]. It is this variability in ice‐cover duration that makes multiyear and interlake prediction of MABT challenging without explicitly accounting for lake surface conditions (i.e., ice covered versus open water). Figure 2 Open in figure viewer PowerPoint The processes controlling lake bed temperature and ice‐cover duration are described with increasing precision by comparison to (a) mean annual air temperature, (b) lake depth, (c) the effective depth ratio during ice‐covered periods, and (d) air temperature during open‐water periods. Partitioning MABT by ice‐covered and open‐water periods allowed us to develop simple and robust relations between bed temperature and the dominant controlling variables (Figures 2c and 2d). The ice‐cover model based on EDR shows a rapid increase in T b‐ice from −12°C to ~0°C in bedfast ice lakes up to an EDR of 0.75, after which T b‐ice increases much more slowly to an EDR of 1.94, with a very slight decrease in T b‐ice at higher EDR values (Figure 2c). Though these breaks in the relation between T b‐ice and EDR were determined statistically, we suggest that they represent important physical transitions in lake behavior. The associated temperatures at EDR thresholds occur just below 0°C (freezing point of fresh water) and ~4°C (maximum density of water). In floating ice lakes with liquid water present throughout the winter, temperature stratification often exceeds 2–3°C, promoting progressive sublake permafrost thaw even during the winter. Cold winter air temperatures drive lake heat loss and ice growth, but it is the amount of snow on lake ice that regulates ice growth and generally exerts a stronger influence on interannual and regional variations in MIT [Zhang and Jeffries, 2000]. Over the 3 year period of our study, MIT on floating ice lakes averaged 1.6 m and ranged from 1.3 m on foothills lakes in 2014 to 1.9 m on coastal plain lakes in 2013. Over longer periods in Barrow from 1947 to 1997, for example, MIT ranged from 1.3 to 2.5 m [Zhang and Jeffries, 2000]. Lake depth prior to freezeup can also vary interannually, but in the absence of geomorphic change, such as catastrophic lake drainage events [Jones and Arp, 2015], this is typically less than 10 cm for Arctic Coastal Plain (ACP) lakes [Jones et al., 2009]. Thus, our results suggest that changes in MABT in shallow Arctic lakes should primarily follow variation in MIT and ice‐cover duration. Long‐term trends toward thinner MIT in regions of the Arctic [Arp et al., 2012] imply similar trajectories for MABT in shallow thermokarst lakes and corresponding impacts on sublake permafrost stability.

3.2 Threshold Sensitivity of Shallow Arctic Lakes Sensitivity analysis showed that variation in winter process and changes in ice thickness and duration caused more variation in MABT than in changing summer conditions. An increase in summer (open‐water period) air temperatures of 3.2°C produced only a 0.5°C change in MABT on average with little variation by depth (Figure 3a). Variation in open‐water duration (OWD) increased MABT by 1.7°C in 0.5 m deep lakes, 1.1°C in 1.0 m deep lakes, and 0.4°C in 3.0 m deep lakes (Figure 3b). The changing sensitivity observed with depth is because temperatures in bedfast ice lakes prior to ice out are <0°C and increase greatly after ice out. For floating ice lakes, bed temperatures are above freezing when ice is present or absent, thus lessening its impact on MABT changes substantially. In fact, the beds of floating ice lakes often experience considerable warming during the period prior to ice out as increasing solar radiation penetrates the ice column as seen in Figure 1b and has been observed more widely in other Arctic lakes [Hinkel et al., 2012; Boike et al., 2015]. Figure 3 Open in figure viewer PowerPoint Results of mean annual bed temperature (MABT) sensitivity to air temperature during the (a) open‐water period (T a‐open ), (b) open‐water duration, and (c) maximum ice thickness (MIT) based on climate and ice thickness records in Barrow, AK, from 1969 to 2015. Lake depths of 0.5, 1.0, and 3.0 m were evaluated for each variable. Changes in MABT were much more sensitive to winter conditions in the form of ice thickness dynamics, both because this is the dominant season in the Arctic and because winter season climate is changing most rapidly [Wendler et al., 2014]. The median MIT for Barrow floating ice lakes from 1969 to 2015 was 1.9 m and ranged from 1.4 m to 2.3 m. In very shallow lakes (0.5 m) with consistent bedfast ice regimes over this full range of MIT, MABT increased by 1.7°C but with maximum values still well below freezing, −4.8°C (Figure 3c). In deeper bedfast ice lakes (1 m), MABT increased by 3.5°C and is projected to rise above freezing at MIT < 1.6 m (EDR > 0.65). Taliks should begin to form in lakes of this depth according to both field studies [Burn, 2002] and heat transfer models [Ling and Zhang, 2003]. Extending this analysis beyond the range of the Barrow record to even thinner MIT, 1.15 m or an EDR of 0.9 for a 1 m deep lake, helps outline the threshold response identified in our model, indicating that MABT stabilizes at temperatures above 0°C. This is a bedfast ice lake that just freezes solid by the end of winter and right at the threshold above which floating ice regimes occur. In deeper lakes (3 m) with consistent floating ice regimes, the impact of MIT is considerably smaller, causing a MABT change of 1.0°C (Figure 3c). Smaller bed temperature responses to climate warming have also been noted for deeper lakes in temperate regions due to summer stratification [Kraemer et al., 2015; Winslow et al., 2015].