Climate and basin drivers of seasonal river water temperature dynamics

. Stream water temperature is a key control of many river processes (e.g. ecology, biogeochemistry, hydraulics) and services (e.g. power plant cooling, recreational use). Consequently, the effect of climate change and variability on stream temperature is a major scientiﬁc and practical concern. This paper aims (1) to improve the understanding of large-scale spatial and temporal variability in climate– water temperature associations, and (2) to assess explicitly the inﬂuence of basin properties as modiﬁers of these relationships. A dataset was assembled including six distinct modelled climatic variables (air temperature, downward short-wave and long-wave radiation, wind speed, speciﬁc humidity, and precipitation) and observed stream temperatures for the period 1984–2007 at 35 sites located on 21 rivers within 16 basins (Great Britain geographical extent); the study focuses on broad spatio-temporal patterns, and hence was based on 3-month-averaged data (i.e. seasonal). A wide range of basin properties was derived. Five models were ﬁtted (all seasons, winter, spring, summer, and autumn). Both site and national spatial scales were investigated at once by using multi-level modelling with linear multiple regressions. Model selection used multi-model inference, which provides more robust models, based on sets of good models, rather than a single best model. Broad climate–water temperature associations common to all sites were obtained from the analysis of the ﬁxed coefﬁcients, while site-speciﬁc responses, i.e. random coefﬁcients, were assessed against basin properties with analysis of variance (ANOVA). All six climate predictors investigated play a role as a control of water temperature. Air temperature and short-wave radiation are important for all models/seasons, while the other predictors are important for some models/seasons only. The form and strength of the climate–stream temperature association vary depending on season and on water temperature. The dominating climate drivers and physical processes may change across seasons and across the stream temperature range. The role of basin permeability, size, and elevation as modiﬁers of the climate– water temperature associations was conﬁrmed; permeability has the primary inﬂuence, followed by size and elevation. Smaller, upland, and/or impermeable basins are the most inﬂuenced by atmospheric heat exchanges, while larger, lowland and permeable basins are the least inﬂuenced. The study showed the importance of accounting properly for the spatial and temporal variability of climate–stream temperature associations and their modiﬁcation by basin properties.

