Study area

The Triangle region of the North Carolina Piedmont is framed by the cities of Raleigh, Durham, and Chapel Hill. Historically, the region was largely agricultural with a few industrial cities. After broad-scale abandonment of agriculture in the 1930s, the area underwent widespread reforestation to pines, which are now succeeding to hardwoods (Kirby 2006). The area is undergoing explosive population growth and now has a range of development intensity from heavily urbanized Raleigh in the east to more agricultural areas in the west (Fig. 1). State parks and large tracts of lands owned and protected by institutions provide remnant areas of natural vegetation with minimal impacts from recent development. The area is a macrocosm in which to study urbanization and its various trajectories and serves as a model landscape for other regions experiencing similar patterns of development.

Fig.1. 

Map of land cover across study area with study watersheds delineated. County names are shown in all capitals.

Site selection

We conducted a synoptic survey of ∼70 low-order streams distributed throughout the area and across a range of land covers (Figs 1, 2A–D) as part of a larger study of the urban stream syndrome. We compiled a list of potential sites that included streams previously monitored by the US Geological Survey (USGS) as part of the Effects of Urbanization on Stream Ecosystems (EUSE) program, North Carolina's Department of Water Quality (DWQ), and Durham Storm Water (DSW), a series of ongoing research sites (Sudduth et al. 2011, Violin et al. 2011), and sites that had not been previously studied. We chose sites on the basis of land cover and ease of access. We worked with a time series of land-cover data created by Sexton et al. (in press), who used signature-extension methods to develop a classification scheme based on the 2001 National Land Cover Dataset (Homer et al. 2004) that could be extended with the same thematic resolution and accuracy to the entire Landsat Thematic Mapper archive. We extracted land-cover data (raster images at 30-m resolution) for 5-y increments (1985–2005). To select study sites, we analyzed 2005 land-cover data with moving-window averages to approximate proportions of development, agriculture, and forest in each candidate watershed to ensure that our sites would represent the full range of land-cover variation in the region. We used a geographic information system (GIS) (ArcGIS, version 9.3; Environmental Systems Research Institute, Redlands, California) to search for areas where streams crossed roads or ran beside them and used 2008 aerial photographs (NAIP 2008) to air-truth these sites. We eliminated sites that appeared to be ponds or lakes rather than streams, did not appear to have any body of water in the area, or appeared inaccessible.

Fig.2. 

Comparison of Goose Creek (83% development, 5% forest) and Stony Creek (7% development, 73% forest) showing the intensity of stormwater piping and stream burial that occurs in developed watersheds (A), land cover from aerial photographs (air photos) (B), watershed-scale thermal regime, as shown by skin temperature from satellite data (C), and ground-level photographs of the study reaches (D).

We visited 118 potential sites and eliminated those that were on private property, extremely difficult to access, too deep to sample in chest waders, or so small that they were likely to dry early in the season. We deployed samplers at 74 sites on 19 May 2009 for 30 d. Each sampler consisted of temperature loggers anchored to large concrete blocks that were secured by rebar inserted into the stream bed. We stationed all samplers in runs or pools, so that they would remain covered by water as long as possible in dry conditions. Some samplers were lost or buried by floods, and some data loggers malfunctioned, leaving a final data set consisting of detailed channel measurements, watershed land-cover analysis, and continuous temperature data at 60 sites.

Temperature data

We used Onset HOBO® Temperature/Alarm (waterproof) Pendant® Data Loggers (UA-001-08; Onset Computer Corporation, Pocasset, Massachusetts) in the samplers to record temperature every 10 min from 20 May to 10 June. These loggers are accurate to ±0.54°C, so we will not describe trends and differences of a magnitude <1.08°C. We focused on 2 biologically relevant temperature variables: 1) maximum temperature surge during stormflow and 2) maximum temperature during baseflow. We began by isolating a 24-h period surrounding a major storm (1100 h 28 May–1100 h 29 May) and calculating the maximum temperature change that occurred over 10 min during stormflow. This change indicates the sudden and severe thermal surges that can occur in urban streams. To analyze differences in baseflow temperature among streams, we focused on the week surrounding this storm (24 May–1 June) to ensure that baseflow and storm-flow temperatures were comparable and to minimize the likelihood of major disturbances caused by sediment burial or low-to-nonexistent baseflow. We calculated mean daily minimum, maximum, and mean temperatures during baseflow, mean observed change over any 24-h period (diel amplitude), and mean degree-days using the double-triangle method and a base temperature of 0°C (Sevacherian et al. 1977).

