Regional, multi-decadal analysis reveals that stream temperature increases faster than air temperature

Stream temperature appears to be increasing globally, but its rate remains poorly constrained due to a paucity of longterm data and difficulty in parsing effects of hydroclimate and landscape variability. Here, we address these issues using the physically-based thermal model T-NET (Temperature-NETwork) coupled with the EROS semi-distributed hydrological model to reconstruct past daily stream temperature and streamflow at the scale of the entire Loire River basin in France (10 km 5 with 52278 reaches). Stream temperature increased for almost all reaches in all seasons (mean = +0.38 °C/decade) over the 1963–2019 period. Increases were greatest in spring and summer with a median increase of +0.38 °C (range=+0.11– +0.76 °C) and +0.44 °C (+0.08– +1.02 °C) per decade, respectively. Rates of stream temperature increases were greater than for air temperature across seasons for 50–86 % of reaches. Spring and summer increases were typically the greatest in the southern headwaters (up to +1 °C/decade) and in the largest rivers (Strahler order >5). Importantly, air temperature and streamflow 10 exerted joint influence on stream temperature trends, where the greatest stream temperature increases were accompanied by similar trends in air temperature (up to +0.71 °C/decade) and the greatest decreases in streamflow (up to -16 %/decade). Indeed, for the majority of reaches, positive stream temperature anomalies exhibited synchrony with positive anomalies in air temperature and negative anomalies in streamflow, highlighting the dual control exerted by these hydroclimatic drivers. Moreover, spring and summer stream temperature, air temperature, and streamflow time series exhibited common change-points 15 occurring in the late 1980s, suggesting a temporal coherence between changes in the hydroclimatic drivers and a rapid stream temperature response. Critically, riparian vegetation shading mitigated stream temperature increases by up to 16 % in smaller streams (i.e., <30 km from the source). Our results provide strong support for basin-wide increases in stream temperature due to joint effects of rising air temperature and reduced streamflow. We suggest that some of these climate change-induced effects can be mitigated through the restoration and maintenance of riparian forests, and call for continued high-resolution monitoring 20 of stream temperature at large scales. 1 https://doi.org/10.5194/hess-2021-450 Preprint. Discussion started: 6 September 2021 c © Author(s) 2021. CC BY 4.0 License.

and Tw stations, and (bottom) altitude (IGN, 2011) and hydrometric stations. The lithology map was adapted from Moatar and Gailhard (2006), based on original data from BRGM (French Geological Survey). All Tw stations were used to assess the performance of the T-NET thermal model, but only the ones with long-term data (in red) were used to assess the model accuracy for long-term trends. Complete information on Tw stations with long-term continuous data are provided in Table 1. All hydrometric stations were used for calibrating the EROS hydrological model, but only stations with long-term data (in red) were used to assess the model performance on long-term trends.

Method and data
We assessed how Tw responded to climate change over the past 57 years in the Loire River basin in four steps. First, we applied physically based Q and Tw models for the period 1963-2019. Second, we assessed model performance by comparing simulated daily Q and Tw to data from observation stations. Third, we assessed Q and Tw long-term trends. Fourth, we analyzed how 90 hydro-climatic changes and landscape features could explain reconstructed trends in Tw. We performed all analyses on both seasonal and annual basis, where winter is December-February (DJF), spring is March-May (MAM), summer is June-August (JJA), and fall is September-November (SON).