Consequently, the effect of climate change and variability on stream temperature is a major scientific and practical concern (Garner et al., 2014).River thermal sensitivity to climate change and variability is controlled by complex drivers that need to be unravelled to better understand (a) patterns of spatio-temporal variability and (b) the relative importance of different controls to inform water and land management, especially climate change mitigation and adaptations strategies (Hannah and Garner, 2015).There is a growing body of river temperature research but there is still limited understanding of large-scale spatial and temporal variability in climate-WT associations, and of the influence of basin properties as modifiers of these relationships (Garner et al., 2014).Due to the focus on large scales, this paper is not investigating higher frequency temporal patterns (e.g.heat waves) or smaller spatial variability (e.g.thermal diversity and refugia).This paper extends on Laizé (2015).
River thermal regimes are complex because they involve many interacting drivers (Hannah et al., 2004(Hannah et al., , 2008)).Caissie (2006) identified atmospheric conditions as the primary group of controls, with hydrology linked to basin physical properties (e.g.topography, geology) as secondary influencing factors.
The main climate variables (Fig. 1), which constitute an "atmospheric conditions" group, can be identified by analysing the theoretical heat budget for a stream reach without tributary inflow, which may be expressed as (adapted from Hannah and Garner, 2015) where Q n is the total net heat exchange, Q * the heat flux due to net radiation, Q h the heat flux due to sensible transfer between air and water (sensible heat), Q e the heat flux due to evaporation and condensation (latent heat), Q bhf the heat flux to and from the river bed, Q f the heat flux due to friction at the bed and banks, and Q a the heat flux due to advective transfer by precipitation and groundwater.The different components of Eq. ( 1) correspond to different processes, related to climatic and hydrological conditions.Q * corresponds to short-wave radiation (insolation from the sun) and long-wave radiation (emitted towards the stream by clouds and overhanging surfaces such as vegetation, and re-emitted back to space (lost) at water surface temperature).Q h corresponds to convective energy exchanges between air and water (at the surface) causing heat loss or gain.Q e represents heat loss by evaporation or gain by condensation.Q bhf and Q f do not relate directly to climate processes but rather local hydrological conditions.Q f can be assumed to be negligible in many systems (e.g.Hannah et al., 2008).Q a corresponds to advective heat exchanges, e.g.inflow or outflow into the river reach, hyporheic exchange, groundwater.A direct, climatic component of Q a is precipitation inputs, which is thought to have a limited contribution (Caissie, 2006).Caissie, 2006, andHannah et al., 2008).
These variables are not independent.Figure 1 features a schematic representation of the interactions between these variables.Downward short-and long-wave radiation increase not only WT but also air temperature, and then there are exchanges between air and water, to influence sensible heating.Additionally, wind plays a significant role by increasing evaporative cooling and in modifying the air-water exchanges by increasing mixing (Hannah et al., 2008).The physical equations underpinning the role of wind can be found in Caissie et al. (2007).
A review of recent international water temperature research can be found in Hannah and Garner (2015).To date, most UK-focused studies (Table 1) tend to be either specific to a few monitoring sites, to have a limited geographical extent (i.e.focused with specific region of the country), and/or to consider few climate drivers.In addition, seasonality, which has huge ecological relevance with regard to phenology, is only explored formally in a small number of papers (e.g.Langan et al., 2001;Hrachowitz et al., 2010).A major difficulty is to pair WT and climate monitoring sites, as monitoring is rarely coordinated, then to identify time series with long enough common periods of record.For example, Garner et al. (2014) undertook an England-Wales-scale study and matched water temperature monitoring sites with climate and hydrological monitoring sites for 38 temperature sites out of ∼ 3000 sites in the Environment Agency's Freshwater Temperature Archive (Orr et al., 2014).Garner et al. (2014) is one of the few studies internationally (e.g.Hrachowitz et al., 2010, in the UK;Isaak and Hubert, 2001, Nelitz et al., 2007, or Isaak et al., 2010, in North America)   and empirically the role of a limited number of basin properties with regard to stream temperature.
In most of these studies, analyses are done on a site-bysite basis, which limits the extent to which broad patterns can be inferred (statistical results for a given site are only valid for that site); Caissie (2006) emphasized this as a limitation when having to work across different spatial scales.In contrast, studies like Garner et al. (2014) group sites together using classification techniques to identify regional patterns.However, doing so causes a loss of information since data points of all sites within a class are summarized and intra-class differences lost, and inferences at group level are not necessarily valid at site level.An alternative analytical/ statistical method, which can characterize broad patterns while preserving individual site information, should be investigated.
The following research gaps are identified (above): (a) climate-WT studies in the UK used a limited number of WT sites or climate explanatory variables (focus on air temperature links to WT) and/or are limited in geographical extent; (b) limited formal analysis of seasonality; (c) limited knowledge of the role of basin properties as modifiers of climate-WT associations; and (d) need for an alternative analysis method to optimize data utility.
Given this context, the aims of this study are (1) to improve the understanding of large-scale spatial and temporal variability in climate-WT associations, and (2) to assess explicitly the influence of basin properties as modifiers of these relationships.This paper resolves the issue of driving data availability by using a comprehensive and consistent set of modelled climate data (see Table 2 below).With a period of record from 1984 to 2007 (24 years), for a total of 35 sites located on 21 rivers within 16 basins (providing a Great Britain-wide geographical extent), six distinct modelled cli-matic variables were taken within 1 km of the sites.The study focuses on broad spatio-temporal patterns; hence, it is based on 3-month-averaged data (i.e.seasonal).Such a temporalscale limits issues of temporal auto-correlation often found in water temperature time series (Caissie, 2006).The study also investigates a wider range of basin properties than previous studies.
Innovatively, this paper investigates both site and national spatial scales at once.Multi-level (ML) modelling with linear multiple regressions is applied as an alternative to sitespecific or to classification-based analyses because it allows for pooling of all site data together while taking into account data structure (i.e.observations at site, sites within same basin) as well as not losing any information due to classlevel data averaging (Zuur et al., 2009).With this modelling technique, it is possible to investigate both study aims (i.e. the broad climate-WT associations common to all sites, and the site-specific responses that may be related to basin properties) within the same analysis framework.

Data
With regard to research aim 1 of this paper, observed river temperature data were assembled with a view to maximize spatial and temporal coverage as much as practically possible.To address the issue of mismatching monitoring networks, climate variables were obtained from a modelled dataset.The paired climate-WT dataset used in this paper has been published online via an open-access data repository (Laizé and Bruna Meredith, 2015).With regard to aim 2, a comprehensive and consistent set of basin properties were derived for all study sites.et al., 2006).Whether recording was done manually or with a logger, measures were instantaneous.Because these original projects were focused on natural rivers, the temperature data used herein may be considered as largely free from artificial influences (e.g.no industrial use for cooling or heated effluent discharges).

Climate data
The Climate Hydrology and Ecology research Support System (CHESS) dataset features six climate variables (Table 2).CHESS is the forcing dataset for the Joint UK Land Environment Simulator model (JULES; Best et al., 2011).CHESS is a UK-wide 1 km grid dataset derived by downscaling the UK Meteorological Office Rainfall and Evaporation Calculation System (MORECS) 40 km grids (Hough and Jones, 1997), except for precipitation that were derived from observed rain gauge data by using the natural neighbour interpolation method, which is a development of the Thiessen approach (Keller et al., 2006).For each 1 km cell, modelled daily time series of all variables are available for the pe- riod 1971-2007.The processes linked to AT, LWR, P, and SWR are given in the stream heat budget overview (see Introduction) and summarized in Table 2. Specific humidity (SH) gives a measure of evaporation potential (i.e. the more humidity, the less evaporation due to reduced vapour pressure gradients; e.g.Hannah et al., 2008).Wind speed (WS) captures the various effects of wind in increasing evaporation (cooling) and convective air-water exchanges (cooling or warming) Each CHESS cell was matched to the study temperature site(s) it contained.

