Longitudinal, Lateral, Vertical, and Temporal Thermal Heterogeneity in a Large Impounded River: Implications for Cold-Water Refuges

: Dam operations can a ﬀ ect mixing of the water column, thereby inﬂuencing thermal heterogeneity spatially and temporally. This occurs by restricting or eliminating connectivity in longitudinal, lateral, vertical, and temporal dimensions. We examined thermal heterogeneity across space and time and identiﬁed potential cold-water refuges for salmonids in a large impounded river in inland northwestern USA. To describe these patterns, we used thermal infrared (TIR) imagery, in situ thermographs, and high-resolution, 3-D hydraulic mapping. We explained the median water temperature and probability of occurrence of cool-water areas using generalized additive models (GAMs) at reach and subcatchment scales, and we evaluated potential cold-water refuge occurrence in relation to these patterns. We demonstrated that (1) lateral contributions from tributaries dominated thermal heterogeneity, (2) thermal variability at conﬂuences was approximately an order of magnitude greater than of the main stem, (3) potential cold-water refuges were mostly found at conﬂuences, and (4) the probability of occurrence of cool areas and median water temperature were associated with channel geomorphology and distance from dam. These ﬁndings highlight the importance of using multiple approaches to describe thermal heterogeneity in large, impounded rivers and the need to incorporate these types of rivers in the understanding of thermal riverscapes because of their limited representation in the literature.


Introduction
In rivers and streams where humans have altered the thermal regime, not only does summer water temperature frequently exceed the thermal tolerance of native salmonids, thus influencing their distribution [1], but these human modifications can amplify or dampen spatial and temporal variation [2]. Specifically, dams and diversions can greatly modify riverine thermal regimes by restricting or eliminating connectivity in longitudinal, lateral, vertical, and temporal dimensions [3,4].
(1) Compare longitudinal patterns of water temperature derived from thermal infrared imagery to instream measurements in the water column. (2) Assess spatial (lateral, longitudinal, and vertical) and temporal thermal heterogeneity and potential cold-water refuges for salmonids at multiple scales. (3) Evaluate statistical associations between hydrologic and geomorphic variables and occurrence of cool-water areas and median water temperature.
We hypothesized that (1) water temperature in the main stem of the Pend Oreille River increases as water moves downstream, but lateral contributions from tributaries create discontinuities in the longitudinal profile; (2) the confluence area formed by the main stem Pend Oreille River and cold tributaries provide potential cold-water refuges for cold-water fishes; and (3) in the summer, thermal stratification from various sources (i.e., seeps, groundwater inputs, low water velocities, and submerged aquatic vegetation) creates vertical temperature gradients that isolate cool water near the river bottom to potentially provide refuge for fish.

Study Area
The study area encompassed the Box Canyon reservoir section of the Pend Oreille River, an 88.5-km portion of the river between the town of Ione, Washington, and Albeni Falls Dam in Idaho (USA) (Figure 1).
The reservoir's surface area varies between 2800 and 3600 ha, depending on water elevation and flow. Flows are maintained for maximum power production below approximately 1980 m 3 s −1 . If flows exceed this level, the reservoir operates in free-flow mode without power production. The Pend Oreille River, a sixth order stream, is the fourth-largest tributary of the Columbia River and drains an area of 62,678 km 2 with a mean annual discharge of 709 m 3 s −1 for the period from 1904 to 2015. Seasonal weather patterns throughout the study area are typical of the Northern Rockies ecoregion with hot, dry summers and cold, relatively wet winters. During the Ice Age, the Pend Oreille Lobe of the Glacial Lake Missoula formed the Pend Oreille River. The glacial Lake Missoula, part of the Cordilleran Ice Sheet, extended south and covered the valley and the surrounding low hills with glacial and alluvial material deposited during the retreat of ice creating an alluvial valley. Alluvial valleys can range from confined (floodplain width/channel width of less than 2) to unconfined (floodplain width/channel width greater than 4). Because landscape metrics such as confinement and tributary confluences are strongly associated with the density and occurrence of cold-water refuges [22], describing these landscape metrics in the Pend Oreille River is important for understanding watershed-scale relationships. The Pend Oreille River has two large, moderately confined sections: (1) Above Box Canyon dam, the first section between river kilometer (rkm) 53 to rkm 93, and (2) below Albeni Falls Dam between rkm 131 to rkm 144. In contrast, the most unconfined section is between rkm 94 and rkm 121 [43]. Calispel Creek and Le Clerc Creek (C6) are the largest tributaries within the study area ( Figure 1).
The Pend Oreille River within the study area is constrained by two main roads and levees where most of the original forests have been converted to agricultural or residential land. Additionally, cottonwood galleries have declined because of limited regeneration and replacement due to flow management. In contrast, there are numerous seasonally flooded wetlands located near the Kalispel Tribe reservation. In addition, flow is managed by the U.S. Army Corps of Engineers (ACOE) for flood control, and comparisons between 1998 and 2015 aerial imagery revealed that there have not been significant changes in channel morphology. The entire Pend Oreille River is considered impaired for water temperature in the states of Idaho and Washington because water temperatures exceed the water quality standards criterion of 20 • C at numerous locations in the river. However, cold-water refuges present within the Pend Oreille have been observed to be used by salmonids during the summer when water temperatures are above their optimum (E.K. Berntsen, J.R. Maroney, and J. Connor, unpublished data). The reservoir's surface area varies between 2800 and 3600 ha, depending on water elevation and flow. Flows are maintained for maximum power production below approximately 1980 m 3 s -1 . If flows exceed this level, the reservoir operates in free-flow mode without power production. The Pend Oreille River, a sixth order stream, is the fourth-largest tributary of the Columbia River and drains an area of 62,678 km 2 with a mean annual discharge of 709 m 3 s -1 for the period from 1904 to 2015. Seasonal weather patterns throughout the study area are typical of the Northern Rockies ecoregion with hot, dry summers and cold, relatively wet winters. During the Ice Age, the Pend Oreille Lobe of Figure 1. Pend Oreille River study area in eastern Washington and northern Idaho, USA. Orange circles are main stem sites labeled as either P (peaks) or T (troughs), and hollow squares are tributary confluences (C) included in the study. Gray tributaries were sampled but not included in analysis because of missing data. Black arrow describes direction of flow from south to north.

