Increasing heatwave frequency in streams and rivers of the United States

Heatwaves are increasing in frequency, duration, and intensity in ocean, coastal, and lake ecosystems. While positive water temperature trends have been documented in many rivers, heatwaves have not been analyzed. This study examined heatwaves in rivers throughout the United States between 1996 and 2021. Riverine heatwaves increased in frequency over the study period, with the most robust increases occurring in summer and fall, in mid‐ to high‐order streams, and at free‐flowing sites and sites above a reservoir. The increase in heatwave frequency was accompanied by an increase in moderate strength heatwaves as well as a doubling of the annual mean total number of heatwave days at a site. Riverine heatwaves were often associated with normal or below‐normal discharge conditions and at sites with a mean annual discharge ≤ 250 m3 s−1. These results provide the first assessment of heatwaves in rivers for a large geographic area in the United States.

Increasing climate variability and rising water temperatures have spurred interest in anomalous extreme hightemperature events called heatwaves due to their potentially disproportionate impact relative to a gradual, long-term change in mean temperature (Jentsch et al. 2007;Vasseur et al. 2014). Positive trends in heatwave characteristics, such as frequency (i.e., total events per unit time), duration (i.e., total days an event lasts), and intensity (i.e., the magnitude of an event) have been documented in ocean, coastal, and lake ecosystems (Lima and Wethey 2012;Oliver et al. 2021;Woolway et al. 2021) and have been linked to shifts in ecosystem metabolism (Berger et al. 2020), organism mortality (Till et al. 2019), and poor water quality (Tassone et al. 2022). Lotic water temperatures have been increasing since the mid-20th century (Webb and Nobilis 1995;Kaushal et al. 2010). However, heatwaves in riverine systems have not been considered, providing uncertainty in riverine heatwave distribution and trends.
Aquatic heatwave development is dependent upon vertical diffusive heat flux across the air-water boundary and horizontal advective heat flux (Holbrook et al. 2019;Oliver et al. 2021). In running waters, the primary control on diffusive heat flux is seasonal and spatial variability in atmospheric temperature (Caissie 2006). There are numerous controls on advective heat fluxes in lotic systems and include precipitation and discharge regimes, flow modification by reservoirs, groundwater supply, and stream size (Vannote et al. 1980;Caissie 2006;Haddeland et al. 2014;Hare et al. 2021). With projected climate warming and increasing demand for freshwater, diffusive and advective heat fluxes are expected to intensify throughout the 21st century and increase the frequency, duration, and intensity of extreme high water temperature conditions (Meehl et al. 2007;Wuebbles et al. 2017).
In this study, we document the frequency, duration, and intensity of riverine heatwaves using in-situ water temperature measurements from 1996 to 2021 for 70 riverine sites located throughout the United States. We also compare regional and seasonal patterns in heatwaves and test relationships between heatwave characteristics, atmospheric temperature, precipitation, discharge, stream order, and position relative to reservoirs. We hypothesize that heatwave characteristics will be positively affected by increasing atmospheric temperature as well as declining precipitation and discharge. We also hypothesize that heatwaves will be more likely at sites above reservoirs where dam operations do not ameliorate temperature extremes and in mid to higher stream orders where riparian cover and groundwater contributions that moderate temperatures are reduced.

