Low snowpack reduces thermal response diversity among streams across a landscape

Spatial and temporal variation in thermal conditions are important dimensions of aquatic landscapes, yet we do not understand how this heterogeneity will respond to climate change. Snowpack in many regions is declining but impacts on aquatic thermal regimes remain poorly understood. Our analyses of summer stream temperatures across a complex river basin in southwest Alaska show that loss of snowpack has disproportionate effects on water temperatures in streams that are cold under typical conditions, thereby homogenizing stream temperatures across the landscape during low snow conditions. Streams draining steep watersheds warmed from a summer average of 4–8°C between high and low snow years and became more responsive to variation in regional air temperature (as proxy for solar radiation), while streams draining flat watersheds changed negligibly. Our results suggest that conservation strategies for aquatic landscapes that focus on snow‐dominated cold‐water refugia may not be robust to future climate warming and changes in snowpack.

Ecosystem responses to ongoing climate change remain difficult to predict but will ultimately determine their capacity to support biodiversity, and the variety of ecosystem functions and services that humans depend on (Lawler et al. 2010). In aquatic systems, water temperature exerts particularly important control over most ecosystem functions, including determination of habitat quality for aquatic organisms (Torgersen et al. 1999). Many of the mechanisms that control thermal regimes within aquatic ecosystems are well understood (Caissie 2006), but the complexities of how ongoing climate change interacts with topography to influence landscape-level thermal regimes is challenging for the development of effective management and conservation strategies at landscape and regional scales (Lawler et al. 2010;Steel et al. 2017).
Landscapes are complex, and local topography filters the effects of regional air temperature, solar radiation, and precipitation on hydrology and aquatic thermal regimes across the landscape (Poole and Berman 2001;Lisi et al. 2013Lisi et al. , 2015Luce et al. 2014). This landscape filtering results in thermal heterogeneity expressed at multiple scales, which plays a vital role for organisms and the functioning of ecosystems Schindler et al. 2013). For example, many terrestrial and aquatic consumers exploit the seasonal progression of thermal variability across the landscape to access food resources for a longer duration (Armstrong et al. 2016). Climate change has the potential to erode this heterogeneity in abiotic conditions; thus, it is important to understand how future climate changes will either support or diminish the thermal heterogeneity that underpins these resource waves and the reliable production of ecosystems services (Armstrong et al. 2016;Brennan et al. 2019).
Conservation efforts have put significant emphasis on identifying habitats that may be resistant to climate warming and act as refugia for cold-adapted species in a warming world (Isaak et al. 2015(Isaak et al. , 2016. Habitats that are less sensitive relative to changes in air temperature are often selected as climate refugia, under the assumption that they will remain cold in the future (Null et al. 2013;Isaak et al. 2015Isaak et al. , 2016. Surface waters, like air temperature, are primarily heated through solar radiation (Mohseni and Stefan 1999). Sensitivity of water temperature to variation in air temperature is commonly used as a reasonable proxy for understanding the influence of solar heating at the air-stream interface and is valuable for understanding the physical processes that regulate stream temperature (Mohseni and Stefan 1999;Lisi et al. 2015). Thermal sensitivities of streams are also controlled by channel and landscape characteristics, which mediate climate warming and thermal habitat shifts (Poole and Berman 2001;Kelleher et al. 2012;Briggs et al. 2018). However, a significant buffer for these systems is the contribution of snowmelt to stream flows (Leach and Moore 2014; Luce et al. 2014;Lisi et al. 2015), and snowpack is changing in a warming climate (Barnett et al. 2005;Fyfe et al. 2017). Altered winter precipitation coinciding with substantially warmer winter temperatures could reduce winter snow accumulation, depending on the rate of warming and the elevation of a site (Barnett et al. 2005, Pepin et al. 2015, Fyfe et al. 2017. The complexities associated with how winter precipitation is delivered and stored on the landscape are a key source of uncertainty in our understanding of how the thermal regimes of aquatic ecosystems will respond to ongoing climate change. In this study, we used stream temperature monitoring along a geomorphic gradient coupled with satellite-based snow coverage to show that loss of snowpack disproportionately affects temperatures in cold streams that drain steep higherelevation watersheds thus homogenizing stream temperatures across the river network in low snow years. Strong contrast in snowpack among years, including two of the snowiest (2012,2013) and two of the least snowy years (2015, 2016) on record, allowed us to explore how future changes in snowpack might affect stream temperatures and temperature sensitivity (τ), the change in daily water temperature relative to variation in air temperature. We evaluated how spatial variation in thermal conditions among streams changes in response to changing climatic conditions and how that may have consequences for the functioning of aquatic landscapes.

Study site
The Wood River basin, located in southwest Alaska, contains a diverse network of lakes, small tributaries, and intermediately sized rivers. Glacial activity in the region created strong contrast in stream geomorphic conditions across the basin. The east side is characterized by flat, meandering streams whose hydrology is rain-dominated, while the west and north sides have fast moving tributaries draining a mountainous landscape with hydrology dominated by water originating from snowmelt (Lisi et al. 2015). These systems have short water residence times and thus minimal deep groundwater inputs. All major tributaries within the Wood River basin drain into a series of connected lakes that all have elevations less than 40 m above sea level. An important aspect of having all tributaries terminate in lakes of similar elevation is that mean watershed slope and mean watershed elevation are highly correlated for each tributary (r = 0.88). Average elevation among watersheds included in this study varied from 45 to 520 m above sea level, and several watersheds have portions above 1500 m elevation.
The Wood River drains into Bristol Bay and is one of several major rivers in the region supporting commercially and culturally important Pacific salmon populations. Dozens of discrete sockeye salmon (Oncorhynchus nerka) populations spawn in tributaries, rivers, and along lake beaches across the basin . The geomorphic diversity within the basin drives differential temperature and stream flow regimes among tributaries, conditions to which local sockeye salmon populations have adapted their life-history strategies (e.g., reproductive timing and strategies, Lisi et al. 2013, Larson et al. 2016).

Data
We monitored water temperatures in streams across the Wood River basin during the summers of 2011 through 2016 (number of streams monitored per year: 34, 34, 25, 29, 23, and 35, respectively). Tributaries were selected to span the geomorphic variation observed for 2 nd -4 th order streams across the river basin. Due to occasional equipment loss or failure, the final set of streams with temperature monitoring was not identical from year to year. In our analyses, we only included streams with at least 3 yr of summer temperature data.
Stream water temperatures were recorded from approximately June 1 st through September 1 st at 60-min intervals with data loggers (i.e., iButton, Onset Tidbit, Onset Hobo Pro V1 temperature loggers, and Onset Hobo U20 stage gauges). The accuracy of these loggers ranged from AE 0.21 C for tidbits to AE 1 C for iButtons. Temperature loggers were placed a few centimeters above the substrate in each stream within about 100-500 m of the stream mouth, upstream of any lake water influence. Reference temperatures were taken manually at deployment and retrieval to ensure logger accuracy. Air temperatures were recorded at 60-min intervals at the University of Washington-Alaska Salmon Program's Lake Aleknagik field station on the Wood River (southern end of the river basin). Because stream temperatures were taken near stream mouths, there was little difference in elevation between stream sites and the weather station. The maximum distance and difference in elevation between the air temperature site and our stream temperature loggers was 57.7 km and 33.8 m, respectively. We averaged across diel temperature cycles by calculating mean daily stream temperatures for each tributary and mean daily air temperatures. All daily stream and air temperature time-series were z-score standardized within each stream-year prior to analysis.
Stream discharge plays an important role in regulating stream temperature and temperature sensitivity by increasing the volume of water and thus the amount of heat required to warm it. We included cumulative summer precipitation (01 June-01 September) as a proxy for contributions to midsummer stream flow, while recognizing that the effects of rainfall on stream discharge vary across basins. Precipitation was calculated from daily observations of a rain gauge (location noted on the map).
In the Wood River basin, tributary elevation and watershed slope are highly correlated. We use watershed slope for all subsequent analyses as it captures more aspects of stream geomorphology that may influence thermal sensitivity, such as channel confinement, soil depth, and catchment aspect. Watershed slopes were calculated as the average of all slope raster pixels within the boundaries of each individual stream's watershed using ArcGIS (v10.0, Environmental Systems Research Institute, Redlands, CA). Catchment aspect statistics (percent of the watershed facing North) were calculated by dividing the area of north facing cells (315-45 ) by the total catchment area for each study watershed. The percent North facing was logit-transformed prior to regression analysis.
We quantified changes in snowpack using satellite measurements of snow cover from 500-m grid resolution MODIS satellite records provided by the National Snow and Ice Center; there are no direct observations of snow depth and coverage for this remote region. These data were further processed by the Geographic Information Network of Alaska to correct for cloud cover interference and to estimate the date of the end of the continuous snow cover season for each grid cell (Lindsay et al. 2015). End of continuous snow cover season is defined by Lindsay et al. (2015) as the first of 14 consecutive days classified as no snow (less than 50% fractional cover and less than 30% albedo) for a given pixel. We then calculated the mean end-of-snow season date across pixels within each tributary to be used as an index of local snowpack conditions.

Statistical analyses
We used dynamic factor analysis (DFA) to identify common trends in stream temperature time-series and to quantify the influence of air temperature (Zuur et al. 2003). DFA is a dimension reduction technique specifically designed for timeseries analysis, allowing characterization of common trends shared among streams, thereby reducing the number of parameters that would need to be estimated if each stream were analyzed individually. The DFA model used in this study can be written as: where y is a matrix of stream temperature time series, t is the day in the time series, x is the latent common trend, d is a vector of explanatory variables (air temperature), and w and v are vectors of process and observation errors, respectively, which are multivariate normal in distribution, with mean 0 and covariance matrices Q and R. Z and D are matrices of stream-specific loadings on the latent trend and explanatory variables, respectively. To make the model identifiable process, error covariance matrix Q is fixed as an identity matrix (Zuur et al. 2003). R is the covariance matrix for the observation errors. We used an unconstrained error structure where each stream has its own variance and deviations among streams are allowed to covary (Lisi et al. 2015). The model was programmed and parameters were estimated using "TMB" package in R (R Core Team 2019). Model code was validated with simulated data sets and accurately recovered "true" parameter values.
To evaluate changes in thermal sensitivity, we calculated stream temperature sensitivity (τ), the per-degree change in daily stream temperatures per daily air temperatures. We computed τ by transforming stream-specific parameter estimates for the influence of the air temperature covariate (D; Eq. 1) from z-score standardized space (σ stream /σ air ) to C stream / C air using stream-year specific means and standard deviations.
We applied linear mixed effects regressions of average summer (June 15 th -September 1 st ) stream temperatures (or stream temperature sensitivity) against watershed slope, snow index, summer precipitation, and percent of the watershed facing North to determine how watershed geomorphology and snowpack affect stream thermal conditions among watersheds. As watershed slope is highly correlated with snow cover (due to elevation being highly correlated with watershed slope), we used the residuals from a linear mixed effects model of snow cover as a function of watershed slope as an index of snowpack (SI): where ESS i,y is the end-of-snow season for stream i in year y, α y is a year-specific intercept, and Slp i is a stream-specific watershed slope. Our snow index SI i,y is calculated from the residuals of Eq. 5 with the year-specific intercept included We estimated the effect of watershed slope, percent of the watershed facing North, cumulative precipitation, and end-ofsnow season on stream temperature and τ using linear mixed effects models. Mixed effects models were estimated using the "lme4" package and candidate models were compared using "MuMIn" package in R (R Core Team 2019). All models included random intercepts to account for repeated measures within individual streams across years. Results from the best models are presented in Fig. 3 and model selection for each is included in Supporting Information Tables S1, S2. Model selection and model weights (Wagenmakers and Farrell 2004) were calculated using Akaike's Information Criterion corrected for small sample sizes, which penalizes models with additional complexity relative to available data (AICc). Model weights were calculated on all models with a ΔAICc < 6 and weighted predictions are presented in Fig. 3.

Results and discussion
Quantifying thermal regimes across the landscape The observed range in average temperatures among streams across the landscape varied among years (Fig. 1A-F). In 2012, average daily stream temperatures ranged from 3 C to 15 C across the river basin, whereas in 2016 average daily stream temperatures ranged from 7 C to 16 C (Fig. 1). In 2011 and 2012, few streams had maximum average daily temperatures above 15 C (Fig. 1A,B). In the four other years, many streams had prolonged periods of temperature in excess of 15 C (Fig. 1E,F). Overall, temporal variation in stream temperatures appears strongly associated with variation in summer air temperatures (a proxy for intensity of solar radiation) but the degree of association differed widely among streams (Fig. 1).
air temperature s shared a common trend of steadily increasing temperatures as summer progressed, and this trend was not correlated with daily changes in air temperature ( Fig. 1A-F, black lines). Overall, this common trend was more prevalent in high-gradient cold-water streams (Supporting Information Fig. S1). There was substantial variation in τ among streams both within and among years. Much of this variation was explained by a negative association with stream watershed slope ( Fig. 1G-L, R 2 for 2011-2016 are 0.54, 0.38, 0.46, 0.21, 0.46, and 0.05) where streams draining steeper watersheds were distinctly less sensitive to changes in air temperature than were streams in flat watersheds. Sensitivity values generally range from near 0 C to 0.7 C, similar to values found for other mountain rivers systems (Isaak et al. 2016), but the magnitude of τ and the relationship between stream geomorphology varied among years.

Changes in snow at the landscape scale
Over the 6 yr of the current study, this region experienced two of the snowiest (2012, 2013) and two of the least snowy (2015, 2016) years on record (Fig. 2), providing a unique opportunity to assess the effects of winter snow conditions on summer thermal regimes across the landscape. A comparison of snowmelt dates between 2012 (high snow) and 2015 (low snow) shows how extreme the differences in early spring snow cover were in these 2 yr ( Fig. 2A,B). In 2012, most of the stream watersheds retained snow into late May (> day of year (DOY) 140), and the steepest watersheds had snow remaining into mid-July. In contrast, in 2015, most of the river basin was devoid of snow by May 1 st , and only small relict snow patches persisted only at high elevations into July. When comparing watershed level changes in summer snow coverage across years, streams draining steep watersheds (slopes > 15 ) lost more than 3 weeks of snow coverage between high-snow and low-snow conditions (Fig. 2D). In many flat streams, snow had disappeared well before May 1 st (DOY 140) in 2015 and is unlikely to figure into summer stream thermal characteristics. The average date at which snow was no longer present in a watershed, the end-of-snow season, was strongly related to the watershed slope (Fig. 2C), which in this system is highly correlated to watershed elevation (r = 0.88). Latitude of the watershed, the percent of the watershed facing north, and interannual variation in winter snowfall all explained additional variation in end-of-snow season dates (Supporting Information Table S3).

The role of snow and geomorphology on stream thermal regimes
Average summer stream temperatures were strongly influenced by stream catchment geomorphology (slope and the percent of the watershed facing north), summer climate conditions (air temperature and cumulative precipitation), and by winter snowpack (Fig. 3A, Supporting Information Table S1). As expected, stream temperatures were strongly related to air temperature and summer precipitation reduced temperatures, likely by increasing stream discharge and thus thermal capacity. Overall, streams draining steep watersheds were cooler and this effect was more pronounced in high snow years. Surprisingly, streams draining more southerly watersheds had cooler temperatures, which may reflect increased snowmelt contributions to streamflow during periods of high solar radiation.
Across the top performing models, these geomorphic and climate variables explained 53% of the observed variation in average daily stream temperatures. Conditionally (including random stream intercepts) the top models explained about 96% of the variation in stream temperature, which suggests strong individual watershed control on stream temperature that is beyond watershed slope, aspect, or interannual variations such as changes in snow. air temperature may also be influenced by channel sinuosity, the presence of small lakes, shading, and subsurface flows or other geomorphic attributes that we have not explicitly accounted for here (Poole and Berman 2001;Caissie 2006). During high snow conditions, streams draining steep watersheds averaged 5 C, and flat watersheds averaged 10 C creating significant heterogeneity in average summer stream temperatures across the river basin (Fig. 3A). In contrast, under low snow conditions, average summer stream temperatures were more homogenous across the basin. Streams draining steep watersheds warmed significantly to 8 C from high snow to low snow conditions whereas flat streams showed only modest warming to 10 C. These changes combined to reduce the thermal variation present on the landscape by more than 20% from high snow to low snow years (Fig. 4).
Similar to average stream temperatures, τ was strongly influenced by stream catchment geomorphology (slope and the percent of the watershed facing north), summer climate conditions (cumulative precipitation), and by winter  Table S2). Together these variables only explained 35% of the variation in temperature sensitivity across streams, but the strength of the effects has important consequences for stream thermal conditions. There was a negative association between steam watershed slope and τ under all snow conditions, but this effect was weakened under reduced snow coverage (Fig. 3B), thus reducing the variation in τ across streams within the basin. Streams draining flat watersheds were equally correlated with air temperature across variable snow conditions, yet streams draining steep watersheds went from being essentially uncorrelated with variation in air temperature in high snow years to having strong positive association with air temperature in low snow years (Fig. 3B). A major implication of these changes in τ is that streams with high values of τ had stronger day-to-day variability in stream temperature, indicative of having lost the significant buffer from snowmelt (Fig. 1). Thus, loss of snow cover increases both average stream temperatures and daily temperature variability in streams draining steep watersheds.

Thermal heterogeneity across landscapes support ecosystem function
An important aspect of climate impacts on aquatic ecosystems is how changing climate affects the spatial and temporal expression of thermal conditions across landscapes. Our results show that loss of snow results in the homogenization of thermal conditions across the landscape, likely leading to important changes in ecological functions and ecosystem services supported by this thermal heterogeneity (Brennan et al. 2019). For example, thermal heterogeneity in the watersheds interannual variation in snowpack. Thermal heterogeneity was measured by computing the coefficient of variation (CV) across the study streams for each day of summer and then averaged from July 15 th to August 15 th . The mean snow index is the average residuals from the landscape snow season linear model ( Fig. 2 and Supporting Information Table S3). The black line indicates the least squares line of best fit (R 2 = 0.79, p < 0.05). in this study is critical for the stability of sockeye salmon (O. nerka) populations and for consumers who rely on salmon subsidies for the majority of their annual growth . Climate-induced changes in spatial variation in thermal heterogeneity will have direct implications for this expression of life-history diversity in salmon, which we know is an important dimension of their biology that makes them particularly sustainable to exploitation by both fisheries and wildlife.
Water temperatures in streams draining steep watersheds were substantially warmer in low snow years than high snow years. While the overall water temperatures in these watersheds were still colder than in flat watersheds, the observed changes in average summer temperature of 3-5 C are ecologically significant for aquatic organisms whose life histories are adapted to cold-water temperatures. These changes also outpace current projections of changes in air temperature and are biologically relevant for spawn timing, summer degree-day accumulation for fishes. Also, populations of salmon adapted to cold conditions lack genetic diversity in the major histocompatibility complex (a gene complex related to immunological function); therefore, they may be more sensitive to moderate changes in temperature and associated effects on disease prevalence (Larson et al. 2016). Ultimately, future change and variation in snowfall may favor organisms that are more flexible in their reliance on cold snow-dominated streams.
Topography, climate change, and the future of snow High latitude and montane ecosystems are experiencing substantial changes in temperature and precipitation regimes with ongoing climate change (Fyfe et al. 2017). We have demonstrated a strong mechanism through which watershed topography interacts with climate conditions to regulate summer stream water temperatures. Across the range of conditions observed in this study, streams draining steep watersheds were cooler and less sensitive to variation in air temperature compared to streams draining flat watersheds, but carryover effects from winter snow conditions control the strength of these relationships. Additionally, our empirical approach to understanding the physical controls on stream temperatures reveals a latent trend of increasing temperature throughout the summer, that likely represents accumulated heat captured by the landscape, a response thought only to be captured by process-based models of stream heat budgets (Leach and Moore 2019). This stored heat most strongly affects streams with large contributions of shallow groundwater or precipitation that slowly percolates through rock and soil in the watershed (steep snow-dominated watersheds in our system), thus contributing to disproportionate effects of climate warming on cold-water systems.
Thermal regimes in snow-dominated systems may show rapid change as small shifts in winter temperature can substantially affect the accumulation of snow (Littell et al. 2018). Snow affects stream thermal regimes by buffering their sensitivity to day-to-day changes in solar radiation, thus regulating the temperatures they attain during the hottest part of the year. Snow also provides an important source of discharge increasing the stream volume and reducing the sensitivity to incoming solar radiation. Snowpack (snow water equivalent) is projected to decline by 48% for the western United States (Gergel et al. 2017) and between 17% and 58% for most of Alaska (Littell et al. 2018) by 2100. SNOTEL sites near our study system recorded~50% difference in maximum snow depth and April 1 st snow depth between 2012 (our highest snow year) and 2015 (our lowest snow year). Therefore, our study reflects changes that may occur before the end of the century and as snowpack is lost or melts earlier in the year, these systems will become more vulnerable to extreme heat.

Ramifications for habitat and climate refugia
Previous research suggests that cold-water refugia have slow climate velocities due to their relative insensitivity to air temperature (Null et al. 2013;Luce et al. 2014;Isaak et al. 2015Isaak et al. , 2016. These refugia are buffered from warming by topographical and hydrological controls, yet these processes are often largely supported by snowpack (Null et al. 2013;Luce et al. 2014;Lisi et al. 2015), or groundwater (Jefferson et al. 2008;Kelleher et al. 2012;Null et al. 2013). In our study, loss of snow disproportionately affects cold streams draining steep watersheds, resulting in both warmer temperatures and increased τ in these streams. Thermal conditions in streams whose hydrological regimes are dominated by snowmelt are likely to change quite rapidly with small changes in winter temperature and precipitation, and may not serve as climate refugia into the future. Also, these streams are likely to be more strongly affected by landscape heating which affects shallow groundwater contributions (Leach and Moore 2019). The mechanisms we show here are probably common to many landscapes, but how they play out in specific ecosystems will be determined by the importance of groundwater contributions (Jefferson et al. 2008;Briggs et al. 2018) and the potential for watersheds to continue accumulating snow during winter months as determined by their elevation and latitude (Winski et al. 2017). Our results suggest that across landscapes there are numerous stream-specific watershed controls on temperature beyond the watershed slope and snow mechanisms demonstrated herein. Nonetheless, our results provide a straightforward empirical framework for developing scenarios of stream temperatures across landscapes in a warmer future.
Future changes in the physical controls on ecosystems such as changing temperature and precipitation are likely to have important effects on ecosystems, but these impacts remain difficult to predict. Our results show that, due to the complex interactions between climate, topography, and hydrology, an approach to conservation that focuses on climate refugia may be vulnerable to future change through losses in snowpack. Approaches that spread the risk and are robust to uncertainty about the future influence of climate change on ecosystems may create better outcomes.