Airborne Thermal Infrared (TIR) Imagery
To characterize longitudinal thermal heterogeneity of the Pend Oreille River, we used airborne thermal infrared imagery (TIR) collected in 2001. TIR surveys were conducted on 16 August 2001 between the hours of 14:00 and 17:00 [44]. The data were collected using a TIR imager with a wavelength range of 8-12 microns [45]. The sensor was mounted to the underside of a helicopter. The helicopter was flown longitudinally along the river channel with the sensor and camera in a vertical or near vertical position at 1372 m above ground level. Flight altitudes were determined based on the river width and floodplain characteristics. The altitude selected allowed for an image with a ground width of 483 m and a pixel size of 0.8 m. Ground-truth measurements of instream temperature at the time of the survey indicated that radiant temperatures were accurate to within 0.5 • C.
To compare spatial patterns of instream water temperatures to TIR data, we selected a day in 2017 that was similar in air temperature and discharge to the conditions when the thermal imagery was collected in 2001. Air temperatures on 22 August 2017 ranged from 11.0 • C at midnight to near 33.

Instream Measurement
We deployed instream Hobo pendant data loggers (Onset Computer Corporation, Bourne, Massachusetts, USA; ± 0.5 • C) and recorded water temperatures at 1-h intervals at 17 sites in the main stem (labeled as either a peak, P, or a trough, T) and at 22 confluences (labeled as C) from 17 July to 14 September 2017 to describe the longitudinal and lateral thermal heterogeneity of the Pend Oreille River from rkm 53 to rkm 144 (i.e., Box Canyon Reservoir) ( Figure 1 and Table 1). At each site, loggers were positioned approximately 30 cm from the bed and 30 cm from the surface of the river to characterize vertical thermal heterogeneity. Each main stem logger was placed at a major temperature deviation (i.e., a peak [n = 7] or a trough [n = 8]) in the 2001 thermal infrared longitudinal profile based on the smoothed trend from the 10-km moving average [26]. We also deployed loggers near Albeni Falls and Box Canyon dams. Characterization of the longitudinal thermal profile is important because not all rivers warm as their water moves downstream from the headwaters to the mouth, and their profile can inform how to implement long-term monitoring programs and climate-adaptation strategies [26]. Additionally, in slow-moving rivers like the Pend Oreille in the summer months, the surface temperature may not be indicative of the overall water temperature or of cold-water areas near the river bottom due to vertical stratification.

Geographical Information System
We used publicly available spatial data layers including land use [46], geology [22,47], and tributaries [17] as explanatory variables in our statistical analysis to assess associations with cool-water areas in the Pend Oreille River. The specific variables in our analyses were land use (% developed land), aquifer permeability (high versus low), degree of confinement (m m −1 ), and tributaries (area in km 2 and number).

