Headwater stream temperature response to forest harvesting: Do lower flows cause greater warming?

This study addressed two hypotheses regarding the relationship between stream temperature response to shade removal and streamflow: (a) that temperature response increases as flow declines and (b) that the relationship can be complicated by shifts in dominant streamflow sources and pathways during low‐flow periods. The study was based on a paired‐catchment design in rain‐dominated headwater catchments in the southern Coast Mountains of British Columbia, Canada, and focused on the effect of a 50% thinning treatment on daily maximum temperature in June, July and August for three sites within the harvest treatment. At the two upstream sites, the treatment response exhibited a negative relationship with daily mean streamflow, especially for days with high incident solar radiation. This result suggests that the effectiveness of forest practice rules for protecting cold‐water habitat may be reduced under future climatic conditions characterized by more frequent extended drought. However, stream temperature response at the most downstream site exhibited a pronounced inverted‐U‐shaped relationship with streamflow measured at a weir. It is hypothesized that the response at the most downstream temperature logger was controlled by the existence of a stable, relatively cool inflow just upstream of the logger, which represented an increasing fraction of flow as streamflow generated higher up in the catchment declined through time. There was a lack of convergent surface topography upstream of the logger, and it is hypothesized that localized inflow may have been controlled by the topography of the soil‐till interface and/or originated as hyporheic discharge.


