Land-Cover and Climatic Controls on Water Temperature, Flow Permanence, and Fragmentation of Great Basin Stream Networks

: The seasonal and inter-annual variability of flow presence and water temperature within headwater streams of the Great Basin of the western United States limit the occurrence and distribution of coldwater fish and other aquatic species. To evaluate changes in flow presence and water temperature during seasonal dry periods, we developed spatial stream network (SSN) models from remotely sensed land-cover and climatic data that account for autocovariance within stream networks to predict the May to August flow presence and water temperature between 2015 and 2017 in two arid watersheds within the Great Basin: Willow and Whitehorse Creeks in southeastern Oregon and Willow and Rock Creeks in northern Nevada. The inclusion of spatial autocovariance structures improved the predictive performance of the May water temperature model when the stream networks were most connected, but only marginally improved the August water temperature model when the stream networks were most fragmented. As stream network fragmentation increased from the spring to the summer, the SSN models revealed a shift in the scale of processes affecting flow presence and water temperature from watershed-scale processes like snowmelt during high-runoff seasons to local processes like groundwater discharge during sustained seasonal dry periods.


Introduction
Formerly intact and connected stream networks are becoming increasingly fragmented, which threatens aquatic ecosystems [1][2][3]. Most studies of fragmentation in stream networks have focused on relatively fixed barriers to movement such as culverts [4] and dams [5]. Other dynamic factors can also fragment populations of stream biota, including interactions with other species [6,7], temperature [8][9][10], and loss of surface flows [11][12][13]. Dynamic spatial patterns of fragmentation in stream networks can have variable consequences for stream biota such as increased mortality and

Data Collection
Field data collection focused on quantifying the stream temperatures and the presence of wetted channels (as inferred by stream temperatures) [58,59]. Sampled sites in both stream networks followed Schultz et al. [28] and were distributed across accessible locations (e.g., public lands or private lands with authorized access) to include a range of environmental gradients (e.g., perennial to non-perennial streams, elevation gradients). Sampling was conducted in 2015, 2016, and 2017 at 106 sites in the WW watershed ( Figure 1a) and in 2016 and 2017 (water years; 1 October to 30 September) at 109 sites in the WR watershed ( Figure 1b). All sites were georeferenced using global positioning system units. Due to deployment times exceeding a year and losses of instruments placed within stream channels, we could not completely control the resulting spatial distribution of recovered loggers. Methods of statistical analysis used to predict stream temperatures and flow presence can accommodate a variety of sampling configurations [60].
Continuous records of temperature were measured at each location using data-logging thermistors (HOBO U-22 water temperature logger; Onset Computer Corporation, Bourne, MA, USA). These instruments were programmed and calibrated using protocols described by Heck et al. [61]. Temperatures were recorded at hourly intervals throughout the water year. In the field, thermistors were shielded with polyvinyl chloride (PVC) housings and secured in the deepest parts of stream channels [61] to ensure that if the thermistor was exposed to air that we could determine whether the stream dried [59]. In some locations, air temperature loggers were also deployed. These instruments were attached to nearby trees and enclosed within a radiation shield (Onset model RS1, Onset Computer Corporation, Bourne, MA, USA).