Habitat measurements

We measured canopy cover, stream channel depth, width, and incision, and streambed habitat diversity in a 100-m reach upstream of 42 of the 60 samplers between June and August 2009 and at the remaining 18 sites between June and August 2010 (Table 1). We counted the number of habitat transitions between pools, riffles, and runs in this 100-m reach. Reach-level description of habitat transitions may reflect watershed urbanization because streams in more urbanized environments typically have fewer habitat transitions consequent to scouring and incision caused by storms (Walsh et al. 2005b, Violin et al. 2011). We measured the wetted width and the depth in the center of the channel at each transition and every 10 m and calculated the minimum, maximum, and mean wetted width and depth for the reach. We surveyed cross-sections at 3 randomly chosen transects in each reach. We used a string level and measuring tape across the top of opposite banks to record bank-to-bank width and measured the height between the tape and a minimum of 7 points: top of left bank, bottom of left bank, left edge of water, thalweg, right edge of water, right bottom of bank, and right top of bank. We added points as needed to account for sandbars and other anomalies. We calculated the incision of the stream as channel depth at thalweg divided by bankfull width at the 3 cross-sections. We measured canopy closure from the ground at the thalweg of the stream with concave forest densiometers (Forest Densiometers, Bartlesville, Oklahoma) (Lemmon 1957) at each cross-section to provide an estimate of canopy closure 100 m upstream of the temperature logger. Densiometer readings are subjective, so field canopy-closure measurements were made by 2 technicians who underwent extensive calibration to ensure their interpretations were consistent. We also calculated canopy closure from aerial photographs taken in 2008 (NAIP 2008). We created a 10 × 10-m grid in ArcGIS that covered the entire study area, overlaid satellite images of each stream reach with the grid, and visually analyzed cover 100 m upstream of the temperature loggers. We counted the grid cells in which the stream was not visible and divided this number by 10. For example, if the stream was clearly visible in 2 of the grid cells, we estimated canopy closure as 80%. These estimates were only slightly correlated with densiometer readings (adjusted R²  =  0.25, p < 0.05), probably because of differences in sampling and photography dates and in resolution between densiometer readings and 30-m grids.

Table 1. 

Habitat- and watershed-scale variables used in our study. Variables with significant relationships with thermal variables and those used in the structural equation model (SEM) are noted. *  =  not significant, **  =  significant, but not selected for model, ***  =  used in SEM, Min  =  minimum, Max  =  maximum, N/A  =  not applicable.

We were unable to measure habitat at 18 sites in 2009, so we revisited these sites in 2010 and used the same procedures to measure habitat indices to increase our sample size of sites with temperature, habitat, and landscape data to 60. We also revisited 13 sites (for which we had 2009 habitat data) that spanned the range of urban development in our study region to test whether habitat variables differed between years. Only variables related to depth differed significantly between years (paired t-tests), primarily because of differences in weather between years. Therefore, we did not use minimum, maximum, and mean depth measurements in the model.

Land cover

We used ArcGIS tools to delineate the watershed upstream of each sampler based on 30-m light detection and ranging (LiDAR)-derived digital elevation model (DEM) data (NCDOT 2007). We imported global positioning system (GPS) coordinates for each sampler and calculated the elevation (range 60–197 m asl) and orientation of each study reach because these variables affect light availability in rivers (Julian et al. 2008). We used a flow-accumulation raster layer to hand-edit the sites by shifting their locations into the nearest cells of highest flow accumulation. This manipulation created a set of pour points from which the GIS could delineate the watershed for the area upstream of each sampler. We checked each watershed with the DEM and aerial photographs and rectified obvious errors.

For each watershed, we computed a set of indices to represent land cover (Table 1). We calculated % developed land in 1985, 1995, and 2005 and % forest and % field in 2005. We addressed potential differences between older development and newer development, such as differing age and types of stormwater infrastructure and best management practices (Kaushal and Belt 2012) by calculating the proportion of development observed in 2005 that occurred between 1985 and 2005 and between 1995 and 2005. We did not use the amount of watershed classified as “field” in our analyses because this classification is difficult to assess. In agricultural areas, field indicated crops and pasture, but in urban areas, field represented golf courses, lawns, or cemeteries.

