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 difﬁculty 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 streamﬂow at the scale of the entire Loire River basin in France (10 5 km 2 with 52 278 reaches). Stream temperature increased for almost all reaches in all seasons (mean = + 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


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 Published by Copernicus Publications on behalf of the European Geosciences Union.(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., 2013Kurylyk et al., , 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., 2017Jackson et al., , 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.

Study area
The Loire River basin is one of the largest in Europe (10 5 km 2 ), encompassing an area with starkly contrasting HydroEco Regions (HERs), land use/land cover, and climatic conditions (Moatar and Dupont, 2016), providing an ideal Hydrol.Earth Syst.Sci., 26, 2583Sci., 26, -2603Sci., 26, , 2022 https://doi.org/10.5194/hess-26-2583-2022 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 %).

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.

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 andMoutzopoulos (1995), andThiéry (2018) for EROS and Beaufort et al. (2016) and Loicq et al. (2018) for T-NET.

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 2 ; mean drainage area = 300 km 2 ), 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 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 Météo-France over the 1958-2019 period (Quintana-Segui et al., 2008;Vidal et al., 2010).ET 0 was computed from Safran variables with the Penman-Monteith equation (Allen et al., 1998).

T-NET thermal model principles and input data
The T-NET thermal model is a fully mechanistic, 1D model that simulates hourly T wi,j at distance, i, along reach, j , by solving the local heat budget.The heat budget of each reach includes six fluxes (W m −2 ): 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 (dT w / dx) for steady-state conditions in order to achieve thermal equilibrium (i.e.H i,j = 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.firstorder 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.https://doi.org/10.5194/hess-26-2583-2022 Hydrol.Earth Syst.Sci., 26, 2583-2603, 2022  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).
-Meteorological variables: hourly T a ( • C), specific humidity (g kg −1 ), wind velocity (m s −1 ), shortwave radiation (W m −2 ), and longwave radiation (W m −2 ) 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: 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: 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 final = Hydrol.Earth Syst.Sci., 26, 2583Sci., 26, -2603Sci., 26, , 2022 https://doi.org/10.5194/hess-26-2583-2022SF × (1 − 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): 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.

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 = mNSE − w × 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)(1972)(1973)(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)(1959)(1960)(1961)(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 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 2 (range = 7-1342 km 2 ) for small/medium streams and 18 926 km 2 (range = 1931-57 043 km 2 ) 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).

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 3 s −1 ).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 −1 for T a and T w and in % decade −1 for Q (percentage of changes over the whole period) to help for comparisons across the basin.

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 nonparametric 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 nonparametric 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 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 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.

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) = −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) = −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 √ 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 −1 ) 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 https://doi.org/10.5194/hess-26-2583-2022 Hydrol.Earth Syst.Sci., 26,2022 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 −1 in observations (simulations) -and summer -+0.62 (+0.58) • C decade −1 in observations (simulations).The smallest increase was found in winter -+0.22 (+0.28) • C decade −1 in observations (simulations).

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 −1 (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 −1 (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.
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 −1 (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 −1 .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).

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 −1 ) (Fig. 4).The greatest increase in median T w over the basin was found in summer (+0.44 • C decade −1 ).
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 −1 (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.

Temporal synchronicity
We observed strong positive correlations between seasonal and annual averages of T w and T a across seasons (r = +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 = −0.40).
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 Hydrol.Earth Syst.Sci., 26, 2583Sci., 26, -2603Sci., 26, , 2022 https://doi.org/10.5194/hess-26-2583-2022(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).

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 = +0.32)and summer (r = +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 −1 (from +0.49 to +0.64 • C decade −1 ) from small streams (Strahler order 1) to large rivers (Strahler order ≥ 5).

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 = −0.56 to −0.37 depending on the season).Unsurprisingly, the miti-gating 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 −1 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 −1 from the lowest shaded reaches to the highest shaded reaches.

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., 26, 2583-2603, 2022et al., 2011;Isaak et al., 2012Isaak et al., , 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.

Comparison with
Global-scale stream temperature modelling suggests trends in annual averages ranging from +0.2 to +0.5 • C decade −1 over France (Wanders et al., 2019), which is consistent with our findings (mostly in the range of +0.2 to +0.4 • C decade −1 , 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 .

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 +0.33 ± 0.03 • C decade −1 in T w , resulting from the joint effects of an increase in T a (+0.39 ± 0.14 • C decade −1 ) and a decrease in Q (−10.1 ± 4.6 % decade −1 ) 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 −1 ) 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 −1 ) and a greater decrease in Q (up to −16 % decade −1 ) 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 −1 ) (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.

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 vari-  ety 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 https://doi.org/10.5194/hess-26-2583-2022 Hydrol.Earth Syst.Sci., 26, 2583-2603, 2022 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.

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 −1 (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 −1 (Fig. S20), probably due to a decrease in Q (up to −2 % decade −1 ; 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).

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 processbased 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.
Data availability.The seasonal and annual T w and Q time series at 67 stations (see Tables 1 and S5) can be downloaded from https://doi.org/10.15454/PTY9R7(Seyedhashemi et al., 2022).Data for other reaches in the basin and codes are available from the corresponding author, Hanieh Seyedhashemi, upon reasonable request.

Figure 1 .
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 TableS5), 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 Table1.The numbers in red in the middle panel correspond to the ID of long-term T w stations in Table1.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 TableS2), 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 TableS4).

Figure 2 .
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 Table1for 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.

Figure 3 .
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.

Figure 4 .
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 .
Figure 5. Difference between T w and T a trends at each reach in • C per decade for all 52 278 reaches.
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 −1 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 https://doi.org/10.5194/hess-26-2583-2022Hydrol.Earth Syst.Sci., 26, 2583-2603, 2022

Figure 6 .
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.

Figure 7 .
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 .
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 (> +1 • C) and corresponding anomalies in T a and Q are identified in red.

Figure 9 .
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.
Validation of the T-NET thermal modelT-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 https://doi.org/10.5194/hess-26-2583-2022Hydrol.Earth Syst.Sci., 26, 2583-2603, 2022 3.2.2

Table 1 .
Characteristics of the 14 long-term T w stations.See Fig.S4for 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.

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.

Table 3 .
Arora et al. (2016) w trends in Europe.A comprehensive review of the relevant literature published prior to 2016 can be found inArora 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 .