Development of spatial regression models for predicting summer river temperatures from landscape characteristics: Implications for land and fisheries management

There is increasing demand for models that can accurately predict river temperature at the large spatial scales appropriate to river management. This paper combined summer water temperature data from a strategically designed, quality controlled network of 25 sites, with recently developed flexible spatial regression models, to understand and predict river temperature across a 3,000 km2 river catchment. Minimum, mean and maximum temperatures were modelled as a function of nine potential landscape covariates that represented proxies for heat and water exchange processes. Generalised additive models were used to allow for flexible responses. Spatial structure in the river network data (local spatial variation) was accounted for by including river network smoothers. Minimum and mean temperatures decreased with increasing elevation, riparian woodland and channel gradient. Maximum temperatures increased with channel width. There was greater between‐river and between‐reach variability in all temperature metrics in lower‐order rivers indicating that increased monitoring effort should be focussed at these smaller scales. The combination of strategic network design and recently developed spatial statistical approaches employed in this study have not been used in previous studies of river temperature. The resulting catchment scale temperature models provide a valuable quantitative tool for understanding and predicting river temperature variability at the catchment scales relevant to land use planning and fisheries management and provide a template for future studies.


| INTRODUCTION
Water temperature (Tw) is an important control on the survival and growth of freshwater fish, with consequences for species distribution, abundance, demographic characteristics and production. This is particularly the case for cold water adapted species including salmonids (Comte, Buisson, Daufresne, & Grenouillet, 2013;Jonsson & Jonsson, 2009;McCullough et al., 2009;Ruesch et al., 2012). Rising Tw has the potential to alter the thermal suitability of rivers for salmonids (Isaak et al., 2010;Isaak, Wollrab, Horan, & Chandler, 2012;Mohseni, Stefan, & Eaton, 2003), which are often the focus of environmental protection and management action. In addition to their ecological importance, salmonid species are frequently associated with a high cultural, social, conservation and economic value. For example, in Scotland, there are 17 rivers designated as special areas of conservation (SAC) for Atlantic salmon (Salmo salar) under the European Union Habitats Directive (Anon, 2009), and Atlantic salmon and sea trout (Salmo trutta) anglers are estimated to contribute over £73 million per year to the Scottish economy (Radford, Riddington, & Anderson, 2004). There is therefore increasing interest in understanding and predicting the spatio-temporal variability of river thermal regime (Webb, Hannah, Moore, Brown, & Nobilis, 2008), the likely effects of climate change and opportunities for mitigation and management, such as riparian tree planting, to inform fisheries management (Hrachowitz, Soulsby, Imholt, Malcolm, & Tetzlaff, 2010;Malcolm et al., 2008).
Spatial statistical models represent the most promising approach for understanding and predicting river temperature at larger spatial scales and for informing and refining sampling network design. However, there have been relatively few studies at the catchment (Chang & Psaris, 2013;Hrachowitz et al., 2010;Imholt et al., 2011Imholt et al., , 2013McNyset, Volk, & Jordan, 2015;Steel, Sowder, & Peterson, 2016) or regional (Hill, Hawkins, & Carlisle, 2013;Isaak et al., 2012;Ruesch et al., 2012) scales necessary for management and fewer still have provided a holistic assessment of thermal regime necessary for improved understanding of ecological processes (Imholt et al., 2011;Malcolm et al., 2008;Steel et al., 2016). In addition, most previous studies have considered only linear responses between Tw and catchment covariates, even though more complex asymptotic or modal responses (i.e. smooth terms) may be required. This paper uses temperature data collected from the summer of 2015 to produce spatial regression temperature models for the River Spey catchment, an SAC for Atlantic salmon in the North East of Scotland.
The paper aims to improve understanding of the spatial variability in temperature across the catchment, highlight areas of potential risk for Atlantic salmon under climate change (i.e. hotter areas in the catchment) and identify optimal locations for riparian tree planting based on model outputs and a conceptual understanding of energy exchange processes.
Secondly, the paper aims to characterise thermal heterogeneity for different river orders to inform future temperature sampling strategies.
These aims were addressed through the following specific objectives: 1.1 | Objectives 1. Develop regression-based models to predict temperature metrics (minimum, mean and maximum) from landscape characteristics representative of energy exchange processes.
2. Understand how catchment scale variability in landscape covariates influences Tw and infer underlying processes.
3. Predict temperatures for unmonitored locations across the Spey catchment and highlight the areas at greatest risk of high temperatures and identify areas suitable for riparian planting.
4. Characterise thermal heterogeneity at different spatial scales as indicated by river order and assess consequences for monitoring.
The network length is around 36,500 km with the main stem constituting 180 km (Hastie et al., 2004). Heather moorland, rough pasture, pastoral farming, both conifer and deciduous woodland and wetlands are the dominant land covers (Hastie et al., 2004). The catchment is relatively un-impacted by artificial flow regulation (Scottish Natural Heritage, 2001).
The Spey catchment is designated as a SAC for Atlantic salmon, Otter, Sea Lamprey and Freshwater Pearl Mussel (Butler et al., 2009), and the main stem is a Site of Special Scientific Interest (Scottish Natural Heritage, 2001). The River Spey has been described as one of Europe's premier salmon and trout rod fisheries (Sandison, 2001) with a study from 2003 estimating the value to be around £10.1 million annum −1 to the local economy (Butler et al., 2009).

