Stream thermalscape scenarios for British Columbia, Canada

Abstract Water temperature is a key feature of freshwater ecosystems but comprehensive datasets are severely lacking, a limiting factor in research and management of freshwater species and habitats. An existing statistical stream temperature model developed for British Columbia (BC), Canada, was refit to predict August mean stream temperatures, a common index of stream thermal regime also used in thermalscapes developed for the western United States (US). Thermalscapes of predicted August mean stream temperature were produced for 680,000 km of stream network at approximately 400 m intervals. Temperature predictions were averaged for 20-year periods from 1981-2100 to produce 86 scenarios: one for each historical period (i.e. 1981-2000, 2001-2020), and 21 for each future period (i.e. six global climate models and an ensemble average under three representative concentration pathways). The final model performance was consistent with other published regional-scale statistical models (R2 = 0.79, RMSE = 1.53 °C, MAE = 1.18 °C), performing well given the relative paucity of data, large geographic extent, and range of climatic and physiographic conditions. Model results suggested an average increase of August mean stream temperature of 2.9 ± 1.0 °C (RCP 4.5 ensemble mean ± SD) by end of century, with significant heterogeneity in predicted temperatures and warming rates across the province. Compared to stream temperature predictions from the western US, the predictions for BC showed good agreement at cross-border streams (Pearson’s r = 0.91), suggesting the possible integration of both products for a thermalscape covering much of western North America. These stream thermalscapes for BC address a major data deficiency in freshwater ecosystems and have potential applications to stream ecology, species distribution modelling, and evaluation of climate change impacts.


