Extreme hydrometeorological events induce abrupt and widespread freshwater temperature changes across the Pacific Northwest of North America

Study area

The Pacific Northwest spans a vast and diverse region of western North America⁴⁷. Here, we define the PNWNA to encompass southeastern Alaska, BC, Washington, Oregon, and parts of Idaho, Montana, and Wyoming that lie within the Columbia River Basin (Supplementary Fig. 9). This region spans from 42 to 60°N and from 110 to 140°W. The area includes major mountain chains like the Cascades, Coast, Columbia, Rocky, and St. Elias Mountains that are generally oriented north-south. Principal waterways draining to the Pacific Ocean include the Columbia, Fraser, Nass, Skeena, and Stikine Rivers, while the upper Peace River is a major system draining to the Arctic Ocean⁴⁸. While coastal rivers typically exhibit pluvial or hybrid (nival and pluvial) regimes, waterways draining mountainous and inland regions exhibit nival or glacial regimes⁴⁹. Several waterways, particularly in the Columbia River Basin, are regulated by hydroelectric dams and reservoirs; however, this is less common in the northern half of the study domain where the largest hydroelectric facilities and control structures lie within the Bridge-Seton, Nechako, and Peace River basins⁵⁰.

The PNWNA includes a wide range of climates from maritime (relatively wet and moderated temperature variability) along the Pacific Coast to continental (relatively dry, particularly to the lee of mountain chains, and amplified temperature variability) in interior regions. Sub-freezing air temperatures and abundant snowfall mark winters in high-elevation and northern regions. As such, freshwater temperatures exhibit pronounced seasonality, with peak temperatures during summer and minimum values in winter⁵¹. Nival and glacial rivers often develop an ice cover during winter and near-surface, unfrozen water remains near the freezing point. Water temperatures generally start rising once the ice cover melts and thereafter exhibit pronounced diurnal cycles.

Data

Time series of sub-daily and daily freshwater temperature (T_w) were assembled from federal and provincial agencies including the United States Geological Survey (USGS), Water Survey of Canada (WSC), Fisheries and Oceans Canada (DFO), and BC provincial monitoring sites. These were supplemented with data from private industry (Rio Tinto, Palmer Environmental Consulting Group) plus data from UNBC networks collected by the authors in the Nechako⁵², Parsnip⁵³, and Quesnel⁵⁴ river basins. Metadata compiled for each site included collection agency, site name and identification number, coordinates, catchment area, and period of record.

Near surface meteorological conditions during the 2021 extreme hydrometeorological events were extracted from the fifth generation of the European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis (ERA5⁵⁵). Variables of interest included the hourly surface air temperatures and precipitation. To illustrate the unusual atmospheric conditions observed in the four events of interest, additional daily air temperature and precipitation data were sourced from Environment and Climate Change Canada for meteorological stations at Prince George, BC (station ID 1096453) and Hope, BC (station ID 1113542).

Site selection

Three primary criteria were used to select study sites: (1) freshwater body (river, lake, or reservoir) within the PNWNA as defined in the study area description; (2) availability of (sub-)hourly water temperature data in °C for part of, if not all of, 2021 and 2022; and (3) availability of corresponding metadata (coordinates, basin area, etc.) for the study sites. Based on these criteria, a total of 554 sites across the PNWNA were chosen for analysis (Supplementary Fig. 9 and Supplementary Data 1). Study sites comprise 536 rivers or creeks, 16 lakes, one reservoir, and one pond. A subset of 193 sites from the full database with long-term (≥10 years of data up to 31 December 2022) records of mean daily or hourly water temperatures were retained for additional analyses. Furthermore, hourly or daily discharge data were extracted for 270 USGS and WSC sites to explore relationships between hydrological and water temperature conditions.

Time series construction

Hourly water temperature data from 00:00 local standard time on 1 January 2021 to 23:00 local standard time on 31 December 2021 were extracted for analysis. The original time stamps for all time series were retained to reflect local standard times (Mountain, Pacific, and Alaska time zones, noting most (91%) study sites fall within the Pacific time zone). Sites where data were available at a sub-hourly scale were degraded to hourly starting at 00:00 local standard time each day. Depending on site data availability, time series of hourly water temperature data were also constructed for 2022 to assess post-E4:LT weekly changes in hydrothermal conditions. While instrument specifications were not available for most sites, water temperatures recorded within UNBC monitoring networks have an accuracy of ±0.50–0.53 °C and a resolution of 0.04–0.14 °C⁵², which we expect to be representative of data accuracy for other networks.