We calculated several variables that incorporated the effects of land cover on streams. We began by assuming that land-cover effects were inversely proportional to their distance from the stream (King et al. 2005). This approach took into account surface connectivity between the stream and watershed, but did not account for the effects of subsurface connectivity created by stormwater infrastructure. Many of the watersheds in our study were outside of municipal boundaries in areas where storm-water infrastructure or information on the location of storm-water infrastructure do not exist. Attempts to estimate these connections would have been extremely haphazard. We chose to estimate connections between development and streams via natural flow paths because storm-water infrastructure tends to follow these flow paths. Thus, our methods provide a first-cut understanding of the effects of land cover in urban watersheds.

We used developed land cover as a proxy to calculate effective development for each watershed. We used the DEM to calculate the distance along a hydrologic flow path from each grid cell with developed land cover to the nearest grid cell containing the stream. This distance was an estimate of the actual distance to the stream along the flow path. We weighted this distance by the relative infiltration capacity of land covers encountered along the flow path based on the curve number method (USDASCS 1986). We calculated curve numbers based on hydrologic soil groups from SSURGO soils data (Soil Survey Staff 2010) and land cover (natural vegetation has high infiltration rates, developed land has low infiltration rates). We used calculated curve numbers to estimate effective resistance to infiltration. Thus, water flowing over land with natural vegetation with high infiltration capacity effectively traversed a longer distance than water flowing over developed land because it had more time to infiltrate. In contrast, impervious surfaces led to shorter effective distances along the flow path. We expected sites farther from the stream along these flow paths to have less effect on the stream than nearby sites, so we modeled the relationship between the distance from the stream and the effect on the stream as a negative exponential decay function. The rate at which this decay should occur and the distance at which the effect of a cell on a stream should be effectively null were not known. Therefore, we selected a range of distances (50–2000 m) over which we could expect the effect of the cell on the stream to become arbitrarily small (1%) and calculated an exponential decay rate based on each distance. We computed infiltration-weighted distances from the stream along surface flow paths for each cell with developed land cover to produce a range of effective development variables with different decay rates. In this approach, a grid cell with developed land cover was down-weighted if it was far from the stream, not in the flowpath, or if water would pass over cells with land-cover with high infiltration capacity before it reached the stream. We calculated the average of these values in each watershed. To calculate effective weighted development, we weighted the land-cover classes, low-density developed, medium-density developed, and high-density developed (corresponding to NLCD classes 22, 23, and 24, respectively) by relative impervious indices of 33, 66, and 100% (Homer et al. 2004) and used the procedures described above for effective development.

We calculated effective forest as for effective development but with forested land cover. The result was an estimate of buffering land cover, weighted by hydrologic proximity to the stream. We expected effective forest to be most important over short distances (e.g., 50–100 m), whereas effective development would be important at longer distances (≥1000 m). Because we computed both indices over a range of distances, we were able to assess these scaling relations directly in exploratory data analysis.

We calculated % connected development as the percentage of the watershed with developed land cover that was closely connected to a stream (based on distance to stream). We calculated % developed land within 100 m of a stream and divided this area by the total watershed area, so that this variable was comparable between watersheds. Dividing by total area in the watershed also placed % connected development in the context of the larger watershed, while weighting it higher than land cover further from the stream.

We indexed the effects of roads and traffic volume on streams by methods described in previous studies of the relationship between roads, traffic volumes, and runoff. We began by computing the Euclidean distance from each raster grid cell in each study watershed to the nearest road. We assumed that road effects decreased rapidly with increasing distance. We log(x)-transformed these distances with a distance–decay constant based on the results of a study of the dissipation of heavy metals with increasing distance from roads (Lygren et al. 1984). We weighted this distance–decay relationship based on estimates of traffic volume provided by North Carolina Department of Transportation (NCDOT 2007). Inverse-distance-weighted traffic volume estimated road effects mediated by traffic and is correlated with road size (Mayer et al. 2011). For example, the heavy metals deposited by car exhaust are directly proportional to traffic volume on nearby roads. Roads that are more heavily trafficked are generally larger. We also computed road density in each watershed to gauge whether cumulative effects might be more important than the influence of the nearest roads (Kratzer et al. 2006). We summed road effects in the watershed as total m of road/ha of watershed. We calculated the mean number of road/stream intersections in each watershed/km of stream to estimate the effect of routing storm water from roads directly into streams (Jones et al. 2000).

