(C) PLOS One
This story was originally published by PLOS One and is unaltered.
. . . . . . . . . .



Daily stream temperature predictions for free-flowing streams in the Pacific Northwest, USA [1]

['Jared E. Siegel', 'Ocean Associates Inc.', 'Contracted To Northwest Fisheries Science Center', 'National Marine Fisheries Service', 'National Oceanic', 'Atmospheric Administration', 'Seattle', 'Wa', 'United States Of America', 'Aimee H. Fullerton']

Date: 2023-08

Supporting sustainable lotic ecosystems and thermal habitats requires estimates of stream temperature that are high in scope and resolution across space and time. We combined and enhanced elements of existing stream temperature models to produce a new statistical model to address this need. Contrasting with previous models that estimated coarser metrics such as monthly or seasonal stream temperature or focused on individual watersheds, we modeled daily stream temperature across the entire calendar year for a broad geographic region. This model reflects mechanistic processes using publicly available climate and landscape covariates in a Generalized Additive Model framework. We allowed covariates to interact while accounting for nonlinear relationships between temporal and spatial covariates to better capture seasonal patterns. To represent variation in sensitivity to climate, we used a moving average of antecedent air temperatures over a variable duration linked to area-standardized streamflow. The moving average window size was longer for reaches having snow-dominated hydrology, especially at higher flows, whereas window size was relatively constant and low for reaches having rain-dominated hydrology. Our model’s ability to capture the temporally-variable impact of snowmelt improved its capacity to predict stream temperature across diverse geography for multiple years. We fit the model to stream temperatures from 1993–2013 and predicted daily stream temperatures for ~261,200 free-flowing stream reaches across the Pacific Northwest USA from 1990–2021. Our daily model fit well (RMSE = 1.76; MAE = 1.32°C). Cross-validation suggested that the model produced useful predictions at unsampled locations across diverse landscapes and climate conditions. These stream temperature predictions will be useful to natural resource practitioners for effective conservation planning in lotic ecosystems and for managing species such as Pacific salmon. Our approach is straightforward and can be adapted to new spatial regions, time periods, or scenarios such as the anticipated decline in snowmelt with climate change.

This is an open access article, free of all copyright, and may be freely reproduced, distributed, transmitted, modified, built upon, or otherwise used by anyone for any lawful purpose. The work is made available under the Creative Commons CC0 public domain dedication.

Here, we model and predict daily stream temperature across the entire calendar year for a broad geographic region, contrasting with previous models that estimated coarser metrics such as monthly or seasonal stream temperature [ 5 , 21 , 38 ] or focused on individual watersheds [ 25 ].Specifically, we improve on previous approaches by incorporating the antecedent air temperature covariate introduced by Siegel et al. [ 37 ] (referred therein as “flow-dependent flexible-lag air temperature”) into the spatiotemporal modeling approach of Siegel and Volk [ 25 ] using publicly available spatially and temporally continuous climate covariates in conjunction with spatially continuous landscape covariates. We present a daily stream temperature model applicable across the entire calendar year for the geographically diverse PNW, a region critical for cold-water salmonid production. We fit the model to stream temperature data from 1993–2013 and used the model to predict daily stream temperatures for all regional stream reaches from 1990–2021. Our results illustrate spatial and temporal variation in stream temperature that can be used directly in management applications now and that may also generate insights about how stream temperatures will be affected by climate change.

Siegel et al. [ 37 ] represented temporal variation in the sensitivity of stream temperature to air temperature by developing a covariate of antecedent air temperature averaged over a variable duration that was determined by flow conditions. They used this covariate (and others) to estimate stream temperatures at ~400 free-flowing sites (i.e., uninfluenced by impoundment) throughout the Pacific Northwest (PNW), USA. They found evidence that snowpack buffered stream temperature long after snow had melted, presumably due to its movement into groundwater that later flowed into streams. Thus, sites with higher snowpack had longer thermal lags than sites uninfluenced by snow, making the snow-influenced sites less responsive to short-term changes in air temperature. Asarian et al. [ 38 ] used a similar approach that included time-varying and nonlinear effects of flow to predict stream temperature for a snowmelt-driven watershed in California. They also found that stream temperatures remained lower well into summer in high flow years. Both studies fit models to single sites, thus their approach did not capture interactions between climatic and landscape influences on stream temperature.

Many statistical models use air temperature as a proxy to represent physical energy exchanges driven by solar radiation and thermal conduction (e.g., [ 33 ]). However, as with the example of groundwater above, the relationship between air temperature and stream temperature can change across seasons and years as it is mediated by snowmelt and river discharge [ 26 , 34 , 35 ]. For example, snow-influenced streams may be less sensitive to air temperatures due to the direct cooling influence of snowpack melt and associated higher flows creating increased thermal inertia [ 36 ]. Stream temperature becomes more correlated with air temperature in low snowpack years and later in the season as the snowpack diminishes and associated flows decline.

Previous statistical stream temperature models have generally assumed non-interacting climate (e.g., air temperature, precipitation) and landscape covariates (e.g., elevation, drainage area, slope) [ 11 , 21 , 30 , 31 ]. In a recent advancement, Siegel and Volk [ 25 ] allowed time-varying climate covariates (e.g., air temperature) and spatially-explicit landscape covariates (e.g., elevation) to interact, thereby producing a model that predicted stream temperatures based on putative mechanistic principles (defined a priori), rather than allowing a covariance function to absorb error that was otherwise unexplained. They argued that these interactions were necessary for models to fit across seasonal and interannual climate variability as the influence of many variables is not unidirectional. Instead, most spatial variables either mitigate or amplify seasonal extremes in stream temperatures, and thus the directional influence of variables changes between cold and warm seasons. For example, groundwater influence can mitigate temperature extremes, cooling river water in the summer but warming it in the winter. Spatiotemporal thermal variability across the entire calendar year was well approximated by these simple interaction terms between spatial landscape covariates and temporal climatic covariates [ 24 , 25 , 32 ].