As the quality of data retrieved from various sources was largely unknown, we developed an automated strategy to flag suspicious and potentially erroneous water temperature data. Quality assurance and quality control (QA/QC) thus followed a method adapted from Gilbert et al.⁵² and Sowder and Steel⁵⁶. Criteria used to flag potentially spurious hourly T_w data at timestep i were: (1) abrupt thermal variations of |T_w,i − T_w,i-1| > 2.5 °C; (2) |T_w,i – \(\bar{{T}_{{{{{{\rm{w}}}}}},i5}}\)| > 1.5 °C where \(\bar{{T}_{{{{{{\rm{w}}}}}},i5}}\) is the rolling mean of five successive T_w data points centered on the current timestep i; and (3) an excessive T_w diurnal range (DR) > 10 °C. All data flagged by the QA/QC criteria were then considered as missing. For criterion 3, all available data during the 24-h period (starting at 00:00 local standard time) were set as missing when DR > 10 °C. All sub-freezing T_w data that passed the QA/QC procedure were set to 0.0 °C. Temporal gaps <24 h were in-filled through linear interpolation; otherwise, data were excluded from analysis if >10% of values were missing for a specific period of interest. Supplementary Fig. 10 demonstrates an application of the automated QA/QC procedure in flagging suspicious hourly water temperatures for the Fraser River at Marguerite, BC during 2021.

Analyses

Analyses focused on four extreme hydrometeorological events across the PNWNA during 2021: E1:LT the cold wave of 7-13 February 2021; E2:HT the heat dome from 25 June to 1 July 2021; E3:HF the catastrophic flooding in southwestern BC and northwestern Washington during 11–15 November 2021 caused by a pair of landfalling atmospheric rivers; and E4:LT the early winter cold snap of 25–31 December 2021. To facilitate interpretation of the results, the four events of interest were assigned a dominant hydrometeorological feature of either low air temperatures (LT), high air temperatures (HT) or high flows (HF). General meteorological characteristics of the four extreme hydrometeorological events are provided in Supplementary Figs. 1–3.

Each site’s 7-day (E1:LT, E2:HT, and E4:LT) or 5-day (E3:HF) average freshwater temperature was first calculated for the four events of interest when <10% of data were missing. Similar averages were computed for the 7-day or 5-day periods prior to and post events. Then, freshwater temperature differences (ΔT_w) between days during and prior to the extreme hydrometeorological events were computed. For E1:LT, E2:HT, and E4:LT, we subtracted the average hourly water temperatures of the week prior to the event from those during the event. Given the more transient nature of the shorter duration E3:HF, a similar process was taken but here focusing on the 5 days prior to and 5 days during the event. A similar approach was also employed for values of standard deviations in hourly water temperatures. Statistics (minimum, median, mean, maximum, standard deviation and interquartile range) on the event-based water temperature changes were tabulated and depicted as box and whisker plots. The Shapiro-Wilk test was applied to assess whether water temperatures were normally distributed using a probability p = 0.05 as the threshold⁵⁷. As all assessments indicated departures from the normal distribution (Table 2), the non-parametric Wilcoxon rank sum test was then applied to the paired values of average water temperatures prior to and during the weekly or 5-day events to assess if changes in the average were significantly (p < 0.05) different from zero⁵⁸. Changes in the means and standard deviations of water temperatures were also plotted spatially to reveal regional patterns during each event. An additional analysis illustrates spatial differences in freshwater temperatures between periods during and post events.

To explore the spatial synchrony of water temperature changes during Events 1–4, an additional metric was assessed. We first standardized all time series of 2021 hourly water temperature data using each site’s overall mean and standard deviation. Thereafter, all possible pair-wise Spearman rank cross correlations (ρ) between standardized hourly water temperatures during the four events were computed. The median value (ρ_m) of all pair-wise cross correlations was then taken as a measure of overall spatial synchrony during Events 1–4.

To place into a longer-term context observed thermal conditions during four extreme hydrometeorological events in 2021, time series of median values and interquartile ranges of all hourly water temperatures plus weekly changes across all sites were then compiled. As well, for a subset of 193 sites with long-term water temperature records (Supplementary Data 1), standardized anomalies from daily average conditions were computed and averaged for the four events of interest. Here, standardization employed daily means and standard deviations based on each site’s period of record. This provided seven values of standardized anomalies for E1:LT, E2:HT and E4:LT, and five values for E3:HF at these sites. These were then averaged for each event and reported as single values in Fig. 5. Thus, these results were restricted to days when the four extreme hydrometeorological events unfolded and excluded comparisons to other times of the year. In addition, we tracked the number of daily minimum and maximum water temperature records during each of the four events of interest for sites with records spanning at least 10 years.