Modeling daily Q and Tw
3.1.1 EROS hydrological model 95 The EROS hydrological model uses daily Ta (°C), solid and liquid precipitation (mm), and reference evapotranspiration (ET 0 , mm) to produce mean daily Q and groundwater flows (Thiéry, 1988;Thiéry and Moutzopoulos, 1995;Thiéry, 2018). Meteorological inputs were provided by the 8 km gridded Safran atmospheric reanalysis data released by Meteo-France over each grid cell.
-Riparian vegetation shading: riparian vegetation is one of the major regulators of shortwave radiation. In the current study, patches of wooded area provided by the BD TOPO® (IGN, 2008) database were used as a proxy of vegetation.
The vegetation species and length of each wooded patch in a buffer of 10 m were extracted for both right and left river banks (van Looy and Tormos, 2013). The vegetation density (vc) was then calculated as the ratio of patch length to 130 reach length for both right and left river banks. In case of multiple wooded patches in any side of a river bank, the average vegetation density of the patches was considered. Then, the model proposed by Li et al. (2012) was used for the calculation of the dynamic shading factor (SF ). The required average tree height for both right and left river banks was estimated based on vegetation species (see Table S1).
In the presence of different vegetation species, the average tree height (m) for each side of a river bank was calculated 135 as follows: with H i and L i the height and length of the tree patch i, respectively, and L the reach length. Next we calculated the proportion of the river width that was shaded (W shaded ) and the dynamic shading factor (SF ) as follows: where H is the average tree height (see Eq. 1), W the river width (see Eq. 7), Ψ the solar altitude angle, δ the angle between solar azimuth and the mean azimuth of T-NET reach, and vc the vegetation density. To take into account 145 the phenology and stages of leaf growth, a coefficient corresponding to each season and transmissivity was applied to SF to calculate the final shading factor: SF f inal = SF × (1 − transmissivity). The transmissivity in leafless months (Jan, Feb, Nov and Dec), months of leaf growth (Mar and Apr), and full-leaf months (May-Sep) was fixed to 0.3, 0.2 and 0, respectively, following Hutchison and Matt (1977). The shortwave radiation was lastly regulated by a factor of

150
-River hydraulic geometry: stream width and depth were calculated using a hydraulic model assuming a rectangular river section being constant over 24 hours (Morel et al., 2020): with f and b being at-a-reach exponents previously predicted by climate, hydrological, topographic and land use descrip-155 tors (see Morel et al., 2020, for more details). Q(t) is the daily mean streamflow provided by the EROS hydrological model. The Q50, W50, D50 (the median of Q, width and height, respectively) and the exponents were available on the Theoretical Hydrographic Network for France (RHT, Pella et al., 2012). There was about 50 % correspondence between the reaches of the T-NET and RHT networks. For the rest of them, required hydraulic geometry variables for the T-NET reaches were extrapolated from the nearest RHT reach. These river hydraulic geometries allowed us to calculate the water velocity by the ratio of Q(t) to rectangular wetted cross-section. The travel time was also defined by the ratio of water velocity to reach length.

T-NET thermal model
The T-NET model does not consider the influence of impoundments on thermal energy balance, and thus produces "natural" thermal regimes. Therefore, model performance was assessed on natural Tw stations, which are weakly influenced by impoundments. 72 near-natural stations with continuous daily data over the 2010-2014 period (see Fig. 1) were identified using the thermal signatures approach that allows distinguishing between natural and altered thermal regimes (see Seyedhashemi   175 et al., 2020, for more details). Of these identified natural reaches, 58 were located on small/medium streams (with distance from the source <100 km), while the remaining are located on large rivers. The mean catchment area of natural stations was 153 km 2 (±225 km 2 ) for small/medium streams, and 16,804 km 2 (±18,833 km 2 ) for large rivers. Seasonal and annual absolute biases were assessed at these 72 stations.
Long-term continuous data was available at 14 of the 72 near-natural stations, including 9 stations with 8-13 years data and 180 5 stations with 20-40 years data. These 14 stations are represented as red points in Fig. 1, middle), and listed in Table 1. The long-term evolution of annual mean Tw at these stations is shown in Fig. S1. These 14 near-natural stations with long-term continuous data comprise the validation dataset for the seasonal and annual trend assessment (see Table 1).

Time series reconstruction and assessment of long-term trends
Daily Q and Tw were reconstructed over the 57-yr period 1963-2019 using the EROS hydrological model and the T-NET and Tw (from T-NET) were reconstructed. Seasonal and annual averages of these 3 variables were considered in the trend assessment.
We estimated the magnitude of trends in these time series with the non-parametric Theil-Sen estimator (Sen, 1968), and evaluated their significance with the Mann-Kendall test (Mann, 1945), commonly used in hydrological analyses (see e.g.

195
Distributions of trends in Tw and Ta were first compared for the whole basin at the seasonal and annual scales using the nonparametric Wilcoxon signed rank test (Bauer, 1972) to determine whether Tw trends were greater than Ta trends. To assess how Ta and Q influenced Tw, we first identified which reaches exhibited 1) jointly increasing trends in Tw and Ta, and 2) jointly increasing trends in Tw and Ta and decreasing trends in Q. Then, for each identified reach, we evaluated the strength and direction of these joint trends using Pearson correlation between 1) Tw and Ta, and 2) Tw and log(Q).

