Greater vulnerability of snowmelt-fed river thermal regimes to a warming climate

The Snoqualmie River Basin (SRB), located on the west side of the Cascade Mountains, Washington, has three main forks (north, middle, and south) and a drainage size of about 1800 km² (figure 1). The SRB has a Mediterranean climate with precipitation occurring predominantly from October to March. Elevation ranges from 7 m to 2138 m, and the strong elevational temperature gradient controls the phase of winter precipitation, leading to a mix of rain-dominated, rain-snow transient, and snow-dominated zones (Jefferson 2011). In the alpine area, ground soil surface is directly underlain by bedrock units that do not contain significant fracture systems, indicating that rain and snowmelt do not infiltrate deeply but are transmitted to the stream system rapidly (Turney et al 1995, Bethel 2004). Snoqualmie Falls, downstream on the mainstem after the three forks converge, signals a transition from forest toward anthropogenic impacts including residential, commercial, and agricultural land use. The watershed provides habitat to five anadromous Pacific salmon: Chinook (Oncorhynchus tschawytcha), chum (O. keta), coho (O. kisutch), pink (O. gorbuscha), and winter steelhead (O. mykiss) (Steel et al 2016, Lee et al 2020).

Zoom In Zoom Out Reset image size

2.2.1.  In situ data for model calibration and validation

2.2.2. Meteorological forcing

A modified version of DHSVM-RBM, which integrates the physics-based hydrological model DHSVM and a vector-based semi-Lagrangian stream temperature model RBM, was used for modeling hydrology and stream temperature over the SRB. The original model neglected the thermal effect of cool snowmelt runoff, which can introduce nontrivial errors in temperature simulations for snow-influenced mountain watersheds like SRB. To address this issue, DHSVM-RBM was modified to account for snowmelt thermal influence from upstream boundary conditions along longitudinal flow paths by advection across the river network. Given a considerable body of literature detailing the physics, formulations and parameterization of DHSVM and RBM (Wigmosta et al 1994, Yearsley 2009, 2012, Sun et al 2015, Yearsley et al 2019), we describe the model components briefly below.

The hydrological model DHSVM (Wigmosta et al 1994) solves the water and energy balance at every pixel for simulating key surface and subsurface hydrological processes. It contains a two-layer energy balance snow model for simulating snow accumulation and melt, a two-layer canopy model for modeling interception and evapotranspiration, a multi-layer unsaturated soil model for modeling vertical soil water distribution, and a saturated lateral subsurface flow model. Surface and subsurface flow is routed downslope driven by topography until reaching the channel network, where it is then routed using a cascade of linear reservoirs (Wigmosta et al 2002). For the SRB application, DHSVM was implemented at a three-hourly temporal scale and 150 m spatial scale, forced by meteorological variables including precipitation, air temperature, wind speed, relative humidity, downward shortwave and longwave radiation, which were disaggregated using the mountain microclimate simulation model (Hungerford et al 1989) from the daily meteorological data (both historical and GCM downscaled). The channel network represented by connected river segments generated from a digital elevation model is used by RBM to establish the topology for the semi-Lagrangian particle tracking scheme in advection-dominated river systems.

Originally, DHSVM provided to RBM the segment-specific time series of forcing data consisting of air temperature, shortwave and longwave radiation, wind speed and vapor pressure, and inflow and outflow simulations for each river segment throughout the river network. Here, we modified DHSVM to provide to RBM additional time series input of segment-specific snow melt to identify melt influenced locations and the magnitude of melt.

RBM simulates water temperature dynamics by solving the one-dimensional, time-dependent equations for the thermal energy budget:

$$\begin{align} & T(t + \Delta t,x) \nonumber\\ & \quad = T(t,x) \!+\! \left\{ {\frac{{F[T(t,x)] \!+\! h[t,x] \!+\! B[T(t,x)]u(t,x)}}{{\rho {C_{\text{p}}}D(t,x)}}} \right\}\!\Delta t\end{align} \tag{ 1 }$$

where T is water temperature (°C) at time t+ Δt and distance x; Δt is the time increment for a Lagrangian parcel to travel from the upstream starting point to the downstream location x determined by flow velocity u(t, x) (m s⁻¹); D(t, x) is water depth (m); ρ is the density of water (kg m⁻³); C_p is the specific heat capacity of water (J kg⁻¹ °C⁻¹); F[] represents thermal energy budget components that are a function of T, given by:

$$\begin{equation}F = {H_{{\text{con}}}} - {H_{{\text{back}}}} - {H_{{\text{evap}}}}\end{equation} \tag{ 2 }$$