Variable selection

Our goal was to select variables that would best represent paths by which urbanization might affect baseflow and storm-flow stream water temperatures (Fig. 3). We used exploratory data analysis to select predictor variables and to assess possible interaction pathways. We classified the variables as reach or catchment scale and grouped them within scales into categories that might describe mechanisms affecting stream temperature (Fig. 3). Catchment-scale variables described forest, road, and development land cover, and reach-scale variables described physical habitat or light. We used simple correlations to test for relationships between each habitat and land-cover variable and the 2 thermal variables, maximum temperature at baseflow and maximum temperature surge during stormflow. Reach orientation was the only categorical variable, so we used a t-test to compare each thermal variable between north–south oriented streams and east–west oriented streams and found no difference between orientations. Except for orientation, we retained the variable in each group with the strongest simple correlations with the response variables to reduce the potential for confounding. For example, we retained a single road-effect variable (road density) and a single light-effect variable (canopy closure from aerial photographs). We also used these correlations to select the most appropriate decay distances for effective development and effective forest. We selected effective forest calculated with a decay distance of 50 m because it was most strongly correlated with the temperature variables. Percent developed area in 2005 was more strongly correlated with the temperature variables than effective development, so we removed effective development.

Fig.3. 

Conceptual metamodel (Grace et al. 2010) showing influence of watershed- and reach-level variables on stream temperature. Each box includes a list of potential variables to describe the given category. This diagram reflects our understanding of the system and was used to select variables to include in the structural equation model (SEM). Arrows show direction of effects. See Table 1 and Methods for descriptions of variables.

We used regression tree analysis (Breiman et al. 1984) to explore the possibility that response variables might be affected by alternative causal pathways expressed in different branches of a regression tree (Urban 2002). This procedure also allowed us to screen for confounding variables that could be highly correlated with both explanatory and response variables in the model. We included these relationships in structural equation models (SEMs) via indirect pathways and intermediate variables, which take into account confounding relationships. For example, the final SEM showed that canopy closure was an intermediate variable between development and maximum temperature, revealing that canopy closure affected maximum temperature but was affected by development. We used the reduced set of predictors and the pathways suggested by regression trees and our knowledge of urban streams to develop an SEM.

SEM

We used SEM (AMOS, version 16.0; IBM SPSS, Chicago, Illinois; Grace 2006, Grace et al. 2010) to evaluate hypotheses about the mechanisms connecting predictors and responses of interest. SEM can be used in a variety of ways (Grace et al. 2012), but our focus was on evaluating and discovering network connections in the system. The conceptual metamodel used to explore the influences of catchment- and reach-level alterations on temperature changes in urban streams (Fig. 3) represented many possible SEMs that could be evaluated with our data. These models included both watershed- (road density, % development) and local-scale variables (stream width, canopy closure, distance-weighted forested area). This conceptual model also suggested possible indirect pathways. For example, % development could affect maximum temperature directly or by association with reduced canopy closure, which also affects maximum temperature. We developed SEMs, evaluated multivariate expectations against our data, and made necessary modifications to pose alternative models. Thus, our use of SEM can be considered a model-building application (Jöreskog 1982).

In model fitting, we used goodness-of-fit variables for the overall model and tests of significance for individual path coefficients to revise the initial model sequentially. We considered removing links in the model if the data suggested they lacked explanatory power. This pruning continued until all path coefficients retained in the model were significant (p < 0.05) or our biophysical understanding of the system required that the variable be retained (Grace 2006). At the same time, we assessed overall model fit after each iteration with the goal of achieving a model that produced results that did not differ from the data (i.e., model–data consistency). We assessed potential violations of the stringent assumptions of maximum-likelihood estimation by fitting this model using the Bayesian estimation procedures (Arbuckle 2007). The methods produced nearly identical standardized regression weights, indicating that maximum-likelihood-estimation methods produced a robust solution.

Longitudinal survey