Seasonal time series
First, sub-daily water temperature data were averaged at a daily time step (Frome, Great Ouse, Tadnoll, LOCAR) while spot measurements (Plynlimon, UKAWMN) were assumed representative of the day on which they were taken, although it is worth keeping in mind that they are only representative of daylight conditions.Second, daily water temperature data were matched to the daily climate data.Third, seasonal averages were computed from these daily data for all variables.Seasons were defined as December-February (winter), March-May (spring), June-August (summer), and September-November (autumn).For winter, these seasonal data for year y were based on data from December of year y − 1 to February of year y (e.g. for 1976, December 1975, January and February 1976).Lastly, five time series were derived from these data: one series per season at an annual time step (i.e. winter 2000, winter 2001, winter 2002, etc.), and one series with all seasons at a seasonal time step (i.e. autumn 2000, winter 2000, spring 2000, etc).These series and their related models are referred to hereafter as "autumn", "winter", "spring", "summer", and "all seasons".

Basin properties
Basin properties were derived from the UK Flood Estimation Handbook (FEH), the UK "industry standard" for flood regionalization studies, which includes 19 basin descriptors (Bayliss, 1999).A subset of descriptors was used.First, the 19 catchment descriptors were derived for each site.Many basin properties co-vary, often substantially, and they are best interpreted as groups of properties ("meta-properties") rather than on their own.Descriptor specifications (Bayliss, 1999), pair plots, and correlation matrices were checked to identify likely groups of descriptors (for example, all FEH rainfall descriptors capturing basin wetness).Three groups were identified, which relate to basin elevation, permeability (i.e.responsive impermeable vs. groundwater-fed basins), and size.These have been found to modify climate-hydrology associations in UK basins (e.g.Bower et al., 2004;Laizé and Hannah, 2010;Garner et al., 2014).Then, a test run of the basin property analysis outlined in Sect.3.3 (analysis of variance, ANOVA) was performed in order to check that all FEH descriptors from a given group of properties had consistent associations (positive or negative) with each model predictor (considering basin properties significantly associated with site-specific coefficients only), while one FEH descriptor was retained to represent each meta-property.
The following meta-properties and their corresponding FEH descriptors were thus selected for the final analysis: -Elevation/wetness ("elevation" hereafter): as noted in Laizé and Hannah (2010), basin elevation and wet-ness are very strongly correlated in the UK; the metaproperty "elevation" is represented by the "mean basin elevation above sea level" (m; FEH descriptor named "ALTBAR"), and, for the winter model only, by the proportion of time basin soils that are wet (%; FEH descriptor named "PROPWET"), based on soil moisture time series classified as wet/dry days (highly correlated to rainfall); elevation is also related to air temperature.
-Size: basin area (km 2 ; "AREA") using its natural log; area is a proxy for discharge, thus for thermal capacity, and is also linked to elevation.
-Permeability: base flow index from hydrology of soil type (BFIHOST; dimensionless); ranging from 0 (less permeable basin) to 1 (more permeable); temperature regimes in groundwater-fed (permeable) basins are expected to be more influenced by groundwater inputs than in impermeable basins.

Methods
This section describes the analytical methods used.First, as stated in the introduction, linear multiple regressions fitted with the multi-level (ML) modelling technique was chosen as the core method because it allowed one to analyse the multiple-site data in terms of both overall climate-WT associations (linked to research aim 1) and site-specific responses (linked to research aim 2; role of basins as modifiers of those associations).Although linear regressions are only approximating climate-WT associations (e.g.AT-WT associations are better described with logistic models; Mohseni et al., 1998), they were considered a sensible compromise.Second, with regards to overall climate-WT associations, ML model selection was done with multi-model inference (MMI), a state-of-the-art technique that selects sets of good models rather than a single best model (Grueber et al., 2011), to yield more robust models than with standard single model selection, especially given the number of climate predictors used.Lastly, any relation between site-specific climate-WT responses and basin properties were tested formally using an ANOVA.
The study work flow is summarized in Fig. 3: (a) WT observed data linked with (b) modelled climate variables, then (c) all are converted to seasonal (3-month) average series used within (d) a ML modelling/MMI framework producing (e) five output models (individual seasons and all seasons; aim 1), and (f) sets of basin properties (aim 2).