Data sources and processing
The United States Geological Survey (USGS) conducts surface water monitoring throughout the United States as part of its national water information system (https://waterdata.usgs. gov/nwis/; U.S. Geological Survey 2016). All sites with daily mean water temperature records available for the 26-year period of 1996-2021 were identified using the R package "dataRetrieval" version 2.7.6 (De Cicco et al. 2022). Tidally influenced and lake stations were removed prior to analysis. Periods of ice cover and flagged data other than those "Approved," "Approved Revised." "Approved Edited," or "Provisional" were considered missing values. Sites with < 90% of their daily records were also excluded. Linear interpolation was applied to water temperature gaps ≤ 2 d. For larger gaps, multiple linear regression models were developed using 1-km 2 resolution, daily climate data provided by Daymet version 4 ( Thornton et al. 2020) using the R package "daymetr" version 1.6 (Hufkens et al. 2018). Only those sites with regressions where R 2 ≥ 0.80 were used in this analysis (mean AE SD R 2 = 0.91 AE 0.04). Of the 70 long-term water temperature sites, 51 had concurrent daily mean discharge that were ≥ 90% complete after removing flagged data. This resulted in a total of 1820 station years of air temperature, precipitation (daymetr), and water temperature (USGS); and 1326 station years of discharge data (USGS) available for analysis.
Site specific variables included region, Strahler stream order, and position in landscape relative to a reservoir. Regional assignment was classified according to historical climatically consistent regions of the United States ( Fig. 1; Karl and Koss 1984). Strahler stream order was determined using the USGS NHDPlus High Resolution geospatial database using ESRI ArcMap version 10.8. Categorical assignment of site position relative to a reservoir (i.e., above, below, none) was determined from aerial photographic visual inspection of each site's location relative to the U.S. Army Corps of Engineers National Inventory of Dams (NID) in ArcMap (https://nid.sec. usace.army.mil, see Supporting Information S1).

Heatwave detection
Riverine heatwaves were identified following the Hobday et al. (2016) definition for marine heatwaves as those periods when daily mean water temperature exceeds a local, seasonally varying 90th percentile threshold for a period ≥ 5 d. Similarly, heatwaves were classified based on their peak intensity relative to multiples of the 90th percentile difference from the long-term climatology (Hobday et al. 2018). Heatwave identification and classification were obtained using the "heatwaveR" R package version 0.4.4 (Schlegel and Smit 2018). For each heatwave, we quantified event duration and cumulative intensity above the 90th percentile. For each site, heatwave frequency was determined as the annual number of heatwave events. High and low discharge was determined as the difference between each site's observed daily mean discharge and an expected daily mean discharge (hereafter referred to as "residual discharge"). The expected daily mean discharge was determined using the "heatwaveR" R package which calculates a centered 11-day mean discharge climatology based on day of year (i.e., 1 January = 1, 31 December = 365) which is then smoothed using a 31-day moving average (Schlegel and Smit 2018). Residual discharge values > 0 and < 0 indicated above and below normal discharge respectively.

Statistical analyses
Long-term trends in monthly mean air and water temperature, discharge, and total precipitation (all forms converted to water equivalent), were determined using Seasonal-Kendall test with Sen's slope estimates ("wql" R package version 0.4.9; Jassby and Cloern 2017). Trends in heatwave characteristics were conducted on annual time series as heatwaves are rare events and often span monthly boundaries. Annual trends were determined using Mann-Kendall tests ("Kendall" R package version 2.2; McLeod 2011) with Sen's slope estimates ("trend" R package version 1.1.4; Pohlert 2020). These nonparametric trend analyses are robust to outliers and have been used in similar studies of long-term change in water temperature (Hirsch et al. 1982;Kaushal et al. 2010) and heatwaves (Perkins and Alexander 2013;Tassone et al. 2022). To reduce Type 1 errors that could arise from multiple comparisons, a 15% false discovery rate was applied to the heatwave analysis prior to determining statistical significance (Benjamini and Hochberg 1995;McDonald 2014). Annual mean discharge was normalized by site to examine relative trends. Data collection and analyses were performed in the R statistical environment (R Core Team 2020).

Season
There were a total of 3984 riverine heatwave events (35,360 d) between 1996 and 2021 for the 70 sites (664,300 daily observations) considered in this analysis. The mean frequency, duration, and maximum intensity above the 90th percentile threshold for heatwaves were 2 events yr À1 , 9 d, and 1.7 C, respectively, and ranged up to 12 events yr À1 , 103 d, and 9.0 C (Supporting Information Fig. S1). Summer heatwaves made up 28% of all heatwaves, followed by spring (26%), fall (24%), and winter (22%). Overall, riverine heatwave frequency significantly increased over time (0.06 events yr À1 , p = 0.004) and for the summer (0.03 events yr À1 , p = 0.004) and fall seasons (0.02 events yr À1 , p = 0.029, Table 1). The increased frequency of heatwaves was accompanied by significant increases in moderate strength events (3.4 events yr À1 , p = 0.003; Fig. 2a). Similarly, moderate strength events significantly increased in duration (30.1 d yr À1 , p = 0.002; Fig. 2b) and maximum intensity (10.5 C yr À1 , p = 0.01; Fig. 2c). There was a significant increase in and approximate doubling of the annual mean total number of heatwave days at a site from 11 in 1996 to 25 in 2021 (slope = 0.53 d yr À1 , p = 0.015).
During riverine heatwaves, discharge generally was at or below normal discharge conditions (Supporting Information Fig. S2a). Some heatwaves occurred when discharge was at or above normal conditions, particularly between November and February. However, those heatwaves that occurred between March and October, which includes North American summer, were associated with periods of low discharge. The annual total number of heatwave events was greatest in rivers that had an annual daily mean discharge ≤ 250 m 3 s À1 (Supporting Information Fig. S2b). The cumulative intensity above the 90th percentile threshold of riverine heatwaves was lowest between July and September and greatest between October and June (Supporting Information Fig. S3).

Region
The South had the largest median (61 events site À1 ) and maximum (76 events) number of riverine heatwave events over the study period. Conversely, the southwest had the lowest median (44 events site À1 ) and minimum number of events (39 events). There were significantly more riverine heatwaves in the northeast (mean AE SD = 60 AE 7), southeast (59 AE 6), and south (65 AE 10) regions compared to the southwest (47 AE 9; p ≤ 0.038). Heatwave frequency trends were positive at 49% of sites, with significant positive trends (n = 19) ranging from 0.06 events yr À1 in the southwest to 0.20 events yr À1 in the northwest (Fig. 3a). Regional heatwave duration trends were predominantly positive (61%), with 11 statistically significant positive trends ranging from 0.13 to 0.46 d yr À1 (Fig. 3b). Similarly, cumulative intensity trends were positive for 60% of sites, with four statistically significant trends ranging from 0.15 C to 0.22 C days yr À1 (Fig. 3c). Heatwave frequency, duration, and cumulative intensity trends were negative at 6% (n = 4), 13% (n = 9), and 20% (n = 20) of sites, respectively.
There was a positive relationship between annual mean air and water temperature trends (slope AE SE = 0.43 AE 0.19, R 2 = 0.07, p = 0.024), with the southeast region having the greatest mean increase in air temperature trend (0.04 C yr À1 ; Supporting Information Fig. S4a) and greatest range of water temperature trends (À0.06-0.06 C yr À1 ; Supporting Information Fig. S4b). Precipitation and discharge trends were positively correlated (slope AE SE = 0.73 AE 0.20, R 2 = 0.21, p < 0.001) with all significant precipitation declines being limited to the northwest (Supporting Information Fig. S4c). The northwest, west, and southwest all had significant declines in discharge trends, ranging from À4.28 to 0.01 m 3 s À1 yr À1 (Supporting Information Fig. S4d). Annual total precipitation and normalized annual mean discharge were significantly positively correlated in all regions (slope range = 0.01-0.27%, p ≤ 0.01) except for the west and west north central. Normalized annual mean discharge was significantly negatively correlated with annual mean water temperature in the south (slope = À0.005 C, p = 0.026), southeast (slope = À0.007 C, p < 0.001), and southwest (slope = À0.004 C, p < 0.001).

Stream order
Riverine heatwave frequency significantly increased in mid (order 4-5) to high (order 8-9) order streams (range = 0.06-0.10 events yr À1 , all p ≤ 0.018; Table 1). The mean duration of 5th order riverine heatwave events increased (slope = 0.08 d yr À1 , p = 0.038, Table 1) however, this result should be interpreted with caution as the p value exceeds the 15% false discovery rate threshold for statistical significance.

Reservoir position
The frequency of riverine heatwaves significantly increased at free-flowing sites with no reservoir (0.06 events yr À1 , p = 0.008) and at sites above a reservoir (0.09 events yr À1 , p = 0.012, Table 1). Similarly, heatwave duration increased at sites above a reservoir (0.06 d yr À1 , p = 0.017). There was not a statistically significant difference in the total number of heatwave events per site with reservoir position, although the median number of riverine heatwaves at sites above a reservoir was greater (median Above = 62) than those below or with no reservoir present (each 56).

Heatwaves are increasing
Riverine heatwaves increased in frequency (0.06 events yr À1 ) over the study period, with significant increases during summer (June-August; 0.03 events yr À1 ) and fall (September-November; 0.02 events yr À1 ). Thus, the frequency of discrete, extended (≥ 5 d) periods of extremely high water temperatures is increasing when water temperature is at its annual peak. However, summer heatwaves generally had the lowest cumulative intensities relative to other times of the year. We did not observe a significant difference in riverine heatwave duration among months, suggesting that the departure from seasonally adjusted water temperature is greatest during fall, winter, and spring. Of the top 1% most extreme and severe heatwaves observed, 92% occurred between December and April, with 8% occurring in July and August. Furthermore, heatwave frequency increased in mid-to high-order streams (range = 0.06-0.10 events yr À1 ) and doubled the average annual total heatwave days from 11 d in 1996 to 25 d in 2021. Sites with an annual daily mean discharge ≤ 250 m 3 s À1 experienced the greatest annual number of heatwave events, indicating smaller gaged streams maybe regularly exposed to heatwaves. The increasing frequency of riverine heatwaves represents an emergent climate change signal outside of linear warming trends, as has been observed in other climate-related variables (Trenberth et al. 2014;Lee 2022).
Riverine heatwaves are increasing in frequency above reservoirs (0.09 events yr À1 ) as well as in free-flowing rivers (0.06 events yr À1 ) but not below reservoirs, suggesting that reservoirs are mitigating increases in heatwave frequency in downstream reaches. This may reflect the release of cool hypolimnetic waters that develop during summer thermal stratification and the release of relatively warm waters during winter that develop due to flow reduction and lack of reservoir shading (Ren et al. 2020). Likewise, the mean duration of riverine heatwaves increased at sites above reservoirs (0.06 d yr À1 ). This increasing duration of heatwaves above reservoirs is likely, in part, due to reduced reservoir discharge and longer hydraulic residence times during periods of drought which have increased in frequency and severity due to climate change (Mosley 2015). Nonetheless, dams provide a management option for reducing extreme high water temperatures in tailwaters by increasing discharge and hypolimnetic release (to the limit of hypolimnetic volume).
We did not observe any significant increasing trends in the frequency, duration, or intensity of riverine heatwaves among regions of the United States. Nevertheless, riverine heatwaves regularly occurred in all regions, with sites  experiencing an average of 2 AE 2 heatwaves yr À1 . Furthermore, riverine water temperature is increasing significantly throughout the United States at a rate of 0.01-0.06 C yr À1 , similar to previous riverine (Kaushal et al. 2010), estuarine (Tassone et al. 2022), and lake studies (O'Reilly et al. 2015), yet exceeding the global ocean (Garcia-Soto et al. 2021).
In inland systems, elevated water temperature increases mobilization of contaminants (i.e., heavy metals, pesticides) within the aquatic food web (Patra et al. 2015) that subsidize riparian and terrestrial consumers (Walters et al. 2008;Moy et al. 2016). As a consequence, fish are likely to carry a larger contaminant burden as freshwater temperatures increase (Patra et al. 2015), leading to stricter advisories on human fish consumption and a greater burden on subsistence fishing communities. Furthermore, extreme high water temperatures interrupt migration (Keefer et al. 2009) and contribute to mass-mortality events (Lynch and Risley 2003;Keefer et al. 2010). Riverine heatwaves represent an emerging climate change threat with potential to impact adjacent ecosystems, exacerbate known stressors, and induce lethal and sublethal damage.

Discharge and precipitation influence heatwaves
Daily mean discharge was often at or below the expected daily mean discharge during riverine heatwaves. Furthermore, the annual total number of heatwaves was greatest at sites with an annual daily mean discharge ≤ 250 m 3 s À1 . These results support model projections and meta-analyses of synchrony between higher water temperature and periods of low flow (Arismendi et al. 2013;van Vliet et al. 2013) and highlight the important effect of discharge on riverine thermal conditions, including extremes. Climate change and anthropogenic activities have altered discharge by: influencing precipitation and low flow timing and intensity, altering snowmelt timing (Mosley 2015;Dudley et al. 2017), increasing water withdrawal, changing seasonal patterns due to impoundment, and altering land cover (Caissie 2006;Jame and Bowling 2020). These collective impacts favor the development of riverine heatwaves.
Normalized annual mean discharge and annual total precipitation were positively correlated across all regions with the northeast (R 2 = 0.49), east north central (R 2 = 0.47), and south (R 2 = 0.43) having the strongest relationships (all other R 2 ≤ 0.40). Annual mean water temperature and normalized annual mean discharge were negatively correlated across all regions, except for Alaska, and were generally weak (R 2 ≤ 0.22) yet significant in the southeast, south, and southwest. This suggests that during years of low precipitation, discharge is reduced, resulting in a lower volume of water that is easier to warm. Climate warming has increased extremes in precipitation such that droughts have become more frequent (Arismendi et al. 2013;Chiang et al. 2021). While we did not consider intra-annual precipitation variability in this analysis, increasing summer synchrony in drought, low discharge, and peak air temperature would increase the risk for riverine heatwave development by promoting warm water conditions.
Increasing riverine water temperature is driven, in part, by rising atmospheric temperature however, water temperature trends at 39% of the sites exceeded the magnitude of atmospheric temperature trends suggesting additional drivers of rising water temperature. Similar results have been reported for some regions of the United States (Isaak et al. 2012;Jastram and Rice 2015), Europe (Michel et al. 2020), Great Britain (Orr et al. 2015), and Japan (Ye and Kameyama 2021) which have implicated hydrology, watershed landcover change, and largescale climatic oscillations as additional drivers of rising river water temperatures. Furthermore, Hare et al. (2021) suggested that United States streams with shallow groundwater contributions, approximately 40% of streams, are warming faster than deep groundwater streams. Our results provide evidence of an indirect relationship between precipitation and water temperature via discharge, and that riverine heatwaves were often associated with normal or below-normal discharge conditions. Interannual variability in precipitation can enhance or reduce heatwave development during wet and dry years respectively with relation to discharge. Similarly, intra-annual variability in precipitation may have contributed to the increasing frequency of summer and fall heatwaves. While atmospheric warming and landcover change will continue to increase riverine water temperature, actions that reduce surface water discharge, such as reduced precipitation, surface water abstraction, and groundwater withdrawal are likely to exacerbate riverine warming trends and heatwave frequency.

Ecosystem management
Water temperature influences many water quality variables which in the United States are regulated by the Clean Water Act. Restoration and antidegradation requirements for water quality criteria will be detrimentally affected by heatwaves potentially resulting in increased violations and failure to reach improvement goals. Broader spatial coverage, greater representation of stream types, and continuous sub-daily long-term data collection will be essential to track heatwave conditions. Furthermore, these improved monitoring efforts will help identify and prioritize restoration of waterways that are exceeding water temperature standards and provide protection for those currently meeting standards (Palmer et al. 2009).

Conclusions
Based on the period covered by the data synthesized in this study, riverine heatwave frequency is increasing throughout the United States, doubling the average annual number of heatwave days between 1996 and 2021. These heatwaves are driven by increases in atmospheric temperature and exacerbated by declines in discharge. Sites most susceptible to riverine heatwaves include those either above a dam or free-flowing, and with an annual daily mean discharge ≤ 250 m 3 s À1 . Similarly, heatwaves were often associated with normal or below-normal discharge conditions. Heatwave frequency significantly increased in summer and fall however, heatwaves were most intense in the winter and spring. Climate models project rising atmospheric temperature throughout the 21st century (Meehl et al. 2007;Wuebbles et al. 2017) along with low discharge conditions (Milly et al. 2005;van Vliet et al. 2013) indicating that riverine heatwaves will continue proliferating. Managing river temperatures will be of increased importance in the future and expanded temperature monitoring should be a priority.