Statistical Models of Flow Presence and Stream Temperature
Spatial stream network (SSN) models, a class of generalized linear models that account for spatial autocovariance within stream networks, have been used to predict a variety of hydrologic, geochemical, and biological attributes of stream networks. Many SSN modeling efforts have focused on water temperature (e.g., [62][63][64]), in part because it is a master variable that determines the occurrence, distribution, and abundance of ectothermic species, but also because of its sensitivity to climate. SSN models applying a Gaussian distribution have been used to predict other continuous data like water chemistry [65,66] and species distribution [67]. SSN models can also be applied to non-normal data such as counts and occurrences through Poisson and binomial distributions, respectively. For example, Isaak et al. [68] developed an SSN model to predict a binary variable of bull trout (Salvelinus confluentus) presence. We are not aware, however, of any efforts to predict a binary variable of flow presence using an SSN logistic regression model. Indeed, the focus of most SSN modeling efforts has been on perennial streams, whereas flow permanence and consequences for stream network fragmentation have not been explicitly considered.
Statistical stream network (SSN) models developed by Ver Hoef and Peterson [69] were adapted for the WW and WR study areas to predict temperature and flow presence as a function of covariates that describe spatial and temporal properties of the stream network and climate. SSN models expand a standard linear model (Equation (1) where is a vector of response variables, is a design matrix of covariate values, is a vector of regression coefficients, and is a vector of random errors to include spatially autocorrelated random effects. Inclusion of spatial structures was necessary to address spatial autocorrelation in the data that would invalidate assumptions of independence and introduce bias into parameter estimates. Ver Hoef and Peterson [69] used a moving average approach to create autocovariance models based on hydrologic distances within the stream network and propose combining them with autocovariance models based on Euclidean distance. Autocovariance models based on hydrologic distance included tail-up and tail-down models, named for the shape of the moving average function with respect to the stream network, that allow for autocorrelation based on flow-connected and flow-unconnected relationships, respectively. A number of moving average functions are available within the SSN, but in the absence of adequate data to differentiate between different spatial autocorrelation models, we limited our models to exponential moving average functions, which describe the relative influence of the moving average function over hydrologic distance. A pair of sites was considered "flow connected" if water originating from the upstream site passed by the downstream sites; if this condition did not occur, the sites were considered "flow unconnected". Vectors of zero-mean random effects based on the Euclidean (zEUC), tail-up (zTU), and tail-down (zTD) models were added to the standard linear model (Equation (1)) such that: The stream autocovariance models were weighted based on upstream contributing area. Two types of SSN models were developed to model different classes of response variables: continuous measures of mean monthly water temperature between May and August were modeled with a Gaussian distribution; binary indicators of flow presence between May and August were modeled with a binomial distribution through a logistic regression model. Spatial autocorrelation within the data was evaluated using Torgegrams that quantify spatial dependence between observations that were hydrologically connected (i.e., "flow-connected") versus observations that were hydrologically disconnected ("flow-unconnected"). Pairs of hydrologically connected observations reside on the same flow path such that water flows from the upstream site to the downstream site, whereas no common flow path links hydrologically disconnected pairs. In nonfluvial settings where a dendritic network structure and directional flow are absent, semivariograms are typically developed using Euclidean distance, but distance plotted on a Torgegram is the hydrologic distance equal to the length of the stream separating two observations. Spatial SSN models were fit to each continuous and binary response variable using a backward stepwise regression technique [63]. Each model of mean monthly water temperature from May to August initially included all five candidate covariates, and the worst fitting non-significant (p-value > 0.10) covariates were sequentially removed until all parameters had a p-value < 0.10. A mixture of exponential tail-up, tail-down, and Euclidean autocorrelation models were included with each spatial SSN model. To evaluate the contribution of spatial autocorrelation functions, nonspatial multiple linear regression models, which excluded spatial autocorrelation functions, were developed using the same backward stepwise regression technique. Spatial and nonspatial models of flow presence were similarly developed using the backward stepwise regression technique. Random effects for year and site were included within each of the models to account for errors associated with repeated measurements at individual sites and across multiple years. Scripts and input data to run SSN models presented within this study were archived and are available for download [70]. The performance of stream temperature models, which used a Gaussian distribution, were assessed with three metrics: mean absolute prediction error (MAPE) between observed temperatures and leave-one-out crossvalidation (LOOCV) predictions; root mean square prediction error (RMSPE) from LOOCV predictions; and predictive r 2 calculated as the square of the correlation coefficient between observed temperatures and LOOCV predictions. The performance of logistic regression models of flow presence, which used a binomial distribution, could not be estimated with these metrics; instead, the discrimination ability of each model was estimated by calculating the area under the curve (AUC) of the receiver operating characteristic (ROC; [71]). Comparisons between spatial and nonspatial models for the same response variables were made by calculating the spatial Akaike information criterion (AIC; [72]), which penalizes models for both the number of covariates and autocovariance parameters.

Diagnosis of Flow Presence/Absence
Differences in the range of daily temperature variability that occur due to the contrast between the heat capacity of air and water have previously been exploited to classify wet and dry conditions within streams [59,73,74]. We used complementary graphical and analytical techniques to determine if observed temperatures recorded by in-channel thermistors indicated wet or dry conditions in the stream channel at a daily interval. Prior to analysis, temperature measurements recorded before and after thermistors were deployed in the stream channel were removed from each record. For each site, hourly temperature measurements, daily temperature range, and daily standard deviation temperature measurements were summarized from May to August of each year and visually inspected to identify abrupt transitions in the daily range and standard deviation that could potentially indicate a change in state (e.g., wet to dry, dry to wet). In addition, temperatures below 0 °C or above 30 °C were classified as dry in most instances.
The graphical determination of wet and dry conditions at each temperature measurement site was complemented by an analytical approach based on comparison of the daily standard deviation of hourly temperature measurements recorded by in-channel thermistors to that recorded by paired out-of-channel thermistors that recorded nearby air temperature. Most in-channel thermistors were compared to a nearby air temperature thermistor, but in the absence of a nearby air temperature thermistor, air temperature data were obtained from nearby weather stations. To smooth the variability in the daily time series, we calculated weekly moving averages of daily standard deviations for stream channel and air temperatures.
The ratios of the daily standard deviation of stream channel temperature (σTw) to the daily standard deviation of air (local or regional data) temperature (σTa) were used to identify a specific cutoff value for ratios of σTw:σTa and diagnose dry conditions. To select an efficient cutoff, electrical resistance data indicating known drying (or wetting) cases from WW were used [59]. A numerical search was performed on the values of potential cutoffs of σTw:σTa between 0.5 and 1.0, with a step size of 0.01. We found that, in WW, a cutoff ratio of 0.8 or greater produced the highest classification rate, as measured by the proportion of correct wet/dry classifications out of the total number of daysite combinations classified ( Figure 2). The classification rate of the analytical technique using a 0.8 cutoff ratio was about 90 percent based on electrical resistivity data in the WW watershed between April and August, which was slightly greater than the classification rate of the graphical technique over the same time period. Therefore, where a discrepancy existed between the wet/dry classification assigned by the analytical and graphical techniques, the analytical technique superseded the graphical technique. By adopting a range of subjective visual and more objective numerical methods for diagnosing flow presence/stream drying, we were more confident in identifying locations that stayed dry or wet during the season of interest (May-August), or that changed in state. These methods identified the presence or absence of water within the channel and did not differentiate between moving and standing water.  [59]. Wet and dry conditions were verified by using electrical resistors. The dashed vertical line represents a cutoff of σTw:σTa = 0.8.

Summary of Temperatures and Flow Presence/Stream Drying
Continuous response variables for water temperatures, including mean monthly water temperatures for May (MayTw), June (JunTw), July (JulTw), and August (AugTw) were calculated at sites that included data for at least 90 percent of the month and at which water was inferred to be present for all days of the temperature record. Sites that included any dry days were excluded. The total number, for all years, of temperature measurement observations included for each of the continuous response variables was 244 (MayTw), 251 (JunTw), 217 (JulTw), and 194 (AugTw), indicating a decrease in the number of stream channels with flow present from May through August. At each temperature measurement site, a binary response variable of flow presence was calculated on May 15 (May15Wet), June 15 (Jun15Wet), July 15 (Jul15Wet), and August 15 (Aug15Wet), where zero denotes a dry channel and one denotes a wet channel. The condition of flow presence was determined on a single day of each month to remove ambiguity regarding the duration of drying at individual sites.

Model Covariates and Hydrography
Heat is transferred into and out of a stream by radiative, conductive, convective, and advective processes, which differ in their relative influence depending on the climatic, hydrologic, and geomorphic regime of a stream and its contributing watershed. The temperature of water at a point in the stream network indicates the accumulated gains and losses of heat through these processes [75]. Similarly, stream discharge, as well as the presence of water itself, results from the cumulative water input to and outflow from the stream. Five covariates related to physical processes hypothesized to affect stream temperature and flow presence (snowmelt, groundwater inflow, evapotranspiration, riparian vegetation abundance, interannual variability in discharge, and interannual variability in air temperature) were introduced as covariates into predictive models of stream temperature and flow presence (Table 1). Candidate covariates represented both temporal and spatial variability in the processes affecting water temperature and flow presence. Multicollinearity within model covariates was assessed by calculating variance inflation factors (VIF) for each model. VIF for model covariates were generally low (<5) and none exceeded 10, above which the strong dependence between covariates impedes model parameterization and interpretation [76]. Landscape, hydrologic, and climatic variables were considered for inclusion as covariates in statistical models of water temperature and presence. Landscape variables such as elevation and drainage area emerged as significant covariates in previous statistical models of stream temperature (e.g., [63,64,77,78]) but are related to multiple physical processes that confound interpretation of their mechanistic link to water temperature. For example, elevation is related to several factors that may influence water temperature, including air temperature, precipitation, snowpack accumulation, and the type and distribution of riparian vegetation. Although variables like elevation and drainage area may explain much of the variance in water temperature throughout a drainage network, they cannot be influenced themselves by management actions. Consequently, elevation and drainage area were excluded from consideration as covariates within models presented here in favor of variables that can be more directly linked to specific physical processes affecting stream temperature.
Three covariates, which varied both spatially and temporally, were calculated annually for each year of observed data within each of the study watersheds: presence of snow in contributing watershed on 1 May (SWEws); mean riparian normalized difference vegetation index (NDVI; NDVIrip); and mean riparian evapotranspiration (ET; ETrip). Gridded snow water equivalent on 1 May was estimated by SNODAS [52] at a 1-km resolution for each year of observations in the model. SWEws was defined as a categorical covariate of the presence or absence of snow on 1 May in the contributing watershed upstream of each temperature and flow presence observation. During the spring, snowmelt lowers the stream temperature and the potential for a loss of surface flow by contributing cool water to the stream network. The relation between vegetation and evapotranspiration with water temperature and presence of surface flow, however, is more complex and feedback mechanisms between the presence of flow, the growth of riparian vegetation, and evapotranspiration exist. Vegetation reduces the input of solar radiation, which typically comprises the largest input of heat to smaller streams [79]. The presence of vegetation, however, may also indicate the presence of shallow groundwater levels [80] that may be associated with cool groundwater inputs and reduced stream temperatures. The evapotranspiration of water from the riparian zone reduces streamflow [81], thus decreasing the potential for surface flow presence, but also contributes to evaporative heat loss thus cooling the stream. ET was estimated using the Operational Simplified Surface Energy Balance (SSEBop) model [82] with 30-m resolution Landsat 8 scenes with less than five percent cloud cover. NDVI, an index of vegetation cover over a landscape derived from the ratio of reflectance of red and near infrared radiation measured by satellites, was estimated from the same 30-m resolution Landsat 8 scenes used by SSEBop to estimate ET. NDVI values range between −1 and 1 with negative values indicating a lack of vegetation. NDVI and ET were estimated over the four-month period from 1 May to 1 September for each year from 2015 to 2017 and assigned to temperature and flow presence observations of the corresponding year. The three-month estimation period for the NDVI and ET increased the reliability of model estimates over shorter (e.g., monthly) time periods [83]. NDVI and ET were then spatially averaged for a 200-m wide riparian corridor for the stream network upstream of each temperature and flow presence observation to generate the riparian-averaged covariates NDVIrip and ETrip. The extent of riparian vegetation in both the WW and WR watersheds varied substantially from 10 m or less in some headwater drainages to nearly 200 m in parts of the mainstem streams ( Figure 3). The 200-m wide riparian corridor width to calculate NDVIrip and ETrip was selected to provide a consistent area around the WW and WR drainage networks that encompassed the potential zone of influence of vegetation on riparian shading and evapotranspiration. To account for interannual climatic variability, discharge and air temperature were introduced as spatially invariant, temporal covariates. The mean monthly temperature (MonthTa) from May to August was averaged over each study watershed from PRISM climate data [51] to account for the expected positive relationship between air and water temperature. We expected discharge to have a negative relationship with stream temperature and was accounted for by introducing monthly averaged unit discharge (discharge per unit contributing area; MonthQ) as a covariate. Because no stream gages existed within either of the study basins during the period of water temperature and presence observations, unit discharges were calculated for neighboring watersheds with contemporaneous discharge data. The mean monthly unit discharges calculated for USGS (U.S. Geological Survey) gage 10352500 (McDermitt Creek near McDermitt, NV, USA) and USGS gage 10321940 (Maggie Creek above canyon near Carlin, NV, USA) [84] were assigned to temperature and flow presence observations for the WW and WR watersheds, respectively.
The hydrography of the WW and WR watersheds was imported from the National Stream Internet (NSI; [85]) data set, which was processed from medium-resolution NHDPlusV2 [86] to meet SSN model requirements for topological consistency by removing braided and diverging channels and confluences with more than two tributaries. Spatial Tools for the Analysis of River Systems (STARS; [87]) was used to further pre-process the spatial and topological attributes for SSN model hydrography, and the stream temperature and flow presence measurement and prediction sites. STARS was then used to create an SSN object for the SSN package version 1.1.13. [88] in R [89].

Flow Presence and Temperature
The annual recovery rates of deployed water temperature loggers ranged from 70 to 97 percent in WR (2017) and WW (2015), respectively. The observed mean monthly water temperatures across all site-years were lowest in May (10.3 °C) and highest in July (15.7 °C); whereas the percentage of sites where flow was present was highest on 15 June (64 percent) and lowest on 15 August (52 percent). The highest observed mean monthly water temperatures occurred in 2015 for all months between May and August, but the lowest observed mean monthly temperatures occurred in 2017 during May and June and in 2016 during July and August. The observed mean monthly temperatures from all months were lowest in the WR watershed and highest in the WW watershed, which was lower in elevation but was more northerly in aspect and at a higher latitude.
Flow presence also varied across watersheds and years (Table 2) (Table 2).

Seasonal Changes in Spatial Dependence of Stream Temperature
Torgegrams of mean monthly stream temperature for the four time periods (MayTw, JunTw, JulTw, and AugTw) were evaluated for individual years from 2015 through 2017 (Figure 4). For flowconnected sites, semivariance of MayTw, JunTw, JulTw, and AugTw was low at short hydrologic separation distances and increased with distance, indicating a broad-scale trend in water temperature. Conversely, spatial autocorrelation was evident among flow-unconnected sites in May and June and leveled off at relatively short (~15 km) separation distances, indicating smaller-scale patchiness in water temperature. Spatial dependence among flow-unconnected sites was more pronounced in May and June ( Figure 4A,B) and became progressively less apparent in July and August ( Figure 4C,D). These relations were largely consistent across years. Semivariance is depicted as a function of separation distance between temperature observations along the stream channel. Water temperature data in 2015 were only available for WW, so Torgegrams for that year were calculated using data from WW, whereas the Torgegrams for 2016 and 2017 incorporated data from both WW and WR. Size of circles is proportional to the number of site pairs averaged to calculate semivariance for each separation distance. The minimum number of site pairs used to calculate semivariance for each separation-distance interval was 30.

Models of Stream Temperature
Spatial stream network models outperformed nonspatial multiple linear regression models of mean monthly water temperature for all four months modeled as measured by higher LOOCV r 2 pred, lower RMSPE, and lower MAPE ( Table 3). The AIC scores of spatial models were 14 to 217 points lower than nonspatial models even after accounting for the additional autocorrelation function parameters included in the spatial models. Predictive performance of both spatial and nonspatial models was highest for the MayTw models (spatial model LOOCV r 2 pred: 0.89) and lowest for the AugTw models (spatial model LOOCV r 2 pred: 0.66); similarly, the relative improvement of AIC scores from nonspatial to spatial models (i.e., by including spatial autocovariance structures) was highest for the MayTw models (ΔAIC: 217) and lowest for the AugTw models (ΔAIC: 14). The relative proportion of variance explained by different spatial structures changed correspondingly between models of MayTw and AugTw. Although the Euclidean autocovariance structure explained 67 percent of variance within the MayTw model, it did not explain any variance within the AugTw model. The tail-up and tail-down autocovariance structures based on the stream network explained the increasing proportions of variance from the models of MayTw to AugTw. Table 3. Summary of spatial stream network (SSN) models predicting mean monthly stream temperature data for the Willow-Whitehorse and Willow-Rock watersheds (significance levels of covariates: p < 0.001 (***), p < 0.01 (**), p < 0.05 (*), p ≥ 0.05 (-)). The percentage of variance explained by covariates was highest for the spatial model of JunTw (57 percent) and lowest for the spatial model of AugTw (16 percent). The relation between covariates and mean monthly water temperatures in fitted models was both consistent and inconsistent with a priori hypotheses. For models of all mean monthly water temperatures, NDVIrip was consistently among the most significant covariates and was associated with cooling water temperature for each month (Table 1). SWEws was significant in the spatial models of MayTw and JunTw, but not for any of the other months, indicating the timing of peak snowmelt in both study watersheds. The mean monthly temperature was significant in all spatial models of mean monthly water temperature except for the model of MayTw, but its effect was opposite in sign of our hypothesis for the models of JulTw and AugTw. Similarly, mean monthly discharge was significant in spatial models of MayTw, JulTw, and AugTw, but its effect was opposite in sign compared to our hypothesis for the spatial models of JulTw and AugTw.

Response Variable
The mean monthly water temperatures from May through August were predicted at unsampled locations at 100-m intervals to develop a continuous map of temperature across the entire stream network of the WW ( Figure 5) and the WR watersheds ( Figure 6). The mean monthly water temperatures from May through August were generally higher in the WR watershed than in the WW watershed. In general, headwater tributaries where NDVI was higher were cooler than mainstem streams where NDVI was lower for all modeled months in both watersheds. The relative difference in temperature between headwater and mainstem streams increased from May through August, and longitudinal heterogeneity in predicted water temperatures also increased from May through August. The mean monthly water temperature was calculated for all stream segments contained within the NSI hydrography, but the presence of water within these segments varied from May to August, as estimated by models of flow presence.

Models of Flow Presence
The SSN logistic regression models fitted to flow presence data are summarized in Table 4. These models predict the probability of flow on the days that we selected for this analysis: 15 May, 15 June, 15 July, and 15 August. The predictive performance of nonspatial models of flow presence was high within the area under curve (AUC) of the receiver operating function (ROC) between 85 percent (May15Wet model) and 91 percent (Jun15Wet model) indicating that flow presence models discriminate well between wet and dry conditions. The inclusion of spatial autocovariance structures increased the predictive accuracy of models of flow presence (Table 4), but only marginally. For example, the AUC of the ROC calculated for the Aug15Wet model increased by seven percent, which was the maximum for all four models of flow presence. NDVIrip was consistently the most significant covariate for all models of flow presence and was positively related to flow presence, as hypothesized. SWEws was also significant in the May15Wet model and was positively related to flow presence. The mean monthly discharge was significant in the Aug15Wet model, but not in any of the other models. Table 4. Summary of SSN logistic regression models predicting flow presence data for the Willow-Whitehorse and Willow-Rock watersheds (significance levels of covariates: p < 0.001 (***), p < 0.01 (**), p < 0.05 (*), p ≥ 0.05 (−); area under curve of receiver operating characteristic (AUC)). The probability of flow presence was predicted at 100-m intervals across the stream networks of both watersheds on 15 May, 15 June, 15 July, and 15 August (Figures 7 and 8). The most spatially extensive and contiguous wetted stream network was predicted by the 15 May model (Figures 7a  and 8a), whereas the least spatially extensive and most fragmented stream network was predicted by the 15 August model (Figures 7d and 8d). The flow presence and fragmentation were not consistent between the WW and WR watersheds but were more pronounced in the WR watershed ( Figure 8). The probability of flow presence predicted by the August 15 model exceeded 50 percent throughout most of the mainstem streams of the WW watershed including Willow, Whitehorse, and Little Whitehorse Creeks and their major tributaries, with limited stream drying and fragmentation predicted. In the WR watershed, flow was predicted to be present throughout mainstem streams, including Willow and Rock Creeks, as well as smaller tributaries like Frazer Creek by the 15 May model, whereas a predicted flow presence exceeding 50 percent was limited to Willow Creek and several headwater streams by the 15 August model.

Discussion
In this study, we sought to quantify annual, seasonal, and spatial variability in stream temperatures and flow permanence (i.e., the maintenance of flow presence) in drought-sensitive stream networks, and to quantify how different landscape processes (climate and land cover) drive stream-network fragmentation. Our results show that stream temperature and flow permanence varied seasonally and spatially at local scales ( Figures 5-8). These findings are consistent with previous work in the region and in similar environments elsewhere, including studies in the WW basin [28,59] and other stream networks (e.g., [90][91][92]), where marked patterns of variability occurred at even the finest level of observation [26,77,93,94]. Although quantifying such patterns of variability can require enormous effort, a larger challenge lies in understanding the processes that drive them. Here, we considered two classes of landscape processes ( Table 1) that we were able to quantify consistently with available spatial data across the two networks studied.
The major drivers of stream temperatures and streamflow permanence are climatic in unmodified rivers. Influences of climate are typically evident at broad regional scales [94]. Broad-scale models of stream temperatures [64,95,96] and streamflow permanence [97] have identified the importance of air temperatures and precipitation, respectively, as strongly associated with these responses. We focused on a smaller spatial extent, however, and associations between stream temperatures and streamflow permanence were not always in agreement with predictions based on hypothesized effects of climatic drivers (Table 1). For stream temperatures, the May snow water equivalent was important for temperatures in May and June, but not for later months (Table 3). This may be expected as the snowpack is largely depleted after June in these watersheds, and other processes can play a greater role in driving hydrologic responses. Unlike snow water equivalent, associations with monthly air temperatures and stream discharge were not consistent across months. Although climate is an important higher-level driver of stream temperatures and streamflow permanence at broader scales [94], variation in the influences of these factors on hydrologic processes within each of the networks studied herein may not have been great enough to detect associations.
We considered two indicators of the influence of land cover on stream temperature and streamflow permanence: a measure of evapotranspiration, and vegetation cover as indicated by NDVI. Although evaporation represents a net export of heat from streams [75], rates of evapotranspiration increase as temperature increases and therefore may be associated with warmer stream temperatures. Evapotranspiration can also be strongly linked to water that is potentially available to streams, particularly in arid regions [98,99]. In this study, however, we did not detect an association between evapotranspiration and water temperature or streamflow permanence. Overall, these findings lead us to conclude that the measure of evapotranspiration we employed does not have a direct, strong influence on stream temperatures or streamflow permanence, although evapotranspiration can be a significant driver of both responses. In contrast to evapotranspiration, riparian NDVI, a variable that could indicate multiple hydrologic processes (Table 1), was consistently associated with both temperature and streamflow permanence across all seasons.
As hypothesized, riparian NDVI was related to cooler water temperatures and an increased probability of surface streamflow. NDVI, which correlates with vegetation productivity, has been widely used to remotely assess vegetation distribution and its interactions with other parts of ecosystems [100]. In small streams such as those studied here, riparian vegetation plays a primary role in regulating temperatures by reducing shortwave radiation from solar insolation through shading [101][102][103]. Although NDVI does not explicitly resolve changes in shortwave insolation due to riparian shading, such as cumulative riparian solar exposure, it does indicate the presence or absence of riparian vegetation that regulates stream temperature in a similar manner to the percentage of riparian canopy cover. In an arid landscape, NDVI may not only differentiate between areas with and without riparian vegetation, but may also reveal areas of shallow groundwater [57,104]. Shallow groundwater, which sustains vegetation and increases NDVI, may discharge into streams where it both increases the probability of flow permanence and buffers summer stream temperatures by contributing cool groundwater discharge to streams.
More broadly, NDVI has been directly applied across the Great Basin to quantify habitats for terrestrial species, such as the greater sage-grouse (Centrocercus urophasianus; [105,106]), and as an indicator of (1) riparian recovery in response to improved grazing management, influences of beaver, and wildfire [56,57,80,107,108]; and (2) longer-term climate-related changes in hydrologic and vegetative conditions [104,109]. Interestingly, direct associations between NDVI and instream responses of macroinvertebrates or fish can be less apparent [110,111]. Based on our findings and those of many prior studies, however, NDVI has broad potential as an indicator of riparian and instream conditions in the Great Basin, though additional work is needed to refine spatial and temporal variability in such relationships [112]. For example, in WW, field-based data on pre-and post-burn responses of stream temperatures to a large and severe wildfire showed a high degree of variability [28]. In our study, unexplained variability in stream temperature (Table 3) may indicate the importance of local processes not captured by NDVI or other model covariates.
The extent to which stream networks are fragmented depends on the specific ecological responses of interest. The streams that we studied are designated critical habitats for Lahontan cutthroat trout, a threatened species [111]. The distribution of this species is strongly influenced by water temperature and perennial stream flows [39,41]. Other species may respond less to temperature and more to the presence of available surface water. We were generally successful in developing models to predict stream temperatures and streamflow permanence. Although we collected highresolution stream-temperature data within the study networks, predictions of stream temperature were statistically less precise than those for a regional model of stream temperature [64]. In other words, the finer-grained variability of temperature within local watersheds can be very complex and more difficult to predict than variability at larger spatial extents. This result underscores the importance of considering scale in applying results from regional [64] and local (our study) models of in-stream hydrologic responses [21] to understand stream-network fragmentation. Measures of model performance for predicting stream temperatures declined from spring (May) through late summer (August), perhaps indicating that as the hydrologic stream network became more fragmented (intermittent), local processes not captured in our spatial covariates were driving localized heat budgets at scales of 10-100 m. In contrast to stream temperature, our predictions of streamflow permanence (probability of drying) were generally much more precise than a regional model of flow permanence [97], particularly in June, July, and August (Tables 1 and 2).
As stream networks become progressively fragmented during the summer dry season, local processes may become increasingly important influences on water temperature and flow permanence as the contribution of processes operating at the extent of watersheds decreases [92,113,114]. Because precipitation is limited in the summer and the snowpack has largely melted, groundwater discharge to streams likely accounted for a larger proportion of stream discharge as the contribution of surface runoff decreased from spring through summer. As a result, covariates related to snowmelt such as SWEws were significant in May and June, but not in July and August water temperature and flow permanence models. Unlike surface runoff, which is broadly distributed across the landscape, groundwater discharge may be very localized within stream networks.
In our study of stream networks, the effect of progressive fragmentation on water temperatures was apparent across seasons (Figures 7 and 8). The spatial autocorrelation among flow-connected sites, indicated by increasing semivariance at greater separation distances, was present between May and August ( Figure 4). In contrast, spatial autocorrelation among flow-unconnected sites was present at distances of less than 15 km in May but became less pronounced in July and August. This suggests that watershed-scale processes and flow connectivity may be relatively more important drivers of stream temperature in May and June than in July and August. This seasonal change in the relative roles of watershed versus local processes was apparent not only in the changes in spatial autocorrelation depicted in the Torgegrams, but also in the significance of covariates and the relative predictive power of stream temperature models.
We hypothesize that, as stream networks become increasingly fragmented, coarse-scale covariates, such as those used in our study, may not adequately capture the localized processes contributing to observed temperature patterns [21,90] because drying of surface flows changes the hydrologic connections between points in the stream network. Although hyporheic connections may exist, streamwise autocovariance functions may change when surface flows are fragmented. As a result, the LOOCV r 2 pred of stream temperature models was greatest in May (0.89) and progressively decreased through August (0.65). Concomitant with this decrease in predictive power was a decrease in the percentage of variance explained by autocovariance structures from May to August. For streamflow permanence, the inclusion of spatial structures only marginally increased the predictive performance of models ( Table 4). The probability of flow presence decreased throughout the network in both the WW and WR watershed from May to August, but a greater proportion of the streams in the WW watershed remained wet from May to August. Paradoxically, the snowpack in the WW watershed was consistently less than in the WR watershed, but because runoff from snowmelt was largely depleted by June, other processes (e.g., groundwater flux) sustained late summer streamflows.
Although much progress has been made toward understanding terrestrial responses to changing climates, much less is known about water quality and availability in streams, which are important not only to fish, but also to a host of wildlife species and other valued ecosystem functions [21,37,115]. This is specifically true of the Great Basin, an arid region where surface water is already chronically scarce [44]. As climate-related drought becomes more likely in the region [116], stream networks will be fragmented and water scarcity will increase in space and time [24,117,118].
In our study, we predicted streamflow permanence using widely applied statistical methods with relatively simple landscape covariates such as NDVI. However, it may be difficult to attribute NDVI to specific, localized processes that influence water availability or stream temperature. Nevertheless, NDVI is an important indicator of changes in land management or influences of wildfires and can be used to track the effectiveness of land and water management actions in a changing climate [56,57,80,107,108]. Additional predictors that are similar to NDVI but finer in spatial resolution are needed to explain the high degree of spatial variability in stream temperature and flow presence that we observed. For example, covariates that characterize localized heat budgets based on shade and groundwater/surface water interactions are needed for better predictive models. Moving from publicly available spatial data sources to high-resolution remotely sensed data (e.g., airborne thermal infrared remote sensing, lidar, and aerial photography) and field measurements may be costly and logistically challenging. However, such efforts may be warranted because broad-scale models of stream temperature may not be sufficient for identifying local patterns of warming and multiscale thermal refuges that are critical to many species in fragmented, increasingly droughtprone landscapes [21,28].

Conclusions
This study provides novel insights into the fragmentation of stream networks with respect to patterns of stream temperature and flow permanence. Although SSNs have been used to predict stream temperature in multiple studies, we are unaware of previous efforts to use an SSN logistic regression model to predict flow permanence for a comparison to contemporaneous stream temperature. The high-density water temperature and flow presence datasets used to develop the SSN models presented within this study enabled us to consider the potential drivers of stream temperature and flow permanence at finer spatial scales than has previously been investigated. This study design made it possible to (1) explain fragmentation in terms of drivers that operate at different spatial and temporal scales and (2) evaluate a seasonal shift in processes that drive water temperature and flow permanence-from snowmelt in the late spring to riparian shade and groundwater discharge in the summer. Quantifying the seasonal and interannual variability in stream temperature and streamflow permanence is critical to understanding climatic variability and vulnerability, especially in the context of drought, which is a concern for aquatic ecosystems globally [17] and locally in the Great Basin [44]. Spatial scaling is also important because coarse-grained models of climate exposure (e.g., stream temperature and flow presence) may not reveal granular details, such as localized refuges that can be critical for species survival during stressful periods [21,119,120]. In streams, it is well known that such refuges can be common [25,77,94]. Accordingly, although fragmentation is often described as a threat to species persistence in stream networks [1][2][3], thermal heterogeneity can also create a patchwork of refuges that can carry species through stressful exposures. Our findings highlight the need for catchment-wide, high-density observations of water temperature and flow permanence to quantify thermal heterogeneity and predict the effects of finerscale covariates such as riparian canopy condition, channel morphology, and groundwater/surface water interactions. Modeling these finer-scale drivers of stream temperature and flow is important to address scales where land managers can effect change [26] and scales that are important to species vulnerability to hydrologic change [21].