We conducted intensive temperature mapping in 4 of the 60 streams in August 2009 to document spatial variability in baseflow temperature. These streams had a range of development and % connected development. Portions of the study reaches of all 3 of the urban streams had been recently restored. The air temperatures during this study were some of the hottest in the year, providing an optimal period in which to examine maximum baseflow temperature and to observe extremes. We used a YSI Model 30 handheld conductivity and temperature probe (Yellow Springs Instruments, Yellow Springs, Ohio) to measure the temperature every 50 m above and below any pipe outfalls or tributaries. We walked upstream and continued until we had data for a ≥1-km-long section of stream. We calculated summary statistics for each stream and used t-tests to compare stream temperatures in shaded vs unshaded locations. We used these analyses to corroborate conclusions drawn from the synoptic survey data.

We did all analyses except SEMs in R (version 2.10.1; R Development Core Team, Vienna, Austria) and the packages tree (Ripley 2009) and ecodist (Goslee and Urban 2007).

Watershed land cover

Watersheds had from 0% to 99% developed land (Table 1, Fig. 4A). Some watersheds were almost entirely developed (99% in 2005), whereas others were nearly completely forested (96% in 2005) (Table 1, Fig. 4B). Other aspects of development varied greatly among watersheds (e.g., % fields; Fig. 4C). Road density ranged from 0 to 60 m/ha (Table 1, Fig. 4D). Even among the 10 most-developed watersheds, road density ranged from 12 to 60 m/ha. The number of road/stream intersections varied from 0 to 4 km of stream (Table 1, Fig. 4E), and % connected development ranged from 0 to 100% (Table 1, Fig. 4F).

Fig.4. 

Histograms showing the distributions of % development (A), % forest (B), % field (C), road density (D), road–stream intersection density (E), and % connected development (F) across watersheds. See Table 1 and Methods for descriptions of variables.

Inverse-distance-weighted land-cover variables

Only 2 inverse-distance-weighted land-cover variables were related to the temperature variables. Effective weighted development had the 4^th-strongest correlation with minimum temperature. Percent development was the best descriptor of temperature variables in the SEM. Effective forest performed better than % forest and was included in the SEM. These variables describe the surface connectivity to streams of both development and forest (King et al. 2005).

Stream thermal regimes

The wide range in development was accompanied by wide variation in stream thermal regimes (Fig. 5). Baseflow temperatures ranged from an absolute minimum of 12.4°C to an absolute maximum of 33.2°C, with an average across all sites of 17.2°C (Table 2). Storm-flow changes ranged from decreases of 3°C to increases of 4.2°C, but mean decreases and increases in temperature were <1°C across all sites. Temperature variables describing minimum, mean, and maximum temperatures, maximum temperature surge during stormflow, and degree-days were all significantly higher in the 10 most-developed watersheds than in the 10 most-forested watersheds (1-tailed unpaired t-test, α  =  0.05). The only thermal variable not significantly different between the 2 groups was the mean diel change in temperature.

Fig.5. 

Temperature data from synoptic survey. Lines are coded by the primary land-cover category in the watershed. Temperature accuracy is ±0.54°C.

Table 2. 

Mean minimum (min) and maximum (max) values for thermal variables for all streams, 10 most-developed streams, and 10 most-forested streams. Average min and max temperature (temp), mean temp, mean diel range, and degree-days describe the period from May 24 to June 1. Max temp increases and decreases were calculated over 10 min during a 24-h period surrounding a storm. * indicates changes of a magnitude <1.08°C, reflecting the instrument loggers' temperature accuracy of ± 0.54°C.

Correlates of stream baseflow temperatures

Baseflow minimum, maximum, and mean temperatures were each strongly correlated with different explanatory factors (Table 3). At the reach scale, warmer baseflow temperatures were associated with wider streams, whereas cooler baseflow temperatures were associated with greater riparian canopy cover and more habitat transitions. At the watershed scale, warmer temperatures occurred in watersheds with higher road densities, and cooler temperatures occurred in watersheds with higher % forest. Reach-scale factors were often stronger correlates with baseflow temperature than were watershed land-cover attributes. Many of these watershed and reach-scale factors were correlated. For example, % development was highly correlated with road density (r  =  0.51) and with canopy closure from aerial photographs (r  =  −0.47). Other factors had opposing effects on temperature, e.g., % development and % forest.

Table 3. 

Correlations between thermal variables and landscape and habitat variables of synoptic survey sites, ordered by correlation with mean temperature (temp). Maximum positive and negative correlations for both reach- and catchment-scales are bolded. Minimum (min) and maximum (max) temp describe the period from May 24 to June 1. Max temp surge over 10 min was calculated during a 24-h period surrounding a precipitation event. Landscape variables were calculated for the watershed of the study site. Habitat variables were calculated for a 100-m reach above the samplers at each site. See Table 1 for a description of variables.