| Water temperature data and metrics
Water temperature was monitored at 25 sites, using cross-calibrated Gemini TinyTag Aquatic 2 (TG-4100) dataloggers, during the summer of 2015 ( Figure 1). The sites are part of the Scotland River Temperature Monitoring Network (SRTMN) and designed to cover a wide environmental range of potential landscape controls (Jackson, Malcolm, & Hannah, 2016). The measurement precision of the temperature dataloggers is 0.01°C. Further details on this network, its design and quality control procedures, can be found in Jackson et al. (2016). Data were collected at 15-min intervals between 22/06/2015 and 31/08/ 2015 and used to produce three summary metrics. These metrics were (a) Tmin: the minimum 7-day moving-average of daily minimum temperatures, (b) Tmean: the mean of daily mean temperatures, and (c) Tmax: the maximum 7-day moving-average of daily maximum temperatures. Moving averages of Tmin and Tmax were used in preference to single daily values, to indicate locations associated with sustained high or low temperatures in agreement with previous studies (e.g. Moore, Nelitz, & Parkinson, 2013). The combination of metrics will allow a broader understanding of the hydrological processes influencing Tw variability, extending beyond the more common focus on maximum temperatures. This can be used to infer the likely consequences of catchment scale Tw patterns on salmonid fish (e.g. growth, mortality) using existing information on thermal preferences and thresholds (see Elliott and Elliott (2010) for detailed discussion of the temperature requirements for Atlantic salmon and brown trout).

| Model covariates
The Spey river network was characterised using a digital river network (DRN; Centre of Ecology and Hydrology, 2014) consisting of a series of line features connected at nodes. Nodes exist at all river sources and confluences as a result of the digitisation process. In common with other DRNs, additional "pseudo nodes" (i.e. nodes present on river segments which were not a source or confluence) were also present throughout the catchment . Additional nodes were generated for each SRTMN monitoring site. Covariates were calculated for all SRTMN sites (for model fitting) and nodes (for prediction), except for source nodes, where it was not possible to obtain covariates that required upstream information, for example, channel gradient (see below).
To enable catchment wide prediction of Tw, a suite of landscape characteristics was used to model river temperature. These represented proxies for the physical processes affecting Tw (Jackson et al., 2016) and have often been used in previous regression-based studies (Chang & Psaris, 2013;Hrachowitz et al., 2010;Imholt et al., 2011).
These covariates were as follows: elevation (elevation), upstream catchment area (UCA), percentage riparian woodland (%RW), hillshading/channel illumination (HS), channel width (width), channel gradient (gradient), channel orientation (orientation), distance to coast (DC) and distance to the sea along the river (RDS). Detailed discussion of the thermal processes represented by these proxies can be found in Jackson et al. (2016).