| INTRODUCTION
Water temperature is an important control on many biogeochemical, biological and ecological processes in aquatic systems (Webb et al., 2008).Decades of research in a range of forest ecosystems have demonstrated that reductions in riparian shade by logging, vegetation management or natural disturbance result in summertime warming by increasing the amount of solar radiation reaching the stream (Bladon et al., 2018;Brown, 1969;Isaak et al., 2010;Leach et al., 2022;Leach & Moore, 2010;Lynch et al., 1984;Raulerson et al., 2020).Many studies have quantified stream temperature response to forestry using paired-catchment studies or other forms of before-after/control-impact (BACI) designs (e.g., Bladon et al., 2016;Gravelle & Link, 2007;Groom, Dent, Madsen, & Fleuret, 2011;Harris, 1977;Janisch et al., 2012;Johnson & Jones, 2000;Macdonald Disclaimer: Use of manufacturers' names and models is for information purposes only, and should not be interpreted as an endorsement of the measurement systems used in this study.et al., 2003;Mellina et al., 2002;Rayne et al., 2008;Rex et al., 2012;Swartz et al., 2020), and some have documented changes in response magnitude over several years post-harvest (e.g., Gomi et al., 2006;Harris, 1977;Johnson & Jones, 2000;Leach et al., 2022).However, relatively few have focused on understanding the drivers of withinseason variability at daily or finer time scales.Notable exceptions include Gomi et al. (2006), Guenther et al. (2014) and Oanh et al. (2021), who showed that the magnitude of the temperature response to shade removal tends to be greater for daily maximum temperatures compared to daily mean or minimum temperatures, and to vary from day to day in response to changing weather conditions.
Understanding the drivers of the temporal variability of temperature response to forestry is important when comparing studies, as in the meta-analysis by Martin et al. (2021) to evaluate the effectiveness of riparian buffers.A number of studies have included only one or two post-harvest years (e.g., Guenther et al., 2014;Rayne et al., 2008;Swartz et al., 2020), and the estimated temperature response to forestry may not be comparable across studies if the post-harvest periods at some sites included unusual hydroclimatic conditions.The relationship between temperature response to shade removal and stream discharge is also of major concern in a climate change context given that many regions, such as the Pacific Northwest of North America, are projected to have increasingly severe summer droughts, which would lead to lower late-summer stream discharge.
To the authors' knowledge, only one study quantified the relationship between stream temperature response to forestry and stream discharge (Oanh et al., 2021).That study found that the springsummer temperature response to shade removal was positively related to daytime solar radiation and negatively associated with stream discharge.They attributed the negative association with discharge to three factors: (1) discharge is positively associated with mean water depth, which reduces stream thermal response to energy inputs; (2) increased discharge is associated with higher velocity and thus shorter residence time within channel reaches affected by shade removal, which reduces the opportunity for heating; and (3) hillslope runoff during and following rain events is associated with advective heat transport, which tends to have a cooling effect during summer.
The results of Oanh et al. (2021) indicate that temperature response to reduced shade should increase with decreasing discharge.
However, the relationship may be complicated by shifts in streamflow sources with changing discharge.For example, Leach and Moore (2017) inferred on the basis of field observations and model simulations that, during extended recession periods, streamflow contributions in an unlogged forested headwater catchment became dominated by water following deeper and/or longer flow paths, which was cooler in summer than contributions originating from shallower and/or shorter flow paths.Other studies that have linked stream temperature dynamics to discharge-related shifts in hydrologic flow paths include Kobayashi et al. (1999) and Subehi et al. (2009).
The objective of this study was to explore the relationship

| Study site
The study was conducted in the University of British Columbia Malcolm Knapp Research Forest (MKRF), located approximately 60 km east of Vancouver, Canada (Figure 1).The area has a maritime climate with warm, dry summers and wet mild winters.Mean annual precipitation varies between 2000 and 2500 mm over the study catchments, with about 70% of the total annual precipitation falling between October and April.Snowfall accounts for approximately 15% of annual precipitation.Soils are permeable, shallow podzols typically about 0.5-1 m in depth, underlain by relatively impermeable compacted basal till or granitic bedrock (Haught & van Meerveld, 2011;Hutchinson & Moore, 2000).Much of the baseflow is generated by lateral discharge from hillslopes rather than from deeper groundwater.Prior to harvesting, forest cover was dominated by mature second growth western hemlock (Tsuga heterophylla), western red cedar (Thuja plicata), and This study was part of a broader project focused on the impacts of forest thinning on headwater stream ecology, which included three treatment catchments and three control catchments (Griffith & Kiffney, 2022).This study focused on Griffith Creek, where streamflow is gauged by a v-notch weir; the other two treatment catchments were ungauged.Griffith Creek is a first-order stream with a catchment area of 10.8 ha at the weir location.The headwater portion of Griffith Creek is an incised channel characterized by gradients of approximately 20%, which decrease in the downstream direction to less than 7% slope.
Of the three control streams, Mike Creek yielded the lowest standard error of the residuals for the pre-harvest period (Guenther et al., 2014), and was chosen as the control for the current study.
Mike Creek has an upstream catchment area of 29.7 ha, mean channel gradient of 8%, mean bankfull width of 1.5 m, and has a dominantly step-pool morphology (Gomi et al., 2006).Like Griffith Creek, Mike Creek's channel has a north-south orientation.
Monitoring began in July 2002 and extended through March 2007.Between September and November 2004, 50% of the basal area within the lower portion of the catchment was removed, with retention of the largest stems.The harvest unit covered the lower 280 m of the creek, with no riparian buffer.Further information regarding the study site and treatment can be found in Guenther et al. (2012), Guenther et al. (2014) and Leach and Moore (2014).

| Stream temperature monitoring and processing
Stream temperature was recorded at 10-min intervals using Onset Tidbit temperature loggers, which have a typical accuracy of 0.2 C. Three temperature loggers were located within the harvest unit along Griffith Creek.The most downstream logger was located just upstream of the weir pond near the downstream end of the harvest unit, denoted 0 m, and at distances of 100 and 200 m.Only one temperature logger was located in Mike Creek.Temperature loggers were installed in deeper portions of the channel to reduce the possibility of de-watering, and were positioned within a length of white PVC pipe to reduce heating by solar radiation.Multiple holes were drilled along the length of the pipe to enhance exchange of water between the pipe and the ambient stream water.
After downloading data from the temperature loggers, periods of de-watering or other anomalies were detected by manual inspection of temperature time series using interactive graphs generated by the plotly R package (Sievert, 2020).Based on this inspection, anomalous spikes near solar noon were observed occasionally at the 200-m site during August 2006.For that period, daily maximum temperature was determined based only on temperatures following 3 PM, given the observation that daily maxima at all sites typically occurred between 3 and 7 PM on days when anomalous spikes were not evident.
Following these pre-processing steps, time series of daily minimum, mean and maximum temperatures were extracted for further analysis.
Only daily maxima are considered here, given that they are most sensitive to the effects of shade removal (Gomi et al., 2006;Guenther et al., 2014;Oanh et al., 2021).

| Meteorological data
A weather station was set up in a clearcut approximately 2 km southwest of and 200 m lower in elevation than the Griffith Creek weir (Guenther et al., 2012).Of relevance to the current study, solar radiation was measured by a Kipp and Zonen CM-6B pyranometer that was scanned every 1 s by a Campbell Scientific CR10x data logger.
Rainfall was measured by a tipping bucket rain gauge with a tipping depth of 0.25 mm.The rain gauge was located at ground level; other sensors were mounted about 1.5 m above the ground surface.All observations were averaged (for solar radiation) or totalled (for rainfall) and stored every 10 min.

| Streamflow data
Streamflow was monitored by M. Feller (Faculty of Forestry, University of British Columbia) at a v-notch weir located at the downstream end of the harvest unit.The weir was constructed by excavating the soil/surficial material down to bedrock across the stream and pouring concrete into wooden forms, followed by fitting a sharp-edged metal "v" into the concrete structure (Figure 2).Stage was recorded using a KPSI Model 7 pressure transducer connected to a Unidata Starlogger model 6004D datalogger located in the weir pond.
Streamflow at the weir was measured using the volumetric method, with different sized buckets for different discharge levels.A total of 21 measurements was made between 2003 and 2006, with discharges ranging from 0.07 to 30.5 L s À1 .A rating curve was developed to convert the stage record at the weir into a streamflow time series based on the following model form: where Q i is discharge (L s À1 ) for measurement i; h i is the stage (as pressure measured by a transducer) for measurement i; a, b and h 0 are parameters to be estimated; and e i is the residual for measurement i.
The model was initially fit using generalized nonlinear least squares (GNLS) regression, using the gnls function in the nlme package in R (Pinheiro et al., 2023), with the error variance modelled as where var e i ð Þ is the variance for measurement i, b Q i is the predicted value for measurement i based on Equation (1), and t is a parameter to be estimated as part of the model fitting, which can account for heteroscedastic error variance.Prediction limits were computed by refitting the model using the nls function in the stats R package, with specified weights computed using the fitted value of t by GNLS, followed by application of the predFit function in the investr R package (Greenwell & Kabban, 2014).
Manual streamflow measurements were also made at distances of 80, 107, 160 and 185 m upstream of the weir during summer 2005 based on constant-rate injection of salt solution using a Mariotte bottle (Moore, 2004a(Moore, , 2004b)).The 80-m and 107-m locations, and the 160-and 185-m locations, represent the upper and lower boundaries of two intensively monitored reaches, respectively (Guenther, 2007;Guenther et al., 2014).For each monitoring reach, salt solution was injected a sufficient distance upstream of the upper boundary that lateral mixing was complete at both the upper and lower boundaries.The injection was maintained until steady state conditions were achieved (i.e., plateau at both boundaries).Calibration was conducted in the field following Moore (2004a).When lateral mixing is complete, as was the case at Griffith Creek, discharge measurements based on constant-rate salt injection can be accurate to within about AE5% (Moore, 2004a).

| Paired-catchment analysis
The analysis focused on June, July and August, the period when seasonal maximum temperatures typically occur and when treatment effects due to harvesting tend to be greatest due to the combination of high solar radiation and low stream discharge (Gomi et al., 2006;Guenther et al., 2014).To account for potential interannual variability of model coefficients, mixed-effects linear models using year as a random effect were fit for the pre-harvest period between stream temperature for each site along Griffith Creek and that in the control catchment.The models have the following form: where y t is the temperature at the treatment stream on day t, x t is the temperature at the control stream on day t; β 0 , β 1 , β 2 and β 3 are fixed-effects coefficients to be estimated; b 0i , b 1i , b 2i and b 3i are random effects for year i; j is day of the year (j = 1 on January 1); T = 365.25, the number of days in a year; and ϵ t is an error term, which was modelled as an autoregressive process of order 1.The sine and cosine terms in Equation ( 3) account for seasonal bias in the residuals (Gomi et al., 2006).The models were fit using the lme function in the nlme R package based on maximum likelihood.
The treatment effect at a given time in the post-harvesting period (T e t ð Þ) was estimated as: where y t and b y t are the observed and predicted temperatures on day t of the post-harvest period, where predicted temperatures were computed from the fixed-effects coefficients from the pre-harvest model.
There is no simple way to compute prediction limits for lme fits with autoregressive error structures.Prediction limits (y pl ) were approximated here by computing residuals for the pre-harvest period based on the fixed-effects coefficients, and computing limits as follows: where s r is the standard deviation of the residuals and n s is the number of standard deviations used to construct limits.For independent samples, n s = 2 would produce limits with approximately a 95% confidence level.Here, n s = 3 generated limits that included all but 2 of the observed temperatures, or 99.7%, for the pre-harvest period, and thus should provide a realistic indication of the precision of the predictions for the post-harvest period.

| Analysis of the treatment effect
The variability of the treatment effect was initially examined graphi- inverted-U shape, also called "hump shaped" (Lind & Mehlum, 2010).
For each site, analyses were conducted separately for three classes of mean daily solar radiation: and K# > 240 W m À2 .Given the highly skewed distribution of streamflow, its logarithm was used in the analysis.
The analysis first addressed hypothesis (b)that the relationship has an inverted-U shape.Testing the significance of non-monotonic relations is not straightforward.For example, one common approach is to fit both linear and quadratic regressions and to evaluate the significance of the quadratic term.However, this approach can provide misleading results because the quadratic term can be significant even for a monotonic relationship (Lind & Mehlum, 2010).Another issue is that quadratic relationships are symmetric around the minimum/ maximum and thus may not fit relationships that are asymmetric.
In this analysis, we applied an ad hoc approach that is nonparametric and can accommodate asymmetric relationships.The analysis involved a first step of testing the alternative hypothesis of an inverted-U-shaped relationship versus the null hypothesis that the relationship is either monotonic or not significant.If the test failed to reject the null hypothesis, then a follow-up test evaluated the alternative hypothesis of a monotonic relationship versus the null hypothesis of no relationship.The steps are described below: 1. Fit a loess smooth with a span of 1 and degree = 2, which can generate concave/convex relations that are asymmetric.
2. Check that the minimum value of the fitted loess relationship occurs at one of the end-points (i.e., the minimum or maximum value of the observed x variable).If not, the relationship cannot be an inverted-U shape.
3. If the criterion in step 2 is met, use the fitted loess relationship to predict the y variable (in this case, the treatment effect) for a sequence of x values ranging from the minimum to the maximum value of x with an interval equal to the range of x divided by 1000.
Denote these predicted values as b y.
4. Find the maximum value of b y and the associated value of x, denoted x p .
5. Split the observed data into two subsets, one for x ≤ x p and one for x > x p .
6. Perform Spearman correlation tests between the y and x variables for each subset.
7. Conclude that the relationship is hump-shaped if the following conditions hold: (a) the Spearman correlation coefficient (r s ) for x ≤ x p is positive and statistically significant at a specified significance level (e.g., α = 0.05), and (b) r s for x > x p is negative and also statistically significant at the specified α.Both hypothesis tests were conducted as two-tailed tests using α = 0.05.
8. If hypothesis (b) is not accepted on the basis of step 2, above, or if the test in step 7 is not met, conduct a Spearman correlation test for the full data set to evaluate whether a monotonic relationship exists, that is, hypothesis (a).
The effective significance level (i.e., the probability of accepting the hypothesis if it were actually false) was determined using a randomization test with n = 10 000 iterations.In each iteration, the test steps above were applied after shuffling the treatment effects using the sample function in R. The probability of incorrectly accepting hypothesis (b), denoted p b , was then computed as where n b is the number of iterations in which hypothesis (b) was accepted.The probability was computed for all subsets and then averaged.

| RESULTS
Summer 2005 tended to be wetter than summer 2006, with lower solar radiation, especially in June and early July (Figure 3).Streamflow variability during summer 2006 was dominated by the seasonal recession trend.In contrast, a rain event in early July 2005 resulted in a marked increase in streamflow that interrupted the seasonal flow recession.
Temporal variations in streamflow monitored at the weir were highly correlated with temporal variations at each of the four upstream sites; Spearman correlation coefficients ranged from 0.94 to 1. Manual streamflow measurements generally indicated gaining conditions (Table 1); however, streamflow at the weir was consistently lower than streamflow measured at 80 m upstream (Figure 4).
As can be seen in Figure 5, the pre-harvest models were Day-to-day variability in the treatment effect was similar among the three sites, and varied in response to both rain events and changes in solar radiation (Figure 3).However, the seasonal trend differed among sites.Whereas the treatment effect at the 100-and 200-m sites tended to increase through most of the summer period, the 0-m site exhibited a peak in response that appears to be associated with streamflow dropping below a threshold about 0.2 L s À1 .
Figure 6 demonstrates that the treatment effect generally increased as solar radiation increased from low (<120 W m À2 ) to high (>240 W m À2 ) classes.As seen in Table 2, the analysis identified significant negative relationships between treatment effect and discharge for all but six cases.In three of these six cases, significant inverted-U-shaped relationships were identified: for the high solar The randomization tests indicated that the probability of incorrectly identifying a relationship as an inverted-U shape is 0.001.For the two cases at the 0-m site with inverted-U-shaped relationships, the streamflow associated with the peak treatment effect was about 0.7-0.8L s À1 .For the one case with a significant inverted-U-shaped relationship at the 100-m site, the threshold discharge was 0.4 L s À1 , and the slope of the left-hand segment of the inverted-U was substantially lower than was the case for the 0-m site.Thus, the relationship for that case is dominantly negative.
As seen in Figure 7, daily maximum stream temperature consistently increased from upstream to downstream for stream discharges greater than 0.2 L s À1 , and for the lowest solar radiation class for stream discharge less than 0.2 L s À1 .However, for days with the T A B L E 1 Manual streamflow measurements made using constant-rate injection of salt solution at four locations upstream of the weir, paired with streamflow at the weir (0 m) as estimated from a rating curve, for late May through August, 2005.F I G U R E 3 Time series of (from top to bottom) daily total rainfall (Rain), mean daily solar radiation (K #), mean daily streamflow (Q), and treatment effect for daily maximum stream temperature (T e ) for the three sites along Griffith Creek for June through August for the two postharvest summers.The y-axis for streamflow is logarithmic.The black lines in the plots of treatment effect were computed using the locally estimated scatterplot smoother (LOESS) algorithm with span = 1.
highest solar radiation and stream discharge less than 0.2 L s À1 , daily maximum stream temperature dominantly increased from 200 to 100 m above the weir, then decreased from 100 to 0 m.

| DISCUSSION
The 100-and 200-m sites dominantly exhibited negative relationships between treatment effect and streamflow, consistent with the findings of Oanh et al. (2021).These results are also consistent with the broader literature on the relationship between stream temperature and discharge (Asarian et al., 2023;Gu et al., 1998;Hockey et al., 1982;Webb et al., 2003).In contrast, the 0-m site exhibited a pronounced inverted-U-shaped relationship between treatment effect and streamflow for days with high solar radiation.
As seen in Figure 7, it was not just the treatment effect at 0 m that decreased on days with high solar radiation and low discharge, but also on many days the actual temperature itself, based on the difference in temperature compared to those measured 100 m upstream.Downstream cooling could result from one or a combination of the following mechanisms: (1) negative surface energy exchange, ( 2) inflow of cooler water (e.g., groundwater), (3) downward-directed bed heat conduction, and (4) hyporheic exchange (Story et al., 2003).
Comparison of discharge at the weir, as determined from a rating curve, and discharge measured 80 m upstream using constant-rate injection of salt solution.Error bars represent 95% prediction limits for weir discharge and a nominal AE10% uncertainty for the salt dilution measurements.The dashed line indicates perfect agreement.
F I G U R E 5 Time series of observed (red) and predicted (black) daily maximum stream temperature at three locations on Griffith Creek along with approximate prediction limits (grey ribbon).Predicted values are based on the fixed-effects coefficients, and prediction limits are computed as the predicted value AE3 s r , where s r is the standard deviation of the residuals from the fixed-effects predictions for the pre-harvest period.
Negative surface energy exchange is unlikely in this case because the thinning treatment extended down to a road and associated rightof-way clearing immediately south of the weir, such that the lower segment of the stream was likely even less shaded than the upstream portions of the treatment reach.Furthermore, negative surface heat flux would also drive a trend to downstream cooling for discharges greater than 0.2 L s À1 , which was not observed (Figure 7).However, as seen in Figure 7, the combination of high solar radiation and low discharge was associated with the highest water temperatures at the 100 m site, which would produce greater energy loss via the emission of longwave radiation at the water surface and increased evaporative loss.These surface energy losses would have reduced net surface T A B L E 2 Summary of statistical tests for inverted-U-shaped relationships between treatment effect and the logarithm of discharge.Note: Values for the discharge associated with the peak treatment effect (Q p ) are shown only for relationships that were found to exhibit an inverted-U shape.r s is the Spearman correlation coefficient between treatment effect and discharge for the full data set, p(r s ) is the associated p-value, and n is the sample size.
F I G U R E 6 Plot of observedpredicted stream temperature versus streamflow during June, July and August for the post-harvest summers.Each site is shown in a separate column, and the rows are based on solar radiation class, with high (K # > 240 W m À2 ) at the top.The streamflow axis is logarithmic.The fitted lines were computed using the locally estimated scatterplot smoother (LOESS) algorithm, with span = 1.
energy exchange and thus have worked together with other cooling mechanisms, such as inflow of cooler water, to promote downstream cooling between the 100-and 0-m sites.
Bed heat conduction and hyporheic exchange would both tend not only to reduce the daily maximum, but to increase the daily minimum (e.g., Johnson, 2004).The earlier analysis by Guenther et al. (2014) showed that the daily minimum temperature tended to increase following harvest, with average increases of about 1 C at both the 0-and 100-m locations.While the increase in daily minimum temperature is consistent with the moderating influence of hyporheic and bed heat fluxes, it is also consistent with an increase in shallow groundwater temperature (Guenther et al., 2014), and thus does not support either mechanism over the other.
We hypothesize that the pronounced inverted-U-shaped response at the 0-m site results from the presence of a localized inflow of water somewhere upstream of the 0-m site, which would be fed by longer and/or deeper flow paths and thus tend to be relatively cool during summer days.Under this scenario, streamflow at the 0-m site would comprise two components: a thermally more responsive component that had flowed downstream through the harvesting unit and been warmed by solar radiation, and a thermally more stable contribution from the localized inflow.As streamflow at the 0-m site declined, water passing the site would have been increasingly dominated by the cooler inflow, which would dilute the warming effect associated with upstream shade removal during daytime.If the water were discharging into the stream via the stream bed, then the bed heat flux could have contributed to or even dominated the cooling due to enhanced temperature gradients and associated downward heat conduction (Caissie & Luce, 2017;Guenther et al., 2014;Story et al., 2003).
The source of the hypothesized subsurface discharge could be lateral hillslope drainage or hyporheic flow, or a combination.If the source of the discharge was lateral hillslope drainage, discharge at the weir should have been higher than at 80 m unless there were complex losing/gaining sequences, for example, as reported by Story et al. (2003).The results in Figure 4 suggest that there was a net loss of flow in the lower 80 m of the stream reach, not a gain.However, it must be borne in mind that different approaches were used to generate streamflow estimates at the weir compared to the upstream sites, which could be subject to differential biases that would influence the inferred rate of inflow.In particular, although the weir was constructed from concrete poured into forms directly over the bedrock after excavating the soil, leakage could still have occurred, which may explain, at least in part, the apparent net flow loss.
Examination of a Lidar digital elevation model and field observations revealed no obvious surface topographic features that would generate concentrated inflow above the 0-m temperature logger.It may be that, during baseflow periods, the spatial pattern of subsurface flow was controlled by the topography of the soil-bedrock or soil-till interface, not the surface topography, consistent with findings by Hutchinson and Moore (2000) and van Meerveld et al. (2015).A hyporheic origin for at least a portion of the inferred subsurface discharge would be consistent with the fact that the temperature loggers were located in pools (to minimize risk of de-watering), which are often locations of hyporheic discharge, as demonstrated at other streams in Malcolm Knapp Research Forest (Moore et al., 2005;Scordo & Moore, 2009) and elsewhere (e.g., Kasahara & Wondzell, 2003).
A major focus of previous research on stream temperature response to shade removal was to evaluate the effectiveness of different forest practice rules, particularly in relation to their success in maintaining water quality standards (Groom, Dent, & Madsen, 2011).
A key takeaway from this study and that of Oanh et al. (2021)   more frequent extended drought.However, the results of this study also suggest that shifts in dominant streamflow sources and pathways during low-flow periods can complicate stream temperature response to forestry by producing zones of reduced impact.This inference is consistent with observations by Moore et al. (2005), who observed downstream cooling within a clearcut at locations with relatively steep increases in contributing area with downstream distance.These results add to a body of work that has demonstrated that discrete inflows of subsurface water to streams, although often difficult to locate and map, are common hydrologic features across stream networks that can have profound influences on aquatic ecosystem functioning (Briggs & Hare, 2018;Lupon et al., 2019;Ploum et al., 2021).
It has been argued that hydrologically adapted buffers that retain more forest cover around discrete inflow locations, by using wider buffers at these locations, would better protect stream ecosystems from harvesting (Creed et al., 2008;Erdozain et al., 2020;Kuglerová et al., 2017).More research is needed on accurate mapping of discrete inflows, quantifying the influence of discrete inflows on stream temperature response to harvesting, and assessing the effectiveness of buffer designs that attempt to preserve the function of discrete inflows on stream environments.

| CONCLUSIONS
This study addressed the relationship between temperature response to shade removal and streamflow in relation to two hypotheses: between summer stream temperature response to forest harvesting and stream discharge, and specifically in relation to two hypotheses: (a) that stream temperature response increases in magnitude as stream discharge declines, and (b) that the relationship can be complicated by shifts in dominant streamflow sources and pathways during low-flow periods.The study focused on a headwater catchment in the transient snow zone of the Pacific Northwest of North America, located adjacent to the catchment studied by Leach and Moore (2017).

1
Map of study catchment, showing the catchment and cutblock boundaries and the locations of the monitoring sites.The catchment boundary is based on a field survey conducted as part of the harvest planning.Locations of the temperature loggers on the control stream are indicated by triangles.Douglas-fir (Pseudotsuga menziesii) that had regenerated following harvest in the 1920s.Tree height prior to harvesting was about 30 to 40 m, and crown closure typically exceeded 90%.
cally.Time series of treatment effect at the three sites, rainfall, solar radiation, and streamflow were generated to visualize the covariation of treatment effect in relation to hydrometeorological drivers.The time series of treatment effect were filtered using the locally estimated scatterplot smoother (LOESS) algorithm (Chambers & Hastie, 1992) using the loess function in base R (R Core Team, 2021) with the "span" set equal to one to visualize the seasonal trend.The relationship between treatment effect and stream discharge was analysed in relation to two hypotheses: (a) that the relationship, if significant, was monotonic, and (b) that the relationship exhibited an F I G U R E 2 Photographs of (left) the weir located near the downstream end of the harvest unit and (right) the stream and adjacent slopes taken about 100 m upstream of the weir.Both photographs were taken in 2005, the first summer following the thinning treatment.
reasonably precise; standard deviations of the residuals from the fixed-effects portion of the model were 0.26 C or less.Whereas the predicted temperatures reached seasonal maxima around 15 C in the post-harvest period, the observed temperatures approached or exceeded 20 C during 2006.
radiation class for both years at the 0-m site and for 1 year (2006) at the 100-m site.The other three non-significant cases were for the low and moderate solar radiation classes at the 0-m site in 2006 and the low solar radiation class at the 100-m site in 2006.

F
I G U R E 7 Daily maximum temperatures versus distance above the weir for the postharvest period, stratified by stream discharge (columns) and solar radiation class (rows).
(a) that stream temperature response increases as flow declines and (b) that the relationship can be complicated by shifts in dominant streamflow sources and pathways during low-flow periods.At the two upstream sites on the treatment stream, the treatment response dominantly exhibited the expected negative relationship with daily mean streamflow, consistent with hypothesis (a), especially for days with high incident solar radiation.This result suggests that the effectiveness of forest practice rules for protecting cold-water habitat may be reduced under future climatic conditions characterized by more frequent extended drought.However, stream temperature response at the most downstream site exhibited a pronounced inverted-U-shaped relationship with streamflow for days with high incident solar radiation, consistent with hypothesis (b).It is hypothesized that the response at the lower temperature logger was controlled by the existence of a stable, relatively cool inflow just upstream of the logger, which represented an increasing fraction of total flow as streamflow generated higher up in the catchment declined.There was a lack of convergent surface topography upstream of the logger, and it is hypothesized that localized inflow may have been controlled by the topography of the soil-till interface, and/or that the inflow was at least partially of hyporheic origin.Future studies of stream temperature response to forestry and other changes to riparian shade should sample the spatial variability of treatment response along the disturbed reach, and also include measurements of streamflow and other hydrological observations to assist in diagnosing the drivers of temporal and spatial variability in thermal response.
All values are in L s À1 ."NA" indicates a measurement was not made.