where H_con is thermal energy flux due to conduction, H_back is the flux due to longwave radiation from the water surface, and H_evap is the flux from the latent heat loss due to evaporation; h[] represents the known inputs to the thermal energy balance that are not a function of T, given by:

$$\begin{equation}h = {H_{{\text{sw}}}} + {H_{{\text{lw}}}}\end{equation} \tag{ 3 }$$

where H_sw is the net shortwave solar radiation reaching the water surface after attenuation from topography and riparian vegetation, and H_lw is the downward longwave radiation from the atmosphere; and B[] represents advected sources of thermal energy from tributaries and lateral surface/subsurface flows, given by:

$$\begin{equation}B = {S_{{\text{trib}}}}[{T_{{\text{trib}}}} - T] + {S_{{\text{lat}}}}[{T_{{\text{lat}}}} - T]\end{equation} \tag{ 4 }$$

where S_trib and S_lat is the ratio of the volume of tributary and lateral inflows to the segment volume, respectively; T_trib and T_lat is the water temperature of the tributary and lateral inflows, respectively. With the model modification, when the snowpack melted (i.e. DHSVM simulated local melt for a river segment exceeds a user-specified depth threshold), T_lat were assumed to be snowmelt runoff with a constant temperature of 0.1 °C (Kobayashi 1985, Ficklin et al 2012). Otherwise, the mean annual air temperature was used as a proxy for T_lat.

Estimation of initial conditions (headwater temperatures) was modified to use snowmelt temperature of 0.1 °C during snowmelt, and the nonlinear regression model of Mohseni et al (1998) was used in the absence of melt:

$$\begin{equation}{T_{{\text{head}}}} = \mu + \frac{{\alpha - \mu }}{{1 + {{\text{e}}^{\gamma (\beta - {T_{{\text{smooth}}}})}}}}\end{equation} \tag{ 5 }$$

where T_head is daily headwaters temperature (°C), T_smooth is smoothed air temperature (°C). Mohseni parameters α, β, μ, γ were estimated and tuned based on multi-year observations of air temperature and water temperature records from K1–K5 stations. With no available observations of headwater temperatures in the upper SRB, we used headwater measurements from a nearby watershed (USGS gauge ID: 12116500).

Segment-specific hydraulic characteristics (i.e. U and D) were estimated from DHSVM streamflows using the power law relationships from Leopold and Maddock (1953):

$$\begin{equation*}D(t,x) = {d_a}Q{(t,x)^{{d_b}}}\end{equation*}$$

$$\begin{equation}U(t,x) = {u_a}Q{(t,x)^{{u_b}}}\end{equation} \tag{ 6 }$$

where Q is segment streamflow, and the parameters, d_a, d_b, u_a and u_b were initialized based on stage-discharge relationships derived from five USGS sites (G1−G5 in figure 1). Both Leopold and Mohseni parameters estimated at the individual gauges were extended to their nearby stream segments with no local observations. Sources of key input data to DHSVM-RBM are presented in table 1.

We used simulated water temperatures from DHSVM-RBM to demonstrate how growth potential for native coldwater Chinook salmon, coho salmon, steelhead, pink salmon, and one nonnative warmwater species, largemouth bass, may change in a future climate. For each species, we estimated growth potential for a 1 g fish in each reach accessible to fish twice daily at 6am and 6pm using the Wisconsin bioenergetics model (Hanson et al 1997). Growth is affected by fish mass, water temperature, food availability, and the energy density of predators and prey. In our scenarios, only water temperature varied. We calculated potential food availability in each reach as the product of stream order and distance upstream, which ranged between a maintenance ration and a maximum consumption ration. We used standard bioenergetics modeling options (Hanson et al 1997), and default values for most parameters (table S1 (available online at stacks.iop.org/ERL/16/054006/mmedia)). We used updated temperature-dependent consumption parameters for Chinook salmon provided by Plumb and Moffitt (2015) because these better represent PNW juveniles and we increased the prey energy density parameter from the default value of 2500 J g⁻¹ to 3000 J g⁻¹ to better match salmon prey in the wild (Beauchamp et al 2007). We summarized seasonal growth potential for each species, focusing on longitudinal patterns in the mainstem from the river outlet to Snoqualmie Falls. We did not consider potential fish movement or other ecological factors that influence growth in the wild. Here, our aim was not to predict actual fish growth, but to illustrate relative changes in growth potential across seasons throughout the mainstem that reflect thermal changes associated with different climate scenarios.