| Covariate calculation
Covariates were characterised using R, version 3.2.3 (R Core Team., 2015) except where specified. For each node, points were established 1,000 m upstream, by travelling up the DRN. Where the upstream path encountered tributaries, it remained on the river line associated with the greatest Strahler river order (order). Where the upstream path encountered confluences with the same order, it proceeded up both tributaries. A river "segment" line feature was established between the nodes and the upstream points, and these segments were then snapped to the nearest ordinance survey (OS) MasterMap river features (lines or polygons). Finally, a buffer was extended 50 m either side of the river features. This combination of points, lines and buffers enabled various covariates to be extracted in a meaningful way. In addition, the approach did not require the removal of braids or complicated confluences (i.e. >2 tributaries joining in the same location) from the DRN.
Elevation was extracted as a point value for each node from the OS.
Terrain 10 digital terrain model, (DTM), 10-m resolution, using the "extract" function in the "raster" package (Hijmans, 2015). Gradient was determined as the difference in elevation between nodes and the associated upstream points, divided by the segment length. Orientation was calculated using the x and y locations of the node and associated upstream points using standard trigonometry. An average of the gradients and orientations was obtained where there was more than one upstream location, that is, where the 1-km segment branched up tributaries.
A UCA raster was created in ArcGIS 10.2.1 using Arc Hydro Tools.
The river network was "burned in" to the DTM to improve correspondence between the raster and DRN and thus catchment delineation (Li, 2014;Peterson, Sheldon, Darnell, Bunn, & Harch, 2011). UCA was obtained by returning the maximum observed UCA in a circular buffer around each node. The size of the buffer corresponded to river order, optimised to maximise the chances of including appropriate values where the DTM and DRN were in poor agreement (e.g. flood plains) or where braided rivers occurred.
Hillshading rasters that provided a metric of potential solar radiation were created using the "terrain" and "hillShade" functions in the "raster" package (Hijmans, 2015) and the available DTM. Solar azimuth and altitude input values were obtained from the U.S. Naval Observatory Astronomical Applications Department (Anon, 2001) for the centre point of the catchment, for every hour the sun was above the horizon. The hillshading raster layers for each time period were then FIGURE 1 Site map of the river Spey catchment with Scotland River Temperature Monitoring Network sites shown as black dots and salmon rivers (where salmon are present or likely present) defined by (Gardiner & Egglishaw, 1986) in dark blue. Non-salmon rivers are shown in light blue. Woodland is overlaid in green summed to create a single layer providing a metric of total potential solar exposure. HS values were calculated for each node by averaging the hillshading raster values in the OS MasterMap river polygons that extended 1-km upstream. If the river was characterised only by a line in the OS MasterMap dataset (meaning they are smaller than the minimum width required for polygons), a 0.5 m buffer was extended either side of river line features, to create a river polygon. Raster cells were weighted by the proportion of the cell within the buffer.
The %RW was extracted by calculating the percentage of the upstream buffers containing MasterMap polygons classified as woodland. The 1-km buffer distances were chosen based on distances used in previous regression-based studies (Hrachowitz et al., 2010;Imholt et al., 2011). Width was calculated by extracting the area of MasterMap water polygons within the 1-km upstream distance and dividing this by the length of the river segment (which could be >1 km in the case of river confluences). All areas were calculated using "gArea" from the "rgeos" R package (Bivand & Rundel, 2016) and lengths using the "SpatialLinesLengths" tools from the R "sp" package (Pebesma & Bivand, 2005).The shortest DC was calculated for each node, using "gDistance" from the "rgeos" R package (Bivand & Rundel, 2016). RDS was calculated by inputting the DRN into the igraph R package and using the "shortest.paths" function (Csardi & Nepusz, 2006).

| Data preparation
Covariates were transformed as required to reduce skewness; this resulted in logging of both UCA and width, taking the square root of gradient, with all negative gradients (<0% change) being made 0. HS was also centred, by subtracting the median from all values, to improve the numerical stability of model fits. Before model selection, Pearson correlation coefficients were used to assess the degree of correlation between potential explanatory covariates. Where a correlation of >0.75 was observed between covariates, one was removed. Given the relatively linear structure of the River Spey, RDS and DC were removed from the analysis due to high correlations with elevation (0.89 and 0.84, respectively). UCA and width were also highly correlated (0.94); thus, UCA was removed from the potential covariates list.

