Climate warming response of mountain lakes affected by variations in snow

Our objectives were to determine how temperatures in mountain lakes respond to changes in climate and to characterize how their responses are mediated by landscape or lake morphometric factors. Our analysis combines the use of high-frequency climate and lake temperature data from 1983 to 2016 in a high-elevation catchment in the Sierra Nevada of California with summer water temperature data from a set of 18 additional lakes scattered throughout the range. Average annual air temperatures warmed at 0.63 C decade, but variation in lake temperature was driven primarily by amount of precipitation as snow. By regulating the duration of ice cover and volume of inflowing spring snowmelt, variation in snowpack size accounted for 93% of variation in summer epilimnetic temperatures. The effect of snow on lake temperatures was mediated by variation in elevation and lake depth at landscape scales, creating a predictable mosaic of lake sensitivities to climate warming. *Correspondence: ssadro@ucdavis.edu Associate editor: Jordan Read Author Contribution Statement: SS, JMM, and JOS conceived of the study. KS maintained field equipment and conducted quality assurance and data analysis. SS collected field data, conducted quality assurance of data, did analyses, and made figures and tables. SS led and JMM, JOS, and KS contributed to interpretation of the data and writing of the manuscript. Data Availability Statement: All climatological, hydrological, and limnological data from Emerald Lake and the lake temperature and snow data for the 19 lakes used in the scaling analysis are archived with the Environmental Data Initiative (S. Sadro, Climate sensitivity of Sierra Nevada lakes, 2018, doi:10.6073/ pasta/fa5b1426f5f89aed4233d4d607554224) at https://portal.edirepository.org/nis/mapbrowse?packageid=edi.237.1 Additional Supporting Information may be found in the online version of this article. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