The final SEM included 5 explanatory variables and explained 47% of the variation in baseflow maximum temperatures across streams (χ²_5df  =  1.655, p  =  0.976; Fig. 6). Maximum temperature appeared to be most strongly influenced by canopy closure via a direct negative path and by mean width of the channel by a direct positive path. Percent developed land cover and road density significantly influenced maximum temperature via a direct path. A significant indirect path indicated that shading effects of % canopy cover dampened the positive relationship between % developed land and stream baseflow temperatures.

Fig.6. 

Final thermal structural equation model (SEM) showing standardized regression weights. The SEM was fit simultaneously with 2 focal response variables, maximum temperature during baseflow and maximum temperature surge during stormflow. Arrows show direction of effects. Arrow line weight indicates strength of the path. See Table 1 and Methods for descriptions of variables.

Baseflow temperature varied by as much as 10°C along a 1-km stretch of urban Goose Creek and as little as 2°C along a 1-km stretch of forested Stony Creek (Table 4). Temperature was spatially variable in the 2 narrower urban streams (Goose Creek and Rocky Branch) but was more uniform (and warmer) in a wider urban stream (North Gate) (Fig. 7A–C). The highest temperatures in the narrower urban streams were associated with canopy gaps (1-tailed unpaired t-test, α  =  0.05). Much of the spatial heterogeneity in Goose Creek and Rocky Branch was a result of variation in riparian shading. Temperatures were significantly lower in sections of the stream with closed canopies than in sections with canopy gaps (1-tailed unpaired t-test, α  =  0.05). Stormwater outlets had inconsistent effects on baseflow temperature. The highest temperature (29.5°C) observed during the longitudinal survey was in a well shaded section of Rocky Branch immediately below a stormwater outlet, whereas the coldest temperature (22.7°C) was in an unshaded section of Goose Creek below piped tributary.

Fig.7. 

Longitudinal thermal profiles of 3 urban streams: Goose Creek (A), Rocky Branch (B), and North Gate (C) compared to a forested stream (Stony Creek). Distance at 0 m represents the upstream beginning of the study reach and the measurements move from upstream to downstream along the x-axis. Stony Creek has no outfalls in the study reach.

Table 4. 

Results of longitudinal survey of 4 stream reaches. Mean air temperature is air temperature over the time of the survey. The number of observations (n) and total distance surveyed differs among streams. Min  =  minimum, max  =  maximum, temp  =  temperature, SD  =  standard deviation, CV  =  coefficient of variation.

Controls of storm-flow temperature surges

Thermal responses of streams to storms were related to watershed land cover but also varied greatly among streams in watersheds with similar land cover (Fig. 8). Maximum temperature surge during stormflow was significantly greater in the 10 most-urban streams than in the 10 most-forested streams (1-tailed unpaired t-test, α  =  0.05). Urban stream temperatures increased intensely and suddenly during storms, whereas forested stream temperatures responded gradually to stormflow and often showed only a small increase or a cooling effect in response to summer storms. Maximum temperature surge during stormflow was more strongly correlated with catchment- than with reach-scale variables. Percent development was positively correlated and % forest was negatively correlated with maximum temperature surge during stormflow (Table 3). Maximum temperature surge during stormflow was greater in reaches with greater channel incision and less canopy cover, factors that were highly correlated with % development.

Fig.8. 

Thermal responses of different streams to the same storm. Lines are coded by the primary land-cover category in the watershed. Bar graphs show precipitation, measured on the alternate y-axis. Temperature accuracy is ±0.54°C.

The final SEM explained 44.8% of the variance in maximum storm-flow temperature surge (Fig. 6). The strongest effect was the positive direct path between % development and the maximum temperature surge during stormflow. Maximum temperature surge during stormflow also was influenced by stream width via a direct negative path. The model also identified a significant indirect path by which road density and % development negated some of the cooling effect of % canopy closure on maximum temperature surge during stormflow.

How much difference do a few degrees make?