| Model fitting and river network smoothers
All analyses were undertaken using R, version 3.2.3 (R Core Team. 2015). Models were fitted using generalised additive models with Gaussian errors where the amount of smoothing was estimated from the data. The models were fitted by maximum likelihood.
Spatial structure in the data, which was not accounted for by the covariates, was considered by fitting a smoother over the river network, using a modified version of the methods described by O'Donnell, Rushworth, Bowman, Scott, and Hallard (2014); herein referred to as the river network smoother (RNS). O'Donnell decomposes a river network into small stream segments and models changes in the response variable between connected segments using penalized regression splines. In particular, where two stream segments meet, smoothness across the confluence is controlled by a penalty that depends on the relative flows of the two joining units. Two modifications were made for the Spey analysis. First, in the absence of spatially distributed discharge data, smoothness across a confluence was controlled by the proportional influence of upstream tributaries, weighted by Strahler river order (Strahler, 1957). A unit increase in river order corresponded to a doubling of flow (Hughes, Kaufmann, & Weber, 2011). River order provided a pragmatic and readily derived weighting and was proposed by Ver Hoef, Peterson, and Theobald (2006) and O'Donnell et al. (2014) for application when discharge data are unavailable. Second, dimension reduction techniques were used so that the models could be fitted using the GAM function in the R "mgcv" package (Wood, 2001). Specifically, the full smoother matrix of the RNS (as described in O'Donnell et al., 2014) was rotated, and the leading eigenvectors, scaled by their eigenvalues, were used as a set of "reduced rank" basis functions that span the river network, with their coefficients penalized by a diagonal penalty (see Wood (2006), page 309). This RNS could then be fitted within the "mgcv" package using a purpose built smooth constructor function, with the amount of smoothing estimated by maximum likelihood.
Because the covariates included in the model selection were carefully chosen to represent physical processes and to ensure that the RNS only incorporated variability that could not be explained by the covariates, any RNS basis functions that were strongly correlated with the covariate responses (>0.75) were excluded from modelling. As such, RNS base 1 and 2 were disregarded due to correlation with elevation, and the next 10 basis functions were then used for modelling.

| Model selection and performance
Given the limited size of the dataset and the expected simplicity of the landscape response relationships, each landscape covariate was only allowed up to 2 degrees of freedom, allowing for responses ranging between linear and modal. The RNS was permitted up to 10 degrees of freedom, which was the maximum that allowed the RNS and all covariates to be fitted in a single model. In practice, this was plenty, given that the effective degrees of freedom (those used by the RNS in the fit) were always much less. No interactions were considered, due to the size of the dataset. All possible model combinations were explored, giving 126 possible models. Due to the small sample size, corrected Akaike information criterion (AICc) was used for model selection (Gallice, Schaefli, Lehning, Parlange, & Huwald, 2015;Hurvich & Tsai, 1989) with the "best" model having the lowest AICc value. AICc values were tabulated for the top 10 models for each metric to identify other candidate models (models with similarly low AICc values).Where smoothed terms in the selected models had an effective degrees of freedom of 1, they were replaced with linear terms.  (Table 2). Prediction error was indicted by the root mean square prediction error. Bias was calculated by taking the mean of the prediction errors. Mean and median absolute deviations were also calculated as the sign independent mean and median of the prediction errors.