200
Ta, Tw, and Q seasonal and annual anomalies -with respect to the 1963-2019 interannual mean -were then used to synchronicity of extreme years. Change-points in time series of anomalies at each reach were computed with the non-parametric Pettitt test that considers no change in the central tendency as a null hypothesis (Pettitt, 1979). Change points were reported at the 95 % confidence level.

Landscape drivers of Tw trends
Stream size and HER (see Fig. 1) were selected as the major potential landscape drivers. The Strahler order of each reach was used as a proxy for stream size. Reaches with Strahler order 5-8 were combined into a single group termed "large rivers".
The Spearman correlation was computed between median decadal trends in Tw (i.e., median across all reaches) and Strahler order. Such correlations were computed across HERs and at seasonal and annual scales to evaluate the spatial heterogeneity and seasonality.

Performance of models against observations
The EROS model performed well at the annual scale with a median relative bias of 0.37% (see Fig. S2). It slightly underestimated winter Q (median bias=-7.26%) and spring Q (-6.79%), and overestimated summer Q (+37.7%) and fall Q (+24.7%).

220
The mean NSE criteria for low, mean and high flows were 0.82, 0.85 and 0.79, respectively. No systematic bias was observed for modeled Tw at the stations located on small and medium rivers (see Fig. S2). Median Tw bias ranged from -0.26°C (in fall) to 0.8°C (in winter). Large rivers exhibited a small Tw underestimation, with a median bias ranging from -0.29°C (in fall) to +0.15°C (in winter), and the overall biases were still quite small across seasons (IQR=0.4-0.7°C).

Spatial reconstruction of long-term trends
Tw increased in almost all modeled reaches for all seasons (Fig. 3, left; Fig. S5). Depending on the season considered, 62 % 240 to 80 % of reaches showed trends in the range of +0.2-+0.4°C/decade (i.e 1.14-2.28°C over the whole 1963-2019 period).
Summer Tw trends were more spatially variable than in other seasons, with more than 50 % of reaches showing values higher than +0.4°C/decade (see Fig. S6). Such reaches were mainly located in the southern part of the basin, in HER A (see Fig. 3, left). Spring Tw trends showed a similar spatial pattern, but with lower trend values.
Likewise, Ta exhibited increasing trends for 99 % of all reaches across spring, summer, and the whole year (Fig. 3, middle).

245
Values were mainly in the range of +0.2-+0.4°C/decade (see Fig. S6). The highest Ta trend values were found in summer (resp. spring), when 67 % (resp. 22 %) of reaches showed values higher than +0.4°C/decade. Such reaches were mainly located in HER A, especially in summer. Non-significant trends were found over the whole basin in winter, and in the southern part of the basin in fall.
In contrast to Ta and Tw, trends in Q were highly variable in magnitude and direction across the basin and across seasons 250 (Fig. 3, right), and most were not significant at p=0.05 (Fig. S5). However, significant decreasing trends were found in the southern headwaters (HER A) in spring, summer, fall, and annual scale (Fig. S5). Decreasing trends were observed for the majority of reaches across seasons (66-83% of reaches), with the exception of winter (37%). Decreasing trends could have magnitudes greater than -5 % per decade, implying a -28 % loss in Q over the whole 1963-2019 period.

Hydroclimatic drivers of Tw trends
255

Relationships between trend estimates
The median of Tw trends were higher than that of Ta trends for every season (p<0.001 according to the Wilcoxon signed rank test), except for summer when the median trend values for Tw and Ta were very similar, but more variable for Tw (+0.08-+1.02°C/decade) (Fig. 4). The greatest increase in Tw was found in summer (+0.44°C/decade). Overall, Tw trends were more spatially variable than Ta trends, suggesting the conditional influence of Q trends (Fig. 4). Indeed, where Tw trends