Multi-level modelling
To take into account the hierarchical nature of the water temperature dataset (e.g.data measured at the same site, sites located on the same river), ML modelling was used to build linear models with water temperature as the predicted variable, and the six climate variables as explanatory variables.When analysing multiple-site datasets, there are two common alternatives: (a) performing one regression for individual sites, or (b) one regression on all sites pooled together.On the one hand, site-specific regressions can make results highly uncertain (sites may have few data points; fitting numerous regressions is more prone to identify spurious relationships, i.e. type II errors).Thus, drawing out general patterns (e.g.variation between sites, effect of site characteristics) can be difficult.On the other hand, full pooling of sites ignores the clustering of samples within groups (e.g.measurements from a given site, or sites on the same river, may be more similar), which may hide important differences between groups and may cause problems with statistical inference (e.g.violation of the assumption of independence between samples, sites with large or small numbers of samples equally influencing the model outcome).
To overcome these issues, ML modelling can take into account the hierarchical structure in a dataset, i.e. the different "levels" at which data can be grouped (e.g.data at sites, sites within basins, basins within countries), thus allowing for the pooling of data from multiple sites.A ML model has two components, which correspond to generic patterns (i.e.similar to a regression on fully pooled data) and to level-specific patterns.The generic patterns, which are described by the explanatory variables as in a standard regression, are called the "fixed component" or "fixed effects" of the model.The unexplained variation between levels (e.g.patterns specific to a site) is termed the "random component" or "random effects".The random component captures the fact that levels may respond differently to a given predictor.For example, stream temperature could be very responsive to climate at one site (high slope value) but unresponsive at another (low slope value).In some cases, levels may have the same response to predictors but may have differing averages, i.e. differing with regards to their intercepts (e.g. two sites with same temporal patterns but with one site systematically cooler than another due to local characteristics or recoding procedure); such ML models are commonly known as "random intercept only".
In our analyses, a three-level data structure was applied: individual observations (level 1) nested within monitoring sites (level 2) nested within river stretches (level 3).In addition, a time variable was included as a predictor to take into account any linear trend in the time series.To avoid instability issues when fitting models, the predictors were centred (i.e.predictor values minus their mean).

Model selection with multi-model inference
Following standard ML modelling practice (e.g.Zuur et al., 2009), the model selection was applied in two stages: (a) selection of the random component variables and (b) selection of the fixed component variables.
First, the random component selection was done as follows.With all predictors included in the fixed component, all combinations of predictors in the random component were fitted.The models were then ranked using Akaike's information criterion (AIC; Akaike, 1974).AIC is used to select models offering the best compromise between fit and predictor parsimony; a model with a lower AIC achieves a better ratio of fit vs. number of predictors.Note that a variation of AIC was used: AICc, which is AIC corrected for small-size datasets.Selection was done for the four seasonal series as well as the "all season" series.In each case, the single combination of predictors giving the lowest AICc was retained as the random component.
Second, with the random component selected, the fixed component model selection followed the MMI approach, which selects sets of "good" models rather a single "best" one.Using a traditional model selection technique, like stepwise regression, the single model with the best (i.e.lowest) AICc would be selected.This presents two issues: (a) due In practice, all possible combinations of predictors using from one to six of the climate variables described above were fitted.The resulting models were ranked based their AICc.All models within four points of the lowest AIC were selected (Zuur et al., 2009).Each set of models was then summarized as an "average model" (predictor coefficients over all models in the set are averaged).Akaike weights (Burnham and Anderson, 2002) were then calculated; these are the re-scaled AICc scores of the models included in a MMI selection set.
The weights, which add up to 1, give an indication of how important relatively to each others the models are within a MMI set.For example, results showed that the "all seasons" model is based on two models with Akaike weights 0.74 and 0.26; the former model has more influence on the resulting average model than the latter.The Akaike weights form the basis to calculate the relative importance (RI) of each predictor: RI is how one reports on the role of each explanatory variable in MMI.For a given predictor, RI is calculated as the sum of the Akaike weights (re-scaled AICc) of the models in which that predictor is included.RI ranges from 0 (variable never included) to 1 (included in all models).For example, results showed that the "all seasons" model is based on two models with Akaike weights 0.74 and 0.26; the explanatory variable P is only included in the latter model, and hence its RI is 0.26, while the other five predictors are in both models and have a RI of 1 (see Table 3 below).With MMI, RI is analogous conceptually to predictor significance, assessed with p values, in a standard regression model.This is why p values are not calculated nor given in the Results section, but instead RI values for predictors are featured (a predictor with a higher RI is more significant).Grueber et al. (2011) covered the above points in details and gave a very good example of such an application of MMI in a natural sciences context.

Analysis of basin property influence
For those explanatory variables that were included in the random effects (i.e.different sites can have different coefficients), any relation between site-specific coefficients and basin properties was investigated by using maps and scatter plots of coefficients against basin properties, and by applying ANOVA to confirm observed patterns.For each coefficient and basin property, ANOVA is formally comparing (a) a model assuming there is no difference in coefficient between sites against (b) a model assuming the coefficient is a function of the basin property.A basin property is considered having significant influence on the WT-climate variable relationship when the ANOVA p value is < 0.05.To quantify the influence of these properties, either alone or combined, linear regressions of the site-specific coefficients against these properties were fitted.