Statistical models can better capture dynamic and changing conditions if parameterized covariate relationships are able to adequately represent seasonally shifting and interacting processes that influence stream temperature [ 25 ]. Numerous studies have compared the performance of model types against each other, offering varying conclusions regarding the best modeling approach [ 13 , 26 – 29 ]. However, model performance has as much to do with the ability of covariates to capture variation in drivers as it does with the model type and statistical methods. While autocorrelation in stream temperature is often discussed as a separate and abstract entity, it is largely a result of the high heat capacity of water causing resilience in the propensity of water temperature to change as it moves through stream networks. This resilience is likely to vary temporally with variation in flow and the influence of snowpack and spatially as a consequence of channel characteristics. We argue that there is substantial room for improvement in statistical models for stream temperature by carefully choosing and parameterizing covariates to better reflect these underlying processes. Modeling efforts are likely to benefit from pursuing the maximization of the ability of covariates before depending on more rigid autocorrelation functions.

In contrast, regressions and other statistical models typically rely on fewer data, simpler relationships, and require less expertise to develop and run [ 16 – 18 ]. These models can, however, suffer from spatial autocorrelation problems that occur when trying to predict stream temperature across a directed-flow connected network, stemming from the underlying processes of water and energy flow [ 11 , 19 ]. Spatial statistical network (SSN) models were developed to address this concern. They use autocovariance functions that are assumed to be stationary across the spatial domain over which models are fit [ 12 , 20 ]. To date, SSN models applications have generally limited predictions to coarse temporal resolution like seasonal or monthly means [ 5 , 21 ] or to small spatial extents [ 12 , 22 ]. More recently, several spatio-temporal statistical stream temperature models have been developed that account for both spatial and temporal autocorrelation [ 12 , 23 , 24 ]. However, because SSN and spatio-temporal stream temperature models depend on substantial data to inform their spatial autocovariance functions or temporal autoregressions, they can require high computational power and they are limited in their ability to predict outside the spatial or temporal bounds of the data used to fit the model. Therefore, these models may not predict well in scenarios where the correlation structures are poorly informed or likely to shift, e.g., with landscape alteration or climate change. There is, therefore, a strong, immediate need to enhance statistical models to better predict stream temperature under dynamic climate conditions.

Climate interacts with landscape spatial heterogeneity, stream flow rates, and stream network connectivity to produce complex spatiotemporal patterns in stream temperature. This makes stream temperature prediction a multi-faceted, dynamic, space-time problem [ 11 , 12 ]. Different modeling approaches have been developed to address this need, ranging from practical models that are not overly complex (e.g., regression models) to complex models that are more difficult to apply yet better represent complex physical processes (e.g., process-based models or models that merge process-based and statistical methods). Complex process-based stream temperature models often track thermal energy budgets through a riverscape over time and aim to reflect mechanistic links between environmental variables and water temperature [ 13 , 14 ]. However, they tend to be data-hungry, time-consuming to develop and run, and are typically applied over limited spatial extents such as the Sacramento River [ 15 ].

The critical role water temperature plays in governing aquatic biological processes and reflecting fluvial system health and function drives our desire to understand riverine thermal regimes and has motivated the development of stream temperature models for decades [ 1 ]. In recent years the rate of model development has increased substantially [ 2 ]. While water temperatures are increasing and expected to continue to increase across the globe with climate change, empirical observations remain relatively sparse, making stream temperature modeling a pressing need for the water resources community [ 3 , 4 ]. However, riverine thermal regimes will change at different rates in different places, as governed by local climate and landscape characteristics [ 5 ]. Having accurate predictions of water temperature over space and time under diverse climate conditions will be important to inform conservation practitioners’ work, especially for cold-water aquatic species such as salmonids. Because many cold-water species are already living near their maximum thermal tolerances, they are likely to be particularly challenged by changing thermal regimes and thus reflect pressing land and riverscape management needs [ 6 – 10 ].

Finally, we evaluated if characteristics of data used to fit models (i.e., reaches and days with empirical data) were spatially and temporally representative of those across the entire PNW where we made predictions. Because of the large number of covariates used in our model, for this dataset evaluation we compressed the multidimensional covariate space of both fitting and prediction data using principal components analysis [ 50 ] and evaluated spatial and temporal covariates separately. All covariates were standardized by relative rank. We compared the spatial covariate space for all reaches in the fitting dataset (n = 10,032) to that of the entire PNW region we predicted (n = 261,207). For the immense temporal dataset (i.e., daily data at each reach), we analyzed a random 10% subset of reaches and aggregated data by month, resulting in ~350,000 unique observations of reach/month combinations. Multiple random subsets were tested and gave similar results to those presented here.

We validated models using root mean squared error (RMSE) and mean absolute error (MAE) statistics. In addition to providing fit statistics, we performed a leave-one-year-out cross-validation and a leave-one-region-out cross-validation, hereafter referred to as the temporal and spatial cross-validations, respectively. Although cross-validation methods often leave out a random proportion of the dataset, our temporal and spatial cross-validation tests evaluated the model’s ability to predict entirely distinct years and regions left out when fitting the model. For the temporal cross-validation, a single year was withheld from the model fitting dataset and used as a validation dataset, whereas for the spatial cross-validation, a region was withheld. The validation processes were repeated until each year/region had been withheld and predicted. For cross-validations, we report root mean squared prediction error (RMSPE) and mean absolute prediction error (MAPE). In addition, we examined how model performance varied across space and time by fitting the relationship between the absolute value of residuals and all covariates using GAM smoothers (knots limited to 5). We evaluated statistics for the daily model predictions and daily predictions summarized as monthly predictions.

To assess the ability of our methods to produce accurate predictions at high spatiotemporal resolutions, we compared the accuracy and precision of predictions from a single model (i.e., a model fit to all available data) and models fit to seasonal and spatial subsets of the data. Prior investigations demonstrated that prediction accuracy improved if models were fit to the warming and cooling seasons separately [ 25 , 30 ]. We therefore split the dataset by the day of year that represented the average peak in stream temperatures in our study (Julian day 210, ~July 29th). For regional models we split the dataset into 10 large sub-regional watersheds within the PNW based on the “processing units” designated in the NorWeST dataset [ 5 ]. We then examined the effect of season/region specific models (20 separate fits).