Synchrony of annual anomalies 265
We observed strong positive correlations between seasonal and annual averages of Tw and Ta across seasons (see Table 2). We further observed strong negative correlation between summer Tw and Q time series. The strong positive correlation observed between winter Ta and Q results from the particular hydroclimate of winters in the Loire basin, where they are either either warm and wet or cold and dry. Annual anomalies of Tw, Ta, and Q exhibited variable patterns, with Tw and Ta generally increasing and Q remaining 270 relatively stationary (Fig. 6). Tw anomalies were generally more variable than Ta, especially in summer, but both time series appeared to exhibit synchronous behavior. Change-point analysis supported this visual observation, where change-points in seasonal and annual averages were largely coincident across these time series (Fig. 7). Tw and Ta   and 1994 for both Tw and Ta. The significant change-points in seasonal Q time series were detected for a substantially smaller proportion of reaches, e.g. less than 40 % of reaches for spring and summer. The majority of these reaches were located in HER A across seasons (66-86 % of such reaches), with the exception of winter (49 %) (see Fig. S8). In spring and summer, 280 they occurred in the late 1980s, similarly to Tw and Ta. Conversely, the significant change-points detected in other seasons were much more scattered in time, probably due to the high interannual variability of Q.
Critically, the largest summer Tw and Ta positive anomalies over the study period were observed in 1976,2003,2006,2015,2017, and 2019, which corresponded to years with the largest negative anomalies in summer Q (Fig. 8). Note that this signal was much less clear for the other seasons (see Fig. S9).

Stream size and HER
Strahler order was strongly and positively correlated with Tw trends in spring and summer, but this effect was modulated by HER ,with HER A exhibiting the strongest correlations (Fig. 9). In other words, larger rivers tended to exhibit the largest increases in spring and summer Tw, especially for reaches located in the HER A. There, median trends in spring (resp. summer) 290 ranged from +0.37 (resp. +0.49)°C/decade for small streams (Strahler order 1) to +0.55 (resp. +0.64)°C/decade for large streams (Strahler order ≥ 5).

Riparian shading
For small streams, i.e. reaches located closer than 30 km from the source, the shading factor SF and median trends in Tw were significantly and negatively correlated for all three HERs in spring, as well as for HER A in summer and at annual scale in the HER A, B and C decreased by 12 %, 5 % and 4 %, respectively, from sparsely shaded reaches (SF < 15%) to highly shaded reaches (SF > 40%). For summer Tw in HER A, the median trends were 16 % lower between the lowest and highest 300 shaded reaches.

Quality and suitability of simulated Tw and Q
Although some biases were observed for both Q and Tw (Fig. S2), we found high correlations between modeled and observed trends in seasonal and annual Q, with the exception of summer (Fig. S3), and Tw, with the exception of fall (Fig. S4), demon-305 strating the adequacy of the modelling approach for temporal trend analyses. For Q, although largely non-significant (Fig. S5), decreasing trends were observed for the majority of the reaches across seasons (66-83% of reaches, with the exception of winter, Fig. S7) strongly suggested an effect on Tw. Moreover, the spatial pattern in simulated Q trends, with significant decreases knowledge, the present study was one of the very few studies used modeled Tw to investigate past trend at a large scale (but see Isaak et al., 2017Isaak et al., , 2012Van Vliet et al., 2011).  Global-scale stream temperature modelling suggests trends in annual averages ranging from +0.2 to +0.5°C/decade over France (Wanders et al., 2019), which is consistent with our findings (mostly in the range of +0.2-+0.4°C/decade, Fig. S6).
We found more pronounced trends in spring and summer, which was also found in other parts of Europe (e.g Kędra, 2020; Arora et al., 2016;Michel et al., 2020). Considerable discrepancies were also found between Tw and Ta trends across seasons for the majority of the reaches (see Fig. 3, and 5), which is a common finding for other sites around the world (Arora et al.,

325
2016; Wanders et al., 2019). This highlights that changes in Ta may not be the only driver of changes in natural Tw.