Model selection and performance
As described above, selecting the five ML models was done in two stages.First, with all predictors included in the fixed component of the ML model, combinations of predictors as random effects were tested, and the combination yielding the lowest AICc was retained.As a result, the following variables were included as random effects (i.e.variables for which different sites have different coefficients): "all seasons" is AT and SWR; winter is SH; summer is P; autumn is SWR; no predictor was included for spring (random intercept only).Second, all combination of the predictors in the fixed components were tested with MMI.The number of models included in each final set as selected by MMI were: all seasons = 2; winter = 4; spring = 12; summer = 6; autumn = 14.
With ML models, standard R 2 are not appropriate; conditional R 2 (Nakagawa and Schielzeth, 2013), which are analogue to standard R 2 but designed for ML models, were calculated.Conditional R 2 were 0.96 for both all seasons models, 0.88 for all four winter models, within 0.88-0.89(mean 0.88) for the 12 spring models (mean 0.88), within 0.84-0.85(mean 0.84) for the six summer models, and within 0.88-0.89(mean 0.88) for the 14 autumn models.
With MMI, each set of models is summarized as an "average model", for which a given variable coefficient is its average value over all models in the set.The average model coefficients are presented in Table 3.All average models have good fits consistent with conditional R 2 given above, and as evidenced by plots of modelled against observed water temperature data in Fig. 4. Thereafter, if unqualified, the term "model" means the average model for a given set of selected models

Form and strength of associations between climate predictors and water temperature
The section focuses on the fixed effect coefficients of the predictors (i.e.coefficients valid for all sites).Predictors AT, SWR and SH have positive coefficients for all models (i.e.increases of these predictors are associated with a consistent warming effect on water temperature).Predictors LWR, WS, and P have positive or (mostly) negative coefficients (i.e.increases of these predictors are associated with warning or cooling, depending on season; Table 3).
The strength of the association varies with season.Comparing the absolute value of the seasonal coefficients for each variable (not between variables as they have different scales): AT, lowest in winter, highest in autumn; SWR, lowest in autumn, highest in winter; LWR, lowest in winter, highest in summer; WS, lowest in autumn, highest in summer; SH, lowest in autumn, highest in winter; P , lowest in summer, highest in autumn.

Relative predictor contributions
By definition, the predictors may have different units and orders of magnitude.Their coefficients cannot be compared directly to get an indication of their relative contribution to WT predictions.Instead, for each generic average model (see coefficients in Table 3), predicted WT values were generated for the whole period of record, then the percentage contributions of each predictor to these predicted WT values were calculated (i.e. a time series of predicted WT and of percentage contributions for the six predictors).Boxplots of the percentage contributions for the six predictors and the five models are featured on the left-hand side of Fig. 5 (for readability, outliers are not displayed).The thick black central line corresponds to the median percentage contribution.The shorter the boxes and whisker extents are, the more constant are predictor contributions to modelled WT, with longer extents representing more variation.While, the boxplots inform about contribution differences between models, plotting predictor contributions against modelled WT (right-hand side of Fig. 5) shows that the contribution variability, for a given model, is in many cases related to WT rather than random (i.e.some predictors are more or less influential depending on thermal conditions).
AT is the main contributor except in winter (second to SH); its median contribution is around 12 % for winter, and 30-35 % for the other models.In all cases, AT contribution increases as WT increases (AT has more influence at warmer WT).
SWR influence is quite constant for all models (medians ranging from +4.5 to 7.5 %; up to a maximum of +15.8 % in winter) except autumn, for which it is very limited (median +0.13 %).Within each model, SWR contribution is fairly stable across the WT range but showing slightly more variability for colder WT.
LWR is the second contributor for the "all seasons" and the summer models.Its contribution is negative except for spring, but in all cases, the contribution decreases as WT increases (i.e.LWR has more influence on colder WT).
WS has a negative contribution for all models except autumn.WS is most influential for colder WT (e.g. down to a minimum of −13.70 % for all seasons model, −11.74 % for summer); its contribution decreases as WT becomes warmer (e.g.around −1 % for most models).WS contributions are more variable for colder WT (i.e. more scatter right-hand side plots; Fig. 5) than for warmer WT.For autumn, WS has limited influence, with its contribution ranging from +0.17 to +0.90 %.
SH contribution is highest in winter (main contributor with median +27.20 %) and for "all seasons", but otherwise limited for the other seasons (medians ranging +2.10 to +7.23 %).SH contributions are independent from WT.
P has limited influence with its contributions ranging from −1.13 % (minimum, spring) to +0.22 % (maximum, winter).Its contributions show very little variability and no pattern in relation to WT.