Spatial landscape covariates that had no temporal component were modeled with linear effects, as they were expected to have simpler relationships with stream temperature, either amplifying or muting the impact of climate across space. Limiting the number of smoothers simplified the model and lowered the risk of overfitting. We allowed each of these variables to interact with antecedent air temperature (Tl) to help account for seasonal variation in their directional influence on stream temperature. For example, covariates that cool stream temperature in the summer when air temperatures are warm, such as groundwater influence (as represented by baseflow index, BFI) or riparian forest cover (F), tend to warm stream temperature during the winter when air temperatures are cold. We considered only two other interactions between linear covariates: annual April 1 st snowpack depth (SA1) and day of year (D) to help account for the delayed effect of snow melt through groundwater connectivity; and stream reach elevation (E) and watershed area (A) to help distinguish high- from low-elevation headwater streams. The final model formula was:

Our modeling approach centered around our primary air temperature covariate (antecedent air temperature Tl), which we allowed to have 6 knots to capture the logistic-shaped relationship with stream temperature [ 33 ]. Most of the other variables were allowed to interact with Tl to parameterize their theoretical impact of either mitigating or exacerbating climate impacts on stream temperature. However, daily watershed averaged air temperature (Tws), which was allowed 3 knots, helped account for the large influence of air temperature on a given day. Tws was considered a better choice than Tl to interact with snowpack (Sws) as it better represented air temperatures encountered by mountain snowpack. For other temporally continuous climate covariates (i.e., snowpack Sws, area-standardized flows Qr, and daylight hours DL), we allowed 4–6 knots and included interactions with antecedent air temperature (Tl).

We parameterized the effects of climate covariates using a combination of GAM tensor product smoothers and interactions fit with the “ti” function in the R package mgcv [ 42 ]. Although the default ti fitting function using thin-plate regression splines has a penalty term to prevent overfitting, our examination of fitted relationships and cross-validated predictions found it to be too permissive for our purposes. To avoid overfitting, we restricted the number of knots within the climate variable smoothers to limit the complexity of fitted relationships while allowing enough flexibility to capture nonlinear relationships based on theory [ 25 , 31 ]. We used the default fitting function (fREML) to fit models.

We evaluated the correlation between standardized covariates across seasons using a Pearson’s correlation test ( S1 Fig ). Other than D and DL, which were correlated in all seasons except winter as expected, there were only three pairs of highly correlated covariates (> 0.8): Tws and Tl during fall (0.9), winter and spring (both 0.8), R and P during summer (-0.8), and SA1 and Sws during winter (0.9). SA1 and Sws were not significantly correlated during the other seasons, so we included both covariates. Although correlated, Tws and Tl were both included in the model due to distinct hypothesized impacts on stream temperature. An interaction term to help the model distinguish the individual impacts was included and found to be important in informing the model. However, because there is seasonal collinearity between some of the variables and non-linear concurvity [ 49 ] between some of the estimated smoothed parameters, the modeled relationships for these variables predicted by the model should be interpreted cautiously. Still, predictive power was improved by including all variables, and thus some remaining collinearity was considered acceptable as the primary goal of the effort was to produce accurate predictions across time and space.

The air temperature data and SWE data are available as continuous 4-km (16-km 2 ) grid cells whereas the NWM flow predictions are available for each reach in the NHDPlus version 2 stream network. To produce distinct local- and watershed-scale air temperature and snowpack covariates, we averaged gridded data over the local catchment area draining directly into a NHDPlus version 2 stream reach and over the entire upstream watershed area draining to a reach. For our local air temperature covariate, we used the antecedent air temperature covariate (Tl) described above. A daily watershed-area averaged air temperature metric (Tws) was also included to better capture the impact of air temperatures upstream on snowpack melt compared to the local values. In addition to daily SWE values, we used an annual metric of April 1st SWE interacting with day of year in an attempt to capture the delayed cooling impact of snowpack following melting. We also used daylight hours (DL) which vary across the year, and reach latitude which defines solar angle irrespective of day length.