Climate change ranks high among influences on lakes  in part because lake temperature regulates a wide range of ecosystem functions. Processes affecting temperatures in lakes influence mixing dynamics (Kirillin 2010), alter biogeochemical rates (Adrian et al. 2009), influence metabolic rates (Yvon-Durocher et al. 2010), and affect primary productivity and foodweb dynamics (Parker et al. 2008;Woodward et al. 2010;Greig et al. 2012;Preston et al. 2016). Although there is ample evidence that lake temperatures across the globe are affected by climate, temperature responses are variable in space and time, and our ability to make predictions is hindered by an inadequate understanding of the mechanisms involved (Kraemer et al. 2015;O'Reilly et al. 2015).
Uncertainty in lake temperature responses stems in part from the multiple ways in which climate interacts with lake heat budgets. The factors contributing to the energy balance include net solar radiation, sensible heat exchange with the atmosphere and lake sediments, latent heat exchange associated with evaporation, and advection of water from outside the lake. In lakes that experience winter ice cover, processes regulating ice phenology and snowmelt runoff may be particularly important because they influence the duration over which lakes gain heat (O'Reilly et al. 2015). Although heat budgets reflect large-scale climate forcing, effects are often modulated by factors operating on nested spatial scales. For example, variation in landscape characteristics such as latitude, elevation, and catchment characteristics or land cover might modulate climate effects on individual lakes (Schmid et al. 2014). The climate signal might be further modified at the individual lake through variation in morphometric attributes such as lake size and shape or through differences in the source and magnitude of water inputs (Rose et al. 2016). As a result, predicting the temperature response of individual lakes requires an understanding of the dominant processes involved and the extent to which they are mediated by other factors.
Mountain regions are ideal ecosystems in which to explore climate warming effects on lake temperatures. They are sensitive to climate change and have some of the highest recorded rates of warming in air temperature (Pepin and Lundquist 2008;Pepin et al. 2015;Preston et al. 2016). Just as importantly, the proportion of precipitation falling as snow in many mountain regions is declining (Lundquist et al. 2009;Berg and Hall 2017). Because they have distinct hydrological patterns driven by the deposition and subsequent melt of winter snowpack (Lundquist et al. 2009;Sadro et al. 2018), the potential impact on lake heat budgets may be substantial (Strub et al. 1985). Finally, there are large gradients in local landscape and lake morphometric features in mountains that are expected to mediate the effect of climate warming (Livingstone et al. 2005;Sadro et al. 2012). While the cumulative effects of warmer air temperatures and reduced snowpack on lake temperatures are reduced ice-duration and warmer surface waters (Livingstone et al. 2005; Thompson et al. 2005;Preston et al. 2016), how these processes vary at landscape scales because of other factors remains poorly understood.
We analyzed the thermal response in the upper mixed layer of a lake in the Sierra Nevada of California to interannual variation in air temperature, snow deposition, and other climate factors to characterize relative effects on water temperature. We then use summer temperature data from 19 lakes to understand how the relative importance of snow and other factors governing lake temperature vary at broad spatial scales and predict sensitivity to warming in over 1600 lakes across the region. Our study has three specific objectives: (1) to characterize how water temperatures in mountain lakes are responding to variation in air temperature, snow deposition, and other climate factors; (2) to evaluate scaling relationships for water temperature using landscape and lake morphometric attributes; and (3) using those scaling relationships, to identify lakes that are most sensitive to warming from ongoing changes in climate. Our results emphasize the high rate of climate warming taking place in mountain ecosystems and demonstrate the substantial role of snowpack in governing lake temperature. We illustrate the extent to which warming within lakes scales with elevation and lake morphometric attributes and use those empirical relationships to identify lakes most sensitive to ongoing changes in climate.

Site description
This analysis combines the use of high-frequency data from a long-term study site (Emerald Lake) with summer lake temperature data from an additional 18 lakes scattered throughout the Sierra Nevada of California (Supporting Information  Table S1). Emerald Lake is located in Sequoia-Kings Canyon National Park in the southern Sierra Nevada. The lake is 2.7 ha in area, has a maximum depth of 10 m, and a morphometry typical of glacially scoured mountain lakes. The lake is representative of thousands of lakes located throughout the Sierra Nevada (Melack and Stoddard 1991;Tonnessen 1991;Sickman et al. 2001). The 18 additional lakes are part of the U.S. National Park Service Inventory Monitoring Program or lakes we had historically sampled. Elevation, lake area, and maximum depth ranges in these lakes were 2475 to 3770 meters above sea level (m.a.s.l.), 0.7 to 12.6 ha, and 5.1 to 35.0 m, respectively (Supporting Information Table S1). Lakes shallower than 5 m were excluded as heat budget dynamics were expected to differ from deeper stratified lakes. Lakes used in this analysis span subalpine and alpine zones where a majority of precipitation falls as snow between October and April. Snowmelt volume in an average year exceeds total lake volume for many lakes, which means lake residence times can range from a few days to over a year on a seasonal basis (Sadro et al. 2018).

Sadro et al.
Climate sensitivity of mountain lakes Meteorological and hydrological data Meteorological measurements were made at a station located approximately 50 m from the Emerald Lake shore (Supporting Information). Air temperature, relative humidity, short and longwave radiation, and wind speed were collected at 10 s intervals and recorded as hourly averages (Sadro et al. 2018). Snow water equivalent (SWE) was computed by measuring snow depth and density at multiple locations throughout the Emerald Lake watershed at the time of maximum accumulation in late, March or early April (Sadro et al. 2018). The date of ice-off was determined through visual inspection of Landsat 5, 7, and 8 GeoTIFF images containing Emerald Lake (http://earthexplorer.usgs.gov). Given the 16 d interval between successive Landsat images, ice-out date was determined as the mid-point between the dates of two Landsat scenes, one containing ice and the other ice free, resulting in an uncertainty of~8 d or 10% of the range in ice-out dates. In years where cloud cover prevented visual determination of ice-out within a 2 week period, no ice-out date was recorded.
Discharge was measured in the outlet stream of Emerald Lake using a calibrated V-notch weir and vented pressure transducers located approximately 5 m downstream from the lake. A continuous record of stage was computed from highfrequency pressure measurements using linear regressions between pressure and staff gauge measurements. Discharge was computed from the stage record using a rating curve. Groundwater inputs to Emerald Lake are small, and we assume outflows = inflows (Kattelmann and Elder 1991). Detailed descriptions of hydrological methods and uncertainty in measurements of discharge and SWE have been published (Sadro et al. 2018).
Emerald Lake outlet temperature (OT) was measured at 10 s intervals with a thermistor located at the weir and recoded as hourly averages on a Campbell Scientific data logger. Water temperatures in the 18 additional lakes used in the scaling analysis (Supporting Information Table S1) were measured manually at the outlet using a multiparameter sensor (Yellow Springs Instruments). We used OT in our analysis as a proxy for lake surface temperature. It is a good approximation for epilimnetic temperature in Emerald Lake, because water clarity is high, maximum lake depth is only 10 m, and a majority of the lake volume often falls within the active mixing layer (Supporting Information Fig. S1). Daily mean outlet and epilimnetic temperatures in Emerald Lake are highly correlated, especially in the spring and summer. Root-mean-square error is lower in the spring and summer (0.5-0.6 C) than in the autumn (1.8 C), when discharge is low and OT becomes more sensitive to changes in air temperature.

Data processing and statistical analyses
Statistical analyses were done using JMP version 12 (SAS Institute) and MATLAB version R2007b (MathWorks) (Supporting Information). Principal component analysis (PCA) was used to characterize relationships among climatic factors influencing lake temperature. The date of autumn turnover (Dmix) in Emerald Lake was determined as the first day when the temperature gradient throughout the water column was less than 1 C as determined from vertical temperature profiles collected at 1 m intervals with a minimum frequency of every 2 weeks. Sensitivity to warming associated with loss of snowpack in each of the lakes was computed by dividing the range in maximum summer water temperature by the range in SWE across all years for which we had data. We computed the range in SWE across years (Supporting Information Table S1) for each lake using data from all available snow course or pillow sites throughout the Sierra (http:// cdec.water.ca.gov/snow/current/snow/index.html).
Because SWE can vary substantially through both time and space in the Sierra (Girotto et al. 2014), for each lake-year in our analysis, we used a mean SWE value corresponding to the elevation of the lake rather than the value from the closest available snow survey or pillow site. Elevation-specific SWE was computed by binning SWE data by elevation at 100 m intervals and computing mean SWE for each elevation bin in each year. SWE anomaly was computed as departure from the long-term mean. Years when SWE anomaly was < −500 mm were classified as "drought" (27% quantile), > +500 mm anomaly as "wet" (74% quantile), and AE 500 mm anomaly as "average." Where analyses were done on a seasonal basis, winter was defined as December-February, spring as March-May, summer as June-August, and autumn as September-November. Additional details regarding data analysis are provided in the Supporting Information.

Results and discussion
Climate trends and the importance of snow to lake temperatures Trends in air temperature and precipitation have differed in the Emerald watershed over the last 30 yr. Air temperature has shown a clear and consistent warming trend. The average annual rate of increase was 0.63 C decade −1 between 1983 and 2016 ( Fig. 1A; Supporting Information Table S2, Eq. 1), caused largely by summer warming rates that reached 1.0 C decade −1 (Supporting Information Table S2, Eq. 2). These are among the highest rates of warming observed globally and consistent with patterns found in other mountain ranges (Pachauri et al. 2015;Pepin et al. 2015). In contrast, precipitation was highly variable and had only a weak long-term trend. SWE, which characterizes the water content of the snowpack, has decrease by 30 mm yr −1 in the Emerald watershed ( Despite the strong climate warming signal in air temperature, average annual water temperature in Emerald Lake over the last 20 yr was variable and did not have a significant trend ( Fig. 1C; Supporting Information Table S2, Eq. 4). The discrepancy in responses between air and lake temperatures reflects the effect of precipitation as snow. Although the relative importance of snow and duration of ice cover on mountain lake heat budgets has been recognized, variation in snowpack is not often included in predictive models (Strub et al. 1985;Roberts et al. 2017). At Emerald Lake, annual variation in SWE accounted for over 90% of variation in mean annual and mean summer lake temperatures (Supporting Information  Table S2, Eqs. 5 and 6), and a significant warming trend in water temperature was only evident during drought years ( Fig. 1C; Supporting Information Table S2, Eq. 7, which characterizes warming through time during only drought years). The linear relationship between SWE and mean annual water temperature provides a powerful statistical tool for predicting lake temperatures in snowmelt-dominated systems and offers a mechanistic explanation for why many lakes around the world do not have warming trends matching air temperatures in their catchments (O'Reilly et al. 2015).
A PCA indicates that the strength and structure of the SWE-lake temperature relationship varies seasonally in Emerald Lake (Fig. 2). During the spring and summer, interannual variation in SWE explained a majority of variation in water temperature ( Fig. 2A,B,D,E). Although lake temperature was negatively correlated with SWE in both seasons, the response differed ( Fig. 2; Supporting Information Table S3). During spring, there was a threshold in SWE corresponding to average or wet years over which lake temperatures remained uniformly cold and below which water temperatures increased linearly with decreasing SWE (Fig. 2B,C; Supporting Information Table S2, Eq. 8). In contrast, summer water temperatures varied linearly by 13 C across the entire range of SWE ( Fig. 2E; Supporting Information Table S2, Eqs. 10 and 11). For every 10 cm decrease in SWE in the Emerald Lake watershed, there was a 0.56 C increase in mean summer lake temperature (Supporting Information Table S2, Eq. 6) and 0.26 C increase in maximum lake temperature (Supporting Information Table S2, Eq. 13).
The role of snow in structuring spring and summer lake temperature was so large that the effect of other climate factors was only evident during drought years (Fig. 2F). In years with little snow, lake temperature was correlated with air temperature and negatively correlated with relative humidity and wind speed (Fig. 2D,F), which when combined with other climate variables explained 96-99% of the variation in water temperature (Supporting Information Table S4). For every degree increase in air temperature during the spring and summer of drought years, the lake warmed by 0.52-0.53 C ( Fig. 2C; Supporting Information Table S2, Eqs. 9 and 11). As a consequence, since 1990, mean summer water temperature in Emerald Lake has increased by 0.4 C decade −1 . This warming trend is consistent with global averages (O'Reilly et al. 2015) but is only evident during drought years when the effects of large snowpacks are not apparent (Fig. 1C).
Autumn differs from spring and summer in that lake temperatures are governed primarily by rates of heat loss. Consequently, there was little effect of snow on lake temperature (Fig. 2G,H). Instead, water temperature was positively correlated with the date at which the lake began mixing and with shortwave radiation levels and negatively correlated with wind speed and relative humidity, which together accounted for 97% of variation in autumn lake temperatures (Supporting Information Table S4). These relationships reflect the importance of convective and wind-induced turbulence in redistributing heat and regulating its loss from the lake in the autumn.

Mechanisms governing water temperature in snowmeltdominated lakes
The relationship between SWE and water temperature reflects a combination of direct and indirect effects associated with snowpack size. The indirect effect stems from the role of snow in regulating ice-off date (Supporting Information Fig. 1. Long-term patterns in climate and lake temperatures for the Emerald Lake watershed. Panel A illustrates the warming trend in air temperature. Panel B highlights variation in precipitation as snow. Panel C illustrates changes in lake temperature through time in relation to precipitation. Color coding of bars and symbols corresponds to precipitation levels: years < −500 mm SWE anomaly are classified as "drought" (orange); years > +500 mm SWE anomaly are classified as "wet" (blue), and years AE 500 mm SWE anomaly are classified as "average" (white). Linear models are computed using ordinary least squares.  Table S2, Eq. 14). Because lake-ice in the Sierra Nevada is a mixture of snow-slush and ice-lenses, it can range from less than 0.5 m to over 3 m in thickness depending on accumulated snowfall. Large snowpacks contain more water, melt faster because they begin melting later in the season, but take longer to melt overall than thinner snowpacks (Harpold et al. 2012;Musselman et al. 2017a). Consequently, in years with large snowpacks, lakes begin heating later than those with shallow snowpacks. Second, the water content of snowpacks directly affects snowmelt dynamics (Supporting Information Fig. S2B; Supporting Information Table S2, Eq. 15), lake residence time, and the duration of time that cold meltwaters enter lakes (Sadro et al. 2018), all important factors affecting the lake heat budget.
These mechanisms can be confirmed through analysis of discharge and water temperature in Emerald Lake (Fig. 3). Differences in snowmelt were substantial over the course of this study (Sadro et al. 2018). SWE varied by an order of magnitude and discharge by a factor of 5. The spring snowmelt pulse began 24 d earlier and reached completion 52 d sooner Fig. 2. Climatic factors governing lake temperature differ seasonally. Panels A-C correspond to spring, panels D-F to summer, and panels G-I to autumn. Plots A, D, and G are PCA, illustrating relationships among various climate factors and lake temperature. Plots B, E, and H illustrate lake temperature anomaly in response to SWE anomaly. Plots C, F, and I illustrate lake temperature anomaly in response to air temperature anomaly. PCA codes are: LT, lake temperature; SWE, snow water equivalent; AT, air temperature; WS, wind speed; RH, relative humidity; SWrad, shortwave radiation; LWrad, longwave radiation; Dmix, day of the year in autumn when lake water column mixing begins. See Figure 1 caption for color coding. on average in drought years, which in turn affected when the lake began warming. The onset of lake warming in years with a large snowpack was delayed an average of 35 d, causing temperatures to peak an average of 50 d later and 27% lower than during the drought years (Fig. 3). As a consequence, dry years warmed~1 C month −1 faster than wet years (Supporting Information Table S2, Eq. 16). Because snow likely plays an important role in the heat budget of many lakes throughout the world, and that in many of these same regions SWE is measured or modeled to assess water availability (Lakshmi and American Geophysical Union 2014), understanding how the snow-lake temperature relationship scales across lakes will be critical for predicting regional responses of mountain lakes to climate variability.
Scaling warming from loss of snow to predict lake thermal sensitivity To determine the thermal sensitivity of Sierra Nevada lakes, we characterized how warming associated with declining snowpack scaled with landscape and morphometric attributes expected to affect temperature (Luoto and Nevalainen 2013;Kraemer et al. 2015) in a set of 19 lakes that includes Emerald Lake (Supporting Information Table S1). The range in thermal sensitivity of summer water temperature across all lakes was 0.16-0.94 C per 10 cm −1 decline in SWE. Thermal sensitivity scaled with elevation and maximum depth (excluding one outlier lake with vestigial glacier in the watershed), but there was no obvious relationship with lake area (Fig. 4). Elevation explained the largest amount of variation (Supporting Information Table S2, Eq. 17). Thermal sensitivity tended to increase with increasing elevation until 3000-3500 m.a.s.l., after which it declined. This matches patterns in snow deposition and melt rate with elevation in the Sierra Nevada (Musselman et al. 2017b). Unlike the pattern with elevation, thermal sensitivity scaled linearly with maximum lake depth (Supporting Information Table S2, Eq.18). Finally, the large variability in thermal sensitivity among lakes smaller than 6 ha prevented a significant relationship with lake area despite larger lakes having generally high thermal sensitivity. When combined in a multiple regression model, elevation, maximum depth, and lake area together explained 41% of variation in temperature sensitivity (Supporting Information Table S2,Eq. 19).
Despite the uncertainty in models, scaling-up to explore lake thermal sensitivity throughout the Sierra Nevada is a useful exercise to demonstrate the extent to which climate warming in lakes is mediated by landscape and lake morphometric attributes and provide an understanding of lake sensitivity at regional scales. We use our depth-elevation-lake area model Fig. 3. A comparison of discharge (A) and lake temperature (B) across years spanning a wide gradient in snow precipitation. Blue lines correspond to wet years (SWE anomaly > +500 mm), orange lines correspond to drought years (SWE anomaly < −500 mm), and gray lines average years (SWE anomaly AE 500 mm). Fig. 4. Scaling relationships for warming in summer water temperatures associated with loss of snow for 19 lakes. Open circle symbol is Emerald Lake.
Warming from loss of snow is scaled according to variation in: (A) elevation, (B) lake area, and (C) maximum lake depth.
(Supporting Information Table S2, Eq. 19) to estimate thermal sensitivity associated with loss of snow in 1625 Sierra Nevada lakes smaller than 30 ha, where all three factors are known (Fig. 5). We use a multiple linear regression model to characterize thermal sensitivity because it is unclear how interactions among factors might be contributing to the nonlinearity in thermal sensitivity observed with elevation alone (Fig. 4A). Thermal sensitivity ranged from 0.2 C to 1.8 C per 10 cm −1 decline in SWE and averaged 0.59 C per 10 cm −1 decline in SWE. Thermal sensitivity increased slightly moving northward through the range and decreased slightly moving eastward toward the Sierra crest and higher elevations (Fig. 5). Over the next century, snowpack across the Sierra Nevada is predicted to decline by over 50% overall, with the magnitude of decline varying with elevation (Knowles and Cayan 2002;Berg and Hall 2017;Fyfe et al. 2017). Based on the loss of snow predicted by the RCP4.5 climate warming scenario, we estimate that summer lake temperatures in the Sierra Nevada may warm 1.1-10.5 C by the end of the 21st century, with the majority of lakes warming at least 3.1 C. While we currently lack the data necessary to make more precise predictions, in general, we expect larger increase in temperature in lakes at low to mid elevations, which are expected to experience larger declines in SWE overall (Berg and Hall 2017;Fyfe et al. 2017).
Our results demonstrate the importance of snow in lake heat budgets and highlight the spatial variability in lake thermal sensitivity to future declines in snowpack. However, there are a number of important limitations that constrain the predictive power of our analysis for individual lakes. First, there is uncertainty in how accurately summer water temperature data from the 18 additional lakes used in our scaling analysis capture true temperature maxima. Second, the spatial distribution of lakes does not represent the full extent of the latitudinal or elevational range found throughout the Sierra. Finally, our data did not have the capacity to explore the role of several factors that are likely to be important in heat budgets and are likely to be contributing uncertainty (e.g., wind exposure, local topographic relief, and light attenuation coefficients). Given the paucity of water temperature data available for lakes and range of different factors that might be contributing to uncertainty, it is difficult to evaluate the extent to which our predictions may be over or under estimating thermal sensitivity in lakes. Despite these limitations, our results demonstrate the importance of SWE in lake heat budgets and the sensitivity of lakes to warming associated with the loss of snow.

Conclusions
Climate change is impacting the Sierra Nevada of California. What warming water temperatures will mean for the ecosystem function of these lakes is an active area of study. Small changes in water temperature are known to affect biological processes and ecosystem function in mountain lakes (Parker et al. 2008;Miller and McKnight 2015;Preston et al. 2016). In the Sierra, snowpack controls on the timing of nutrient delivery and warming have important implications for phytoplankton biomass (Sadro et al. 2018). Our study suggests that throughout the Sierra Nevada of California. Inset for California identifies the Sierra Nevada along the 2200 m.a.s.l. elevation contour (black area) and the boundaries of the enlarged map area (gray shaded). The enlarged map shows the locations of the 19 lakes used in the scaling analysis (black). Thermal sensitivity ( C warming in water temperature per 10 cm decline in SWE) of lakes is coded according to upper 25% quantile (orange), lower 25% quantile (blue), and middle 50% (yellow). ecosystem processes in snowmelt-dominated lakes are most at risk from climatic shifts that reduce the accumulation of snow during the winter, extend summer-like conditions into the spring, and delay the onset of cooling in the autumn. We have shown that interactions among multiple factors can result in complex lake thermal responses, producing years with cold lake temperatures despite record high air temperatures. Such temporal variation, along with spatial variation in lake sensitivity to loss of snow, may offer important refugia for ecosystems (Sedell et al. 1990;Magoulick and Kobza 2003;Isaak et al. 2010) and possibly mitigate or delay longer term climate warming trends. Understanding how local-and regional-scale processes interact to affect water temperature in individual lakes will be critical for predicting the ecosystem responses at regional scales.