DHSVM-RBM was run from 1 October 1990 to 31 December 2013 with the first two years as a spin-up period. Streamflow in the first six water years (1992–1998) at station G1 were used for model calibration, and the model was validated for the remaining 15 water years (1998–2013) at G1 and all streamflow records collected at G2–G5 from 1992 to 2013. Given shorter durations of stream temperature measurements, full records at the two sites (T1, T2) were used for calibration. Details of the calibration procedure are provided in the supplementary information.

Figure 2 presents DHSVM-RBM calibration and validation results for selected snow, streamflow, and water temperature sites. Full calibration and validation results for all sites are shown in figures S1–S3. For calibration, the daily SWE Kling–Gupta efficiency (KGE, equation S1) was 0.75 and 0.72 for S1 and S2, respectively; the daily streamflow KGE was 0.73 for G1. Model performance in the validation periods was slightly degraded compared to calibration. The daily SWE KGE was 0.93 for S1 and 0.67 for S2; the daily streamflow KGE ranged from 0.60 to 0.65 for G1–G5. Simulated daily water temperatures for the five tributary sites (K1–K5) had Pearson's correlation coefficients (r) of 0.94–0.96 and mean absolute error (MAE) of 0.92 °C–1.43 °C; the r for the two mainstem sites (T1 and T2) was 0.91 and 0.84, and MAE was 0.97 °C and 1.35 °C, respectively. The results suggest that DHSVM-RBM-simulated snow, hydrologic and stream temperature predictions match in situ observations fairly well across the basin.

Zoom In Zoom Out Reset image size

Figure 3 compared stream temperature simulations during 2011–2013 with and without considering snowmelt effects for the N1 site which is in the alpine area with seasonal snow cover. Results suggest that during the snowmelt season from March through August, simulations that did not include snowmelt effects showed substantial overestimation of observed stream temperatures. Accounting for snowmelt effects helped to correct these biases: the MAE of simulated daily stream temperature and the absolute error of 7-DADMax decreased from 2.4 °C and 3.4 °C without snowmelt effects to 2.0 °C and 0.1 °C with snowmelt effects considered, respectively. Although the overall streamflow in SRB has limited contributions from deep groundwater (Leibowitz et al 2016), field surveys confirmed that some locations in the river were potentially influenced by localized hyporheic exchange and/or groundwater inputs (King County 2016, McGill et al 2021). Not representing thermal input from hyporheic exchange or deep groundwater in confined aquifers in DHSVM-RBM may explain the generally overestimated diurnal range of river temperatures, and the underestimated winter stream temperatures.

Zoom In Zoom Out Reset image size

Figure 4 presents the basin-wide monthly mean stream temperature changes between two simulations (with snowmelt minus without snowmelt) from March to August over 2001–2013. As observed in figures 4(a)–(f), beginning in March snow in the low-to-moderate elevations starts to melt. In early snowmelt periods, cool snowmelt water has a major influence on temperatures of headwaters and lower-order tributaries. As air temperatures continue to rise, the spatial extent and rates of melt increase, resulting in amplified effects of cool snowmelt water extending from tributaries to the mainstem from advection. By June, the model suggests that snowmelt from upstream rivers reduces the basin outlet stream temperature by about 0.5 °C. Towards the end of the melt season (July to August), snowmelt cooling effects start to recede and occur mainly in the highest elevations (i.e. reducing Middle Fork temperatures by 2.6 °C to 3.4 °C).

Zoom In Zoom Out Reset image size

Comparisons were also made in simulated basin-wide 7-DADMax with and without considering snowmelt effects, which were summarized by five elevation bands of the SRB shown in figure 1(c) (figures 4(g) and (h)). Results suggest that snowmelt generally has greater cooling effects on the 7-DADMax of up to 5.8 °C for streams at higher elevations such as the Middle Fork drainage, and a cooling effect of >0.5 °C extends downstream until near the confluence of the Tolt River. For example, the differences of 7-DADMax are −1.0 °C at the elevation band >717 m, as compared to −0.1 °C at the elevation band of <85 m (most of which are downstream of the lower Snoqualmie-Tolt River confluence). At elevations higher than 717 m, with snowmelt influence considered, about 33% and 26% of stream segments showed a decrease in 7-DADMax greater than 1.0 °C and 2.0 °C, respectively. This spatial variability of snowmelt influence on 7-DADMax, which occurred during July and August in all river segments, is related to the timing and rates of snowmelt from both the local reach and its upstream contributing reaches. High elevation areas such as the Middle Fork subbasin are characterized by higher snow accumulations and later onset of above-freezing air temperatures. This topographical influence resulted in greater snowmelt late into the summer that was able to reduce 7-DADMax locally and downstream as melt runoff traveled along its flow path.