Role of basin properties
The site-specific coefficients were initially mapped against elevation and permeability to explore basin modification of the WT-climate relationship, and any pattern linked to easting/northing.While there was no clear easting/northing pattern, the maps showed potential associations between coefficients and basin properties.
Then, ANOVA was run on those descriptors to identify the ones significantly associated with the model site-specific coefficients.Associations between meta-properties/descriptors and site-specific coefficients are showed in Table 4.Note that no property was found to be associated with P coefficients in summer.
To quantify the influence of the properties, either alone, or combined, simple linear regressions of the site-specific coefficients were fitted and ranked with AICc following the MMI technique used above.Models are featured in Table 5.The best models are the ones with the lowest AICc (displayed in bold characters); while all models featured are within four AICc points, and hence are considered equally good (Zuur et al., 2009).Depending on the site-specific coefficient, the R 2 range from 0.125 (autumn SWR) to 0.411 ("all seasons" AT).In each case, a single regression (on BFIHOST or AL-TBAR) is the best model AICc-wise, although most of the multiple regressions are within 4 AICc points, and therefore equally valid models.In the UK context, these metaproperties are themselves not independent: (i) high upland basins are more often impermeable because permeable geology predominantly occurs in the UK lowlands; (ii) there are comparatively more larger basins at lower elevations.Results in Table 5 demonstrate this.For the "all seasons" AT coefficient models, single regressions on BFIHOST, ln(AREA), and ALTBAR achieve a R 2 of 0.370, 0.284, and 0.127, respectively, but the multiple regressions with either two or all of them only achieve R 2 within 0.381-0.411.The comparatively small gain when adding several predictors is due to the three properties co-varying.Similar comments can be made on the other models.

Influence of climate drivers
This section discusses results related to the fixed component of the ML models, which provide information on national-  scale patterns (i.e.patterns valid for every sites used in the analysis).As explained above, these patterns would be analogue conceptually to those sought by using cluster analysis or fully pooled regressions but without their shortcomings (e.g.loss of information, issues with dependent observations).The use of ML modelling addressed one of the limitation of empirical regression-based models, for which temperatures are predicted at specific sites only.Note that the four seasonal models are by definition related to the "all seasons" model, since they are based on subsets of the same original dataset; therefore, seasonal patterns are not independent from the "all seasons" patterns.The six climate predictors investigated were identified as significant within the MMI framework (note that MMI applied to the selection of the fixed component part of the ML models only).Standard model selection techniques (e.g.stepwise) would have most likely excluded the predictors that are not retained in all models of the MMI selected model sets (i.e.predictors with lower RI values).In this regard, this study illustrated how MMI can be useful in picking the effect of secondary controls, otherwise masked by dominant primary drivers.
The models broadly make sense against known physical processes.In interpreting model results, it important to bear in mind that the aim of the study was to assess the relative empirical associations between WT and the set of climate drivers, therefore the models are not explicitly process-based.In addition, the climate variables are inter-related in some extent (e.g.P associated with more cloud cover, hence reduced SWR and greater SH), and the analysis is based on 3-monthaveraged data, which may cause some aspects of the physical processes to be lost by the averaging (e.g.distinction between variable like SWR, only contributing during daylight and others like LWR contributing continuously).
All models flag a close association between AT and WT.This finding is consistent with the literature; it is well doc-umented that AT and WT are both influenced by similar climatic drivers (e.g.incoming radiation), and tend towards thermodynamic equilibrium (Caissie, 2006).Both variables consequently tend to co-vary positively, making AT a very useful predictor (as it has been widely demonstrated in the literature; e.g.Webb and Nobilis, 1997), although the association is only partly causal (Johnson, 2003).SWR (insolation from sun) is physically a positive input of energy; and it is appropriately captured in the models with positive coefficients.In this study, LWR is the downward component of long-wave radiation (see Table 2).From an energy budget perspective, LWR therefore corresponds to a positive flux toward the river water.Consequently, LWR contribution to WT should be positive.Results (Table 3 and Fig. 5) show this is not necessarily the case.LWR corresponds to radiation diffused by clouds, and therefore co-varies positively with cloud cover (in addition, a pairwise plot of the study dataset shows that within a given season LWR inversely co-varies with SWR).Therefore, the negative WT-LWR associations would either be due to LWR acting as a proxy for processes driving colder water temperatures (e.g.cloud cover), or be a model artefact due to the LWR-SWR collinearity.SH represents the mass of water vapour in moist air.The rate of evaporation at the water surface is directly proportional to the SH gradient (the more humid the air, the lower the evaporation rate).All models give a positive association between SH and WT.As SH increases, the evaporation rate decreases, and consequently, cooling due energy loss as latent heat decreases as well.WS has a cooling effect by increasing evaporation at the water surface, which would be captured by a negative contribution to WT.In addition, WS plays a significant role in air-water energy exchanges by increasing mixing, which would manifest as increased cooling or warming depending on the AT-WT gradient.For all models but autumn, WS has an overall negative contribution (cooling).For the autumn model, the variable RI and its percentage contribution are both low, so the positive association has to be considered with caution.P have positive or negative coefficients depending on model.When rainfall occurs, its temperature may be higher or lower than that of the river depending on season.In addition, P can also act as a proxy for cloud cover, thus for reduced SWR and increased LWR; in some cases it might also capture the effect of increased streamflow and thermal inertia.P has limited importance and percentage contribution in all the models, which is probably due to precipitations being event-based, whereas other variables are continuous (e.g.AT).
The form and strength of the climate-WT association vary depending on season and on WT range, as showed by the variability in predictor coefficients and contributions.This most likely captures that the dominating climate drivers and physical processes (e.g.evaporation/condensation, radiative fluxes; see energy budget above) may change from one season to another, or within the same season, from colder to warmer weather conditions.As a consequence, the impact of short (e.g.seasonal climatic drought) and long-term climate variability or change, as well as of mitigation schemes (e.g.increasing riparian tree shading) on stream temperature, may not be uniform across time (e.g. higher long-term temperature increases in winter and spring; Langan et al., 2001).
Probably because AT performs very well as a predictor (e.g.Webb and Nobilis, 1997), most empirical models have been based on single AT-WT regressions (Caissie, 2006) with very few using other climate predictors (e.g.AT and solar radiation; Jeppesen and Iversen, 1987).The present study demonstrated the potential of several other climate variables to contribute explanatory power (even if they are weaker predictors than AT), which can be beneficial when trying to tease out the relative influences of the various interconnected processes controlling water temperature regimes.Although this was not the primary objective of the study, the models could be used to generate seasonal water temperatures for the whole spatial and temporal extent of the CHESS datasets (whole country, 1971-2007 period of records), for example allowing one to investigate broader geographical pattern, or the impact of extreme events like drought.

