Regional, multi-decadal analysis on the Loire River basin 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 long-term 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² with 52 278 reaches).

Stream temperature increased for almost all reaches in all seasons (mean $=+\mathrm{0.38}$ ^∘C decade−1) over the 1963–2019 period. Increases were greatest in spring and summer, with a median increase of + 0.38 ^∘C (range $=+\mathrm{0.11}$ to +0.76 ^∘C) and +0.44 ^∘C (+0.08 to +1.02 ^∘C) per decade, respectively. Rates of stream temperature increases were greater than for air temperature across seasons for the majority of reaches. Spring and summer increases were typically greatest in the southern part of the Loire basin (up to +1 ^∘C decade⁻¹) and in the largest rivers (Strahler order ≥5). Importantly, air temperature and streamflow could exert a 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 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 0.16 ^∘C decade⁻¹ 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.

1 Introduction

Stream temperature is a critical water quality parameter affecting the distribution of aquatic communities (Poole and Berman, 2001; Ducharne, 2008), but its future under global change remains uncertain. As air temperature (T_a) increases worldwide due to climate change, stream temperature (T_w) is expected to follow a similar trajectory (Mohseni et al., 1999; Kaushal et al., 2010; Van Vliet et al., 2011; Isaak et al., 2012; Arora et al., 2016). Indeed, there is growing evidence that stream warming is occurring around the world, affecting freshwater ecosystems through structural and functional changes in biological communities throughout the food web (Woodward et al., 2010; O'Gorman et al., 2012; Scheffers et al., 2016). Deleterious warming effects are documented from bottom-dwelling microorganisms (e.g. Romaní et al., 2016; Majdi et al., 2020) up to macroinvertebrates (e.g. Floury et al., 2013; Bruno et al., 2019) and fish communities (e.g. Maire et al., 2019; Stefani et al., 2020). However, the paucity of long-term time series of T_w (Webb and Walling, 1996; Nelson and Palmer, 2007; Webb et al., 2008; Arora et al., 2016) has impaired the larger-scale assessment of such trends, especially in light of confounding factors like hydrological changes and land-use change. Hence, analyses of T_w trends, especially at large spatiotemporal scales, remain scarce (but see Kaushal et al., 2010; Orr et al., 2015; Arora et al., 2016; Michel et al., 2020; Wilby and Johnson, 2020).

To overcome the lack of T_w data, large-scale ecological studies typically use T_a as a proxy for T_w to assess the impact of climate change on the spatial distribution of aquatic organisms (e.g. Buisson et al., 2008; Buisson and Grenouillet, 2009; Tisseuil et al., 2012; Domisch et al., 2013), but T_a can overpredict changes to stream fish assemblages with climate warming compared with T_w (Kirk and Rahel, 2022). Indeed, T_a can be an imprecise surrogate for T_w (Caissie, 2006), since many landscape and basin characteristics (e.g. stream discharge (Q), streambed morphology, karst resurgences, topography, and vegetation cover) contribute to the response of T_w to climate change over time and space (Stefan and Preud'homme, 1993; Webb and Walling, 1996; Webb et al., 2008; Hannah and Garner, 2015). For instance, riparian vegetation can obstruct solar radiation, which is the dominant heat flux at the air–water surface (Hannah et al., 2004; Caissie, 2006), and therefore decrease T_w response to T_a (Johnson, 2004; Loicq et al., 2018). However, while riparian vegetation shading can greatly decrease the temperature of small rivers (Dan Moore et al., 2005; Loicq et al., 2018), it has limited effects on larger rivers since the width of such rivers is large enough that only a small part of it can be shaded. Rising groundwater temperature (Taylor and Stefan, 2009; Kurylyk et al., 2013, 2014) and reduced groundwater flows (Kurylyk et al., 2014) due to climate change may further contribute to T_w trends (Meisner, 1990; Arora et al., 2016), leading to asymmetric controls (vis-à-vis T_a) on T_w (Moatar and Gailhard, 2006), especially in headwaters (Caissie, 2006; Kelleher et al., 2012; Mayer, 2012). Finally, intensification of the water cycle (Huntington, 2006), with more frequent and severe droughts (Mantua et al., 2010; Giuntoli et al., 2013; Prudhomme et al., 2014), as well as more intense and sudden floods (Blöschl et al., 2019) may decouple T_a–T_w trends, exacerbating T_w increases that will most likely be evident during low summer flows when thermal inertia and flow velocity are at their minima (Webb, 1996; Webb et al., 2008).

There is thus a clear need to improve our estimates of T_w trends to assess how stream ecosystems will respond to the climate change. Unfortunately, extrapolating trend estimates derived from short time series may lead to paradoxical results, e.g. cooling streams in a warming world (Arismendi et al., 2012). This discrepancy in short- and long-term dynamics is likely due to confounding influences of T_a and hydrology, with implications for the persistence of specialized aquatic organisms (e.g. for cold-water biota, Arismendi et al., 2013b) and the completion of their life cycle (e.g. for diadromous fish, Arevalo et al., 2020). Hence, from an ecological perspective, it will be critical to understand and deconvolve the joint influences of changing T_w and Q regimes. In the absence of more robust data sources, modelling is thus an indispensable tool in meeting these goals.

T_w models output data sets that can then be used to assess the magnitudes of long-term trends, but model selection entails important considerations. For example, T_w can be estimated by developing a statistical, or stochastic, model based on multiple independent drivers (Benyahya et al., 2007), which is a common practice for large-scale studies (e.g. Mantua et al., 2010; Jackson et al., 2017, 2018).

However, these statistical models lack mechanisms; they cannot reveal the specific energy transfer mechanisms responsible for the spatiotemporal patterns of T_w (Dugdale et al., 2017). They are also unable to predict T_w for periods other than those used for their calibration due to, for instance, the non-stationary relationship between T_a and T_w over time (Arismendi et al., 2014). Alternatively, physically based, or deterministic, models are entirely mechanistic: they predict T_w dynamics through a heat budget, accounting for energy exchanges and effects of landscape characteristics on energy transfer (Sinokrot et al., 1995; Webb and Walling, 1997; Yearsley, 2009; van Vliet et al., 2013). Critically, such process-based models can be used not only to reconstruct past time series, but they can also be used in forecasting or in predicting T_w response to projected climate or land-use changes (Caissie et al., 2007; Dugdale et al., 2017).

Here, we used a physical process-based thermal model coupled with a semi-distributed hydrological model to understand how T_w has responded to recent climate change at a large scale. To do so, we first assessed the performance of the models against field observations over the Loire basin, France, and then reconstructed daily Q and T_w over the past 57 years over the whole hydrographic network. We then used reconstructed time series to compute the magnitude of decadal trends in seasonal and annual T_w, T_a, and Q. To understand the relative influences of T_a and Q (as the main hydroclimatic drivers of T_w) on T_w, we compared their trends and assessed their spatial and temporal links. Finally, we sought to understand variation in T_w trends as a function of stream size, landscape attributes, and riparian shading.

2 Study area

The Loire River basin is one of the largest in Europe (10⁵ km²), encompassing an area with starkly contrasting HydroEco Regions (HERs), land use/land cover, and climatic conditions (Moatar and Dupont, 2016), providing an ideal case study to disentangle the drivers of the spatial heterogeneity of trends in T_w. Mean annual precipitation (549–2130 mm), mean annual potential evapotranspiration (PET) (550–850 mm), mean annual T_a (6.0–12.5 ^∘C) (see Fig. S1), and altitude (0–1850 m) (see Fig. 1, right) provide spatially variable controls on stream thermal regimes. Wasson et al. (2002) identified several HERs over France, representing homogeneous areas in terms of land use/land cover, geology, and climate conditions. A grouping over the Loire basin was done to identify three contrasted HERs (identified as A, B, and C in Fig. 1) that can be used to describe the spatial heterogeneity in T_w at the basin scale. Granite and basalt dominate the southern headwaters of the basin (mostly in the Massif Central, HER A), whereas sedimentary rocks occupy the middle reaches with a potential for groundwater input (HER B), followed by granite and schist in the lower reaches (HER C) (see Fig. 1, left). The percentage of riparian vegetation cover (mean over both sides of a river bank at a buffer of 10 m, Valette et al., 2012) is more important in HER A (median = 73 %) and in HER B (median = 68 %) (Fig. 1, middle). In HER C, the presence of riparian vegetation is quite heterogeneous (median = 50 %).

Figure 1 Maps of the Loire basin. Left panel: basin lithology adapted from Moatar and Dupont (2016), based on original data from BRGM (French Geological Survey). The figure in the top right corner of the left panel shows the position of the Loire River basin in France. Middle panel: riparian vegetation cover (provided by Valette et al., 2012) and T_w stations. All T_w stations were used to assess the performance of the T-NET thermal model (see Table S5), but only the ones with long-term data (in red) were used to assess the T-NET thermal model accuracy for long-term trends. Complete information on T_w stations with long-term continuous data is provided in Table 1. The numbers in red in the middle panel correspond to the ID of long-term T_w stations in Table 1. Right panel: altitude (IGN, 2011) and hydrometric stations (extracted from the French national Banque Hydro database: http://www.hydro.eaufrance.fr/, last access: 12 July 2021). All hydrometric stations were used for calibrating the EROS hydrological model (see Table S2), but only stations on the French Reference Hydrometric Network (RHN) with long-term data (in red) were used to assess the EROS hydrological model performance on long-term trends (see Table S4).

3 Method and data

We assessed how T_w responded to climate change over the past 57 years in the Loire River basin in four steps. First, we applied physically based Q and T_w models for the 1963–2019 period. Second, we assessed both hydrological and thermal model performance by comparing simulated seasonal and annual Q and T_w to data from observation stations. Third, we assessed Q and T_w long-term trends. Fourth, we analysed how hydroclimatic changes and landscape features could explain reconstructed trends in T_w. We performed all analyses on both seasonal and annual bases, where winter is December–February (DJF), spring is March–May (MAM), summer is June–August (JJA), and autumn is September–November (SON). Analyses were also performed by HER (see Fig. 1 and Sect. 2) to look at the large-scale influence of landscape characteristics.

3.1 Modelling daily Q and T_w

We used two models to calculate T_w in the Loire River basin according to the method developed by Beaufort et al. (2016). The first model is the EROS semi-distributed hydrologic model, which estimates daily Q at sub-basin outlets. The second model is the fully distributed, mechanistic temperature model called Temperature-NETwork (T-NET; Beaufort et al., 2016; Loicq et al., 2018) that uses Q from EROS and meteorological reanalysis data from Safran atmospheric reanalysis data (Quintana-Segui et al., 2008; Vidal et al., 2010) to estimate T_w at each reach in the Loire River basin. These models are briefly described below and are detailed in Thiéry (1988), Thiéry and Moutzopoulos (1995), and Thiéry (2018) for EROS and Beaufort et al. (2016) and Loicq et al. (2018) for T-NET.

3.1.1 EROS hydrological model principles and input data

The EROS semi-distributed hydrological model simulates daily Q at the outlet of 368 sub-basins (ranges between 40 and 1600 km²; mean drainage area = 300 km²), designed to be as homogeneous as possible with respect to land use and geology. On each of these sub-basins, the water balance is modelled by a lumped model using three reservoirs (see Fig. S2 in the Supplement) and a routing function for propagation across sub-basins (Thiéry, 1988; Thiéry and Moutzopoulos, 1995; Thiéry, 2018). Water abstractions are not considered in EROS, and it produces natural Q. This hydrological model has already been used in several other studies, including climate change impact studies (Ducharne et al., 2011; Habets et al., 2013; Bustillo et al., 2014).

EROS used daily T_a (^∘C), solid and liquid precipitation (mm), and reference evapotranspiration (ET₀, 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 Météo-France over the 1958–2019 period (Quintana-Segui et al., 2008; Vidal et al., 2010). ET₀ was computed from Safran variables with the Penman–Monteith equation (Allen et al., 1998).

3.1.2 T-NET thermal model principles and input data

The T-NET thermal model is a fully mechanistic, 1D model that simulates hourly ${T_{\mathrm{w}}}_{i,j}$ at distance, i, along reach, j, by solving the local heat budget. The heat budget of each reach includes six fluxes (W m⁻²): net solar radiation, atmospheric longwave radiation, longwave radiation emitted from the surface water, evaporative heat flux, convective heat flux, and groundwater flux. Briefly, the model calculates the longitudinal change in T_w at time t (d$T_{\mathrm{w}}/\mathrm{d}x$) for steady-state conditions in order to achieve thermal equilibrium (i.e. $\mathrm{\Sigma }{H}_{i,j}=\mathrm{0}$, where H is a heat flux) while accounting for confluence mixing. Detailed information about the T-NET model principles and calculation of the six heat fluxes at the water–air and water–stream bed interfaces and thermal propagation was provided in Beaufort et al. (2016) and Loicq et al. (2018).

The hydrographic network of the model over the Loire basin consists of 52 278 reaches delimited either by confluences of two streams or a headwater source (i.e. first-order reaches) (Beaufort et al., 2016; Loicq et al., 2018). The mean reach length is 1.7 km, and 74 % of the reaches have a Strahler order lower than 3 (see Fig. S3). To compute the six heat fluxes and the water travel time for each reach, the following input data were used.

• Meteorological variables: hourly T_a (^∘C), specific humidity (g kg⁻¹), wind velocity (m s⁻¹), shortwave radiation (W m⁻²), and longwave radiation (W m⁻²) were provided by the 8 km gridded Safran atmospheric reanalysis (Vidal et al., 2010). All reaches within a grid cell were assigned the values of the meteorological variables in that grid cell. For reaches flowing through more than one grid cell, meteorological variables were weighted by the relative length of the reach within 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 the right and left river banks (van Looy and Tormos, 2013). The vegetation density (vc) was then calculated as the ratio of patch length to reach length for both the right and left river banks. In the case of multiple wooded patches on 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) at an hourly time step. The required average tree height for both the 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 as follows:

  $$\begin{array}{}\text{(1)}& H={\displaystyle \frac{\mathrm{1}}{n}}{\mathrm{\Sigma }}_{i=\mathrm{1}}^n{H}_i{\displaystyle \frac{L_i}{L}},\end{array}$$

  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 SF as follows:

  $$\begin{array}{}\text{(2)}& {\displaystyle }{W}_{\mathrm{shaded}}\left(t\right)={\displaystyle \frac{H_{\mathrm{left}/\mathrm{right}}\times \cot \mathrm{\Psi }\left(t\right)\times \sin \mathit{\delta }\left(t\right)}{W\left(t\right)}},\text{(3)}& {\displaystyle }{\mathrm{SF}}_{\mathrm{right}}\left(t\right)=\left(W_{\mathrm{shaded}}{)}_{\mathrm{right}}\right(t)\times (\mathrm{vc}{)}_{\mathrm{right}},\text{(4)}& {\displaystyle }{\mathrm{SF}}_{\mathrm{left}}\left(t\right)=\left(W_{\mathrm{shaded}}{)}_{\mathrm{left}}\right(t)\times (\mathrm{vc}{)}_{\mathrm{left}},\text{(5)}& {\displaystyle }\mathrm{SF}\left(t\right)=\mathrm{Max}\left({\mathrm{SF}}_{\mathrm{right}}\right(t),{\mathrm{SF}}_{\mathrm{left}}(t\left)\right),\end{array}$$

  where H is the average tree height (see Eq. 1), W is the river width (see Eq. 7), Ψ is the solar altitude angle, δ is the angle between solar azimuth and the mean azimuth of T-NET reach, and vc is the vegetation density. Loicq et al. (2018) already tested this method and showed that T_w simulated by the above method was close to the observed T_w. To take into account the phenology and stages of leaf growth, a coefficient corresponding to each season and transmissivity was applied to the SF to calculate the final shading factor: SF${}_{\mathrm{final}}=\mathrm{SF}\times (\mathrm{1}-\mathrm{transmissivity}$). The transmissivity in leafless months (January, February, November, and December), months of leaf growth (March and April), and full-leaf months (May–September) 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 1−SF_final.

• Reach streamflow: the daily Q simulated at the outlet of 368 homogeneous sub-basins by the EROS model (see Sect. 3.1.1) was redistributed along the river network inside each sub-basin according to the reach drainage area for informing the T-NET model at the reach scale. The ratio of the sum of the lengths of all reaches upstream of a reach to the sum of the lengths of all reaches located in a sub-basin was used as a proxy for the drainage area of a reach. To have Q at an hourly time step, Q was assumed to be constant over 24 h.

• River hydraulic geometry: stream depth and width were calculated using a hydraulic model assuming a rectangular river section being constant over 24 h (Morel et al., 2020):

  $$\begin{array}{}\text{(6)}& {\displaystyle }D\left(t\right)=D\mathrm{50}\times {\left[{\displaystyle \frac{Q\left(t\right)}{Q\mathrm{50}}}\right]}^f,\text{(7)}& {\displaystyle }W\left(t\right)=W\mathrm{50}\times {\left[{\displaystyle \frac{Q\left(t\right)}{Q\mathrm{50}}}\right]}^b,\end{array}$$

  with f and b being at-a-reach exponents previously predicted by climate, hydrological, topographic, and land use descriptors (see Morel et al., 2020, for more details). Q(t) is the daily mean streamflow provided by the EROS hydrological model. 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 reach length to water velocity.

3.2 Model assessment and validation

3.2.1 Calibration and assessment of the EROS hydrological model

For calibration, 352 hydrometric stations with observed Q data (see Table S2) extracted from the French national Banque Hydro database (http://www.hydro.eaufrance.fr/, last access: 12 July 2021) were used (Fig. 1, right). Stations along the main Loire (14 stations) and Allier (11 stations) rivers were influenced by operations of four large dams for hydropower, flood control, and low-flow management, notably through summer releases (see Table S3). Time series at these stations were first naturalized by EDF (Électricité de France) based on variations in water storage in reservoirs due to operations based on target storage curves (Naussac, Villerest) and optimization of the hydropower energy generation for the electric grid needs (Grangent, Montpezat). EROS was then calibrated over the 1971–2018 period (which maximizes the number of available streamflow observations) against daily Q at all 352 sub-basins with at least 10 years of daily observations. Note that, in the current study, calibrating EROS by performing a classical split-sample test was not possible due to the numerous and diverse gaps in streamflow observations at this large scale. Therefore, calibration was done on the complete set of available information over the past decades, which is a standard approach for calibrating a hydrological model in climate change impact studies (e.g. Habets et al., 2013; Vidal et al., 2016).

The calibration aimed at optimizing all unknown parameters (soil capacity, recession times, and propagation times) by maximizing the Nash–Sutcliffe efficiency (NSE) criterion on the square root of streamflow and minimizing the overall bias. Considering the NSE criterion on the square root of the flows provided an estimate of model performance without favouring either high flows or low flows. The overall calibration criterion, C, was $C=\mathrm{mNSE}-w\times \mathrm{mRB}$, where mNSE is the mean NSE, mRB is the mean relative bias, and w is a weighting factor fixed at 0.05. A 3-year warm-up period (1971–1974) was discarded from the overall calibration period for the next assessments.

The calibrated model was then used to simulate streamflow over all 368 homogenous sub-basins in the Loire basin over the whole 1963–2019 period. Note that, although meteorological variables were available from 1958 onwards, the first years (1958–1962) were discarded from the analysis to ensure the EROS model convergence.

Of the 352 calibration stations considered by EROS, 44 were part of the French Reference Hydrometric Network (RHN) described by Giuntoli et al. (2013) with long-term continuous high-quality data over the 1968–2019 period, especially for low flows (see Table S4). These stations are shown with red points in Fig. 1, right. Seasonal and annual relative biases, together with Nash–Sutcliffe efficiency on Q, ln (Q), and $\sqrt{Q}$ were computed on both 352 calibration (over the 1975–2018 period) and 44 RHN stations (over the 1968–2019 period) to provide an overview of the performance of the EROS model. Moreover, to assess EROS model ability to capture long-term trends, decadal trends (in percent per decade) over the 1963–2019 period on seasonal and annual averages from the EROS simulation were compared with corresponding observed trends at each of the 44 RHN stations with long-term continuous daily data.

3.2.2 Validation of the T-NET thermal model

T-NET does not consider the influence of impoundments on thermal energy balance and thus produces “natural” thermal regimes. Therefore, model performance was assessed on near-natural T_w stations, which are also weakly influenced by impoundments, i.e. at 67 stations with continuous daily data over the 2010–2014 period (see Fig. 1, middle, and Table S5). These stations were identified using the thermal signatures approach that allows one to distinguish between natural and altered thermal regimes (see Seyedhashemi et al., 2020, for more details). Of these identified natural stations, 55 were located on small/medium streams (with distance from the source < 100 km), while the rest were located on large rivers. The mean catchment area of natural stations was 151 km² (range = 7–1342 km²) for small/medium streams and 18 926 km² (range = 1931–57 043 km²) for large rivers. T_w data at these 67 stations were provided by different organizations: EDF, OFB (Office français de la biodiversité), DREAL (Direction régionale de l'environnement, de l'aménagement et du logement), and FD (Fédération Nationale de la Pêche) (see Tables 1 and S5). T_w data provided by OFB can be downloaded from http://www.naiades.eaufrance.fr/ (last access: 1 August 2020). The rest of the T_w data are not archived for public use. Seasonal and annual absolute biases and root mean square error (RMSE) were assessed at these 67 stations. Note that, unlike the EROS hydrological model, the T-NET thermal model does not have any free parameter, and hence it did not require calibration.

Long-term continuous data were available at 14 of the 67 near-natural stations, including 9 stations with 8–13 years of data and 5 stations with 20–40 years of data. These 14 stations are represented as red points in Fig. 1, middle panel, and listed in Table 1. The long-term evolution of annual mean T_w at these stations was shown in Fig. S4. These 14 near-natural stations with long-term continuous data comprised the validation data set for the seasonal and annual trend assessment (see Table 1).

Table 1 Characteristics of the 14 long-term T_w stations. See Fig. S4 for the evolution of observed annual T_w at these stations. The IDs of stations correspond to the numbers (in red) shown in Fig. 1, middle. The coordinate system is Lambert 93.

3.3 Time series reconstruction and assessment of long-term trends

Daily Q and T_w were reconstructed over the 57-year period 1963–2019 using the EROS hydrological model and the T-NET thermal model. For each of the 52 278 reaches, daily time series of T_a (from the Safran reanalysis), Q (from EROS), and T_w (from T-NET) were reconstructed. Seasonal and annual averages of these three variables were considered in the trend assessment. Maps of annual mean T_w, T_a, and Q over the 1963–2019 period were presented in Fig. S5 with annual T_w > 11.5 ^∘C for large rivers (annual T_a > 10.5 ^∘C and log (annualQ)>1 m³ s⁻¹).

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. Giuntoli et al., 2013) but also thermal analyses (e.g. Kaushal et al., 2010; Arismendi et al., 2013a; Arevalo et al., 2020). This test is robust to non-normal data, non-linear trends, and series with outliers and missing values. Trend magnitudes were reported in ^∘C decade⁻¹ for T_a and T_w and in % decade⁻¹ for Q (percentage of changes over the whole period) to help for comparisons across the basin.

3.4 Hydroclimatic drivers of T_w trends

Many factors affect the spatiotemporal variability of T_w. In the current study, we considered T_a as a proxy for heat fluxes and meteorological variables and Q as a proxy for thermal inertia and hydraulic geometries (which depends on Q; see Eqs. 6 and 7). To understand variability in T_w trends, we assessed spatial coherence and temporal synchronicity between trends in T_w and trends in T_a and Q as the main hydroclimate drivers influencing T_w. To do so, we first assessed the spatial coherence between these variables. In this regard, distributions of trends in T_w and T_a were first compared for the whole basin at the seasonal and annual scales using the non-parametric Wilcoxon signed rank test (Bauer, 1972) to determine whether median T_w trends were greater than median T_a trends. Then, the spatial coherence across reaches in terms of difference in trends between T_w and T_a on the one hand and sign of trend in Q on the other hand was assessed to explain the discrepancy between T_w and T_a found in the previous step with respect to Q.

Then, we assessed the temporal link between these three variables (i.e. T_w, T_a, and Q). To do so, for each identified reach, we first evaluated the strength and direction of joint trends using Pearson correlation between (i) T_w and T_a and between (ii) T_w and Q. T_a, T_w, and Q seasonal and annual anomalies – with respect to the 1963–2019 interannual mean – were then used to assess the synchronicity of extreme years and change points in time series. 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.

3.5 Landscape drivers of T_w trends

Stream size, within individual large-scale homogeneous HERs, was selected as the first major potential landscape driver. 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 decadal trends in T_w (i.e. 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. Finally, in order to better illustrate the relationship between trends in T_w and Strahler order, median T_w trends of each group of rivers with respect to Strahler order were presented.

Second, the influence of riparian vegetation shading on trends in T_w was assessed using the daily average of the riparian vegetation shading (SF) simulated by T-NET. Seasonal shading was computed as the average of the daily SF over each season. For this analysis, only low-order reaches – distance from the source <30 km – were considered, based on previous observations that riparian shading primarily influences T_w at this scale (Dan Moore et al., 2005; Loicq et al., 2018). Then, as for the previous analysis for the influence of stream size, the correlation between decadal T_w trends and five levels of riparian shading (<15 %; 15 %–25 %; 25 %–40 %; 40 %–60 %; >60 %) was computed across HERs and seasons. Finally, median T_w trends were compared for each level of riparian shading.

4 Results

4.1 Performance of models against observations

The EROS model performed well at the 352 calibration stations at the annual scale with a median relative bias close to 0 % (see Fig. S6, left). It slightly underestimated winter Q (median relative bias (across stations) $=-\mathrm{6.27}$ %) and spring Q (−3.47 %) and overestimated summer Q (+34.7 %) and autumn Q (+20.9 %). The EROS model also performed well at 44 RHN stations with long-term continuous daily data at the annual scale with a median relative bias of 0.37 % (see Fig. S6, right). It slightly underestimated winter Q (median relative bias (across stations) $=-\mathrm{7.26}$ %) and spring Q (−6.79 %) and overestimated summer Q (+37.7 %) and autumn Q (+24.7 %). The mean NSE criteria for Q, ln (Q), and $\sqrt{Q}$ were > 0.7 for at least 75 % of both calibration and RHN stations (see Fig. S7). A good performance of the EROS model in reconstructing daily Q was also seen at three different hydrometric stations located in the upstream, middle, and downstream parts of the Loire River basin (see Fig. S8).

No systematic bias was found for T_w modelled by T-NET at the stations located on small and medium rivers (see Fig. S9, top left). Median T_w bias ranged from −0.26 ^∘C (in autumn) to 0.8 ^∘C (in winter). Large rivers exhibited a small T_w underestimation (see Fig. S9, top right), with a median bias ranging from −0.29 ^∘C (in autumn) to +0.15 ^∘C (in winter), and the overall biases were still quite small across seasons (interquartile range (IQR) = 0.4–0.7 ^∘C depending on season). On the other hand, the median RMSE of the T-NET thermal model, for small and medium rivers, ranged between 0.52 ^∘C (in annual) and 0.91 ^∘C (in DJF and JJA) across seasons (see Fig. S9, bottom left). For large rivers, the median RMSE of the T-NET thermal model ranged between 0.38 ^∘C (annual) and 1.11 ^∘C (SON) across seasons (see Fig. S9, bottom right). Moreover, 53 %–83 % stations (50 %–100 %) on small and medium (large) rivers had a RMSE < 1 ^∘C across seasons.

Trends in observed and modelled Q were significantly correlated for all seasons (Fig. S10), with the highest correlation across stations found in spring and autumn (r=0.69 and 0.71, p< 0.05) and the lowest correlation found in summer, which was also non-significant (r=0.17, p=0.26). The trends of both modelled and observed Q were slightly decreasing (up to −11 % decade⁻¹) for the majority of stations across all seasons, but the trend was significant for a very few of them (and mostly at the annual scale), all located in the southern part of the basin in HER A (red points of Fig. S10). Moreover, there were only a few discrepancies between estimates of trend significance in modelled and observed Q across seasons (11 %–18 % of stations).

Modelled and observed T_w trends also correlated significantly (see Fig. S11) across seasons, with the highest correlation in summer (r=0.94, p<0.001) and the lowest correlation in autumn, which was also non-significant (r=0.29, p=0.32). There was also good agreement between estimates of trend significance in modelled and observed T_w across seasons. Indeed, the trends in observed and simulated T_w were either both significant or neither was significant. Contrasting with trends in Q, trends for T_w increased for most stations across seasons, but the short period of record led to mostly non-significant trends. However, stations with long-term data showed significant increasing trends for all seasons, with the exception of winter (Fig. S11).

A visual comparison of observed and modelled T_w time series and trends at stations with long-term data (> 20 years) indeed suggested a strong coherence and agreement between observations and simulations provided by T-NET for all seasons (Fig. 2). For the four stations along the main stem of the Loire River (Fig. 2), the greatest increase occurred in spring – +0.61 (+0.71) ^∘C decade⁻¹ in observations (simulations) – and summer – +0.62 (+0.58) ^∘C decade⁻¹ in observations (simulations). The smallest increase was found in winter – +0.22 (+0.28) ^∘C decade⁻¹ in observations (simulations).

Figure 2 Seasonal and annual time series of observed and simulated T_w at stations with at least 20 years of long-term continuous data between 1977 and 2019 (see Table 1 for more information). The red and blue solid lines show the trend in observed and simulated T_w. Numbers in red and blue in the top left corner of each graph show trend magnitudes (Sen's slope) in observed and simulated T_w. Numbers in black in the bottom right corner of each graph show the mean bias of the reconstruction.

4.2 Spatial reconstruction of long-term trends

T_w significantly increased in almost all modelled reaches for all seasons (Figs. 3, left, and S12, left). Depending on the season considered, 62 % to 80 % of reaches showed trends in the range of +0.2 to +0.4 ^∘C decade⁻¹ (i.e. +1.14 to +2.28 ^∘C over the whole 1963–2019 period; see Fig. S13). Summer T_w 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. S13). Such reaches were mainly located in the southern part of the basin, in HER A (see Fig. 3, left). Spring T_w trends showed a similar spatial pattern but with lower trend values.

Figure 3 Spatial variability of trends in seasonal and annual T_w, T_a, and Q over the 1963–2019 period based on Sen's slope estimator. Solid black lines show the Hydro-Ecoregion (HER) delineation (see Fig. 1). The statistical significance of these trends is given in Fig. S12.

Likewise, T_a exhibited increasing trends for 99 % of all reaches across spring, summer, and the whole year (Figs. 3, middle, and S12, middle). Values were mainly in the range of +0.2 to +0.4 ^∘C decade⁻¹ (see Fig. S13). The highest T_a trend values were found in summer and spring, when 67 % and 22 % of reaches, respectively, showed values higher than +0.4 ^∘C decade⁻¹. Such reaches were mainly located in HER A. Non-significant trends were found over the whole basin in winter and in the southern part of the basin in autumn.

In contrast to T_a and T_w, trends in Q were highly variable in magnitude and direction across the basin and across seasons (Fig. 3, right), and most were not significant at p=0.05 (Fig. S12, right). However, significant decreasing trends were found in the southern part of the basin (HER A) in spring, summer, and autumn and at an annual timescale (Fig. S12). Decreasing trends were observed for the majority of reaches across seasons (66 %–83 % of reaches), with the exception of winter (37 %) (see Fig. S14). Decreasing trends could have magnitudes greater than −5 % per decade, implying a −28 % loss in Q over the whole 1963–2019 period (see Fig. S14).

4.3 Hydroclimatic drivers of T_w trends

Spatial coherence

The medians of T_w trends were higher than those of T_a trends for every season (p<0.001 according to the Wilcoxon signed rank test), except for summer, when the median trend values for T_w and T_a were very similar, but more variable for T_w (+0.08 to +1.02 ^∘C decade⁻¹) (Fig. 4). The greatest increase in median T_w over the basin was found in summer (+0.44 ^∘C decade⁻¹).

Overall, T_w trends were more spatially variable than T_a trends (see Figs. 3 and 4), suggesting the conditional influence of other factors like Q trends. Indeed, for the majority of reaches, T_w trends exceeded T_a trends (see Fig. 5), with the exception of summer, when it was the case for half of the reaches. The difference between T_w and T_a trends could go up to 0.5 ^∘C decade⁻¹ (i.e. up to 2.85 ^∘C over the whole 1963–2019 period) (see Fig. 5). At the majority of such reaches (where T_w trend > T_a trend), decreasing Q trends occurred coincidentally (43 %–94 % of all reaches) for all seasons (with the exception of winter), irrespective of whether all significant and non-significant trends were considered (Fig. 6, top) or only significant trends of all three variables (i.e. T_w, T_a, and Q) were considered (Fig. 6, bottom). Of these specific reaches where T_w trend > T_a trend and Q trend < 0, most were located in HER A (see Figs. S15 and S16). Moreover, wherever T_w trend < T_a trend, increasing Q trends occurred coincidentally (see the dark blue bars in Fig. 6), suggesting again the influence of Q trends as a possible explanation for discrepancies found between T_w and T_a trends.

Figure 4 Distributions of seasonal and annual trends in T_w and T_a for all 52 278 reaches over the 1963–2019 period. Sen's slope is used as the trend value estimate (see text). This representation includes both significant and non-significant trends.

Figure 5 Difference between T_w and T_a trends at each reach in ^∘C per decade for all 52 278 reaches.

Figure 6 Percentage of reaches with consistent 1963–2019 trends in T_w, T_a, and Q, categorized with respect to two criteria: (1) comparison between the T_w and T_a trends and (2) sign of the Q trend. Sen's slope is used as the trend value estimate.

Temporal synchronicity

We observed strong positive correlations between seasonal and annual averages of T_w and T_a across seasons ($r=+\mathrm{0.72}$ to +0.82 depending on the season; see Table 2). We further observed a strong negative correlation between summer T_w and Q time series ($r=-\mathrm{0.40}$).

Table 2 Pearson correlation over the 1963–2019 period between seasonal and annual T_w and T_a and/or Q time series, averaged over all reaches. Percentages in brackets show the proportion of reaches with a significant correlation at the 95 % confidence level.

Annual anomalies of T_w, T_a, and Q exhibited variable patterns, with T_w and T_a generally increasing and Q remaining relatively stationary (Fig. S17). T_w anomalies were generally more variable than T_a, especially in summer, but both time series appeared to exhibit synchronous behaviour. Change-point analysis supported this visual observation, where change points in seasonal and annual averages were largely coincident across these time series (Fig. 7). T_w and T_a anomalies exhibited clear negative-to-positive change points in the late 1980s at nearly all reaches, with median values increasing by approximately +2 ^∘C after the change point. These change points were observed in all seasons but were most pronounced and synchronous around 1988 in spring and summer. The change points detected in winter T_w and T_a time series were less concomitant, occurring mostly in the early 1990s (1992 and 1993) for T_w and in the late 1980s (1986–1989) for T_a. The autumn change points were distributed between 1980 and 1994 for both T_w and T_a. The significant change points in seasonal Q time series were detected for a substantially smaller proportion of reaches, e.g. fewer 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. S18). In spring and summer, they occurred in the late 1980s, similarly to T_w and T_a. 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 T_w and T_a positive anomalies over the study period were observed in 1976, 1994, 1995, 2003, 2005, 2006, 2015, 2017, 2018, and 2019, which corresponded to years with the largest negative anomalies in summer Q, with the exception of 1994 and 2018 (Fig. 8). Note that this signal was much less clear for the other seasons (see Fig. S19).

Figure 7 Change points in T_w, T_a, and Q time series at the seasonal and annual scales, plotted as a proportion of reaches experiencing a shift in a given year. Only the first change point detected at the 95 % confidence level is considered, and non-significant change points were removed, leading to curves not reaching 100 %.

Figure 8 Relationship between summer anomalies in T_w and T_a on the one hand and summer T_w and Q anomalies on the other hand. Individual years are identified from their median values across all reaches. Years with the highest anomalies in T_w ($>+\mathrm{1}$ ^∘C) and corresponding anomalies in T_a and Q are identified in red.

4.4 Landscape drivers of T_w trends

4.4.1 Stream size and HER

Strahler order was strongly (p<0.001) and positively correlated with T_w trends for all HERs in spring and for HER A in summer and autumn and at the annual scale. HER A exhibited the highest positive correlations in spring ($r=+\mathrm{0.32}$) and summer ($r=+\mathrm{0.15}$) (Fig. S20). In other words, larger rivers tended to exhibit the largest increases in spring and summer T_w, especially for reaches located in HER A. There, median trends in spring (summer) ranged from +0.37 to +0.55 ^∘C decade⁻¹ (from +0.49 to +0.64 ^∘C decade⁻¹) from small streams (Strahler order 1) to large rivers (Strahler order ≥5).

4.4.2 Riparian shading and HER

For small streams, i.e. reaches located closer than 30 km from the source, the SF and trends in T_w were significantly (p<0.001) and negatively correlated in HER A in all seasons as well as in HER B and C in spring and in HER B in autumn (Fig. 9). Indeed, across seasons, the highest negative correlation was found in HER A ($r=-\mathrm{0.56}$ to −0.37 depending on the season). Unsurprisingly, the mitigating effect of shading on trends in T_w for small streams was observed for all HERs in spring and only for HER A in summer and, to a lesser extent, in autumn and winter. The median spring T_w trend was mitigated by 0.12 ^∘C decade⁻¹ from sparsely shaded reaches (SF <15 %) to highly shaded reaches (SF >40 %). For summer T_w in HER A, the median trends were mitigated by 0.16 ^∘C decade⁻¹ from the lowest shaded reaches to the highest shaded reaches.

Figure 9 Relationships between shading factor and median trends in T_w over the 1963–2019 period for small streams, by HER and by season. Note that some shading factor classes are not observed in autumn and winter. Correlations and associated p values are shown in the top right corner of each graph, and significant relationships at the 95 % confidence level are identified by full solid lines.

5 Discussion

5.1 Quality and suitability of simulated T_w and Q

Although some biases were observed for both Q (Fig. S6) and T_w (Fig. S9), we found significant correlations between modelled and observed trends in seasonal and annual Q, with the exception of summer (Fig. S10), and T_w, with the exception of autumn (Fig. S11). The low correlation value found in summer Q (Fig. S10) originated from poor simulation at very few stations, all located in HER B. Two of these stations gauged catchments where numerous small ponds were found, and the highly decreasing observed trends might be due to the increasing evaporation from these ponds, which were obviously not included in the EROS hydrological model. This was also true to a lesser extent for the other hydrometric station, in which a canal followed a large part of the course of the river and might play a similar role with respect to summer evaporation trends. Apart from these specific stations, in summer, a coherence as good as in other seasons was found between trends in simulated and observed Q. Moreover, the spatial pattern in simulated Q trends, with significant decreases in the southern headwaters, was consistent with observations at the high-quality reference hydrometric stations (Giuntoli et al., 2013).

The low correlation between simulated and observed T_w trends found in autumn (Fig. S11) originated from two stations with 8–13 years of T_w data, while such correlation was really good at stations on the Loire and Vienne rivers with longer (20–40 years) T_w data. Therefore, poor correlation in autumn could be due to insufficient T_w data at these two stations. Moreover, consistent with Moatar and Gailhard (2006) and Arevalo et al. (2020), we found no trend (p>0.05) in both observation and simulation on the Loire (Dampierre) in winter.

5.2 Comparison with T_w trends in other European basins

T-NET simulations over the 1963–2019 period showed significantly increasing trends in T_w for almost all reaches over the Loire basin across seasons (see Figs. 3, left, and S12, left), with an increase of +0.38 ^∘C decade⁻¹ on average at the annual scale. To the best of the authors' knowledge, the present study was one of the few studies using modelled T_w to investigate past trends at a large scale (but see Van Vliet et al., 2011; Isaak et al., 2012, 2017). Table 3, summarizing recently published findings based on observations, demonstrates that the present results are consistent with trends observed in other European basins, with clear increases in T_w over the recent decades. It also shows that the much larger scale and finer spatial resolution of the current study clearly stands out as unique. Although start year, end year, and length of the study period can have a significant influence on trend estimates and trend detection (Arora et al., 2016), comparing trends with other studies conducted over different periods gives us a comprehensive view of the overall magnitude of changes in T_w and possible related drivers.

Niedrist and Füreder (2021)Wilby and Johnson (2020)Michel et al. (2020)Kędra (2020)Maire et al. (2019)Ptak et al. (2019a)Zhu et al. (2019)Zobrist et al. (2018)Arora et al. (2016)Orr et al. (2015)Ptak et al. (2016)Moatar and Gailhard (2006)

Table 3 Recent studies on T_w trends in Europe. A comprehensive review of the relevant literature published prior to 2016 can be found in Arora et al. (2016). Magnitudes with unspecified seasons are related to the annual scale. Note that, except for the present study, the others used observed T_w.

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 to +0.4 ^∘C decade⁻¹, Fig. S13). 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 T_w and T_a trends across seasons for the majority of the reaches (see Figs. 3 and 5), which is a common finding for other sites around the world (Arora et al., 2016; Wanders et al., 2019). This highlights that changes in T_a may not be the only driver of changes in natural T_w.

5.3 Drivers and spatial patterns of trends in T_w

Consistent with our findings (see Fig. 5), Moatar and Gailhard (2006) found T_w increased more quickly than T_a in spring and summer and at the annual scale for all four stations on the Loire River. Arora et al. (2016) also found T_w trends > T_a trends in summer. In Switzerland, Michel et al. (2020) described an increase of $+\mathrm{0.33}\pm \mathrm{0.03}$ ^∘C decade⁻¹ in T_w, resulting from the joint effects of an increase in T_a ($+\mathrm{0.39}\pm \mathrm{0.14}$ ^∘C decade⁻¹) and a decrease in Q ($-\mathrm{10.1}\pm \mathrm{4.6}$ % decade⁻¹) over the 1979–2018 period. In contrast with our results, they found T_w trends lower than T_a trends due to the influence of snowmelt and glacier melt in Alpine catchments. Consistent with their findings, Orr et al. (2015) also found T_w trends < T_a trends in the UK. They suggested that such differences between T_w and T_a trends could be a result of different processes driving T_w. In the current study, we found spatial coherence between trends in T_w and trends in T_a and Q. Indeed, the greatest increases in T_w (up to +1 ^∘C decade⁻¹) were predominately located in the southern part of the basin, in HER A (Massif Central), where a greater increase in T_a (up to +0.71 ^∘C decade⁻¹) and a greater decrease in Q (up to −16  % decade⁻¹) occurred jointly. We also found, at the majority of reaches where T_w trend > T_a trend, decreasing Q trends occurred coincidentally for all seasons (with the exception of winter) (Fig. 6).

The decrease in Q could itself be due to a significant increase in PET (up to +10 %) over the whole basin and a decrease in total precipitation (P) (up to −5 % decade⁻¹) (Figs. S21 and S22). 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. We also observed, for the majority of reaches where T_w increased less than T_a, an increase in Q that occurred jointly, suggesting that an increase in Q or in other words an increase in thermal inertia could also explain the discrepancy between T_w and T_a trends at these reaches (see Fig. 6).

A strong synchronicity between T_a and T_w 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 T_w could be due to co-occurrence with the increase in summer T_a (average correlation: +0.82) and with the decrease in summer Q (average correlation: −0.40). These findings are consistent with those of Michel et al. (2020): average T_w–T_a correlation +0.61 and average T_w–Q correlation −0.66. For the middle Loire River, Moatar and Gailhard (2006) found that the increase in T_a (decrease in Q) explains 60 % (40 %) of the increase in T_w. Moreover, the significant change point in T_w, T_a, 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 T_w time series like the Loire at Dampierre also displayed a similar change point.

Trends in T_w might also be explained by trends in additional drivers (not considered directly in the current study), like shortwave radiation (Wanders et al., 2019), which is the dominant flux at the air–water surface and is notably increasing over Europe (Sanchez-Lorenzo et al., 2015). This might explain why T_w increased more quickly than T_a in spring, summer, and autumn, when no decreasing trends in Q were found (see Fig. 6).

The current study suggests that T_a and Q could exert a joint influence on T_w, based on an analysis of the spatial coherence and temporal synchronicity of these variables. Assessing the causal influence of these factors on T_w trends goes beyond the scope of this paper and is left for future research. In this regard, one could devise a formal attribution framework where one may e.g. remove trends in Q and trends in T_a alternatively in T-NET inputs.

5.4 Natural trends and anthropogenic influence on T_w

Natural T_w time series were used in the current study for detecting trends, as both the 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 T_w regimes in a variety 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 lower downstream T_w (Olden and Naiman, 2010) and mitigate the increasing trend in T_w (Cheng et al., 2020). Nevertheless, it is anticipated that a considerable proportion of streams regulated by large dams 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 T_w trends and gathers most of the existing large dams.

On the other hand, ponds and shallow reservoirs, by releasing warm water, can increase downstream T_w (Chandesris et al., 2019; Seyedhashemi et al., 2020; Zaidel et al., 2020) and exacerbate increasing trends in T_w (Wanders et al., 2019; Michel et al., 2020). The warming effect of such surface waters in the current study can be significant for streams located in lowlands in the middle and north of the Loire River basin, where most of the shallow reservoirs and ponds are located. In these streams, anthropogenically induced trends in T_w may be greater than natural ones, and the warming process can get worse through the increasing demand for storing water in small reservoirs for irrigation. Nevertheless, the warming effect can also be local, and streams located closer to such regulated streams may show limited to no warming (Wanders et al., 2019).

Note that, although there are nuclear power plants in the Loire basin, their impacts on T_w are considered negligible according to Moatar and Gailhard (2006) and Bustillo et al. (2014). Moreover, it was observed in the current study that the T_w trend at Belleville located upstream of power plants and, consequently, not influenced by them had the same trend magnitude as the other three stations located downstream (see Fig. 2), showing a negligible influence of nuclear power plants on T_w trends.

5.5 Implications for river management and aquatic biota

The removal of riparian vegetation can increase T_w (Caissie, 2006), and changes in T_w can be even more sensitive to changes in riparian vegetation than to changes in T_a or Q (Wondzell et al., 2019). We showed that, in small streams, an increase of > 25 % of riparian shading (from 15 % to 40 %) could mitigate the median trend in spring and summer T_w by up to 0.16 ^∘C decade⁻¹ (Fig. 9). Spring and summer T_w trends were more pronounced in large rivers, especially in the south of the basin, with a difference in median T_w trends of up to +0.18 ^∘C decade⁻¹ (Fig. S20), probably due to a decrease in Q (up to −2 % decade⁻¹; see Fig. S23), a greater thermal sensitivity, and the absence of mitigating factors like riparian vegetation shading or groundwater inputs (Kelleher et al., 2012; Beaufort et al., 2020).

Restoring riparian vegetation and shading can therefore substantially mitigate future increases in T_w. 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 trees, 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 (Edmonds et al., 2000; Caissie, 2006). 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 availability and freshwater quality for organisms (Lennox et al., 2019), and shifts in species distribution (Comte et al., 2013) are already observed as consequences of the long-term increase in T_w. Some major changes in fish density and community structure have already been reported in large rivers over France (Maire et al., 2019), for which we also found greater trends in T_w compared with 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 T_w (Morales-Marín et al., 2019).

6 Conclusions

Regional trends in T_w 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 three variables (T_a, Q, and T_w) and five timescales (seasons plus annual), we found consistent increasing T_w trends at the scale of the entire Loire River basin, regardless of the season. Critically, the rate of warming for stream temperature was in the majority of reaches higher than the rate of atmospheric warming, suggesting a decoupling of thermal trajectories linked to decreasing Q, especially in the southern part of the basin, supported by observed coherent spatial and temporal patterns and well-understood physical processes. Moreover, in the southern part of the basin, T_w 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 small rivers.

The synchronicity of extreme events of low flows and high stream temperature in the southern part of the basin will likely generate a double penalty for existing cold-water aquatic communities. However, riparian shading in small streams in the southern part of the basin may mitigate such warming. These findings underscore that air temperature alone is likely not an adequate proxy for explaining 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 T_w sensor networks. It also highlights that the physical process-based thermal models like T-NET can be used to assess the various stresses on freshwater habitat sustainability due to changes in Q and T_w. This knowledge can help develop appropriate management strategies to conserve thermal refugia and mitigate extreme thermal events induced by climate change.