For consistency with the historical stream temperature simulations, we also used a 13 year period to estimate the potential changes of stream temperature under a warming climate. We ran DHSVM-RBM using ten GCM ensembles for a historical period (1993–2005) and a future period at the end of century (2087–2099).

Figure 5 presents the changes of monthly precipitation and air temperature between historical and future projections based on ten GCMs averaged over the basin. It is projected that the mean annual temperature over the basin will increase by 4 °C–7 °C by the end of the century with the largest increases occurring during the summer months. Most GCMs project increased precipitation of up to 58% during the winter and early spring and substantially decreased precipitation of up to 70% during the summer months. Similar to previous studies (Tohver et al 2014, Mote et al 2018, Yan et al 2019b, Lee et al 2020), the model suggests that the SRB shifts from a snow-influenced basin to a rain-dominated basin in future climates. Over the basin, the mean snow accumulation decreases substantially (i.e. 92% reduction in peak SWE) as a result of warming, which prompts earlier onset of snowmelt by about 2 months and earlier disappearance of snow cover (figure 6). Alterations in snowpack along with changing climate result in higher streamflows in winter and early spring, and significantly lower streamflows from late spring through fall (figure S4). For example, summer streamflow in the high-altitude Middle Fork is projected to decrease by 52%–82%. Overall, river segments at higher elevations, where snow regimes exert a greater influence on streamflow, exhibit greater streamflow seasonal variability in response to changing climates (figure S5).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 7 shows the changes in 7-DADMax, comparing future and historical periods for each stream segment across the basin. Consistent with Steel et al (2019), simulations suggest that 7-DADMax will not uniformly increase in the future over the SRB basin; increases are projected to be most prevalent in the snow-dominated high elevation areas (e.g. Middle Fork) and the downstream mainstem until the convergence with the Tolt River. This is due to the compound effect of increasing heatwave events both in duration and intensity, diminishing snowpack and thus reduced snowmelt cooling effects in summer, as well as diminished streamflow from decreasing precipitation and snowpack. It should be noted that our findings are not consistent with previous studies, that analyzed multi-year stream temperature observations for interannual correlations with air temperature and found that cold streams are less sensitive to air temperature variation mainly due to snowmelt input (Isaak et al 2010, Isaak and Rieman 2013, Luce et al 2014). The changing sensitivity of snow-dominated streams to a warming climate reflects a 'tipping point' that the long-term model projection captures but short data records may not, especially during extreme hydrometeorological events. For example, the recent 2015 drought in Washington State provided a useful harbinger of the future (Marlier et al 2017). In 2015, Washington State experienced extreme drought attributed to high winter air temperature that resulted in exceptionally low winter snowpack in the Cascades Range. Steel et al (2019) compared 2015 stream temperature observations to more typical conditions (2012–2014) throughout the Snoqualmie River and found that the largest changes for summer temperature occurred in the snowmelt-dominated Middle Fork.

Zoom In Zoom Out Reset image size

Based on the Washington State Legislature spawning rearing and migration criteria (figure S6), we found that river temperature exceeding the state standards dominated in the summer season, and about 65% of the river segments (70% of anadromous-accessible reaches) over the basin fail to meet the water quality standard in the historical summer period and this number will increase to 92% by the end of the 21st century averaged over the ten GCMs; the basin average number of days exceeding the state standards annually will increase substantially from 77 (86) to 135 (143) d (figure S7). Our bioenergetics modeling suggested that warmer future stream temperatures could enhance growing opportunities during winter; however, the increasing air temperature coupled with loss of snow-cooled water during late spring and summer will substantially lower growth potential for coldwater salmonids (figure 8). High future summer temperatures may be more detrimental for some species than others. For instance, pink salmon smolts enter marine waters in spring, whereas coho, steelhead, and some runs of Chinook remain in streams during summer. Warmwater species such as largemouth bass are likely to experience minor increases in habitat suitability year-round, which could compound risks to juvenile salmon if warmer temperatures increase largemouth bass's energetic needs and if newly thermally suitable habitat expands their range and encroaches on salmon habitat (Hawkins et al 2020). We noted a longitudinal decrease in summer growth potential progressing downstream for Chinook and steelhead, suggesting the substantial loss of coldwater habitats in summer.

Zoom In Zoom Out Reset image size