Role of basin properties
The analysis of the random component of the models (i.e.site specific) identified permeability, elevation, and basin size as the main modifiers of the climate-WT response (note that unlike for the fixed component, the random predictors were selected using standard AIC, i.e. there is only one random component formulation for each of the five models).The use of ML modelling addressed the limitations of empirical regression-based models to work across different spatial scales (see above; Caissie, 2006).The basin properties are first reviewed individually, then together to assess how their respective influences may combine within a basin (i.e. are all influences cumulating, or one property dominating?)For all models and for all predictors (all seasons AT, autumn SWR, winter SH), the more (less) permeable the basin, the lower (higher) the coefficients.Thus, water temperature in impermeable basins appears to be more sensitive to seasonal climate data than in permeable basins.Indeed, in permeable basins, the temperature regime is comparatively more influenced by the groundwater input to the river; groundwater temperature tends to have more inertia and to have a dampening effect on river WT (groundwater warmer than river in winter, cooler in summer) -see, for example, Webb andZhang (1999), Hannah et al. (2004), Caissie (2006), andKelleher et al. (2012).This pattern is consistent with Garner et al. (2014), who used different temperature monitoring sites and basin properties to investigate air-water temperature associations only.
With regard to basin size, with the "all seasons" model, WT in smaller basins is more sensitive to AT but less sensitive to SWR than in larger basins.With the autumn model, WT in smaller basins is more sensitive to SWR.With the winter model, WT in smaller basins is more sensitive to SH.
Although, there are seemingly contradictory patterns for SWR, this can be explained by the modelling.Where studies typically use only one variable to represent the whole climate (e.g.AT, Garner et al., 2014), several climate predictors are considered herein.As noted in the Introduction, AT and SWR co-vary to some extent.In the "all seasons" model, AT and SWR were both selected to capture the between-site variability of the climate-WT response, while in the autumn model, only SWR was retained.As a consequence, in the autumn model, SWR represents climate control, most probably capturing part of the WT variability explained by AT when both variables are included as in the "all seasons" model.Overall, WT is more sensitive to seasonal climate data in smaller basins.Then, the inclusion of both AT and SWR in "all seasons" allows one to refine the assessment of river thermal sensitivity beyond climate as a whole, to different types of energy processes: smaller streams are more sensitive to air-water heat exchanges but less sensitive to radiative fluxes than larger streams.One can hypothesize that smaller streams have a lower volume of water to heat up than larger streams but also are likely to experience greater relative shading by riparian trees than wider rivers downstream.
This finding, at first, looks partly inconsistent with Garner et al. (2014), who concluded that larger basins were more sensitive to climate than smaller ones because (i) headwater stream being located at the start of the network have less time than larger streams to reach equilibrium with AT further downstream, and (ii) headwater streams are more likely to be shaded (riparian woodlands, topography).However, Garner et al. ( 2014) was based on cluster analysis; small basins were included in one cluster only, which also included permeable basins.As a consequence, it is likely that permeability and size influences were in some extent confounded.In contrast, the sites used in this paper cover all combinations of size/permeability basin types.Second, as noted by Kelleher et al. (2012), within the small stream type, one needs to distinguish between shaded (i.e.due to with riparian woodland or topography) and exposed streams, with shaded streams behaving more like permeable streams.Only basin-wide land cover information was available for 29 out of 35 sites: 27 basins are under 20 % woodland.While one cannot exclude woodland being concentrated on the riparian corridor of each site, it is sensible to assume the 35 sites have a mix of shaded and exposed streams.Although it would explain the pattern with "all seasons" SWR (more shading, less incoming sun), the shaded headwater argument has to be considered inconclusive in relation to the wider climate controls.
With regard to basin elevation, results can be summarized as follows: (i) "all seasons" model, WT in higher elevation basins is more sensitive to AT but less sensitive to SWR; (ii) winter model, WT in higher elevation basins is more sensitive to SH.These patterns can be explained partly by elevation, partly by the fact that permeability, size and eleva-tion are not strictly independent in the UK.As noted above, elevation and rainfall co-vary greatly in the UK, so that upland basins are wetter than lowland basins, and hence associated with greater precipitation (i.e. with more cloud cover and consequently, less influenced by SWR).In terms of basin types, the study sites have no upland permeable basins (the UK geology is such that this type hardly occurs in any case), plus high elevation basins tend to be smaller basins.The patterns observed with elevation, which are consistent with those for permeability and size, are most likely partly reflecting the upland basins are also largely impermeable and smaller.
Although each property has been statistically identified as having an influence, the latter point leads to investigating how these influences may combine.The regression models of site-specific coefficients against permeability, size, and elevation presented in Table 5 provide some quantification of the influence of basin properties, both on their own, and combined.In each case, the best model uses a single basin property, although the retention of other properties in the MMI sets confirms the role of all three.In three cases out of four ("all seasons" AT, autumn SWR, winter SH), permeability (BFIHOST) is dominant.Therefore, the patterns described above would be primarily set by basin permeability, then by size and elevation.At one end of the spectrum, small, upland, and/or impermeable basins are the most exposed to atmospheric heat exchanges, at the other end, large, lowland, and permeable basins are the least exposed.