Our sources for climate data were: daily mean air temperatures from the PRISM Climate Group at Oregon State University ([ 45 ], http://prism.oregonstate.edu ); daily modeled snowpack data represented as snow water equivalent (SWE) provided by the National Snow and Ice Data Center [ 46 ]; and daily flow predicted by the National Water Model Version 2 [ 47 ]. Siegel et al. [ 37 ] performed comparisons between the National Water Model predictions and corresponding measured flow data from free-flowing streams in the region. They found that the model generally captured hydrologic patterns and flow magnitudes, but that its ability to capture the magnitude and timing of peak flows was less consistent.

We selected commonly used model covariates that represent well-known drivers of stream temperature; see [ 43 ] for a summary of energy flux processes. Most of these covariates have previously been shown to be useful in other stream temperature modeling efforts and aligned with our pre-hypothesized effects. We performed iterative backwards stepwise variable selection and retained only variables that aligned with hypothesized effects and that significantly contributed to model performance (p < 0.05). We used the R package visreg [ 44 ] to visualize conditional effect relationships for each potential variable. As found by Siegel and Volk [ 25 ], we found that information criteria such as AIC or BIC were not adequate for model selection due to the large size of our datasets because they tended to support inclusion of covariates that did not align with hypothesized effects.

Next, we developed a metric representing the best moving average window size for each reach and day. To do this, we calculated daily log-scale area-standardized flows (Qr) by subtracting estimated average annual flows based on a reach’s watershed area (A) from estimated daily mean log-scale flows from the National Water Model (see below) using the relationship developed by Siegel et al. [ 37 ]: Qr = -4.10 + 0.93*A. We then split each of the snow-dominated, transitional, and rain-dominated data subsets into 200 groups separated by 0.5% quantiles of Qr. We correlated stream temperature in each Qr quantile with local air temperature averaged over antecedent moving windows of different sizes (1, 3, 6, 9, 12, 15, 20, 25, 30, 40, 50, and 60 days including the day of prediction). To identify the optimal window size (number of days) over which air temperature would be averaged for each hydrographic classification, we fit Generalized Additive Model (GAMS, [ 41 ]) smoothers with <6 knots using the R package mgcv [ 42 ] for rain, transitional, and snow data independently, where the response variable was the window size with the highest correlation from the previous step, and the explanatory variable was mean Qr in each quantile. Finally, we used these models to predict the optimal window size for each daily value of Qr for all reaches (rounded to the nearest tested window size). The antecedent air temperature covariate (Tl) covariate was then created by selecting the air temperature associated with the predicted window size.

Siegel et al. [ 37 ] found that water temperature in snowpack-influenced streams demonstrated slower responses to changes in air temperature than rain-dominated streams, even after the snowpack had melted. Specifically, they described that the average window size over which antecedent air temperature was most highly correlated with stream temperature tended to be larger in snowpack-fed streams. Based on their results, we categorized each stream reach by the average snow accumulated within its drainage area (Sws) for each calendar year and defined each reach and year combination as having “rain-dominated” (<20 mm), “transitional” (20–100 mm), or “snow-dominated” (>100 mm) hydrology. Although technically these represent low, moderate, and high snow classes, we retain the terminology defined above to be consistent with previous related research. The average annual snow-water equivalent (SWE) was calculated by averaging all daily values across the year, including days with no snowpack, and thus the metric reflects the duration of snowpack as well as the magnitude of accumulation. In contrast to Siegel et al. [ 37 ], we used annual SWE accumulation rather than the multi-year average, which allowed each reach to transition between classes among years to better capture interannual variability. These classifications were used only to develop the antecedent air temperature covariate Tl; data were not subset by these classes when fitting the stream temperature model.

To reflect variation in the sensitivity of stream temperature to changes in air temperature, we used a covariate of air temperature averaged over an antecedent period, the length of which varied in conjunction with flow conditions. The rationale for including this covariate was to account for temporal variation in thermal inertia caused by dynamic processes that alter the rate at which streams warm (e.g., precipitation stored as snow and subsequent snowmelt-driven increases in flow) and processes that add direct inputs of cooler water that suppress stream temperatures in comparison to air temperatures (e.g., groundwater and snowmelt). Our antecedent air temperature covariate (Tl) was adapted from the approach described in Siegel et al. [ 37 ], as outlined below.

We used ~20 years of daily stream temperatures from 1993 through 2011–2013 (with the end year varying by region) that were contributed by various organizations to the NorWeST project; raw data were curated, quality controlled, summarized, and made publicly available by the USDA Forest Service [ 5 ]. We downloaded “All Days Daily Summary” files within the PNW study region, which included individual stream temperature measurements typically taken at 15–30 min intervals that had been summarized into daily metrics ( https://www.fs.fed.us/rm/boise/AWAE/projects/NorWeST/StreamTemperatureDataSummaries.shtml ). We used ~5.16 million daily stream temperature estimates for free-flowing reaches (87% of ~5.91 million available, after excluding data from reaches heavily influenced by dams) ( Fig 1 ).

For our spatial framework we used the National Hydrography Dataset 1:100k (NHDPlus version 2) network of stream reaches created by the US Geological Survey [ 39 ]. The PNW portion of the dataset (Region 17) is primarily composed of Columbia River Basin streams in the USA, but also includes coastal areas of Oregon and Washington and Puget Sound drainages. It includes 270,941 stream reaches of 1.98 km in length on average (SD = 1.74 km), >96% of which were defined as “free-flowing”. As in Siegel et al. [ 37 ], we classified free-flowing streams as those upstream of the influence of major dams, i.e., those with capacity >0.1 km 3 from the Global Reservoir and Dam Database (GRanD v1.3; [ 40 ]) and those with height >15.24 m (50 ft) from the National Inventory of Dams ( https://nid.usace.army.mil ). In contrast to Siegel et al. [ 37 ], we also included streams with minimal dam influence in our definition of free-flowing reaches. To calculate dam influence, we calculated the percentage of the upstream drainage area (PDA) for each reach that lies upstream of a major dam. For fitting models, we implemented a conservative cutoff, only including reaches where PDA < 5%. We used the resulting model to predict stream temperatures for reaches with a PDA < 25%. When more than 25% of the flow in a reach comes through a large dam, the temperature in that reach is frequently uncoupled from the physical processes we modeled (e.g., when water is released from below the thermocline in deep reservoirs, or when warm surface water is spilled over dams at unnatural times). The free-flowing streams in our dataset comprised small to large streams (stream order mean 1.8, max 9), represented a large range of elevations (mean 1200 m, range 0–3362 m), and included headwater streams with minimal drainage areas to large rivers with drainage areas up to 20,000 km 2 .

Some of the more influential spatial covariates lacking a temporal component included reach elevation (E), the difference between the reach elevation and the average watershed elevation (Ed), watershed area (A), a representation of base flow (BFI), and the percent riparian forest cover in the 100-m buffer (R). The relationship between Tl and E suggests that streams at higher elevations are generally cooler irrespective of air temperature, which is accounted for separately in the model ( Fig 7E ), whereas Ed suggests that at high air temperatures, streams that drain steeper and higher elevation basins are generally cooler ( Fig 7F ). The relationship between Tl and A suggests that headwater streams are less sensitive to air temperature, which may reflect greater relative groundwater or snow influence ( Fig 7G ). Similarly, streams with high BFI (the covariate best representing potential influence of and connection to subsurface flow pathways) are also less sensitive to air temperatures ( Fig 7H ). Finally, streams with higher riparian forest cover (F) tended to be cooler at higher air temperatures, representing shading from solar radiation during growing seasons when trees have foliage ( Fig 7I ). Additional conditional effects plots for covariates not discussed here are provided in S3 Fig .

We note that the influence of Sws is also captured indirectly by the variables Qr and Tl. Consequently, the modeled influence of Sws is more muted than might be expected from snowpack melt. In addition, the two air temperature variables (Tl and Tws) and Dl also have some overlap in their impacts. As a result, estimated concurvity is high amongst the temporal variable smoothers in the model (estimates 0.66 to 0.97). The modeled relationships presented in Fig 7 should be interpreted cautiously with this in mind.

As intended when designing our model, much of the variation in stream temperature predictions was associated with air temperature metrics: the antecedent air temperature at a local reach (Tl) and the mean daily air temperature in the upstream watershed (Tws). Both Tl and Tws had roughly equal impacts on stream temperature as shown through their interaction with each other ( Fig 7A ). In addition, interactions between air temperature covariates and other covariates allowed a shift in directional impact in accordance with theory (e.g., BFI and Ed and A). Daylength (DL) was important at higher values, consistent with the expectation that air temperature is more related to stream temperature drivers when days are longer ( Fig 7B ). Area-standardized flows (Qr) were also important at higher values, which likely reflects effects of snowmelt that was not fully captured in Tl or Sws ( Fig 7C ). Higher values of daily snowpack in the upstream watershed (Sws) led to lower stream temperature predictions when air temperatures in the upstream watershed (Tws) were above 0°C, reflecting the direct cooling influence of snowmelt ( Fig 7D ). In addition to the daily snowpack metric (Sws), annual April 1 st snowpack (SA1) had a slightly negative relationship with stream temperature (~ 1°C per 1m), which was stronger in the late summer and fall after the snowpack had generally melted ( Fig 7K ).

PCA of (A) spatial and (B) temporal covariates used in the stream temperature model, indicating that reaches and days with empirical data that were used when fitting models (magenta) are broadly representative of the parameter space for predicting into locations throughout the full region (gray). The spatial dataset used all data; the temporal dataset used a random 10% of reach-days, aggregated by month.

The available stream temperature data used to fit models were broadly representative of the physical characteristics ( Fig 6A ) and of the temporal variation of climate covariates during the study period ( Fig 6B ) across the Pacific Northwest region where we made predictions of stream temperature. For the spatial PCA, three axes were significant and explained 64.0% of the total variation in the data. The spatial covariates with the highest loadings were canopy cover (C), precipitation (P), forested riparian buffer percentage (F), and annual temperature range (R) ( S1 Table ). PC1 (26.3% of variance) predominantly explained variation associated with forest cover and regional climate characteristics; PC2 (23.1%) was associated with geomorphology, and PC3 (14.5%) included elevation and similar elements as the first two axes. For the temporal PCA, three axes were significant and explained 85.2% of the total variance. The temporal covariates with the highest loadings were watershed temperature (Tws), antecedent air temperature (Tl), and day length (DI) ( S1 Table ). PC1 43.8%) predominantly explained variation associated with air temperature, PC2 (25.0%) was associated with streamflow, and PC3 (16.4%) centered on the snowpack.

Modeled relationships between residuals and covariate values are shown in Fig 5 for a set of covariates discussed here; others are shown in S2 Fig . This analysis demonstrates that prediction error tended to be higher at warmer air temperatures during the summer and during low flow and snowpack conditions. In addition, prediction error was higher in reaches with the highest BFI values (over 75), low local canopy cover, average annual precipitation values below 1000 mm, and that exhibited higher air temperature ranges across the year.

Prediction error in the temporal cross-validation procedure (MAPE = 1.33 (<+0.01%), RMSPE = 1.76) was nearly identical to the single model, demonstrating the ability of the model to expand outside of the temporal parameter space used to fit the model. Prediction error increased slightly more in the spatial cross-validation test (MAPE = 1.44 (+0.09%), RMSPE = 1.88 (+0.07%)), suggesting that model has comparatively more difficulty extrapolating into new watersheds, but that spatial extrapolations may still be useful in many cases. Summarizing daily predictions to monthly values improved prediction error by 10–14% across all metrics.

Statistics are shown for daily predictions and for monthly averages (note: not all month/site combinations had data for the entire month and in these cases the available data were averaged for the month). RMSE = root mean squared error; MAE = mean absolute error; RMSPE = root mean squared prediction error; MAPE = mean absolute prediction error.

We predicted daily stream temperature over 30 years for all stream reaches throughout the PNW in a single model with reasonable accuracy and precision (MAE = 1.32, RMSE = 1.76, Table 2 , total effective degrees of freedom = 98.5). The combined predictions from models fit to seasonal, regional, and seasonal/regional subsets of the data demonstrated slightly lower error than the single model (RMSE = 1.70 (-0.03%), 1.63 (-0.07%), and 1.57 (-0.10%) respectively).

Moving average window sizes increased with positive Qr for snow and transitional data ( Fig 4B ), but snow data had larger window sizes of around 9 to 12 days at intermediate conditions Qr (-1.5 to 0). In contrast, rain data demonstrated relatively smaller window sizes of around 3 to 6 days at all flows. The highest maximum correlations between stream temperature and antecedent air temperature occurred for snow data at mid to low Qr ( Fig 4C ). Maximum correlations were generally lowest for rain data and intermediate for transitional data.

Across hydrology classes (snow, rain, and transitional), the moving average window size used to calculate the antecedent air temperature covariate Tl reflected seasonal flow patterns, with the lowest flows occurring during summer and late fall ( Fig 4A ). Snow reaches demonstrated a strong distinct flow peak during spring snowmelt, whereas rain reaches had elevated flows throughout the late fall and winter that more closely mirrored regional precipitation patterns. The transitional class demonstrated annual flow patterns similar to the snow class, but with an earlier and more modest spring-melt peak.

Maps of the PNW illustrating (A) winter average air temperature and (B) mean annual precipitation from the PRISM dataset, (C) mean accumulated snow water equivalent (SWE) in upstream watersheds from the National Snow Water and Ice Dataset, and (D) regions classified in this paper as having rain-dominated, transitional, and snow-dominated hydrology. The mean accumulated snow water equivalent (SWE) was used to classify the hydrograph classes shown in (D). Here, we show the average snowpack accumulation and most common hydrograph classification across the entire dataset. In contrast, for the analyses, annual values were incorporated (see text). Basemap source: a combination of public domain data customized by DH in ArcMap (Esri). Bathymetry data are from NOAA’s National Geophysical Data Center and the British Columbia Marine Conservation Analysis Project Team. Political boundaries and shoreline are from the United States Geological Survey. Hydrography is from the National Hydrography dataset.

Winter air temperature ( Fig 2A ) and annual precipitation ( Fig 2B ) both influence average snowpack accumulation ( Fig 2C ). Each reach in the PNW was classified as rain-dominated (hereafter, “rain reaches”), snow-dominated (hereafter, “snow reaches”), or intermediate (hereafter, “transitional reaches”) based on the mean annual SWE ( Fig 2D ). Note that Fig 2 illustrates the average values across all years in our study period, but for our analyses, we incorporated annual values from each individual year. We used these hydrology classifications to develop our antecedent air temperature metric Tl from an average of 508 rain reaches (mean across years; maximum 1182), 612 snow reaches (maximum 1094), and 566 transitional reaches (maximum 1148). Note that this subset of reaches used to inform Tl were those for which we had empirical stream temperature data (i.e., 10,032 of >261,200 free-flowing stream reaches). Rain reaches occurred at lower elevations and had warmer air temperatures and lower flows but higher precipitation on average compared to snow reaches ( Fig 3 ). Mean annual watershed snow depths were usually < 500 mm for snow reaches (median = 218 mm, but >2000 mm for some reaches). Environmental characteristics of transitional reaches were generally intermediate to snow and rain reaches, but were more like snow reaches for air temperatures (Tws), elevation (E), and base flow (BFI). However, transitional reaches were drier on average (median annual precipitation of 656 mm, compared to 950 mm and 1103 mm for snow and rain reaches, respectively), suggesting that precipitation represented the primary reason for lower snowpack accumulations in comparison to snow reaches on average.

Discussion

In this work, we present a novel stream temperature modeling framework that extends previous efforts [5, 21, 25, 37, 38]. We applied our model to predict stream temperature at fine temporal (i.e., daily) and spatial (i.e., reach) resolutions across the broad spatial extent of the PNW (Fig 8). Using publicly available data and a straightforward statistical GAM framework, we provide daily predictions for >261,200 free-flowing stream reaches from 1990 through 2021. Our model performed well outside the spatial and temporal extent of the data used to inform the model, suggesting that our model captured many of the complex and interacting drivers that determine stream temperature. These results should be immediately useful in conservation planning for lotic systems of the PNW that support cold-water species such as salmonids. Our approach is straightforward and can be easily adapted to new spatial regions or time periods.

PPT PowerPoint slide

PNG larger image

TIFF original image Download: Fig 8. Stream temperature predictions averaged by season (top panel: summer; bottom panels: spring, fall, winter) to illustrate the spatial extent over which predictions were made and temporal variability. Values were mapped to individual reach-contributing areas to improve visualization, but daily data are linked to each stream reach identifier in the National Hydrography Dataset version 2. Gray indicates reaches where we did not make predictions such as in reservoirs or because reaches were heavily influenced by dams (PDA > 0.25). Basemap source: a combination of public domain data customized by DH in ArcMap (Esri). Bathymetry data are from NOAA’s National Geophysical Data Center and the British Columbia Marine Conservation Analysis Project Team. Political boundaries and shoreline are from the United States Geological Survey. Hydrography is from the National Hydrography dataset. https://doi.org/10.1371/journal.pwat.0000119.g008

Model utility The performance of our model compared favorably to other statistical stream temperature models recently applied in the western USA [5, 21, 25] that produced prediction errors of around 0.5 to 2°C (S4 Fig). However, our objectives were more ambitious than most other models in that we predicted daily rather than coarser metrics such as mean monthly stream temperature [5, 21], across a broad geographic region rather than within individual watersheds [23], and across the entire calendar year in contrast to seasonally [51]. Although our results demonstrated that predictions were improved by fitting models to smaller spatial or temporal extents, average prediction errors were only 0.1 to 0.2°C better than the single model, an improvement that is small enough to be within the error range of most measuring devices. The flexibility of our statistical model structure was essential to its ability to accurately predict across a geographically diverse region characterized by substantial seasonal and interannual climate variability. Most previous models did not explicitly account for hydrograph class (snow, rain, or transitional) despite that previous research has clearly demonstrated that the amount of winter precipitation falling as snow can strongly influence stream temperature and its relationship with air temperature across multiple seasons [26, 34, 37, 52]. Also, previous statistical models typically used linear, non-interacting effects of air temperature on stream temperature, which does not reflect the nonlinear and spatiotemporally variable nature of this relationship [53, 54 but see 38]. Our model overcame these limitations by including (1) an antecedent air temperature metric that was dependent on snowpack and flow, (2) nonlinear relationships, and (3) covariate interactions primarily with air temperature that allowed the influence of covariates to shift seasonally as theoretically expected. For example, streams with high base flow (BFI, which represents groundwater influence) should be (and were) less sensitive to air temperature, and therefore should be (and were) warmer in the winter and cooler in the summer on average compared to streams with low base flow. These adaptations allowed us to better capture the dynamism of the air-stream temperature relationship and overcome a known drawback of linear statistical models [13, 14]. A common concern in using GAMs is that relationships can be over-fit, capturing small patterns in the data that are not caused by covariate influences [41]. Following previous research, we limited model complexity by reducing the number of smoother knots [25, 31]. Automated parameter penalties within the fitting functions of the modeling package alone led to overfit relationships. The success of our model in the cross-validation tests suggests that our model is not over-fit, and other research has supported GAMs as appropriate statistical models for stream temperature [27, 28, 31]. In addition, the relationships parameterized in our statistical model align with theoretical expectations, as demonstrated in conditional effects plots (Fig 7 and S3 Fig). These results support our decisions to limit model complexity by constraining the number of knots in climate variable smoothers and to fit linear relationships to temporally constant variables. Visually examining concordance between expected and observed conditional effects is an important model quality control step when considering covariates to include in a model [13, 25]. Our statistical model was less expensive in data requirements, expertise, execution time, and computing power than a more complex process-based model. The diversity of spatial and temporal scales and physical forcing pathways that drive stream temperature over broad spatial and temporal extents is so vast as to make purely mechanistic physical models impractical [13, 14, 16–18]. The variables we used in this model are publicly available throughout the contiguous United States, although retrospective expansion is limited by the historical availability of the climate metrics. Our results suggest that statistical methods that can generally and robustly model patterns that emerge from physical processes are an appropriate compromise between physical process representation and data and computational requirements. Our model successfully predicted stream temperatures outside of the spatiotemporal bounds of the data used to fit the model, suggesting that our model’s statistical relationships represented mechanistic controls and can therefore be useful in predicting stream temperatures under un-monitored conditions. Some statistical models have shown poor prediction accuracy when extrapolating to novel time periods or regions (e.g., [53, 55]). This likely stems from a combination of extrapolating predictions to streams with dissimilar hydrological classes and not including variables that represent distinct thermal regimes. The ability of our model to extrapolate predictions was improved by fitting the model across a large, geographically diverse region and over multiple decades representing substantial interannual climate variability. Accordingly, data in temporal expansions (e.g., filling gaps or extending time-series) or spatial expansions (e.g. predicting in adjacent basins) are more likely to be represented in the wide-ranging covariate space used to inform the model and thus more accurately predicted. While successful across a broad geographic area and variable climate conditions, model performance was not consistent across space and climate. To better understand how model performance varied, we examined the relationship between the absolute value of residuals and model covariates. This analysis showed that the absolute value of residuals is higher during warm and dry summer conditions. However, as also described in Siegel and Volk [25], there is much more variation in stream temperature across the landscape during the warm summer months. Thus, while absolute prediction error is higher, the model actually captures a higher proportion of the variation in stream temperature during the warm season in comparison to the cooler winter season. In addition, the model performed more poorly at reaches that receive little precipitation and have comparatively low riparian canopy cover. These stream reaches are more likely to consistently exhibit lower flows, be more influenced by variation in solar radiation due to a lack of shading, and be more sensitive to local forcing factors that can be difficult to accurately capture in model covariates (e.g. localized groundwater springs or agricultural runoff/diversions). Finally, the model also performed more poorly in the reaches with the highest level of groundwater influence as measured by BFI. Reaches with substantial groundwater influence are more likely to be decoupled from air temperature.

Potential applications Our model’s ability to predict temperature across a large range of spatiotemporal scales makes it a powerful tool for informing and supporting species and habitat management. Predictions can be used for evaluating different management and restoration scenarios [56, 57]. For example, consulting maps of predicted stream temperature during times of potential ecological stress may indicate where increased riparian shade, floodplain reconnection, or other temperature-moderating actions might be necessary to improve ecosystem function. If our modeling framework were to be applied in other regions, a cross-region comparison could provide information useful for considering whether management approaches that have worked well in one region may be applicable in another. For protected cold-water species like salmonids, a potential application of our results could be to evaluate thermal habitat across life history stages to guide management strategies. Our stream temperature predictions can help determine if management actions should resist change (e.g., conserve functional thermal habitat), accept change (e.g., recognize that some habitats are unlikely to become as functional as they may have been in the past), or direct change (e.g., actively alter stream temperatures or move sensitive species) [10]. For example, our predictions can identify streams that are more stable throughout the year, such as snowpack-influenced or groundwater-fed streams. Streams with stable thermal regimes can be more productive than regulated reaches, and may be more resistant to climate change, making them better targets for species reintroductions [26, 58–61]. As another example, our model could help practitioners envision potential stream temperatures as if a dam did not exist (since we did not model dam effects), which may be useful for assessing dam management and/or removal strategies. Another powerful application of our model would be to predict the effects of climate change on stream temperature. The model successfully predicted stream temperature outside of the spatiotemporal bounds of the data used to fit the model, so it is likely to perform well in unsampled time periods. Our methods provide improvements for parameterizing the seasonal and annual impacts of snowpack and associated flows, which are expected to be strongly influenced by climate change. In addition, fitting the model across such a large and diverse region across numerous climatically distinct years will likely be beneficial in predicting climate impacts. Even in years with annual climate values outside the range of those used to fit the model, few of the daily covariate combinations are likely to fall outside combinations already seen by the model. For example, warmer air temperatures are expected to result in decreased snowpack throughout most of the Northern Hemisphere, shifting many regions from snow-dominated to rain-dominated hydrology [62], but our model already accounts for the predicted range of future snowpack conditions and the influence of different hydrology types on stream temperature. If future air temperatures are warmer and snow melts earlier, spring conditions will likely resemble past summer conditions, and higher elevations and latitudes may mimic historical conditions at lower elevations and latitudes. Therefore, only the most extreme future conditions (e.g., the warmest sites during the warmest time of year) would be outside of our modeled parameter space. The production of accurate predictions of future air temperatures, stream flows, and snowpack accumulations to inform climate covariates may be a larger source of error than the stream temperature model.

Model limitations and future directions While we believe that our results are useful for examining patterns in thermal regimes across landscapes and in conservation planning, our stream temperature predictions should be interpreted within the limitations of the modeling approach and prediction accuracy. Prediction error was low on average, however local model residuals should be consulted before using predictions for finer scale applications. If local error is considered large, our dataset and methodology could be used to fit a more tailored model if enough local stream temperature data exists to inform it. Our methods produced predictions representing well-mixed mainstem temperatures of stream reaches that were generally ~1–2 km in length. Accordingly, the model does not reflect finer scale thermal heterogeneity from features such as localized groundwater inputs, backwater areas, deep pools, or localized riparian shading. Such heterogeneity can have important impacts on biota. Similarly, as with other statistical models, the model is not designed to predict temperatures in stagnant or thermally stratified waters such as slow-moving pools, lakes, or reservoirs [5, 21]. The model does not account for flow regulation and dam operations; we therefore did not predict stream temperature for reaches likely to be heavily influenced by dams. Because regulated watersheds often produce unnatural flow conditions and thermal regimes [61], statistical models have difficulty accurately predicting below dams (e.g., [26, 37, 55] without additional dam-specific covariates [5]. Larger dams often release water from below the thermocline in the reservoir, which can lead to large differences in stream temperature (e.g., reductions in the summer and increases in the winter) that can persist many kilometers downstream of the dam [63]. The strength of these effects depends on the size of the dam and how it is managed, discharge volume, water release temperature, and environmental conditions like flow and air temperature [53, 63]. These factors make accounting for dam impacts in a spatially expansive model like ours particularly difficult. An analysis comparing empirical measurements of stream temperature to predictions from our existing model (that did not account for dams) would provide insights into dam-specific spatiotemporal impacts on stream temperature and potentially could inform the development of covariates to account for the impact of dams across a regional scale. Local geology is one of the main factors controlling groundwater influence on stream temperatures [64–67]. For example, Tague et al. [68] found that water temperature in western Oregon streams underlaid by volcanic rock >7 million years old (Western Cascades) was more synchronous with air temperature compared to streams underlaid by volcanic rock <2 million years old (High Cascades), inferring that groundwater played a larger role in the geologically younger, High Cascades streams. Here, we included the percent of the watershed covered by extrusive volcanic rock (V), but it is unclear how broadly this applies to other regions. Similarly, Dralle et al. [69] found that underlying geology in the Eel River basin, California, influenced water storage capacity, flow dynamics, vegetation type and canopy cover, all of which subsequently impacted water temperature. The depth of groundwater aquifers can also be important in determining groundwater influence on surface stream temperature [70]. Although we represented potential surface-subsurface interactions with an index of base flow (BFI), it is difficult to empirically assess groundwater inputs and losses at a large spatial scale. The model demonstrated increased prediction error at the highest BFI levels, suggesting that the model could be improved with further exploration into the ways to parameterize the impact of groundwater interactions. A recently developed base flow model produced by Lombard et al. [71] seems a promising way to improve representation of the interaction of subsurface and surface water on stream temperature. In addition, planform and longitudinal channel morphology influences rates of surface and hyporheic aquifer exchange, and thus stream temperature, in a complex and not well characterized manner [72]. This is strongly influenced by the degree of riverscape impairment, particularly floodplain disconnection [73]. Thus, more detailed riverscape context could be incorporated as new covariates or used to define more stream types beyond snow-dominated, rain-dominated, and transitional [74]. Vegetation covariates (i.e., canopy cover and riparian shading) in this and other statistical stream temperature models are temporally constant [5, 21, this study], but vegetation cover changes seasonally and spatially [75] and will likely be impacted by climate change. Future improvements should attempt to incorporate vegetation as a temporal covariate. Similarly, the effect of wildfires should be incorporated because wildfires can dramatically decrease riparian vegetation and alter flow paths, impacting stream temperatures, with the strength and direction of the effect varying by season and the severity and location of the wildfire [76, 77]. Furthermore, it is likely that the impacts of wildfires on stream temperature will increase as wildfire severity and frequency intensify with climate change [78–80]. As discussed by Siegel et al. [37], the parameterization of the predicted averaged antecedent air temperature window size could potentially be improved with more predictive information. While we included a third “transitional” category in addition to the two “snow” and “rain” dominated categories previously used, further explorations could examine the impact of further dividing snowpack or other factors to create additional categories. Such factors could include watershed/channel and size, riparian shading, geology, channel morphology, and the characteristics of the watershed’s groundwater hydrology. For instance, antecedent air temperature window size in rain-dominated catchments may be influenced by watershed area. A relatively high level of nonlinear dependency among predictor variables, or concurvity [49], was estimated amongst the smoothed terms of the model. This result is not surprising given that air temperature, snowpack, flow, and day length are related factors. Nevertheless, we found that each of these variables also had a unique impact which was necessary for fitting our model across a broad geographic area and for capturing climate variability. Included interactions were necessary for the model to distinguish the impacts amongst these related variables. However, due to this interdependence estimated variable relationships should be interpreted with this conflation in mind. Further exploring how to better isolate the independent impacts of climate variables to lower concurvity [49] would make parameterized relationships easier to interpret, would likely lower computing time, and could potentially improve predictions. Finally, we recognize that the application of autocorrelation functions in conjunction with the modeling techniques we present in this study may further improve prediction accuracy in some cases. One of our objectives in this study was to explore the replacement of rigid autocorrelation functions that aim to account for processes that limit the temporal and spatial rate of change of stream temperature with covariate relationships that can vary over space and time. In addition, we aimed to produce a model that could be easily used to predict unmonitored time periods and watersheds which may lack the real world stream temperature data to accurately inform autocorrelation functions. Accordingly, we did not include autocorrelation functions in our model to best test the ability of our covariates to account for these processes. The ability of our model to produce accurate cross-validated predictions across space and time suggests that our covariates were largely successful in this effort. However, we recognize that a level of autocorrelation remains in our model residuals and are in support of future modeling efforts exploring combining autocorrelation functions with our methodology. In general, we believe that stream temperature modeling remains a dynamic field with further gains to be made by continuing to examine covariate parameterization and statistical techniques.

[END]
---
[1] Url: https://journals.plos.org/water/article?id=10.1371/journal.pwat.0000119

Published and (C) by PLOS One
Content appears here under this condition or license: Creative Commons - Attribution BY 4.0.

via Magical.Fish Gopher News Feeds:
gopher://magical.fish/1/feeds/news/plosone/

You are viewing proxied material from magical.fish. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.