| Visualising and interpreting spatial patterns in Tw
Maps of predicted Tmin, Tmean and Tmax were used to illustrate the spatial variability of Tw. In particular, the map of Tmax provides an indication of risk that can also be used to explore opportunities for mitigation. To identify the level of thermal heterogeneity observed at different spatial scales (Objective 4), the predictions for each node and metric were summarised by river order (min, median, max 5th and 95th percentiles) and illustrated in box plots.
The effect of individual covariates on Tw was illustrated using plots and maps of partial effects, where Tw predictions were made for each covariate with the others held at median values. Temperature predictions were restricted to river nodes where the environmental characteristics were within the environmental range observed at temperature monitoring sites, thereby preventing extrapolation beyond the range of the data (Jackson et al., 2016). In practice, this approach constrained prediction to salmon rivers (Figure 1) which were the focus of the network design and the likely target of future management action (Isaak et al., 2010(Isaak et al., , 2012Mohseni et al., 2003).

| RESULTS
The summer of 2015 was wetter than average, and this was reflected in observed Tw. Regional data collected by the MET Office for Eastern Scotland suggested that summer rainfall was 124% of the 1981-2010 mean and that maximum, minimum and mean air temperatures were 0.4°C, 0.5°C and 0.5°C lower than average over the same period.

| Temperature models
The top 10 models for each temperature metric are shown in Appendix. The final models (lowest AICcs) were the following: Tmin and Tmean were highly correlated (0.97). It is therefore unsurprising that the final model for both metrics was the same (Equa-  (Table 1). There was only weak evidence of an effect of gradient (Table 1). Cross-validation showed negligible bias, small median absolute deviations (≤1°C) and root mean square prediction error for the two models ( Table 2). Plots of residuals against Given the similarities between the Tmin and Tmean models (Figure 3), the predicted patterns of spatial variability were also very similar and varied by 7.3°C and 7.0°C across the catchment (Figures 4a and 5a). The warmest Tmin and Tmean were predicted in the mainstem reaches, predominantly near the river mouth and coolest in the headwater tributaries (Figures 4a and 5a). This primarily corresponded with lower elevations and channel gradients in the lower mainstem areas of the catchment (Figure 3). Figures 4 and   (Figures 4d and 5d).
The partial effect of %RW was "patchier," than that observed for elevation or gradient, reflecting the local distributions of trees (   Figure 1). For both models, the standard error of predictions was generally small (<0.5°C) across the catchment (Figures 4b and 5b).     across second-order streams and 3.5°C across seventh-order streams.

| Thermal heterogeneity and spatial scale
The range of Tmax was greatest across third-order streams (7.0°C) and lower across higher-order streams (5.6°C and 3.6°C in order 5 and 7, respectively; Figure 7c).

| Modelling approach
In contrast to data intensive process-based models (e.g. Garner, Malcolm, Sadler, Millar, & Hannah, 2015), statistical models of Tw that incorporate landscape proxies for energy exchange processes have the potential to predict Tw at the large spatial scales appropriate to river management with limited field-based data collection. Historically, large-scale spatial regression models have assumed only linear responses between Tw and landscape covariates and have ignored, or have not needed to deal with, spatial correlation in the data (e.g. Hrachowitz et al., 2010;Imholt et al., 2011;Isaak & Hubert, 2001;Isaak et al., 2010;Mayer, 2012;McNyset et al., 2015). More recently, increasingly sophisticated geostatistical modelling approaches have been developed and applied to Tw data (Isaak et al., 2014). These models account for spatial correlation with the associated benefits for model selection and Tw prediction but typically still assume linear responses between Tw and covariates (e.g. Detenbeck, Morrison, Abele, & Kopp, 2016;Roberts, Fausch, Peterson, & Hooten, 2013;Ruesch et al., 2012;Steel et al., 2016). An alternative approach for addressing spatial correlation involves the use of river network smoothers (O'Donnell et al., 2014). The current study extends the RNS approach to also include smoothed responses between Tw and landscape covariates. As far as the authors are aware, this paper is the first to incorporate these recent developments in a study of Tw with consequent benefits for Tw process understanding and prediction.

| Physical interpretation of temperature models and spatial patterns of predictions
Elevation was a significant covariate in many previously published regression-based Tw models (Chang & Psaris, 2013;Hrachowitz et al., 2010;Imholt et al., 2011;Isaak & Hubert, 2001;Ruesch et al., 2012;Steel et al., 2016;Trumbo et al., 2014). In this study, elevation had a significant negative effect on Tmin and Tmean, producing strong and consistent spatial patterns. Elevation was also present in the second best Tmax model (Appendix). Elevation represents adiabatic lapse rates, which reduce air temperature (Ta) and thus Tw with increasing altitude (Hrachowitz et al., 2010). Dry adiabatic lapse rates are ca.
9.8°C/km, and there are potentially Ta differences in the order of 12.8°C across the whole Spey catchment or 5.0°C across the SRTMN sites used in this study. This is remarkably similar to the estimated elevation effect size of 5.8°C and 4.4°C for the Tmin and Tmean models respectively. Although Ta is not the main control on Tw, both are influenced by similar controls and are therefore often strongly correlated (e.g. Krider, Magner, Perry, Vondracek, & Ferrington, 2013).
Percentage riparian woodland was also a significant predictor of Tmin and Tmean. Shading reduces the amount of incident shortwave radiation reaching the river, decreasing rates of warming in shaded reaches during daylight hours (Garner, Malcolm, Sadler, & Hannah, 2014;Hill et al., 2013;Moore, Spittlehouse, & Story, 2005). Shaded reaches also experience reduced wind speeds, longwave and evaporative heat losses relative to more open moorland locations, reducing thermal variability and increasing nocturnal minimum temperatures Moore et al., 2005). Consequently, %RW has also been a significant covariate in a number of previous regression-based Tw models, particularly during summer months (Chang & Psaris, 2013;Hrachowitz et al., 2010;Imholt et al., 2011;Isaak & Hubert, 2001;Trumbo et al., 2014). Importantly, datalogger site selection for this study incorporated a range of %RW values spread across the catchment, thereby allowing this effect to be separated from other spatial covariates (Jackson et al., 2016). At 2.6°C (Tmin) and 2.4°C (Tmean), the effect size for %RW was also plausible and broadly comparable to previous studies (Garner, Hannah, Malcolm, & Sadler, 2012;Hannah et al., 2008;Malcolm et al., 2008;Simmons et al., 2014).
In contrast to previous investigations (e.g. Hrachowitz et al., 2010), %RW did not appear in the best Tmax model, although it was present in two of the top 10 models (Appendix). This could reflect the overriding importance of other processes controlling maximum river temperatures, a failure to precisely characterise %RW from available mapping resources or an inability to allow for appropriate interactions between covariates given the relatively limited availability of data (monitoring sites).
The inclusion of gradient in the final Tmin and Tmean models is again consistent with previous studies (Chang & Psaris, 2013;Hill et al., 2013;Hrachowitz et al., 2010;Imholt et al., 2011;Mayer, 2012;McNyset et al., 2015) although the small δAICc and p-values give only relatively weak support for this effect. Nevertheless, the negative relationship between gradient and Tw in the Tmin and Tmean models is physically interpretable as transit times are longer in lower gradient channels (notably the mainstem), providing greater opportunities for warming (Webb et al., 2008). In contrast to elevation, the effect of gradient (which was moderately correlated with Elevation, 0.54) was more spatially variable locally, reflecting the presence of landscape features such as bedrock outcrops and post-glacial moraines that create spatially heterogeneous gradients.
Width was the only significant landscape covariate in the final model of Tmax. Although previous studies have discussed the potential importance of width, or more precisely width-depth ratios for energy exchange processes (Imholt et al., 2011), they have not typically characterised this potential covariate. This is presumably due to the difficulties associated with obtaining suitable river size datasets and extracting representative river width data. In this study, a spatial river polygon dataset was used to characterise river width. However, it did not allow characterisation of the widths of smaller rivers and thresholds at which rivers were represented as lines was not spatially consistent (e.g. <1 m in urban areas and <2 m in rural areas). Furthermore, there was strong correlation between width and UCA in the Spey catchment which meant that it was not possible to specifically attribute spatial variability in Tmax to spatial variability in width. Rather, the width metric should be considered as a composite measure of width and UCA, which is in turn a proxy for discharge, water volume and thermal capacity Ver Hoef et al., 2006).
Larger water volumes have a greater thermal capacity, taking longer to warm but also retaining heat for longer (Imholt et al., 2011) making larger waterbodies more thermally stable but susceptible to antecedent conditions.

| Tw predictions, prediction accuracy and sources of error
All three models were associated with good measures of fit and low standard error that compare favourably with previous studies (Hrachowitz et al., 2010;Imholt et al., 2011;McNyset et al., 2015).
As expected, higher standard errors were observed at the extremes of the environmental range. The low and consistent bias was also reassuring as some previous studies have found slight bias towards under predictions at high temperatures (Isaak et al., 2010;McNyset et al., 2015).
The lack of a RNS in the final models for Tmin or Tmean indicates that there was no evidence of substantial residual spatial correlation in these datasets (McGuire et al., 2014) and that covariates explained much of the variation in Tmin and Tmean. However, an RNS was included in the final Tmax model indicating substantial spatial variability relating to river network structure that could not be explained by the covariates (Steel et al., 2016). Accounting for spatial structure within Tw models where present substantially improves prediction accuracy where data exist (Isaak et al., 2010(Isaak et al., , 2014Peterson & Urquhart, 2006). This is shown by the low standard errors on predictions for rivers and tributaries containing SRTMN sites, but higher prediction errors for tributaries without data where there was no information to constrain the smoother.
This was especially the case where there were large changes in river order (e.g. 7 to 2), resulting in low RNS weightings, potentially allowing rapid rates of change in Tw.
Covariate characterisation was another potential source of uncertainty in the Tw models (Gallice et al., 2015;Millar, Millidine, Middlemass, & Malcolm, 2015). For example, differences in the scale of characterisation between river lines and land use maps or imprecise DTMs directly influence the covariate values (Gallice et al., 2015).
Detailed discussion of these potential errors and methods to correct them can be found in Millar et al. (2015).
4.4 | Thermal heterogeneity across river orders: implications for monitoring networks There are two potential processes by which between-stream variability in Tw might depend on river order. Firstly, low-order rivers are more numerous and spatially extensive than high-order rivers, thereby covering a greater part of the environmental range of landscape controls on Tw. Secondly, high-order rivers will have a greater thermal capacity, reducing their thermal variability in response to local landscape controls (Hrachowitz et al., 2010;Imholt et al., 2011;Isaak & Hubert, 2001). This study shows that the greatest thermal variability occurred at river order 5 and below. Such information should inform future monitoring strategies where the objective is to develop cost-effective approaches for characterising thermal variability across whole catchments. Specifically, future networks may wish to target monitoring effort in line with the variance observed in different river orders with more loggers deployed in low-order rivers and fewer in high-order rivers.

| Opportunities for management and future research
Riparian tree planting has often been suggested as a mechanism for reducing maximum Tw under climate change. Using the temperature maps produced in this study, it would be possible to identify sites that are associated with the highest temperatures. This information could be used as an initial assessment of where riparian planting should be targeted. However, process-based Tw studies have also shown that riparian shading is most effective in reducing temperatures where gradient is low, channel widths are narrow (Caissie, 2006;Hrachowitz et al., 2010;Malcolm, Hannah, Donaghy, Soulsby, & Youngson, 2004;Reiter, Bilby, Beech, & Heffner, 2015), thermal capacity is small (Hrachowitz et al., 2010;Moore et al., 2005) and channel orientation maximises the effects of bankside shading. Such information could provide a metric of "planting potential" to be considered alongside maps of temperature variability. The decision on which tributaries have the greatest planting potential could then be based on the available temperature models and resulting maps. This would greatly improve on previous generalised suggestions of planting headwaters to reduce Tw (Caissie, 2006;Hrachowitz et al., 2010) and ensure more targeted, effective and thereby cost-efficient management action.
To provide a more detailed quantitative assessment of planting Although the models developed in this study were physically plausible, they still require additional expert knowledge to inform optimal planting strategies. The incorporation of interaction terms offers the potential of predicting optimal locations at large spatial scales within a modelling framework.
Given the logistical and financial costs of temperature data collec- River network smoother was given 10 degrees of freedom, and all other covariates were smooth terms with 2 degrees of freedom, if the effective degrees of freedom was 1 or less, the term was subsequently made linear.