Drivers and spatial patterns of trends in Tw
The greatest increases in Tw (up to +1°C/decade) were predominately located in the southern headwaters of the basin, in HER A (Massif Central) where a greater increase in Ta (up to +0.71°C/decade) and a greater decrease in Q (up to -16 %/decade) occurred jointly. The decrease in Q could be due to a significant increase in potential evapotranspiration (PET) (up to +10 %) 330 over the whole of basin and a decrease in total precipitation (P) (up to -5 %/decade) (Figs. S10 and S11). Such trends, computed here based on variables from the Safran surface meteorological reanalysis (Vidal et al., 2010), are consistent with larger-scale studies (see e.g. Spinoni et al., 2017;Tramblay et al., 2020;Hobeichi et al., 2021). Moreover, Vicente-Serrano et al. (2019) attributed annual streamflow trends in southern France mostly to trends in precipitation and potential evapotranspiration, as opposed to irrigation and land-use changes that have additional strong effects e.g. in the Iberian peninsula. In Switzerland, 335 Michel et al. (2020) described an increase of +0.33 ± 0.03°C/decade in Tw, resulting from the joint effects of increasing in Ta (+0.39 ± 0.14°C/decade) and decreasing in Q (−10.1 ± 4.6%/decade) over the 1979-2018 period. In contrast with our results, they found Tw trends lower than Ta trends due to influence of snow melt and glacier melt. Trends in Tw might also be explained by trends in additional drivers, like shortwave radiation (Wanders et al., 2019), which is the dominant flux at air-water surface, and is notably increasing over Europe (Sanchez-Lorenzo et al., 2015). This might explain discrepancies between Tw and Ta 340 trends in spring and summer, when no decreasing trend in Q are found (see Fig. 3).
Strong synchronicity between Ta and Tw anomalies was observed in the present study in the warmest years, and these years were also among those with the largest negative Q anomalies (see Fig. 8). Indeed, increase in summer Tw could be due to cooccurrence with the increase in summer Ta (average correlation: +0.82), and with decrease in summer Q (average correlation:  Moatar and Gailhard (2006) found that the increase in Ta (resp. decrease in Q) explain 60 % (resp. 40 %) of the increase in Tw. Moreover, the significant change-point in Tw, Ta and Q time series in the late 1980s has also been found in other studies in Europe (Moatar and Gailhard, 2006;Arora et al., 2016;Zobrist et al., 2018;Ptak et al., 2019b;Michel et al., 2020). Long-term observational time series like that of the Loire at Dampierre also display a similar change-point.

Natural trends and anthropogenic influence on Tw
Natural Tw time series were used in the current study for detecting trends, as both EROS and T-NET models are used in a non-influenced set-up (see Sect. 3.1). However, anthropogenic impoundments (e.g. large dams, small reservoirs, and ponds) influence downstream Tw regimes in a diversity of ways that depend on their structure and position along the river continuum (Seyedhashemi et al., 2020). In this regard, on the one hand, large dams, by releasing cold hypolimnetic water in summer, can 355 lower downstream Tw (Olden and Naiman, 2010), and mitigate increasing trend in Tw (Cheng et al., 2020). Nevertheless, it is anticipated that a considerable proportion of streams regulated by large reservoirs may still experience high temperatures and low flows under future climate change (Cheng et al., 2020). The mitigating influence of dams could be of importance for streams in the southern headwaters of the Loire basin since this area both experienced the greatest Tw trends and gathers most of existing large dams. On the other hand, ponds and shallow reservoirs, by releasing warm water can increase downstream 360 Tw (Seyedhashemi et al., 2020;Zaidel et al., 2020;Chandesris et al., 2019) and exacerbate increasing trends in Tw (Wanders et al., 2019;Michel et al., 2020). Nevertheless, the warming effect can be local, and unregulated streams being located closer to such regulated streams may show limited to no warming (Wanders et al., 2019). The warming effect of such surface waters in the current study seem more significant for streams located in lowlands in the middle and north of the Loire River basin where most of the shallow reservoirs are located. In these streams, anthropogenically-induced trends in Tw may be greater than 365 natural ones, and the warming process can get worse through the increasing demand for storing water in small reservoirs for irrigation.

Implications for river management and aquatic biota
The removal of riparian vegetation can increase Tw (Caissie, 2006), and changes in Tw can be even more sensitive to changes in riparian vegetation than to changes in Ta or Q (Wondzell et al., 2019). We showed that in small streams, an increase of > 370 25 % of riparian shading could decrease the median trend in spring and summer Tw by 5-16 % (Fig. 10). Spring and summer Tw trends were more pronounced in large rivers, especially in the south of the basin, with a difference in median Tw trends of up to +0.18°C/decade (Fig. 9), probably due to decrease in Q (up to -2 %/decade, see Fig. S12), greater thermal sensitivity, and the absence of mitigating factors like riparian vegetation shading or groundwater inputs (Kelleher et al., 2012;Beaufort et al., 2020).

375
Restoring riparian vegetation and shading can therefore substantially mitigate future increases in Tw. In addition, riparian restoration may lessen the impacts of climate change on flood damage to human infrastructure, on riparian biodiversity, on ecosystem vulnerability and on changes in Q (Palmer et al., 2009;Seavy et al., 2009;Perry et al., 2015). However, riparian restoration is not an easy task since the survival, persistence, growth rate of planted species as well as required time for thermal regime recovery under possibly severe future conditions should be studied beforehand (Perry et al., 2015). For instance, it may take between 5 and 15 years for rivers to recover their thermal regime following vegetation growth (Caissie, 2006;Edmonds et al., 2000). Moreover, the efficacy of riparian planting is also highly dependent upon the type and structure of forest stands (Dugdale et al., 2018), and this should also be considered in long-term projects.
Stream warming affects cold-water fish populations negatively at the warmer boundaries of their habitat (Hari et al., 2006). Furthermore, changes in spawning and migration timing (McCann et al., 2018;Arevalo et al., 2020), decreases in habitat 385 availability and freshwater quality for organisms (Lennox et al., 2019), and shifts in species distribution (Comte et al., 2013) are already observed consequences of the long-term increase in Tw. Some major changes in fish density and community structure has already been reported in large rivers over France (Maire et al., 2019) for which we also found greater trends in Tw compared to small ones. Therefore, physical process-based thermal models like T-NET can also be used to assess the various stresses on freshwater habitat sustainability due to changes in Q and Tw (Morales-Marín et al., 2019).

Conclusions
Regional trends in Tw at the reach resolution were detected and assessed by using the T-NET physical process-based model coupled with the EROS hydrological model over the Loire basin. Using model outputs across 52 278 reaches over the Loire basin, for 3 variables (Ta, Q, and Tw), and 5 time scales (seasons plus annual), we found consistent increasing Tw trends at the scale of the entire Loire River basin, regardless of the season. Critically, the rate of warming for stream temperature was 395 in the majority of cases higher than the rate of atmospheric warming, suggesting a decoupling of thermal trajectories linked to decreasing Q, especially in the southern headwaters. Moreover, Tw trends in all seasons except winter were greater in rivers with Strahler order > 5, which we attributed to the mitigation effect of riparian shading for large rivers.
The synchronicity of extreme events of low flows and high air temperature in the southern headwaters will likely generate a double penalty for existing cold-water aquatic communities. However, riparian shading in small mountainous streams may 400 mitigate such warming. These findings underscore that Ta alone is likely not an adequate proxy to explain stresses and shifts experienced by aquatic communities over time and space, especially in regions with more pronounced stream warming, and thus there is a need to grow and maintain Tw sensor networks. This knowledge can help develop appropriate management strategies to conserve thermal refugia and mitigate extreme thermal events induced by climate change. Agence de l'eau Loire-Bretagne (AELB) and EDF (Hynes Team). The SAFRAN database was provided by the French national meteorological service (Météo-France). Stream temperature data were provided by the Office français de la biodiversité (OFB), Fédération de Pêche (fishing federation). We are grateful to Electricité de France (French Electricity, EDF) for providing long-term data on observed stream temperature and naturalized streamflows. We also thank André Chandesris for his assistance in the incorporation of vegetation cover into the