Potential associations between weekly (5-day) changes in water temperatures and site metadata, hydrological and meteorological conditions were then assessed. Here, ρ values between site metadata (latitude, longitude, catchment area, and mean elevation where available), average water temperatures for the week (5 days) prior to an event, potential incoming solar radiation (K↓), and water temperature differences were computed. Relationships between changes in observed streamflow and water temperatures were investigated for 270 WSC and USGS sites where both quantities were available (Supplementary Data 1). Changes in streamflow were expressed as a percentage relative to the week (5 days) prior to the event. For air temperature, hourly ERA5 data at 2 m above the surface were extracted for the grid cell overlying the water temperature measurement site and averaged for the week (5 days) prior to and during the event. Weekly (5-day) changes in ERA5 air temperature were then correlated to the corresponding water temperature changes. Potential K↓ for a flat surface was calculated at 1-min intervals using subroutines from the Catchment-based Land Surface Model⁵⁹ and time-integrated to daily radiation totals (in MJ m⁻² day⁻¹) based on each site’s coordinates. All ρ values were considered statistically significant when p < 0.05, noting that any significant correlation does not necessarily imply causation⁶⁰.

Although linear or nonlinear regressions are commonly used to assess the thermal sensitivity (TS, °C °C⁻¹) of stream temperatures^{11,22,61,62,63}, we defined TS here as the ratio of changes in weekly (5-day) water temperature to the corresponding changes in air temperature at each site:

Results for TS were plotted spatially and the median of all values (TS_m) taken as an overall measure of thermal sensitivity across the PNWNA during each event.

A final analysis synthesizes key results by ranking various statistics of T_w responses to the four extreme hydrometeorological events across the PNWNA relative to the remainder of 2021. Statistics used include the means, standard deviations and interquartile ranges of T_w and ΔT_w for all available sites which were then each assigned a nearest-rank percentile relative to all other 7-day (E1:LT, E2:HT and E4:LT) or 5-day (E3:HF) periods during 2021. Here, the T_w standard deviation was computed at each site using all available hourly T_w data for each event, and the standard deviations were then spatially averaged. Conversely, the interquartile range was computed across all available sites at each hour and then averaged over the duration of a period. Statistics for the means and standard deviations in ΔT_w represented differences between values prior to and during an event. All analyses were performed across the entire PNWNA as well as the northern, central, and southern sub-domains.

Data processing, limitations, and uncertainty

Data QA/QC and all analyses were performed using the Fortran computing language with the exception of statistics that were computed using R Studio (version 4.1.3). Spatial plots were generated in QGIS (version 3.28.5) while line, box and whisker plots were created using Grace (version 5.1.25).

Although the automated QA/QC procedure flagged on average 0.35% of the hourly T_w data across all sites during 2021, additional erroneous measurements potentially remained within the database. Similarly, some valid data may also be omitted from the database through this step. Thus, caution must be used in interpreting results for a specific site; as such, results presented herein focus on overall statistics and regional patterns during the four events of interest. Some river systems (e.g., the Columbia and Fraser Rivers) include multiple measurement sites along their main stems and tributaries thereby violating data independence. Another limitation with water temperature data used in this study is that site representativeness, instrument specifications, and set-up details are largely unknown for 522 sites outside of the UNBC networks. Nevertheless, intercomparisons of hourly water temperature data collected by different agencies in close proximity suggest they are entirely consistent with each other (Supplementary Fig. 11 and Supplementary Table 2). Thus, the quality-controlled database we have assembled is relatively homogenous and representative of temperature variations of surface waters.

Another potential shortcoming of this study is the 5-day averaging period for E3:HF but otherwise 7-day averaging periods for E1:LT, E2:HT, and E4:LT. A 5-day averaging period was selected for E3:HF given the shorter duration of this event⁴ relative to the other three extreme hydrometeorological events (Supplementary Figs. 3, 7). A comparison of binned freshwater temperature changes averaged over five and seven days centered on 13 November 2021, however, showed statistically indistinguishable results (Supplementary Fig. 12). Thus, to focus on the direct impacts of the two atmospheric rivers and high flows on freshwater temperatures, a 5-day period of analysis was retained for E3:HF.

A further limitation on this effort was the limited number of records used to assess departures from long-term records. This was most prevalent at Canadian sites where long-term records were particularly sparse with only one site (Fraser River at Hope, BC) having ≥30 years of data for the period covered by E2:HT and none for the other three events (Supplementary Fig. 13). Meanwhile, there were 37–40 sites with at least three decades of data in the US during periods spanning all four events. As such, this analysis was biased toward the southern portion of the study domain. Nevertheless, standardized anomalies for those sites with ≥10 years of data were representative of the subset of sites with ≥30 years of data (Supplementary Fig. 14).