Introduction
Water temperature is a key factor in aquatic environments, affecting ecological, biological, and chemical processes (Angilletta 2009;Kingsolver 2009), but available datasets are sparse, especially for entire stream networks.This situation represents a major data deficit in freshwater ecosystems, where the complexity inherent in freshwater systems (Ward 1989), and the accompanying heterogeneity of habitats (Ward 1998;Poole 2002), requires data at a large enough spatial extent to capture relevant catchment-scale processes and of a sufficiently high resolution for the scale of investigation (Torgersen et al. 2022).Stream temperature models (see reviews in Benyahya et al. 2007;Dugdale, Hannah, and Malcolm 2017) can potentially address this deficit.
Stream temperature models generally fall along a spectrum spanning two distinct approaches: process-based and empirical (statistical/stochastic).Process-based models have the potential to account explicitly for the complex sets of catchment-and reach-scale hydrological and thermal processes that drive stream temperature variability in both time and space, including lateral and longitudinal advection and vertical energy exchanges at the stream and bed surfaces (e.g.Bosmans et al. 2022;Khorsandi, St-Hilaire, and Arsenault 2022;Larabi, Schnorbus, and Zwiers 2022;Leach et al. 2023).However, process-based models require a substantial investment in model set-up and are computationally intensive, especially at regional and larger scales.Further, when applied over large regions, process-based models may lack the resolution for more fine-grained applications, especially in headwater reaches, where lateral advection via hillslope runoff can exert a major influence (Gallice et al. 2016;Leach and Moore 2017).Large-scale models typically parameterize rather than simulate headwater thermal dynamics; for example, Sun et al. (2015) used a logistic curve relation between stream and air temperatures to specify the upstream boundary conditions.However, this approach can lead to inaccurate simulation at those scales (Leach and Moore 2015).This issue is important given that headwater reaches constitute a significant portion of the total channel length in stream networks (e.g.Gomi, Sidle, and Richardson 2002).
Empirical approaches include a range of model types, including site-specific regression models that relate stream temperature time series to air temperature (e.g.Mohseni, Stefan, and Erickson 1998), multiple regression models that relate the spatial patterns of a static stream temperature metric to catchment attributes and climatic variables (e.g.Moore 2006;Wehrly, Brenden, and Wang 2009), and more complex models that simulate spatiotemporal variability at a regional scale (e.g.Jackson et al. 2018).In North America, the Northwest Stream Temperature (NorWeST) project (Isaak et al. 2017) is the gold standard for stream temperature monitoring and empirical modelling (https://www.fs.usda.gov/rm/boise/AWAE/projects/NorWeST.html).It consists of a large, centralized, quality-controlled database of stream temperature summaries collected from organizations and individuals across the contiguous western United States (US) (Chandler et al. 2016).Those data were then used to fit spatial statistical network models for predicting stream temperature that account for the unique spatial autocorrelation structure of hydrologic networks (Ver Hoef et al. 2006).Thermalscapes, which are continuous maps of water temperature across the stream network, were produced for August mean stream temperatures (AugTw) under a range of historical and future time periods and climate scenarios (Isaak et al. 2016a).The model performance (R 2 ¼ 0.91, RMSE ¼ 1.10 � C, MAE ¼ 0.72 � C; Isaak et al. 2017), spatial extent (western US), and resolution of the thermalscapes (1:100,000 scale stream network) has made them a useful tool for a wide range of user groups (Isaak et al. 2017).The broad use of NorWeST products underscores the value of this type of stream temperature dataset and the need for comparable products in other regions.
British Columbia (BC), Canada, does not have a province-wide dataset of continuous historical or future stream temperatures that is comparable to the NorWeST thermalscape scenarios (Isaak et al. 2016a).This lack is a major data deficit given the province's freshwater resources, which includes 2 million km of streams and rivers and over 700,000 lakes and wetlands (BC Freshwater Atlas (FWA) 2020), and the pressures of a warming climate (IPCC 2021), land use changes (e.g.deforestation), and human uses (e.g.water abstraction, hydro power generation) on the freshwater systems and species that they support.In addition to within-province applications, BC stream thermalscapes that are directly comparable to NorWeST could facilitate production of a thermalscape for much of western North America's stream network, overcoming the issues of (1) freshwater datasets restricted by political boundaries (Domisch et al. 2015a;Domisch, Amatulli, and Jetz 2015b), and (2) artificially-bounded ranges used to assess and model species-environment relationships (Barbet-Massin, Thuiller, and Jiguet 2010;Domisch et al. 2015a).
Following the example of the NorWeST project to produce comparable thermalscapes for BC, while ideal, is not an expedient solution to the current lack of stream temperature predictions for the province.Stream temperature data in BC have typically been collected and maintained by individual user-groups (e.g. government agencies, professionals, community groups), and there is currently no central data repository for stream temperature observations in BC.Further, the necessary spatial datasets to run spatial statistical network models are lacking for BC, including both topologically conditioned stream network layers (Isaak et al. 2013) and important model covariates used in NorWeST, like groundwater influence and stream discharge.An alternative approach would be to use existing BC-specific stream temperature models to produce thermalscapes (e.g.Moore 2006;Moore, Nelitz, and Parkinson 2013); however, the thermal metrics predicted by these models are not directly comparable to NorWeST.
In light of current data limitations and the pressing need for stream temperature predictions, a previously developed statistical model for maximum weekly average temperature (MWAT; Moore, Nelitz, and Parkinson 2013) was refit to cover the provincial stream network and provide a full suite of projections under different climate change scenarios.Further, the model was updated to predict a different thermal metric, AugTw, to align with the NorWeST stream temperature predictions (Isaak et al. 2016a) and evaluate the potential of an integrated thermalscape product for much of western North America.The thermalscapes presented here for BC's stream network are the first in Canada and represent a major step forward for application in aquatic research and management.

Study area
The province of BC covers approximately 950,000 km 2 , spanning 12 � of latitude and a diverse range of climates, physiographic features, and ecological zones (Eaton and Moore 2010;Moore et al. 2010) (Figure 1, Table S1).Land elevations range from sea level to over 4,000 m, and include drainages for some of the largest rivers on the continent (e.g.Fraser, Columbia, and Mackenzie).

General overview
The MWAT statistical stream temperature model (Moore, Nelitz, and Parkinson 2013), hereafter referred to as the 'BC MWAT model', was selected as the foundation for developing BC thermalscapes.Full details of the BC MWAT model can be found in Moore, Nelitz, and Parkinson (2013), but key details are briefly summarized here for context.Moore, Nelitz, and Parkinson (2013) fit a multiple linear regression model to the mean MWAT from 418 monitoring stations throughout BC.Stations had up to nine years of in situ stream temperature data; 64% of stations had only one year of data.For stations with multiple years of data, the average of annual MWATs was used.The catchment area for all stations was between 1 and 10,000 km 2 .Model covariates included catchment area, catchment elevation, channel slope, glacier cover, lake cover, normal summer air temperature, annual summer air temperature anomaly, and k 2 .The k 2 coefficient from Eaton, Church, and Ham (2002) is related to flood intensity and bankfull width, functioning as an index of hydroclimatic conditions.The RMSE from a 10-fold cross-validation was 2.1 � C.
The BC MWAT model was selected for a number of reasons.First, the model was regionally appropriate; it was developed for BC using stream temperature observations from stations across much of the province.Second, the in situ daily observations used to calculate MWATs were of sufficient duration and temporal resolution to calculate AugTw values and refit the model, hereafter referred to as the 'BC AugTw model'.Refitting the model to predict AugTw was necessary so modelled stream temperatures could be directly compared to the NorWeST thermalscapes (Isaak et al. 2016a).Finally, the spatial data layers required as inputs for the model and necessary to predict stream temperatures under current and future climate scenarios were available with sufficient temporal coverage and resolution for the study area.
Besides using a different thermal metric, the approach to refitting the BC MWAT model deviated from the original methods (Moore, Nelitz, and Parkinson 2013) in several notable ways (detailed in following sections).First, an expanded dataset of in situ stream temperature observations was used.Second, the model predictions were made at the stream reach scale (i.e. an uninterrupted section of the stream between two stream junctions, or between a junction and an origin or outlet), instead of for specific stations, to better align with the spatial data framework.Finally, a limited set of additional covariates was considered when refitting the model to improve fit and better facilitate projections under future climate scenarios (see Table 1 for a comparison of NorWeST and BC stream temperature model covariates).
All statistical analyses were conducted in R 4.1.2(R Core Team 2021).All spatial analyses, data manipulation, and mapping were done in ArcGIS Pro v 2.7 (ESRI, Redlands, California).

Spatial data framework
The BC Freshwater Atlas (FWA) was used as the spatial framework for this project.The FWA is a geospatial database describing BC's freshwater features (e.g.streams, lakes, watersheds) at a 1:20,000 scale that serves as a standardized dataset for mapping and analyses of freshwater systems within the province (Province of BC 2010).The FWA uses a hierarchical coding scheme that describes the relative position of features within the hydrologic network; this facilitates catchment delineation and analyses of network connectivity.The 'fundamental watershed' is the base sampling unit of the FWA.Fundamental watersheds are the smallest areal unit in the FWA and represent the surface area that drains directly into a stream reach.Within the FWA, there are > 3.2 million fundamental watersheds with a median area of 13 ha.The median length of a stream reach associated with a fundamental watershed is approximately 400 m.Hereafter, the term 'watershed' is used in reference to fundamental watersheds; a 'watershed' is the areal unit that directly corresponds to a linear stream 'reach'.The term 'catchment' is used to refer the entire upstream area that drains into a specific point along a reach.Within the spatial framework of the FWA, a catchment consists of all watersheds upstream of a specific point, inclusive of the watershed that contains the point.Since the FWA does not resolve to a finer scale than a watershed, this means that all catchments are defined relative to the most downstream point of a particular reach/watershed.
Most spatial data layers in the FWA are limited to the BC provincial boundary, so the watersheds dataset does not include areas that drain into BC from neighbouring provinces, territories, or states.To address this deficiency, additional watersheds were delineated to capture all areas outside of BC that drained into the province's stream network and create a revised fundamental watersheds dataset that included the entire catchment area for all of BC (Figure 1).Watersheds that intersected the BC boundary were identified, then their respective catchment areas were delineated using the Hydrology toolset in ArcGIS Pro and a digital elevation model (DEM).The DEM was compiled using data from the Canadian Digital Elevation Model (CDEM) (2015; 3 = 4 arc second Table 1.Covariates from regional-scale stream temperature models developed for BC and Western US (i.e.NorWeST) that were considered as templates for this modelling effort.Listed covariates are representative of an accepted mechanism of influence on stream temperature, but data type, naming conventions, and calculation methods (e.g.data transformations, spatial and temporal resolution) vary between models.Modelled stream temperature metrics were Maximum Weekly Average Temperature (MWAT), August median temperature, and August mean temperature.

Covariate
BC (Moore 2006) BC (Moore, Nelitz, and Parkinson 2013) NorWeST (Isaak et al. 2017) resolution) for all areas in Canada and data from the National Elevation Dataset (NED) (2019; 1 = 3 arc second resolution) for all areas in the US.Both elevation datasets were resampled to a 25-m resolution using a shared coordinate system (BC Albers: https:// epsg.io/3005) then mosaicked together.The FWA's hierarchical coding scheme was used to assign each of the new watersheds a position within the hydrologic network and to maintain the functionality of the original FWA database.Where necessary, the hierarchical codes of existing watersheds within BC were revised to accommodate the new watersheds and maintain the correct drainage patterns.
A notable development in the modelling process was the shift from a station-specific to a reach-or watershed-specific approach to better align with the FWA spatial framework.All data used in this project, including model covariates and in situ stream temperature observations, were initially summarized by watershed (Figure S1).Catchment-scale covariates (e.g.fractional lake cover) were subsequently calculated for each watershed by delineating the associated catchment, then aggregating the contained watershedscale values (Figure S1).For the BC MWAT model development, observed thermal metrics and temporal covariates were averaged by station when more than one year of data were available.This maintained the assumption of independence for each observation while using all available information.That approach was adapted here to watersheds, averaging the AugTw observations and all temporal covariates for each station-year by watershed instead of station (Figure S1).For example, two in situ stations along the same reach (i.e.within the same watershed) with a combined total of three station-years of AugTw observations and associated covariates would be averaged to produce a single AugTw value and set of covariate values for that watershed.

Stream temperature data
Stations included in the BC MWAT model were limited to catchments with areas between 1 -10,000 km 2 , but the upper limit on catchment size was removed for this study to capture a greater portion of the BC stream network.This allowed for inclusion of some stations that had been excluded from Moore, Nelitz, and Parkinson's (2013) original project dataset, as well as more recent observations (� 2010 -2020) available from Water Survey of Canada hydrometric monitoring stations within BC that record water temperature data.All additional data were quality controlled according to the same protocols used for the BC MWAT model (see Appendix A in Lewis et al. 2000).The updated in situ stream temperature dataset used to refit the model included observations from 562 unique stations and 1,544 station-years.Observations spanned 1981-2020, with an average of 2.7 years of observations per station.Approximately 50% of stations had only one year of data, and 74 stations had observations from more than five years.Of stations with at least two years of data, the range of AugTw observations averaged 1.9 � C (SD ¼ 1.4 � C).Stations drained areas between 1 and 217,000 km 2 , with the majority of stations (510) draining areas less than 10,000 km 2 .The type and extent of human impacts in catchments varied (e.g.forestry, flow regulation, agriculture).The in situ observations from all 562 stations were aggregated into 534 unique watersheds within the FWA spatial framework.

Covariate data
The BC MWAT model provided an initial suite of candidate covariates to consider during the model refitting process.Given the high correlation between MWAT and other metrics of summer thermal maxima (Sullivan et al. 2000;Isaak and Hubert 2001;Dunham et al. 2005;Nelitz, MacIsaac, and Peterman 2007), it was expected that the suite of covariates in the BC MWAT model would have been sufficient if the model were simply refit to predict AugTw as an alternate thermal metric.However, given the intent to use the BC AugTw model to model stream temperatures under future climate scenarios, efforts were made to include temporal covariates that could capture interannual changes in key processes.From the BC MWAT model covariates, suitable datasets were available to treat summer air temperature and glacier cover as temporal covariates.Hydroclimatic conditions, as represented by k2, were also expected to change through time, but data for future scenarios were not available.Instead, annual precipitation was substituted for k2 due to readily available datasets of historical and future conditions.Of the potential covariates from other regionally similar stream temperature models (Table 1), most lacked available datasets with sufficient spatial and temporal coverage for BC (e.g.groundwater, discharge, flow regulation).Latitude was included for consideration to represent any geographic effect of solar radiation that was not captured by air temperature and precipitation.The full suite of considered model covariates and the data processing steps are illustrated in Figure S1.
The FWA fundamental watersheds dataset was used to delineate catchment boundaries and derive catchment area for each watershed.Latitude was taken from the centroid of each watershed.Mean catchment elevation was calculated from the 25-m study area DEM.Slope was calculated as the gradient (rise:run) of the stream reach flowing through a watershed; reach length was calculated from the FWA's linear stream network layer and total elevation change was derived from the study area DEM.The FWA contains a provincial lakes dataset that was used to calculate total lake area within each watershed in BC.Lake areas in transboundary drainages within Canada were derived from the National Hydro Network (NHN 2016), and from the National Hydrography Dataset (United States Geological Survey (USGS) 2017) in the US.
ClimateNA v7.30 (Wang et al. 2016) was used to produce monthly gridded climate layers, downscaled to a 1-km resolution, for both August mean air temperature ( � C) and total annual precipitation (mm).Each climate variable was represented by both a baseline value representing the broad-scale geographic pattern and an anomaly value representing the temporal variability.Baseline values were calculated as the 40year average of mean August air temperature or annual precipitation from the 1981-2020 period.The air temperature anomaly was calculated as the difference between the mean August air temperature in a given year and the baseline value.The precipitation anomaly was calculated as the ratio of annual precipitation in a given year to the baseline value.The Randolph Glacier Inventory Version 6.0 (RGI) was used to derive glacier area for the baseline period (RGI Consortium 2017a).The RGI was developed as a global inventory of glacier footprints, approximating conditions at the start of the twenty first century (RGI Consortium 2017b).The RGI footprints were considered to be representative of the entire baseline period (i.e. 1981-2020).The majority of RGI glacier outlines in BC were dated from 2004-2010, but there was no glacier inventory of comparable coverage and resolution that aligned with the earlier half of the baseline period.Glacier loss in western Canada has been documented during the late twentieth century (Menounos et al. 2019;Hugonnet et al. 2021;Bevington and Menounos 2022), so use of the RGI for the entire historical period may have resulted in some underestimation of baseline glacier area.
Datasets for temporal covariates under future conditions (i.e.2021-2100) were driven by six CMIP5-era global climate models (GCMs) forced under three representative concentration pathways (RCPs): RCP 2.6 ('best case'), RCP 4.5 ('middle of the road'), RCP 8.5 ('worst case') (Table S2).Projected air temperature and precipitation datasets are widely available, so the selection of future climate scenarios was dictated by those available in the projected glacier dataset (Clarke et al. 2015).Glacier projections from the Clarke et al. (2015) regional glaciation model (data accessed from https://couplet.unbc.ca/data/RGM_archive/) were used to quantify glacier cover for future climate scenarios.Clarke et al. (2015) developed a high-resolution (200-m grid) model to project changes in glacier volume and area in western Canada under a range of future climate scenarios.Change-factor values (Ho et al. 2012) derived from the projected glacier layers were applied to the baseline glacier values to produce annual estimates of glacier cover for each future scenario (see Supplemental Materials for full details of glacier layer processing).Air temperature and precipitation layers for future years and scenarios were produced using the same downscaling process used for the baseline climate data (ClimateNA; Wang et al. 2016).Projected air temperature and precipitation data were bias-corrected (Ho et al. 2012) to a 2001-2020 reference period to better align with the baseline climate data; the GCM forecasts started in 2005/2006 while the future period for this project started in 2021.Annual anomaly values for air temperature and precipitation under future scenarios were calculated relative to the 1981-2020 baseline value.
Data for temporal covariates (i.e.air temperature anomaly, precipitation anomaly, and fractional glacier cover) were averaged over 20-year increments for the production of thermalscape scenarios (Table S2).Twenty-year periods were used, instead of 30-year normal periods, to provide finer temporal resolution for the thermalscape products that would be more useful in a conservation and resource management context.Thermalscape scenarios were split into two general time periods: a 'historical period' (1981-2020) and a 'future period ' (2021-2100).

Model fitting and evaluation
The full suite of covariates considered for model refitting included all of the covariates from the BC MWAT model except k2, as well as precipitation and latitude.All covariates had well-established relationships with physical processes known to affect stream temperature and have been widely used in other modelling efforts (e.g.Table 1).Multi-collinearity was assessed using variance inflation factors and was not assessed to be a problem: values for all candidate covariates were less than 3.5.
Multiple linear regression was used to fit a model to each unique subset of candidate covariates (MuMIn package in R, Barto� n 2022).The best model was determined using the lowest ranked AIC score.Of the models with DAIC score < 2, the best model was determined by the lowest ranked BIC score that penalises more complex models.A 10-fold cross-validation was used to assess the predictive accuracy of the bestfit model (Moore, Nelitz, and Parkinson 2013;Jackson et al. 2018).Data were randomly split into 10 groups of approximately equal size.Each group was used as a test dataset for the model fit to the other nine groups.Observations and predictions from each test group were merged, and a suite of performance metrics were calculated (i.e.R 2 , RMSE, and MAE).
The degree of extrapolation for modelled stream temperatures was characterized using the dsmextra package in R (Bouchet et al. 2020).The Extrapolation Detection tool was used to quantify the total length of stream network where conditions were analogous (i.e.no extrapolation) to the reference dataset (i.e. the data used to build the BC AugTw model) and identify which covariate most contributed to extrapolation in the non-analogous sections of stream network.The tool identified both Type 1 extrapolation (values outside the range of individual covariates found in the reference dataset) and Type 2 extrapolation (unique combinations of values within the range of individual covariates found in the reference dataset) (Mesgaran, Cousens, and Webber 2014).Extrapolation was also characterized using the Percent Nearby tool that estimates how 'well-informed' a prediction is based on the percent of reference data within a neighbourhood defined by one geometric mean Gower's distance (King and Zeng 2007).Extrapolation assessments were completed for the 2001-2020 historical period, and for each of the three end-of-century (i.e.2081-2100) RCP ensembles to capture a range of environmental conditions.
To evaluate the compatibility between the NorWeST and BC AugTw models, standard performance metrics for each model were assessed and predicted stream temperatures from the respective models were compared along shared stream reaches.All cross-border (i.e.BC-US) stream reaches with associated temperature predictions from each model were identified.Since the FWA and NorWeST stream networks were derived from different data sources (FWA and NHDPlus, respectively) and at different spatial resolutions (FWA ¼ 1:20,000; NHDPlus ¼ 1:100,000), they did not align seamlessly.Therefore, only cross-border reaches where segments on each side of the border could be confidently assigned to the same reach were included; decisions were supported by analysis of basemap imagery and the study area DEM.A total of 76 cross-border streams reaches were identified, two of which were excluded because they were approximately 1-2 km downstream of major dams.Unlike the NorWeST models, there was no covariate in the BC AugTw model that captured the effect of flow regulation structures, so direct comparison of predicted temperatures for those reaches was not appropriate.Predictions from the NorWeST baseline period (1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011) were compared to BC AugTw predictions for the entire historical period (average of 1981-2000 and 2001-2020 periods) and the matching 1993-2011 time period.Both historical comparisons were included to assess whether the production of a separate 1993-2011 thermalscape scenario for all of BC was necessary for integration with NorWeST or if the 20-year historical thermalscapes were sufficient.Comparisons were not made between future temperature projections from the respective models because the time periods, GCMs, and climate forcing scenarios did not match.

Thermalscapes
The BC AugTw model was used to predict August mean stream temperature in each watershed.Predictions were made for all six 20-year periods spanning 1981-2100; this included the two historical periods and each unique future scenario between 2021 and 2100 for a total of 86 unique temperature predictions (Table S2).For all future periods, stream temperature predictions corresponding to each of the six GCMs were used to calculate an ensemble mean for each 20-year period and RCP scenario.All temperature predictions were transferred from the watersheds to the corresponding reach in the FWA's linear stream network that served as the template for the thermalscapes.Modelled stream temperatures below 0 � C were possible, especially in catchments with high glacier cover.Sub-zero predictions were rounded up to 0 � C, but retained in the thermalscape so a baseline comparison was available if stream segments warmed above 0 � C under future scenarios.The final thermalscapes were produced by filtering the network to exclude secondary flow paths (e.g.braided channels, isolated wetlands, connecting lines), lakes, and all reaches with catchment areas less than 1 km 2 (the minimum catchment size of the in situ stations).All analyses were performed on this final version of the thermalscapes.
To further characterize the thermal conditions, we calculated the interannual stream thermal sensitivity (e.g.Isaak et al. 2016b;Luce et al. 2014) for each stream reach from the modelled stream temperature response.Thermal sensitivity describes the expected change in water temperature for a given change in air temperature (Mohseni, Stefan, and Erickson 1998;Kelleher et al. 2012), and can highlight regions more vulnerable to temperature shifts.For each stream reach, AugTw predictions from all scenarios were paired with the corresponding air temperature covariate (74 pairs total -ensemble means were not considered; Table S2) and a simple linear regression was fit for each reach.The coefficient of the slope parameter was taken as the sensitivity value (e.g.Leach and Moore 2019;Isaak et al. 2016b;Luce et al. 2014;Mayer 2012).

Model re-fitting and evaluation
The best model from the re-fitting process is summarized in Table 2 (see Table S3 for comparison of top-ranked models).The directionality of the coefficients was consistent with the assumed effects represented by each covariate.The scatterplot of observed versus predicted values showed good agreement through the majority of observed temperature range (adjusted R 2 ¼ 0.79) (Figure 2).Residuals from the cross-validation appeared normally distributed with a RMSE of 1.53 � C and a MAE of 1.18 � C (Figure 3).There was no apparent spatial pattern in the distribution of the prediction errors throughout BC (Figure 4) and no obvious issues with model fit based on scatterplots of model covariates and prediction errors (not shown).
The range of environmental conditions represented in the reference dataset captured much of the variation across the considered scenarios with extrapolation limited to 4.2 − 23.2% of the total network length (Table 3).High elevation areas in the Columbia and Fraser drainage regions were a persistent source of extrapolation under all considered scenarios (mean catchment elevation >2,100 m; Figures S6-S9).Extrapolation attributed to low mean August air temperatures (<8.8 � C) in northern BC during the 2001-2020 period was largely absent under the   The cross-border comparison of the BC AugTw predictions for the full historical period to NorWeST yielded good agreement (Figure 5).Predictions from the two models were highly correlated (Pearson's r ¼ 0.91), and the slope coefficient from a simple linear regression fit between the two sets of temperature predictions was 1.03 (Figure 5).The BC AugTw predictions averaged 1.1 � C warmer than the predictions from NorWeST.The comparison between NorWeST and the 1993-2011 period for BC (not shown) produced a similar result with a Pearson's r of 0.90 and a slope coefficient of 1.03.Therefore, the thermalscape for the full historical period was considered sufficiently aligned with NorWeST's baseline period to combine the two products, and production of a separate 1993-2011 thermalscape for BC was not necessary.

Thermalscapes
The final thermalscape product consisted of 681,142 km of linear stream network across the entire province (Figure 6).For the 2001-2020 historical period, AugTw averaged 8.6 � C (SD ¼ 3.3 � C), with a maximum value of 21.0 � C and a minimum value of 0 � C. A total of 8,941 km of stream network (�1% of total) were assigned a 0 � C temperature value during the 2001-2020 period; that value dropped to <0.1% of network length under end-of-century climate scenarios, with 0 � C stream reaches only persisting in the most high-elevation, glacier-dominated catchments.The spatial distribution of modelled stream temperatures across BC followed expected trends: warmer temperatures at low elevations and in larger streams and river systems, and cooler temperatures in lowerorder, high-elevation streams (Figure 6).Due to the high resolution of the thermalscape stream network, smaller catchments made up a relatively high proportion of the total network length.For example, third order streams and smaller (as defined using the FWA's stream order attributes) accounted for 76% of the total thermalscape length.Therefore, high-level summaries of thermalscape characteristics were largely driven by the characteristics of those smaller streams.
Between the two historical periods, average reachlevel AugTw in BC increased by 0.3 � C (SD ¼ 0.3 � C).By end-of-century (2081-2100), province-wide increases in stream temperature were projected under all three RCP scenarios (reported values are ensemble Table 3. Percentage of BC stream network where environmental covariates were not extrapolated beyond environmental conditions at the in situ stream temperature stations.Historical and end-of-century scenarios were used to illustrate the degree of extrapolation under a range of environmental conditions.ranging from 0.7 − 10.5 � C. Modelled stream temperature increases under RCP 2.6 averaged 1.9 � C by endof-century, with temperatures generally stabilizing by mid-century.As expected, modelled temperature changes under the RCP 8.5 scenario were most extreme, with warming rates approaching 0.8 � C/decade and average AugTw increasing 6.1 � C by the 2081-2100 period. A finer-scale analysis of the thermalscapes illustrated the effects of the key model covariates and the spatial heterogeneity in AugTw.The thermalscape for the 2001-2020 period in the Fraser Valley region of southwestern BC (Figure 7A) showed the warmest temperatures in the mainstem of the Fraser River and streams in the low-lying valley floor.In the Coast Mountains to the north, the effects of elevation and glacier cover were evident in the cooler temperature predictions in streams draining to the coast or the Fraser mainstem.The variation in stream temperature change was evident in the DAugTw between 2001 and 2020 and the end-of-century RCP 4.5 scenario (Figure 7B).The greatest increases in modelled temperatures (>4.5 � C) occurred in the headwaters of glacier-dominated streams, whereas other high-elevation streams in the Coast Mountains that were not glacier-fed exhibited the least amount of change (<2.5 � C).Modelled temperature changes in the Fraser River and other low-elevation catchments were in an intermediate range (2.5-3.5 � C).
Thermal sensitivity values ranged from 0 to 1.9 with a median value of 0.75; nearly 90% of the network had sensitivity values between 0.7 and 0.8.The spatial distribution of thermal sensitivity values across the province was consistent with the relative influences of the temporal covariates in the BC AugTw model (Figure 8).Specifically, province-wide increases in August mean air temperature were correlated with warming stream temperatures.Projected increases in annual precipitation resulted in lower thermal sensitivity values across much of the province, particularly along the coast where the greatest increases in precipitation were predicted.The geographically broad, but modest, effect of precipitation on thermal sensitivity was evident from the median sensitivity value being slightly lower than the air temperature anomaly coefficient (0.75 and 0.77, respectively).Higher-sensitivity streams were associated with glacier-influenced catchments, where both projected increases in air temperature and declines in glacier cover contributed to higher warming rates.The influence of thermal sensitivity was evident at a finer-scale in the modelled DAugTw in the Fraser Valley region where the most

Discussion
The BC AugTw stream temperature model provides thermalscapes for the entire province under a range of past and future scenarios.This is the first dataset of continuous stream temperature predictions of this extent and resolution in Canada, bridging a gap between site-or catchment-specific monitoring and global-scale water temperature products (e.g.Bosmans et al. 2022).This work addresses a major environmental data deficit in BC's freshwater ecosystems and can benefit myriad research and management applications.As exemplified by the use of thermalscape products developed elsewhere, especially NorWeST (Isaak et al. 2017), the thermalscapes that we developed can contribute to: investigations of thermal heterogeneity across stream networks and projected changes under future climate scenarios (e.g.Ficklin et al. 2014;Isaak et al. 2017;Jackson et al. 2018); projected changes in species distributions (e.g.Al-Chokhachy et al. 2016;Bell et al. 2021) or community assemblages under future climate scenarios (Walters, Mandeville, and Rahel 2018;Kirk and Rahel 2022); and identification and prioritization of habitat for conservation or restoration.The refit of an existing model to allow more direct comparisons of predictions with an existing, neighbouring dataset represents a somewhat novel approach to mosaicking stream temperature datasets to leverage existing products and expand spatial coverage.In this case, the mosaicked thermalscape produced from NorWeST and the BC AugTw thermalscapes covers much of western North America.
The main benchmark for the BC AugTw model's performance was the BC MWAT model that was used as a template.The BC AugTw model showed some improvement over the original with a lower RMSE (1.5 � C vs. 2.1 � C), which was consistent with the modifications to the model.While removing the upper catchment size limit may have generalized some scale-dependent processes (Jones and Schmidt 2017), improved model performance was primarily attributed to the longer averaging period of AugTw relative to MWAT (Stefan and Preud'homme 1993).Similar to the BC MWAT model, the precision of the BC AugTw model was in line with other regionalscale statistical stream temperature models (reviewed in Gallice et al. 2015).All covariates included in the model were consistent with expected physical processes.The precipitation covariates had a cooling effect, as was found by Moore (2006) and Isaak et al. (2017), and were an effective substitute for the k2 variable included in the BC MWAT model.Latitude, a covariate not in the original model, had a cooling effect consistent with the geographic variation in solar radiation.Slope was not included in the final model (Table S3); including slope did not improve model performance and it had a warming effect in alternative models, contrary to its expected cooling effect (Webb et al. 2008).The lack of importance of slope in the BC AugTw model may be attributed to the addition of sites in larger, typically lower-gradient catchments, and the shift to a coarser measurement scale for stream gradient (i.e.reach vs. station).Agreement between the BC AugTw and NorWeST stream temperature predictions was good, especially given the differences in modelling approach, data sources, considered covariates (Table 1), and data processing methods.For example, the BC AugTw model used air temperature values averaged across watersheds, whereas NorWeST used air temperature values averaged across 'processing units' that can cover tens of thousands of km 2 .
The use of the FWA as the spatial framework for this project was a notable departure from the original methods and was a potential source of uncertainty in the model.Whereas covariate values used in the BC MWAT model were calculated for specific stations (i.e.points), data in this project were only resolved to the scale of watersheds.The use of the FWA as a spatial framework for organizing data was similar to Isaak et al.'s (2017)   The effect of this uncertainty is expected to be minor given that the typically small area of the FWA's watersheds (median � 13 ha) and associated reaches (median � 400 m).Each watershed typically constitutes a small fraction of the total catchment area in all but the most headwater streams, so the variation in covariate values within a watershed is usually small.For example, in >90% of the watersheds with in situ temperature data, catchment-scale covariate values differ by less than 5% when calculated at the most upstream point in the watershed instead of the most downstream point.Of the in situ stations in more headwater systems, where position within the watershed would have a greater effect on covariate measures, most were towards the downstream end of watershed and should be reasonably well represented by the watershed-scale covariate data.
A main driver of this work was a need for stream temperature projections under future climate scenarios.This project delivers on Moore, Nelitz, and Parkinson's (2013) recommendation to use the BC MWAT model for first-order estimates of the effects of climate change on stream temperature, albeit with the refit model.The applicability of stream temperatures predicted with statistical models to future periods, as done here, has been questioned due to potential non-stationarity in the statistical relationships (Arismendi et al. 2014;Leach and Moore 2019).S4).
Similar to Isaak et al.'s (2017) defence of the NorWeST scenarios: (1) the BC AugTw model better captures the physical processes involved in stream heat budgets than simple air temperature to water temperature models, and (2) the range in values from the static and temporal covariates in the model dataset covered much of the projected range in future conditions for the majority of the province.This position was supported by the extrapolation assessment (Table 3, Figures S5-S12), which demonstrated an approach to detect and characterize the uncertainty introduced by extrapolation (Isaak and Luce 2023).Further, assessment of climate change impacts will continue to be a priority in both freshwater research and management contexts, necessitating the creation of this type of product.The BC thermalscapes represent an improvement over commonly used alternatives to existing water temperature observations or predictions: direct use of air temperature as a proxy for water temperature (e.g.Rieman et al. 2007;Wenger et al. 2011;2013) or simple air temperaturewater temperature relationships (e.g.Mohseni, Stefan, and Erickson 1998;Mantua, Tohver, and Hamlet 2010;Al-Chokhachy et al. 2013).These approaches often lead to over-prediction of warming, especially in mountainous regions (Kirk and Rahel 2022).
The projected changes in BC stream temperatures under future climate scenarios were consistent with the forecasted changes in the model's temporal covariates, broadly including increasing air temperature, increasing coastal precipitation, and declining glacier cover.The BC AugTw model predicted overall warming rates during the historical periods of approximately 0.14 � C/decade , similar to NorWeST's estimate of 0.17 � C/decade from a similar 40-year period , or Ferrari, Miller, and Russell's (2007) estimate of 0.12 � C/decade in the lower Fraser River, BC.Thermal sensitivities derived from the BC AugTw model were in good agreement with simulated sensitivity values in the Fraser Basin (0.43-1.01;Islam et al. 2019), and western Canada (0.25-0.75;Shrestha and Pesklevits 2023).Ficklin et al.'s (2014) projected warming in the Columbia basin of 4.1-5.0� C for an end-of-century RCP 8.5 scenario was in line with the BC AugTw model, but projections by Morrison, Quick, and Foreman (2002) and Ferrari, Miller, and Russell (2007) of approximately 2 � C warming in the lower Fraser River were markedly lower than those of the BC AugTw model for end-of-century warming of 2.8 � C or 6.5 � C (RCP 4.5 and RCP 8.5, respectively).NorWeST projected changes in stream temperatures of approximately 2 � C by end-of-century under an A1B emission scenario for the processing units bordering southern BC (Isaak et al. 2017), compared to the BC AugTw projection of roughly 3 � C warming under the RCP 4.5 scenarios (A1B roughly corresponds to RCP 6.0).Differences between the projected warming could be attributed to climate forcings from different model eras (CMIP5 vs. CMIP3) and a different set of temporal covariates.For instance, glacier cover and annual precipitation were treated as static variables in NorWeST models, whereas discharge was included as a temporal covariate that could not be considered in the BC AugTw model.
Glacier cover is a common covariate considered in stream temperature models (e.g.Moore, Nelitz, and Parkinson 2013;Isaak et al. 2017), although it has been typically treated as a static variable (Isaak et al. 2017) or climate change projections were not explicitly considered (Moore 2006;Moore, Nelitz, and Parkinson 2013).Declining glacier cover in western North America is expected to result in warmer summer stream temperatures due to reduced contribution of glacier meltwater to total streamflow (Moore et al. 2009), which is consistent with the BC AugTw model's projected warming from reduced glacier cover.Warming-induced increases in glacial discharge, and correspondingly cooler stream temperatures, are expected during the initial phase of a climatic warming period (Moore et al. 2009); however, declining streamflow attributed to glacier retreat has been documented through much of the study region (Stahl and Moore 2006;Moore et al. 2020), suggesting that most of BC has passed this initial warming period.Therefore, the BC AugTw model should provide reasonable first-order estimates of the effects of glacier retreat on stream temperature for most of the province.
The precipitation covariates in the BC AugTw model had a cooling effect as expected, given that more precipitation would be associated with higher discharge and more groundwater (Isaak and Hubert 2001).Precipitation covariates were based on watershed-specific estimates of total annual precipitation (similar to the station-specific approach by Moore 2006); however, potential changes in the type of precipitation (i.e.rain or snow) and its intra-annual distribution should be considered, particularly when interpreting future stream temperature projections.The consensus among considered GCMs was for increased precipitation across the study area, especially along the coast.The more recent generation of GCMs from CMIP6 similarly project increased precipitation, as well as higher likelihood of high intensity precipitation events (IPCC 2021).Islam, D� ery, and Werner (2017) projected increased precipitation for the Fraser basin under RCP 4.5 and RCP 8.5 scenarios, but with the fraction of annual precipitation as snow projected to decline nearly 50% by mid-century.The combination of intense precipitation events and less snowpack is expected to result in lower summer discharge in streams, which would be expected to offset some of cooling projected by the BC AugTw model based on increased annual precipitation.
There are limitations inherent to statistical models of stream temperatures, and required precision is an important consideration depending on the application of the thermalscapes.Regional-scale status and trends appear represented, but more fine-grained applications may warrant further scrutiny.In a more localized setting, validation of stream temperature predictions against external datasets may be a useful approach to quantify, and potentially correct for, model bias.Lakes were removed from the thermalscapes, but stream temperature predictions immediately downstream should be interpreted with caution due to the potential influence of the lake outflow temperature.Similarly, flow regulation structures were not explicitly considered, and the outflow temperature from reservoirs is dependent on the heights of the dam and the outlet structure (Hamblin and McAdam 2003;Olden and Naiman 2010;West and Moore 2020).Stream predictions in the headwaters of glaciated streams also warrant caution where fractional glacier cover exceeds the range of the model development data.The most extreme of these were reaches where sub-freezing predictions of AugTw were truncated to 0 � C. Data limitations highlighted by Moore, Nelitz, and Parkinson (2013) in the BC MWAT model, notably the effects of riparian shade and groundwater influences on stream temperatures, persist and were not addressed in the BC AugTw model.Lastly, landscape activities such as logging can increase stream temperatures (Moore et al. 2005).Anthropogenic land cover or land use changes were not considered due to the lack of a consistent spatial dataset for the entire study area with sufficient spatial and temporal resolution to align with the in situ observations.Future improvements to the BC AugTw model are possible (but not currently planned) should relevant datasets become available, including spatial layers for covariates that could not be considered here (e.g.discharge, riparian shade) and stream temperature observations from underrepresented regions of BC (see Figure 4).

Conclusions
A regression model was developed to predict August mean stream temperatures for British Columbia, Canada from land cover, physiographic, and climatic characteristics.A 10-fold cross-validation suggested the model's performance is in line with similar statistical stream temperature models (R 2 ¼ 0.79, RMSE ¼ 1.53 � C, MAE ¼ 1.18 � C).The model was used to create thermalscapes, continuous predictions of stream temperatures, for the provincial stream network.Thermalscapes were produced for a total of 86 scenarios, spanning the period of 1981-2100 and a range of future climate scenarios.Modelled stream temperature warming averaged 2.9 � C by end-of-century under an intermediate climate forcing scenario, with significant spatial heterogeneity in warming rates across the province.Predicted stream temperatures for BC under baseline conditions were well-correlated with stream temperature predictions from the NorWeST thermalscape product that spans much of the western United States (Pearson's r ¼ 0.91).The thermalscapes provided through this modelling effort address a significant data gap for freshwater systems within BC.More broadly, the potential applications extend beyond BC through the integration of BC and NorWeST thermalscapes, providing continuous, highresolution stream temperature predictions across much of western North America, and marking a major step forward for applications in aquatic research and management.

Figure 1 .
Figure 1.Major drainage regions within British Columbia, Canada (inset: dark gray).Transboundary drainage area was included to capture the entire catchment area for all stream reaches within the province.Source: Data was adapted from the British Columbia Freshwater Atlas under the Open Government Licence -British Columbia.

Figure 3 .
Figure 3. Histogram of prediction errors from the cross-validation.Standard deviation of prediction errors is 1.53 � C.

Figure 2 .
Figure 2. Observed vs. predicted August mean stream temperatures from the cross-validation with a 1:1 reference line (dashed line).

Figure 4 .
Figure 4. Mapped prediction errors from BC AugTw model, binned by standard deviation (SD ¼ 1.53 � C).Spatial offset applied to error points for display purposes.
means).Under RCP 4.5, average AugTw rose by 2.9 � C (SD ¼ 1.0 � C) from 2001-2020, increasing at approximately 0.4 � C/decade through mid-century before stabilizing in the 2061-2080 period.While the average increases in AugTw were generally representative of the degree of warming at the provincial scale, local influences on stream temperature had great effects: end-of-century warming was projected in all stream reaches with DAugTw (i.e.change in AugTw)

Figure 5 .
Figure 5. NorWeST vs. BC predictions of August mean stream temperature for transboundary stream segments (n ¼ 74) along the BC-US border.NorWeST predictions were for a 1993-2011 reference period and BC predictions were averaged from the 1981-2000 and 2001-2020 predictions.Solid line is a simple linear regression between NorWeST and BC temperature predictions.Dashed line is a 1:1 reference.

Figure 6 .
Figure 6.Predicted August mean stream temperatures for all of BC for the 2001-2020 period (for display purposes, temperature predictions are mapped to watersheds instead of stream reaches).
use of the NHDPlus dataset in developing the NorWeST thermalscapes.With respect to the catchment-scale covariates in the FWA framework, this means that the values for a given watershed

Figure 8 .
Figure 8. Modelled thermal sensitivity of BC stream reaches.

Figure 9 .
Figure 9. Predicted change in August mean stream temperature at individual sites (inset map bottom right) with different thermal sensitivities (High, Intermediate, Low).Change in stream temperatures were calculated relative to the 2001-2020 baseline period.Lines are the mean of a 6-member GCM ensemble for each RCP scenario and shaded error bars are ± 1 SD.River name, thermal sensitivity (sens.), and baseline temperature (temp.)included in each panel.Sites correspond with Water Survey of Canada gauge locations for hydrometric monitoring (TableS4).

Table 2 .
Details of the best-fit model.