Land Use
We grouped land use into two categories (percent developed or percent undeveloped) from the 2010 data produced by the Washington Department of Ecology (https://ecology.wa.gov/Research-Data/ Data-resources/Geographic-Information-Systems-GIS/Data, accessed on 16 March 2018). Of the 84 land use categories, seven were described as nondeveloped land, and the remaining 77 categories were classified as developed land. To estimate percent developed land we generated points every 1 km on the centerline of the Pend Oreille River from rkm 53 to the Albeni Falls dam (rkm 144). We buffered the Pend Oreille River shoreline at 1 km and generated Thiessen polygons using every other river kilometer point that then were intersected with the buffered river layer to demarcate river segments. These river segments were bisected to create 1-km river segments within which we estimated the percent contribution of developed land on each bank.

Aquifer Permeability
We used aquifer permeability data for the Columbia River basin (https://catalog.data.gov/dataset/ pacific-northwest-pnw-hydrologic-landscape-hl-polygons-and-hl-code, accessed on 19 October 2017). These data are based on similarities in lithology. They estimate hydraulic conductivity value, identify areas of low or high aquifer permeability, and describe locations where shallow subsurface groundwater mixes with deep groundwater flows and where the loss or gain of water through groundwater fluxes may be important [48].

Valley Width and Confinement
Confinement was calculated as the ratio of floodplain width to bankfull channel width. Confined streams had ratios of less than 2, moderately unconfined streams had ratios between 2 and 4, and streams with ratios greater than 4 were considered unconfined [49]. We mapped confinement using the NetMap floodplain width layer [50] and bankfull channel width (m) data obtained from Beechie and Imaki [43] Tributaries We used the NetMap stream layer to count the number of tributaries per kilometer and to determine the watershed area of each tributary to infer tributary size (km 2 ). The NetMap stream layer is derived from digital elevation models (DEM) ranging in resolution from 10 m to 1 km. This layer was produced using a consistent method to extract a synthetic routed and attributed network including drainage area and drainage area per contour length [50].

Three-Dimensional Hydraulic Mapping and Channel Morphology
To characterize the Pend Oreille River morphology, Freshwater Map, Inc. (https://freshwatermap. com) conducted 3-D bathymetric mapping of 81 km of the Pend Oreille River from the Albeni Falls Dam to the Box Canyon Dam. Eight, one-person pontoon boats deployed Teledyne, RiverPro acoustic Doppler current profilers (ADCP, Poway, California) with fully integrated Hemisphere, Vector V102GPS compass units deployed from catamaran rafts while floating on the river to measure water depth and velocities and current vectors. Mapping was conducted in two separate surveys: On 11-14 September 2017, and on 26 March 26-1 April 2018. Sampling dates in March were chosen because SAV was not present at that time. In shorelines areas, technicians stayed in deeper water to avoid SAV. SAV can also cause incorrect readings of the streambed with ADCP.

Water Depth and Velocity
To measure depth, velocity, and discharge in the September 2017 survey, the ADCP team made at least two separate data collection runs each day. The first run covered one half of the river channel to a midpoint of the day's collection objective (river left). The second half of the river channel (river right) was measured as the team returned to the day's measurement starting point. The team maintained a preset distance apart and semi-staggered positioning. River discharge at the Newport U.S. Geological Survey (USGS) gage 12395500 during the survey ranged from 150 m 3 s −1 to 385 m 3 s −1 . For the March 2018 survey, depth, velocity, and discharge data were collected on the southern 48 km of the river, and the eight-person team made four passes each day over 8 km for a total 32 passes each day for 6 days. River discharge ranged from 694 m 3 s −1 to 745 m 3 s −1 . For both sample periods, the measurement error was within 1% of the actual flow and depths being sampled.
The survey team also measured river discharge at 20 transect locations along the river, once at the beginning of each day and once at the end. Discharge measurements were replicated 5 to 8 times at each location following standard USGS protocols for discharge measurements, resulting in a mean variance in discharge measurements less than 1% for most locations and less than 3% for all, following the methods of Lorang and Tonolla [51] and Marotz and Lorang [52]. A DEM that included bathymetry tied to publicly available lidar (light detection and ranging) data was produced from the water depth data. The DRIFTER routine in River Analyzer was used to create bottom boundary and top velocity maps.

Transect Velocity and Discharge Calculations
The reservoir was divided into 50,000 cross-channel slices using the SLICER routine in River Analyzer. Each slice was 12 m apart, spanning from river left to river right, and contained data from 32 longitudinal transects composed of 10-cm velocity bin ensembles extending from the surface to the bottom minus a 10-cm blanking distance at the surface and a variable bottom blanking distance of 10 cm to 50 cm depending on variable signal attenuation factors. Velocity ensembles covered depths from 30 cm to 8 m resulting in 8000 to 10,000 velocity data points that went into each slice discharge calculation. Wetted perimeter was calculated based on bottom ping data which defined the bathymetry that gave the distance along the bottom as the wetted perimeter for each slice. Hydraulic radius for each slice was derived by summing the cross-sectional area for each transect slice and dividing by the wetted perimeter.
We also calculated the Froude number (Fr), a dimensionless value that describes the energetic state of the water column and values greater than 1 represent supercritical flow: where V is the average slice water velocity (m s −1 ), D is the hydraulic depth (cross-sectional area/top width, in meters), and g is the acceleration of gravity (9.8 m s −2 ). For our case we wanted to quantify the relative energetic state of the water column as a potential proxy for potential water column mixing.
Lastly, continuous transects along the river, each 100-m wide and perpendicular to the river centerline, were created to summarize depth and velocity (i.e., mean, standard deviation, median, 25th quantile, 75th quantile, minimum, and maximum). Each of these transect polygons were associated with Universal Transverse Mercator (UTM) coordinates of the point where the river centerline intercepted the midline of the perpendicular transect, a rectangular polygon with 100 m width.

Submerged Aquatic Vegetation
SAV is a dominant feature of this river and can influence sediment dynamics, water velocity, and thermal gradients [15]. To characterize SAV distribution, we used data mapped by The Pend Oreille County Public Utility District in 1997. We recognize that although these data may not reflect current conditions, they are representative of 2001 conditions when the TIR data were collected. For our analysis, we only used the areas classified as densely vegetated regardless of plant species composition, and we did not consider unvegetated areas. The midpoint of these polygons was then linearly referenced to the nearest river kilometer.

Spatial and Temporal Analysis
To describe spatial and temporal thermal heterogeneity of the Pend Oreille River, we used digital hydrography data (mixed-scale hydrography, version 3.1, or "MSHv3.1") at 1:24,000 scale or finer for the state of Washington from the StreamNet Project, Pacific States Marine Fisheries Commission as the map template in a geographical information system (GIS) (https://www.streamnet.org/gisdata/ map_data_base/Hydrort_MSHv31_September2012.xml, accessed on 13 October 2016). We then reduced the thermal and environmental ancillary data (e.g., flow velocity, depth) to one-dimensional positions, expressing these geographic positions as measurements along the river's path in relation to the confluence with the Columbia River [53]. We calculated minimum, average, maximum, daily minimum, and daily maximum differences between surface and bottom water temperatures, and the consecutive 7-day moving average of daily maximum (7DADM) temperature for all sites. We calculated the 7DADM because this metric incorporates the average of maximum temperatures that fish are exposed to over a week and is used in water quality standards to protect against acute effects such as mortality and migration blockage [54]. Confluences were considered stratified when differences between hourly surface and bottom water temperatures were equal or greater than 2 • C, and we estimated the percentage of time that confluences were stratified.
To determine whether main stem or tributary confluences could potentially be used as cold-water refuges by adult salmonids, we used the period from 9 August to 13 September 2017. Sites were considered potential cold-water refuges based on Greer et al. [55]: (1) Ambient refuges, defined by differences of 1 • C or 2 • C, and (2) physiological refuges, defined as the physiological threshold for migrating adult salmonids (18 • C) [54][55][56][57]. More specifically, bottom water temperatures in the main stem can potentially function as cold-water refuges when the difference between the surface water temperature and bottom water temperature was greater than 1 • C or 2 • C for at least an hour a day, or when bottom 7DADM water temperatures were lower than 18 • C. Confluences can be potential cold-water refuges when their 7DADM temperatures are 1 • C or 2 • C cooler than the 7DADM main stem temperature, or the 7DADM confluences temperature were below 18 • C. In the Pend Oreille River, acoustic-tagged salmonids were observed near shore and at shallow depths, but detection in deeper areas was difficult because of the loss of signal (J. Connor, unpublished data). In contrast, in the main stem Columbia River, tagged bull trout used deep (12.6 m on average), slow habitats [58]. We characterized the temporal variability and persistence (i.e., the number of times a potential cold-water refuge at a given location was observed) [18] for each of the thermal criteria outlined above. Finally, we estimated the spacing between all persistent physiological cold-water refuges and compared them to the spacing of physiological cold-water refuges present at least 50% of the time to determine (1) how far salmonids would have to travel before they experienced thermal relief during their migration, and (2) how often salmonids may use these ambient cold-water refuges as stepping stones because the distances between physiological cold-water refuges may exceed the swimming abilities of trout under these stressful conditions [59].

Statistical Analysis
We calculated z-scores for all the channel morphology metrics measured or estimated during 3-D mapping to normalize measures for temporal variation among the two surveys. Z-score distributions all have the same mean and standard deviation, and individual scores from different distributions can be directly compared. The complete suite of variables for which we calculated z-scores included surface width (m), cross-sectional area (m 2 ), wetted perimeter (m), hydraulic radius (m), water column velocity (ms −1 ), maximum water depth (m), mean water depth (m), and Froude number.

Generalized Additive Modeling
We explained the probability of occurrence of cool water areas (i.e., troughs in the longitudinal profile) and the median water temperature from the TIR data using generalized additive models (GAMs). Generalized additive models were selected because of their ability to compute nonlinear relationships between the response and the set of explanatory variables [60]. The response variables were (1) the presence of troughs, defined as areas of cooler water than the median radiant water temperature of the river, and (2) the median radiant temperature. The GAM analysis uses smoothers to describe the empirical relationships between explanatory variables and response variables and, therefore, does not assume a specific relationship. We used the GAM function in the Mixed GAM Computation Vehicle with Automatic Smoothness Estimation (MGCV) package of the statistical program R [61,62]. For the probability of occurrence model and the median radiant water temperature, we used a binomial distribution with a logit-link function and a Gaussian distribution with a log function, respectively, to determine whether there were significant relationships between the response and explanatory variables.
Because rivers are hierarchical systems, we considered two spatial scales and selected explanatory variables that characterized temperature patterns at each scale: Reach (100 m to 1 km), and subcatchment (1 km to tens of kilometers) ( Table 2). We checked for concurvity using the MGCV package [63]. Concurvity is the existence of nonlinear dependencies among predictor variables [64], which we minimized by removing variables with a Pearson correlation greater than 0.40. Presence of concurvity in the data may lead to poor parameter estimation and loss of interpretation as the effect of a predictor on the response variable may change depending on the order of the predictors in the model if concurvity is present [64].
To automate model selection from the global model for each response, we used the Multi-Model Inference (MuMIn) R package [65]. We used Akaike's Information Criterion (AIC) with the small-sample bias adjustment (AICc) and calculated Akaike weights to assess the relative fit of each candidate model and compared candidate models for the same data with the best-fitting models with the lowest AICc [66]. We then selected the models that had an AICc difference less than 2. Models with a difference less than 2 were considered to have a substantial level of empirical support for the most plausible candidate model [66]. Finally, we selected the most parsimonious model from the models that had AICc difference less than 2.

Comparisons between TIR and Instream Measurements
The TIR data indicated a downstream warming trend, with several peaks and troughs along the longitudinal thermal profile (Figure 2).

Comparisons Between TIR and Instream Measurements
The TIR data indicated a downstream warming trend, with several peaks and troughs along the longitudinal thermal profile (Figure 2). The median radiant temperatures were 23.4 °C for the right bank and 23.6 °C for the left bank, with standard deviations of ±0.50 and ±0.48, respectively. We found that the 22 August 2017 instream surface water temperature lowest trough and highest peak coincided with the TIR lowest trough and highest peak. The lowest TIR trough was 22.3 °C and was located near river kilometer (rkm) 80 in a The median radiant temperatures were 23.4 • C for the right bank and 23.6 • C for the left bank, with standard deviations of ±0.50 and ±0.48, respectively. We found that the 22 August 2017 instream surface water temperature lowest trough and highest peak coincided with the TIR lowest trough and highest peak. The lowest TIR trough was 22.3 • C and was located near river kilometer (rkm) 80 in a narrow and deep section of the river and near a cold tributary, Ruby Creek (C5, 17.9 • C). The troughs were generally associated with either cold tributaries or springs and seepages. In contrast, the largest TIR peak was 24.7 • C. This section of the river is wide and shallow and has warm water inputs from Cusick Creek (C9, 28.3 • C) and a large side channel (25.7 • C).
The TIR sensor only detects thermal radiation emitted from the upper 0.1 mm of the water surface, but there were numerous instances where thermal stratification may have occurred. Examples of midstream stratification were observed intermittently (e.g., rkm 57 to rkm 80, rkm 88 to 102, and rkm 122 to rkm 126) and in tributary and side channel junctions where slow-moving water cooler than the main stem generally sank to the bottom when it mixed with the warmer main stem water. The instream measurements confirm some of this stratification; for example, in midstream near rkm 102 on 22 August 2017, maximum surface water temperature was 25.2 • C and minimum bottom water temperature was 21.1 • C, a 4.1 • C difference between surface and bottom water temperatures.
A total of 25 tributaries and side channels were mapped during the TIR surveys, and of those, only six were cooler than the mean temperature of the Pend Oreille River. Four of these six cooler tributaries were found between rkm 81 and rkm 93, near where the largest TIR trough was observed. However, the bottom instream loggers deployed at confluences in 2017 revealed that most of the tributaries sampled were cooler than the main stem Pend Oreille River and had strong vertical temperature gradients at their confluences ( Figure 3). The confluences of tributaries such as Renshaw Creek (C3) and Tiger Slough (C4) had average temperature differences of 9.5 • C and 11.9 • C between the surface and bottom, respectively. Streams such as Skookum Creek (C13) and Indian Creek (C16) were well mixed.

Longitudinal
We did not observe a warming trend in the 7DADM surface or bottom water temperature along the longitudinal profile, nor when we averaged the entire period of collection (17 July to 14 September), by each month. However, the daily maximum water temperature data showed several peaks and troughs along the longitudinal profile that the daily minimum water temperature did not (Figures 4a  and S1). The range of mean bottom water temperature along the longitudinal profile was 21.2-23.0 • C, and the range of mean surface water temperature was 22.7-23.4 • C. The minimum (19.6 • C) and maximum (28.8 • C) surface water temperature and minimum (19.0 • C) bottom water temperature were all observed at P4, whereas the maximum bottom water temperature (25.4 • C) was observed at P7.

Lateral: Tributary Confluences
Of the 22 confluences we surveyed, 17 were retained for analysis because they did not have any missing data (Figure 3). Contrary to the limited spatial longitudinal variability in the main stem, spatial thermal variability at the confluences was approximately an order of magnitude greater than the main stem. Overall, the range of water temperatures for all confluences during the entire period varied considerably. Surface temperatures had a wider range (7.9-32.3 • C) than bottom temperatures (7.9-27.4 • C). C17 was the most variable confluence at the surface (∆24.4 • C) and C14 was the most variable confluence at the bottom (∆16.9 • C), whereas the least variable confluence was C4 for both surface and bottom (∆2.7 • C). The 7DADM for the surface water temperature ranged from 18.8 • C (C13) to 30.2 • C (C17), and the 7DADM bottom water temperature ranged from 12.4 • C (C4) to 24.5 • C (C14). As expected, the 7DADM for the bottom water temperature was cooler than the surface water 7DADM by at least 1 • C degree at all confluences (except at C13, 18.7 • C). Remote Sens. 2019, 11, x FOR PEER REVIEW 13 of 32

Longitudinal
We did not observe a warming trend in the 7DADM surface or bottom water temperature along the longitudinal profile, nor when we averaged the entire period of collection (17 July to 14 September), by each month. However, the daily maximum water temperature data showed several peaks and troughs along the longitudinal profile that the daily minimum water temperature did not (Figure 4a, S1). The range of mean bottom water temperature along the longitudinal profile was 21.

Vertical: Water Column
Differences between daily minimum bottom and daily minimum surface water temperatures in the main stem were less than 0.5 • C, with daily minimum water temperatures near 20.3 • C (SD ± 0.24 to 0.36 • C) except for P4, where minimum bottom and surface water temperatures were 19.0 • C and 19.6 • C, respectively. Maximum daily surface water temperatures in the main stem were on average 1.2 • C warmer than maximum daily bottom water temperature. Furthermore, hourly differences between surface and bottom water temperature ranged from −0.2 • C to 4.8 • C: 17 sites had differences greater than 1 • C, and 13 sites had differences greater than 2 • C (Figures 4b and S1).

Lateral: Tributary Confluences
Of the 22 confluences we surveyed, 17 were retained for analysis because they did not have any missing data (Figure 3). Contrary to the limited spatial longitudinal variability in the main stem, spatial thermal variability at the confluences was approximately an order of magnitude greater than the main stem. Overall, the range of water temperatures for all confluences during the entire period varied considerably. Surface temperatures had a wider range (7.9-32.3 °C) than bottom temperatures (7.9-27.4 °C). C17 was the most variable confluence at the surface (Δ 24.4 °C) and C14 was the most variable confluence at the bottom (Δ 16.9 °C), whereas the least variable confluence was C4 for both surface and bottom (Δ 2.7 °C). The 7DADM for the surface water temperature ranged from 18.8 °C (C13) to 30.2 °C (C17), and the 7DADM bottom water temperature ranged from 12.4 °C (C4) to 24.5 °C (C14). As expected, the 7DADM for the bottom water temperature was cooler than the surface water 7DADM by at least 1 °C degree at all confluences (except at C13, 18.7 °C).

Vertical: Water Column
Differences between daily minimum bottom and daily minimum surface water temperatures in the main stem were less than 0.5 °C, with daily minimum water temperatures near 20.3 °C (SD ± 0.24 to 0.36 °C) except for P4, where minimum bottom and surface water temperatures were 19.0 °C and 19.6 °C, respectively. Maximum daily surface water temperatures in the main stem were on average 1.2 °C warmer than maximum daily bottom water temperature. Furthermore, hourly differences between surface and bottom water temperature ranged from -0.2 °C to 4.8 °C: 17 sites had differences greater than 1 °C, and 13 sites had differences greater than 2 °C (Figure 4b, S1).
Unlike the main stem where water temperature differences between the surface and bottom were small, temperature differences at the confluences in the tributaries were more pronounced. Hourly differences between surface and bottom water temperature ranged from -11.1 °C to 19.2 °C, Unlike the main stem where water temperature differences between the surface and bottom were small, temperature differences at the confluences in the tributaries were more pronounced. Hourly differences between surface and bottom water temperature ranged from −11.1 • C to 19.2 • C, which was nearly four times greater than the vertical variability in the main stem. Furthermore, we identified five confluences that had daily minimum bottom water temperatures of at least 1 • C cooler than the daily minimum surface water temperature, with the largest hourly differences observed for the entire period at C3 (12.4 • C) and C4 (15.2 • C). We identified eight confluences that had periods where hourly surface water was at least 1 • C cooler than bottom water. These occurrences were most pronounced at C7 (11.1 • C), C8 (4.3 • C), and C14 (4.6 • C) in late August and September when dense SAV patches trapped warm water near the streambed while surface water cooled at night ( Figure S1).

Temporal
We initiated data collection on 8 July 2017, but we restricted our data analysis to the relatively warmer 17 August to 13 September 2017 period for comparison among sites due to missing data. Water temperature in the main stem (surface and bottom) was above 20 • C except at one site (P4) during our last week of sampling in September. Diel thermal patterns in the main stem were relatively uniform along the longitudinal profile ( Figure 5), and typically water was warmest from late morning to late afternoon, with considerable variation from day to day. These diel patterns were also more pronounced at the surface than at the bottom ( Figure 5). late morning to late afternoon, with considerable variation from day to day. These diel patterns were also more pronounced at the surface than at the bottom ( Figure 5).
Diel patterns at confluences, like the main stem, revealed peaks in mid to late afternoon and troughs late at night or early morning for both surface and bottom thermographs. Over the course of the summer, most of the confluences water temperature decreased, but some were highly variable (S2).

Implications for Potential Cold-Water Refuges
The presence of a potential cold-water refuge depended on the definition used. We did not observe main stem water temperatures cooler than 18 °C during our study, so none of the main stem Diel patterns at confluences, like the main stem, revealed peaks in mid to late afternoon and troughs late at night or early morning for both surface and bottom thermographs. Over the course of the summer, most of the confluences water temperature decreased, but some were highly variable (Supplementary Material S2).

Implications for Potential Cold-Water Refuges
The presence of a potential cold-water refuge depended on the definition used. We did not observe main stem water temperatures cooler than 18 • C during our study, so none of the main stem sites were classified as a potential physiological cold-water refuge. The percentage of time during which the main stem sites met the 1 • C ambient criterion ranged from 14% (P5) to 100% (P4). However, when we used the 2 • C ambient criterion, 5 sites (T1, T3, P3, P5, and T53) were not considered potential ambient cold-water refuges, an additional 5 sites were considered potential ambient cold-water refuges less than 10% of the time (BC, T4, P6, T7 and AF), and one site (P4) decreased slightly to 98% (Table 3). Table 3. Seven-day average daily maximum (7DADM) water temperature at main stem sites where instream loggers were deployed at the surface and on the bottom of the river. Cold-water refuge persistence (percentage of time) is presented for two ambient criteria based on whether the difference between the surface water temperature and bottom water temperature was greater than 1 • C or 2 • C for at least an hour a day.

Site
Main Five confluences were classified as potential physiological cold-water refuges over 96% of the time when bottom 7DADM water temperatures were below 18 • C (Table 4), whereas only two confluences never experienced temperatures below this physiological criterion (C2, C9). The remaining 10 confluences were classified as potential cold-water refuges between 10% and 62% of the time. Over twice as many confluences (11) were always classified as potential ambient cold-water refuges in relation to the main stem, and only one confluence, C2, was never a potential cold-water refuge regardless of the ambient criterion used (2 • C versus 1 • C). The remaining six confluences were potential cold-water refuges 52-96% of the time when using the 1 • C criterion or 31-86% of the time when using the more restrictive 2 • C criterion (Table 3). Overall, persistence of potential cold-water refuges was associated with the difference between main stem 7DADM bottom water temperature and confluence 7DADM bottom water temperature regardless of criterion (Spearman r = 0.89, 0.81, 0.82, Figure 6). The spacing along the longitudinal profile among persistent physiological cold-water refuges formed by confluences was greater than 2 km but less than 41 km (median = 7 km). When we considered physiological cold-water refuges that were present at least 50% of the time, the spacing (distance between nearest neighbors) was reduced to less than 18 km, but the median distance increased from 7 to 11 km.

GIS Data
Land use (% developed land), aquifer permeability (high versus low), degree of channel confinement (m m −1 ), and tributary contributions (area in km 2 and number of tributaries) were used to characterize the riverine landscape in the Pend Oreille River (Figure 7). The percent developed land along the Pend Oreille River ranged from 3% near Box Canyon dam to 75% near the city of Cusick at rkm 101. Three areas with a high percentage of developed land (peaks) were observed: One near the town of Ione, Washington; another at the largest thermal peak near the town of Cusick; and a third peak near the town of Oldtown, Idaho, near Albeni Falls dam. Aquifer permeability was low in most of the Pend Oreille River except in two locations: At rkm 115 and between rkm 133 and rkm 144 where numerous seeps and springs were observed. Similarly, most of the Pend Oreille River was considered unconfined (>4) with the section from rkm 97 to rkm 122 having confinement values greater than 10. Additionally, there were no confined sections (<2) and only 9 km along the longitudinal profile were moderately unconfined, with most of these instances located near Box Canyon Dam and from rkm 83 to rkm 85.
The average number of tributaries per kilometer was 1.85 (SD ± 1.40). A maximum of 5 tributaries per kilometer were found at rkms 72, 84, and 87, but the longest section without tributary contributions was located from rkm 117 to 122, followed by the section from rkm 98 to 101. Most tributaries drained areas near or less than 50 km 2 . From those tributaries that drained considerably larger areas, the largest tributary was Calispel Creek, followed by C6 and C11 (Figure 7).

Three-Dimensional Hydraulic Mapping and Channel Morphology
The z-scores for water velocity showed four distinct sections of the river along the longitudinal profile. The section immediately upstream of Box Canyon to rkm 78 and the section from rkm 95 to near rkm 125 were slower than the mean water velocity, and the sections from rkm 78 to 94 and the section nearest the Albeni Falls dam were faster than the mean velocity ( Figure 8). The z-scores of the Froude number also showed like patterns to those observed for water velocity. As expected, the z-scores for wetted perimeter and channel width were similar and highly correlated (0.99, p < 0.00001). The z-scores of the cross-sectional area mirrored those of water velocity (Pearson correlation coefficient = −0.90, p < 0.0001). The z-scores of average depth and maximum depth indicated that the deepest locations in the river were near the Albeni Falls dam (rkm 131 and rkm 142), whereas most of the shallowest locations in the river were found between rkm 71 and rkm 87 with a few locations near Albeni Falls dam (Figure 8). Most of the river had slow, tranquil flow (e.g., 0.05-0.1 Fr), but some portions of the river in the upper reaches were faster and more turbulent (e.g., 0.4 Fr, depth < 1 m, and velocity 1 m s −1 ).
Remote Sens. 2019, 11, x FOR PEER REVIEW 21 of 32 near rkm 125 were slower than the mean water velocity, and the sections from rkm 78 to 94 and the section nearest the Albeni Falls dam were faster than the mean velocity ( Figure 8). The z-scores of the Froude number also showed like patterns to those observed for water velocity. As expected, the zscores for wetted perimeter and channel width were similar and highly correlated (0.99, p<0.00001). The z-scores of the cross-sectional area mirrored those of water velocity (Pearson correlation coefficient =-0.90, p < 0.0001). The z-scores of average depth and maximum depth indicated that the deepest locations in the river were near the Albeni Falls dam (rkm 131 and rkm 142), whereas most of the shallowest locations in the river were found between rkm 71 and rkm 87 with a few locations near Albeni Falls dam (Figure 8). Most of the river had slow, tranquil flow (e.g., 0.05-0.1 Fr), but some portions of the river in the upper reaches were faster and more turbulent (e.g., 0.4 Fr, depth <1 m, and velocity 1 m s -1 ).

Submerged Aquatic Vegetation
The most abundant dense patches of submerged aquatic vegetation were found from rkm 98 to 104. Aquatic vegetation density was low near both dams (rkm 55 to 71 and rkm 139 to 143) ( Figure  7).

Submerged Aquatic Vegetation
The most abundant dense patches of submerged aquatic vegetation were found from rkm 98 to 104. Aquatic vegetation density was low near both dams (rkm 55 to 71 and rkm 139 to 143) (Figure 7).

Associations between Hydrologic and Geomorphic Variables and TIR Imagery
The best reach-scale model within the set of models predicting median water temperature included confinement and the z-scores of cross-sectional area and maximum depth as explanatory variables (AIC c = 60.2, w i = 0.48, Table 5). Based on this model prediction, the z-scores of cross-sectional area, maximum depth, and confinement explained approximately 59% of the variability in the median water temperature at the reach scale ( Table 6). The best reach-scale model within the set of models predicting the probability of occurrence of cool water included confinement and the z-scores of cross-sectional area as explanatory variables (AIC c = 65.0.1, w i = 0.20, Table 4). This model explained 51% of the variability (Table 6). Cross-sectional area had a positive relationship with median water temperature and a negative relationship with probability of occurrence of cool-water areas. Table 5. Candidate models with reach and subcatchment scale explanatory variables predicting the median radiant water temperature and the probability of occurrence of cool water in the Pend Oreille River. For each model, the number of parameters (K), −2 log-likelihood (−2logL), Akaike's information criterion adjusted for small sample size (AIC c ), the difference in AIC c between the given model and the best-performing model (∆AIC c ), and Akaike weight (w i ) are shown. An asterisk (*) indicates a model selected based on parsimony in addition to AIC c scores and weights. At the subcatchment scale, the best models included distance from the dam for the set of models predicting median radiant temperature (AIC c = 34.7, w i = 0.23, Table 5), and distance from the dam and percent land developed for the set of models predicting probability of cool water occurrence (AIC c = 59.6; w i = 0.50, Table 5). These models revealed that distance from Albeni Falls Dam explained approximately 68% of the variability of median water temperature and 63% of the probability of occurrence of cool-water areas. All the models explained at least half of the variability of the response variables, but the subcatchment models performed slightly better (Table 6).

Discussion
We combined the high spatial resolution of airborne thermal infrared (TIR) imagery with the high temporal resolution of instream measurements to describe broad-scale patterns of water temperature in relation to localized drivers of heterogeneity. Our findings demonstrated that (1) lateral contributions from tributaries along the longitudinal profile dominated thermal heterogeneity in the Pend Oreille River, (2) spatial and temporal thermal variability at confluences was approximately an order of magnitude greater than that of the main stem, (3) most of the potential cold-water refuges were found at confluences regardless of the refuge definition used, and (4) the probability of occurrence of cool areas and the median radiant water temperature were closely linked to geomorphic metrics at the reach scale and to the distance from the Albeni Falls Dam at the subcatchment scale. These results highlight the importance of using multiple approaches to describe thermal heterogeneity in large impounded rivers. TIR data alone provided an incomplete picture of the thermal variability found at tributary confluences and in areas of the main stem with vertical stratification. Instream thermographs at the bottom of most confluences recorded considerably cooler temperatures than surface loggers and revealed discontinuities in the longitudinal profile that were not observed in the TIR data. These results have important implications for the identification and restoration of cold-water refuges.
This study also expands understanding of thermal heterogeneity in large impounded rivers by complementing earlier work, mostly in Europe and North America (e.g., [8,32,33]). Overall, we observed a slight downstream increase in water temperature along the longitudinal profile with multiple peaks and troughs. A downstream increase in water temperature is a pattern that has been reported in many rivers and, until recently, was our conceptual model for rivers in general [26]. Like Arscott et al. [30], we found that lateral heterogeneity dominated thermal heterogeneity. Lateral heterogeneity in our study was an order of magnitude greater than longitudinal heterogeneity, but unlike Arscott et al. [30] where this variability was related to floodplain habitat, lateral variability in our study is associated with tributary confluences, as most floodplain habitat has been modified in the Pend Oreille River. Other studies [67][68][69] have investigated lateral heterogeneity at the segment and reach scales but longitudinal heterogeneity was not considered.
Because the rivers surveyed with TIR remote sensing are generally well mixed vertically except at isolated large features such as pools [16,26], previous studies have focused mostly on thermal variability in longitudinal and lateral dimensions [26,67,68]. Persistent thermal stratification in main river channels seldom occurs because streamflow is usually turbulent enough to cause complete mixing [70]. Like previous studies on riverine thermal stratification, we found that stratification in the main stem was only persistent at one site, P4, where the river was shallow, slow moving, unconfined, and had dense patches of aquatic vegetation. However, it is important to consider temporal variability of vertical thermal gradients because of their potential use by fish as ambient cold-water refuges. We found that 10 main stem sites had differences between surface and bottom water greater than 1 • C for at least one hour per day about 50% of the time. This difference appears small, but it may provide salmonids thermal refuge during migration while fish travel to and from cold-water refuges that better meet their physiological needs.
Our results corroborate those of Wawrzyniak et al. [71] who found that TIR-derived temperatures overestimated in situ measurements at most confluences where thermal stratification was regularly observed. In addition, four confluences were stratified persistently, near 100% of the time (C1, C3, C4, and C5), two confluences were hardly ever stratified (C2, C13), and the rest of the confluences were stratified between 13 and 87% of the time. In a tributary of Three Gorges Reservoir in China, Jin et al. [72] found that the main drivers of stratification were reservoir operations (i.e., water level and water fluctuations) and air temperature which were not included in our study but are likely important factors in our study. We also speculate that temperature of the incoming tributary, volume of water, and water velocity at the confluences are important drivers influencing thermal stratification because the four confluences where we observed persistent stratification (95-100%) had low water velocity and temperature differences of at least 4 • C at the bottoms of the incoming tributary and the main stem.
These findings illustrate that the identification of a potential cold-water refuge depends on the definition used (i.e., the criteria) and the temporal scale observed. For instance, we did not observe main stem water temperatures cooler than 18 • C during our study; therefore, no main stem sites were classified as potential physiological cold-water refuges. However, when we assessed the main stem sites using our least restrictive 1 • C ambient criterion, all sites were considered cold-water refuges at depth during some portion of the day, albeit over different periods of time. Consequently, cold-water refuge definitions based solely on physiological thresholds or that exclude fine temporal and spatial scales overlook the importance of thermal heterogeneity for providing thermal refuge when temperatures exceed physiological thresholds even if it is for few hours a day [55]. Because most studies have focused on the spatial variability of cold-water refuges [1,16,45,59], cold-water refuge persistence has not been well studied. Dugdale et al. [18] documented in the Miramichi River that 61% of thermal refuges (≥0.5 • C colder than the ambient river channel temperature) were only identified once during multiple TIR surveys. In our study, we found that the persistence of confluence cold-water refuges was associated with the difference in 7DADM for the main stem and the confluence; larger differences in temperature created refuges that were persistent over longer periods, regardless of how refuges were defined.
We designed our study to describe water temperature patterns in the Pend Oreille River when conditions were warmest to understand how thermal heterogeneity during these stressful conditions influence cold-water refuge distribution for salmonids. Because water temperature was already above 20 • C when we started sampling in 2017, and water temperature in the main stem Pend Oreille River reached 18 • C during the first week of July [73], bull trout and, to a lesser extent, cutthroat trout were under thermal stress for most of the summer. Thus, using physiological thresholds alone for the main stem Pend Oreille was not very informative. However, comparing the spacing of cold-water refuges together with their persistence can help evaluate whether these refuges are adequate for salmonids migrating through the Pend Oreille River. Spacing along the longitudinal profile between persistent cold-water refuges (i.e., those that were present 100% of the time) that are formed by confluences ranged between 2 and 41 km (median = 7 km), but when we considered physiological cold-water refuges that were present at least 50% of the time, the spacing was reduced to less than 18 km, with an increase in the median separation distance (11 km). Additionally, areas that are 1-2 • C cooler than the median main stem bottom water temperature (main stem ambient cold-water refuges) but are still above the physiological thresholds can provide connectivity and can be used as stepping stones during migration [59]. Spatially continuous bottom water temperatures would be needed to determine the spacing of ambient cold-water refuges in the main stem along the longitudinal profile. Furthermore, we did not include every confluence in our analysis, so our numbers for cold-water refuge spacing are most likely overestimated, especially if localized groundwater inputs occur at depth or are too small to be detected at the surface using remotely sensed methods (e.g., TIR). Nevertheless, our spacing estimates for the confluences were within the range of variation of previous studies. Fullerton et al. [59], in a review of rivers in the Pacific Northwest, reported two to three times the median spacing of this study, but the spacing had high variability and ranged from 6 to 49 km, whereas Dzara et al. [28] in the Walker River in Nevada found that median spacing was 0.9 km and ranged from 0.3 to 37 km.
The high resolution of the geomorphic metrics measured along most of the Pend Oreille River provides a unique opportunity to examine how these metrics affect the occurrence of cool areas and median water temperature at the 1-km scale. The strong associations with channel confinement corroborated findings of previous studies [17,22,30,74]. Cross-sectional area, another variable with strong explanatory power for both reach-scale models, was also inversely correlated with water column velocity (Pearson correlation = −0.90). These robust associations support our hypotheses that channel confinement, water velocity, and cross-sectional areas explain different mechanisms that directly affect water temperature patterns. More specifically, at confined sections where the river is narrow and deeper, two alternative mechanisms may be at play depending on water velocity. At high velocities, water may be cooler as a result of mixing of the water column and topographic shading, whereas at low velocities, water may be warmer at the surface due to thermal stratification. In contrast, at unconfined sections, where channel width and cross-sectional area are largest, warmer water temperatures are the result of greater surface area being exposed to solar radiation.
At the subcatchment scale, our models suggest that the distance from Albeni Falls Dam was the most important factor determining the occurrence of cool water areas and explaining median water temperature. This result was unexpected for our study because of the weak longitudinal downstream patterns observed in the TIR data. However, our statistical results corresponded to those found in other studies when location in the network was considered [74,75]. The mechanism for this downstream increase was not entirely clear from our results. A regional study of the impact of dams on water temperature in Eastern Canada found that dams with impounded runoff index (storage capacity divided by median annual runoff) less than 10% had little impact on the thermal regime of medium-size rivers [76]. The Pend Oreille River is at least an order of magnitude larger than the rivers described by Maheu et al. [27], and has an impounded runoff index of less than 1%, which is typical of run-of-the-river systems. However, during most of the 2001 summer when the TIR data were collected, assuming no major withdrawals, the reservoir was storing water. The amount of time that the daily reservoir outflow was less than daily reservoir inflow was 82% in July, 70% in August, and 13% in September. The largest daily difference between outflow and inflow was 18% on August 7. Despite the Pend Oreille River having run-of-the-river characteristics on an annual basis, the reservoir appears to operate as a storage reservoir in the summer and is likely to experience a warming effect. Our statistical results did not point to the size or number of tributaries as explanatory variables. These results, together with TIR imagery from tributary plumes suggest that lateral contributions were operating at the reach scale, except at the lowest trough (rkm 80 to 90), where C6 and C5 were large enough to potentially have an impact on downstream water temperature patterns.
We recognize that our TIR survey was conducted in 2001 and our instream thermographs were deployed in 2017. Although, this mismatch could influence our interpretations, Fullerton et al. [26] found that spatial patterns of water temperature from repeat TIR surveys were consistent among years and, thus, comparisons with broad spatial patterns data collected in 2017 should be consistent over time given that the river morphology has not changed. Similarly, 3-D mapping was conducted at considerably different discharges: Average discharge from the second survey was almost three times the discharge from the first survey (268 m 3 s −1 versus 720 m 3 s −1 ). We used z-scores to normalize these values to account for the temporal variability. The z-scores directly compare the patterns of variability rather than the magnitude of the variables and can be used as explanatory variables by removing any bias potential from sampling during different flow conditions. We also recognized that some confluences were highly spatially variable due to either their low stream bed stability (substrates were mostly sand or silt) or their dense aquatic vegetation, so placement directly below the thermal plume and foreseeing the path of the tributary was not possible. For these types of streams (e.g., C14, C8, C7, and C17), distributed temperature sensing (DTS) cables could have captured the plume path more consistently than instream loggers as the water elevation decreased. However, budget constraints and potential logistical and maintenance issues prevented our use of DTS. Incorporating this variability can be informative to management for potential cold-water refuges at these confluences with fluctuating conditions. Finally, our study was limited to one large, impounded river with specific geomorphic and hydrologic characteristics and flow operations. However, inferences from this study may still be relevant to other run-of-the-river reservoirs with similar physical characteristics.
Many questions remain about the thermal heterogeneity of large impounded rivers, and the implications for cold-water refuges. Because thermal heterogeneity is likely dampened in the main stem river, studies that consider flow management scenarios affecting the size, spacing, and persistence of cold-water refuges may be critical to develop management strategies for conserving, augmenting, and restoring areas suitable for salmonids. Our study also raises questions about the sufficiency of the potential cold-water refuges described here. Understanding the adequacy of cold-water refuges for bull trout and cutthroat trout migration is important for the survival of these species in the Pend Oreille River system. Suitability of cold-water refuges is not limited to water temperature; therefore, other factors such as dissolved oxygen and food availability should also be considered when determining the potential use of a cold-water refuge. For instance, we observed bull trout at C3 and C17 where dense abundant vegetation and night hypoxic conditions were present. Dense SAV may negatively impact these areas by decreasing dissolved oxygen at night near the bottom where it is cooler and may prevent fish from using these areas if SAV becomes so thick that it creates a barrier. Understanding the interactions between submerged aquatic vegetation, water temperature, and dissolved oxygen may provide management strategies at confluences where these conditions occur. Lastly, our findings suggest that the distance from Albeni Falls Dam was an important variable explaining occurrence of cool areas and median water temperature. Hence, understanding statistical associations of hydroclimatic variables related to conditions upstream of the dam, such as upstream water temperature, discharge, ratio of reservoir outflow to inflow, and air temperature to downstream water temperature could be used to examine whether Albeni Falls Dam affects downstream water temperature patterns.
Given that climate change will increase temperatures globally, thermal heterogeneity will become increasingly important for fish to behaviorally thermoregulate during times of thermal stress even as the ability of fish to locate cold-water areas is also diminished because of habitat homogenization [77]. Because unidentified refuges are at risk of being lost through habitat modification and degradation or unsuitable management [78], locating cold-water refuges is the first step to manage, restore, augment, and create thermal heterogeneity to conserve salmonid habitat [79].

Conclusions
Our study described thermal heterogeneity across space and time and identified potential cold-water refuges for salmonids in the Pend Oreille River, a large impounded river with an epilimnetic-release located in inland northwestern USA. Our findings demonstrate that the complexity of thermal patterns in impounded rivers require a combination of approaches, including assessment and modeling to understand the implications of thermal heterogeneity for cold-water fishes. These findings also highlight the importance of considering thermal heterogeneity when restoring and maintaining aquatic landscapes that can support native salmonids.
The results were consistent with our hypotheses that water temperature in the main stem of the Pend Oreille River increased as water moved downstream, but tributaries created discontinuities in the longitudinal profile. These lateral discontinuities (1) constituted most of the thermal heterogeneity in the river, (2) provided most of the potential cold-water refuges for salmonids, and (3) were more persistent in areas where thermal differences between the main stem and confluences were largest. In addition, we found that vertical thermal stratification was generally present at slow-moving confluences where vertical temperature gradients isolated cool water near the river bottom. Modeling results support mechanistic understanding of the interaction between confinement, cross-sectional area, water velocity, and depth to explain median water temperature and the probability of occurrence of cool-water areas. These results can be used to develop management and restoration strategies to enhance thermal heterogeneity.