The hottest temperature observed in our survey was 28.5°C, a temperature above the critical thermal maxima for Salmonidae and Cottidae, among other families (Beitinger et al. 2000). For most North American freshwater fish species, the critical thermal maxima fall within the range of 32 to 40°C (Beitinger et al. 2000). In most cases, even during summer thunderstorms, our urban streams did not enter this range. Thus, we found little evidence to suggest that urban stream heating will lead to fish mortality. However, temperatures that exceed the thermal optima of aquatic organisms can lead to slower rates of development or growth (reviewed by Hester and Doyle 2011). The higher temperatures in the urban streams in our study have the potential to be stressful for in-stream biota and are likely to exacerbate and extend zones of hypoxia and anoxia in benthic sediments and to speed biogeochemical reactions. Quantifying the biological effects of the baseflow temperatures observed in our study streams is difficult because so little is known about this topic.

Thermal surge during stormflow

The effect of temperature surges as high as 4°C on stream biota are poorly understood (Nelson and Palmer 2007). In laboratory studies of thermal tolerance, temperature changes typically are gradual. Thus, results of these studies may have little relationship to outcomes of the temperature surges in urban watersheds. In urban areas, storms result in high flows, intense scouring, and large inputs of pollutants and nutrients in streams. Thus, the effects of heat surges are entangled with a large suite of changes that occur during storms. The extent to which heat shock might contribute to the biological degradation that accompanies flashy urban-stream flow regimes is unclear. However, the sudden and intense temperature surges observed in our study conclusively demonstrate rapid conveyance of heat from impervious surfaces into streams and are an effective indicator of the degree of connectivity between impervious surfaces (and their associated contaminants) and nearby streams. The rapidity of the response highlights the role of storm-water infrastructure in connecting the entire watershed during storm events. Our results are a convincing demonstration of the capacity of the watershed to absorb and dissipate urban heat. Even our most extreme observations of storm-flow temperature increases of 4°C in 10 min involved only a small fraction of the heat stored watershed-wide in impervious cover.

Management implications

Surface-connected development, described here by effective development, was an ineffective predictor of thermal regimes in streams. This result emphasizes the importance of including storm-water infrastructure and % developed land in the entire watershed in management practices, rather than focusing on the state of the riparian zone (see Walsh et al. 2005a and Bernhardt et al. 2008 for further discussion). Baseflow temperatures strongly influenced the magnitude of the temperature surge during stormflow. This result also has important management implications. Local managers could effectively maintain cooler baseflow by increasing canopy cover in riparian areas and by altering stream width, but these efforts are unlikely to be effective mitigators of storm-flow heat surges. Riparian cover can reduce baseflow temperatures, but urban streams with high canopy cover may experience much larger temperature increases during storms than those with low canopy cover because large volumes of hot storm water enter the stream during stormflow. Thus, the difference between baseflow and storm-flow temperature is likely to be greater in a cooler than in a warmer stream. These findings reinforce the importance of considering the entire catchment and its effects when designing local conservation and restoration projects and when considering how changes in local variables will interact with largely unchanged landscape variables.

Future research

Documentation of temporal and spatial variations in temperature of urban streams is an important first step in understanding how urbanization influences the thermal regime of streams. Our study provides a detailed starting place for research designed to identify the extent of thermal changes caused by urbanization in warm-water streams in other areas. Our descriptions of thermal changes in urban streams and the variables that best explain them are likely to hold true for warm-water streams in urban areas across the world. However, the strength of the relationships is likely to change depending on the specific landscape and the local context of the site. Our work also provides a basis for using temperature as a tracer of urban effects on streams, especially for contaminants and nutrients likely to enter the stream during stormflow.

More research is needed to explore the thermal variability observed in watersheds with moderate amounts of development. For example, we observed an almost 5°C difference between the minimum and maximum mean temperature of streams across the landscape. The thermal variability in moderately developed watersheds was not well explained by the model and many mechanisms were left unclear. Within this moderate range of development, storm-water management, specifically subsurface infrastructure and best management practices, are likely to have intense effects on the thermal regimes of streams. Researchers need more accurate descriptions of locations and conditions of stormwater pipes. We suspect that the lack of explanatory power of our hydrologically nuanced indices of surface connectivity reflect the reality that these effects are overridden by subsurface connectivity in urban watersheds. Incorporation of fine-scale effects of canopy cover and imperviousness into models also is essential for understanding urban stream ecosystems. Last, more in-depth longitudinal studies with increased spatial and temporal resolution will increase understanding of the influences of variables like outfalls and canopy openings on baseflow and stormflow and will enable researchers to better understand controls on the thermal regimes of urban streams.