Conclusions
By focusing on a nation-wide set of water temperature sites and extensive climate dataset, this study addressed some of the limits of previous UK papers (limited number of WT sites, climate predictors, and/or geographical extent); it also investigated formally seasonal patterns, and, by using a wide range of basin descriptors, improved knowledge of the role of basin properties as modifiers of climate-WT associations.
With regards to the need to explore alternative modelling techniques to maximize data utility, ML modelling allowed us to model climate-WT responses both at site and at national scales, thereby addressing the limitation of empirical regression-based models compared to deterministic models (Caissie, 2006).While the present ML models took into account discrepancies in temperature sampling (e.g.data from sites with 15 min recording may show different patterns from sites with weekly data), the effect of these discrepancies were not investigated explicitly, and would merit further research.In addition, the model selection based on the MMI approach permitted us to investigate climate variables that would been most likely excluded by standard selection techniques, and identify their influence as secondary controls.
In relation to research aim 1 (improved understanding of large-scale climate-WT associations), the modelling exer-cise showed that all of the six climate predictors investigated in this study play a role as a control of water temperature.AT and SWR are important for all models/seasons, while LWR, SH, and WS are important for some models/seasons only.The form and strength of the seasonal climate data stream temperature association vary depending on season and on water temperature.The dominating climate drivers and physical processes may change across seasons, and across the stream temperature range.The impact of climate variability or change, whether short-or long-term (e.g.seasonal supraseasonal, or inter-annual climatic drought, long-term air temperature increase), and the benefit of mitigation measures (e.g.increasing shading) on stream temperatures need to be assessed accordingly.While this study focused on wider spatial patterns, it is noteworthy that stream temperature could also be influenced by micro-climate effects (as far as metadata could be scrutinized, the study sites were free of such effects), future research could investigate how micro-climate and climate data spatial resolution may influence the models.
In relation to research aim 2 (assessing influence of basin properties as modifiers of climate-WT associations), the study confirmed the role of basin permeability, size, and elevation as modifiers of the climate-WT associations.The primary modifier is basin permeability, then size and elevation.Smaller, upland, and/or impermeable basins are the ones most influenced by atmospheric heat exchanges, while the larger, lowland and permeable basins are least influenced (note that some basin types occur less frequently or hardly in the UK, e.g.upland permeable).This means that, in addition to seasons and temperature range, the impact of seasonal climate data on stream temperatures and the benefits of mitigation schemes may vary with location.This study shows the importance of accounting properly for the spatial and temporal variability of climate-stream temperature associations and their modification by basin properties.

Figure 1 .
Figure 1.Multiple interdependent climate controls of water temperature; the sum of K and L corresponds to Q * (heat flux due to net radiation); Q a corresponds to advective heat exchanges, which include precipitation (direct climatic component) and smaller fluxes due to inflow/outflow into river, hyporheic exchange, or groundwater (not shown on figure); (adapted fromCaissie, 2006, andHannah et al., 2008).

Figure 2 .
Figure 2. Location map of the study sites.

Figure 4 .
Figure 4. Plots of observed and modelled water temperature for the five models.

Table 1 .
Climate-water temperature studies carried out in the UK.
* Includes different measurements of related climatic variables.

Table 3 .
Generic response for the five average models.AICc values have similarly good performance, it is not statistically correct to keep the lowest AICc model only as the best model and discard the others.MMI addresses these issues by selecting sets of good models.
* Tested on natural log.

Table 5 .
Linear regressions of site-specific coefficients as function of basin properties (models ordered by increasing AICc; best model in bold characters, all other models are within four AICc points of best model and hence selected via MMI).