Recent bark beetle outbreaks have little impact on streamflow in the Western United States

In the Western United States (US), the current mountain pine beetle (MPB; Dendroctonus ponderosae) epidemic has affected more than five million hectares since its start in 1996, including headwater catchments that supply water to much of the Western US. There is widespread concern that the hydrologic consequences of the extensive pine tree die-off will impact water supply across the Western US. While forest disturbance studies have shown that streamflow increases in response to tree harvest, the actual effect of bark beetle infestations on water supply remains widely debated. The current study evaluates watershed-level response following bark beetle outbreak for 33 watersheds in seven western states. Streamflow records were investigated to assess whether the timing and amount of stream discharge during bark beetle outbreak and early recovery periods were significantly different to pre-outbreak conditions. Results show no significant modification in peak flows or average daily streamflow following bark beetle infestation, and that climate variability may be a stronger driver of streamflow patterns and snowmelt timing than chronic forest disturbance.


Introduction
Watersheds in the Western United States (US) are undergoing unprecedented levels of tree die-off due to wildfire, drought, and insect outbreaks (Williams et al 2010, Adams et al 2012, Anderegg et al 2012. Bark beetle outbreaks are typically endemic, causing minimal impact to catchment areas. However, the ongoing mountain pine beetle (MPB; Dendroctonus ponderosae,) epidemic has affected conifer forests at historic levels (Raffa et al 2008). Over 6 million forested hectares in the US and British Columbia have been impacted by bark beetles, with more than 5 million hectares affected by the MPB outbreak alone . This includes headwater catchments to the Colorado, Arkansas, Rio Grande, and Missouri Rivers. Many species of bark beetle introduce bluestain fungi into the tree xylem, which inhibits water flow (Paine et al 1997). The interruption of transpiration in MPB-attacked pine trees occurs in weeks to months (Hubbard et al 2013) ultimately leading to tree mortality (Edburg et al 2012). This is termed the 'green phase' of the attack. Frank et al (2014) hypothesize that mortality takes longer in spruce forests because, unlike pine trees, spruce trees can survive without tightly coupling stomatal conductance to plant hydraulics. The needles on the MPB-attacked tree change from green to red during the year following the attack (termed the 'red phase'), eventually falling off within three to four years following attack and the tree appears gray (Wulder et al 2006, Edburg et al 2012. Alterations to watershed vegetation can have a profound impact on the forest hydrology (Brown et al 2005, Adams et al 2012. An extensive body of literature exists on the hydrological consequences of tree harvest and has been used to predict the impact of bark beetle outbreaks on water resources (Adams et al 2012, Edburg et al 2012. Reviews of paired watershed studies by Brown et al (2005) and Stednick (1996) conclude that deforestation leads to an initial increase in streamflow due to reduced interception and transpiration and increased baseflow, with a more intense impact observed in wetter regions. The water yield increase is reduced as regrowth occurs and the Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.
Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. streamflow eventually returns to pre-harvest conditions. Canopy reduction decreases intercepted snow while also increasing through-canopy solar radiation, leading to higher snow accumulation in tree harvest areas and earlier, faster snowmelt (Jost et al 2007, Varhola et al 2010. However, the progressive reduction in forest canopy due to bark beetle outbreaks is different from forest harvest disturbance in that it occurs gradually and does not cause the complete removal of understory vegetation and canopy or directly impact non-host trees, muting the disturbance signal and making it challenging to predict the impact of widespread bark beetle infestation on watershed hydrology (Adams et al 2012, Mikkelson et al 2013.
The literature supporting hydrologic change resulting from bark beetle outbreak is limited and conflicting. While early watershed studies by Potts (1984) and Bethlahmy (1974) report increased streamflow, more recent paired watershed studies found highly variable streamflow modification Jenson 2007, Somor 2010) or no impact to streamflow (Biederman et al 2015) following bark beetle outbreak. Increased baseflow was reported by Bearup et al (2014) and remote sensing based studies report reduced MODIS-based evapotranspiration following bark beetle attack (Bright et al 2013, Maness et al 2013, Vanderhoof and Williams 2015. In contrast, several recent studies using eddy-covariance methods found that increased evapotranspiration by understory, secondary structure, and surviving canopy vegetation and the soil offset reductions due to tree mortality (Biederman et al 2014b, Brown et al 2014, Reed et al 2014. The results of snowpack studies are also inconsistent. While several stand-level studies report increased snow water equivalent in the snowpack under gray phase forests (Boon 2007, Pugh and Small 2012, 2013, Boon (2012) found this effect was diminished during high snow years due to the ability of large snowfall to exceed the interception capacity of the canopy. Pugh and Small (2012) found earlier spring snowmelt in red and gray phase forests in Colorado. Biederman et al (2014a) found that increased snowpack sublimation compensated for decreased canopy sublimation, resulting in no net change to snow water equivalent in gray phase forests. Recent literature suggests that the hydrologic response to bark beetle infestation is complex and dependent not only on tree mortality, but on ecosystem response to tree die-off, physical characteristics of the watershed, and regional climate.
The limited number and scale (stand or hillslope) of recent studies make it difficult to predict streamflow response to beetle infestation at larger, watershed scales, where water resource managers are tasked with critical decision-making. In addition, there is a paucity of literature that considers the hydrologic response of more than a few watersheds. The current study evaluates the streamflow and baseflow response of 33 watersheds in seven western states using a several hydrologic metrics, including timing and amount of peak flows and daily streamflow statistics. Watershed statistics after bark beetle outbreak were compared to pre-infestation conditions to answer the question: to what degree is the timing and amount of watershed discharge across the Western US modified by bark beetle infestation?
2. Data collection and methods 2.1. Study area description Spatial data identifying areas of bark beetle infestation in the contiguous Western US from 1997 to 2014 were obtained from the United States Forest Service Insect and Disease Detection Survey database (USDA Forest Service, Forest Health Protection and its partners 2015). The spatial data is based on aerial detection surveys (ADSs) conducted by the USDA. Accuracy assessments of ADSs data collected in USDA Forest Service Rocky Mountain Region compared to ground reference points show >70% accuracy at a 500 m scale, indicating that ADS data are appropriate for assessing forest disturbance at coarse scales (Johnson and Ross 2008). The tree mortality data contained in the ADS database were not used in this study because they have been shown to be inaccurate . Tree mortality estimates may be obtained using high resolution satellite imagery (e.g., as shown by Meddens et al 2012); however, it was not feasible to obtain these data products over all 33 watersheds included in this study. The ADS data were intersected with watershed boundaries from the GAGES II: Geospatial Attributes of Gages for Evaluating Streamflow dataset (GAGES II; Falcone et al 2010, Falcone 2011) to derive a time series of annual percent watershed area impacted by the MPB and spruce beetle (Dendroctonus rufipennis) for the catchment area of each GAGES II reference gage. The GAGES II reference gages are US Geological Survey (USGS) stream gages considered to have minimal regulation in their contribution catchments.
Watersheds were included in this study if 15% or more of their area was impacted by MPB and/or spruce beetle during at least one year between 1997 and 2014. This resulted in 33 watersheds (shown on figure 1) for the analyses. The cumulative area impacted since 1997 ranged from 21% to 90% of the total watershed area, with five watersheds having greater than 75% area impacted, 13 watersheds having between 50% and 75% area impacted, and 15 watersheds having between 21% and 50% area impacted. Dominant forest types in the bark beetle-impacted area, identified by the National Atlas of the United States (2002), are: lodgepole pine (Pinus contorta; 14 watersheds), ponderosa pine (Pinus ponderosa; eight watersheds), fir-spruce (mixed species; ten watersheds), and Douglas fir (Pseudotsuga menziesii; one watershed). Catchment size ranged from 15.5 to 3355 square kilometers. The portion of precipitation falling as snow ranged from 20.5% to 72.3% (Falcone 2011).

Identification of infection periods and other impacts to the study area
The ADS-derived time series of bark beetle impact were used to delineate three infestation periods for each watershed: (i) a pre-infestation phase representing the period prior to the bark beetle outbreak in the watershed; (ii) a mixed/red phase representing the period of most intense bark beetle infestation; and (iii) an early gray phase representing the early recovery period after the bark beetle outbreak. The mixed/red phase period was defined as the period inclusive of the first and last year where the ADSs indicate that MPB and/or spruce beetle impacted 15% or more of the watershed area. This is considered a mixed outbreak phase because the watershed includes trees in the green, red, and gray phases during this period. Because the ADSs data identify ongoing bark beetle infestation by the red foliage signature, the green phase is missed by the surveys and the gray phase is not recorded (Ciesla 2000). The early gray phase period was inferred from the ADS data as the five years immediately after the last year of the mixed/red phase period. The year immediately prior to the first year of the mixed/red phase period was considered to be the first year of the bark beetle attack (green phase) and was not included in the pre-outbreak period. Outbreak periods with gaps in the mixed/red phase period where the bark beetle-impacted area dropped below 15% were considered separate outbreaks if they included two or more consecutive years. This study did not include data collected during the second outbreak.
This study considers the watershed impact of the MPB and spruce beetle outbreaks together because these infestations tend to occur concurrently within western watersheds, though sometimes at a lagged timescale, based on our analysis of the ADSs data. The ADSs database indicates that some watersheds were impacted by other bark beetles [primarily the Douglas fir beetle (Dendroctonus pseudotsugae), western balsam beetle (Dryocoetes confuses), and fir engraver beetle (Scolytus ventralis)] during the red/mixed phase period, but not at levels above 10% areal impact prior to the outbreak phases or during the early grey phase periods. The maximum area burned in any watershed is less than 7% of the total watershed area.

Streamflow data
Time series of mean daily discharge and peak instantaneous flow data for the USGS gages corresponding to the 33 study area catchments were obtained from the USGS (2015) for water years (WYs) 1970-2014 (the water year begins on 1 October). The streamflow dataset includes a complete hydrologic record for seven or more water years prior to the green phase period of each study area watershed, as well as a complete hydrologic record for the red/mixed phase periods. Streamflow data is available for the complete early gray phase period for 20 of the 33 study area watersheds. A shorter period of record is available for the early gray phase of the 13 remaining watersheds due to late onset of the outbreak or the presence of a second outbreak.

Climate data
Temperature and precipitation data for the contiguous Western US for WYs 1970-2014 were obtained from the Parameter-Elevation Regressions on Independent Slopes Model (PRISM) database at a spatial resolution of approximately 4 km (Daly et al 2008, PRISM Climate Group, Oregon State University 2015. Time series of mean monthly temperature and total monthly precipitation depth for each watershed were derived from the gridded data.

Streamflow and climate statistics
Statistics describing streamflow distribution, magnitude, high flows, low flows, and baseflow were implemented in R (R Core Team 2015) on a WY basis. The magnitude and distribution of streamflow discharge were described by the minimum (min), maximum (max), mean, and median 1 day mean daily discharge statistics and by total annual discharge (discharge). Three statistics were used to describe high flows: timing of the center of mass of the annual flow (CT), peak flow (peak), and water year day of the peak flow (peak day). CT was calculated from mean daily streamflow data as described by Stewart et al (2005) as follows: where t i is time in days corresponding to the WY day (1 October is day one of the WY), and q i is the corresponding mean daily flow measurement. As snowmelt runoff is the largest contribution to streamflow in snow-dominated watersheds, CT may be used as a proxy for the timing of the snowmelt pulse (Stewart et al 2005). Peak flow and peak day were determined from the peak instantaneous flow dataset.
For the purposes of this study, peak day is defined as the water year day during which the instantaneous peak flow occurred. Low flows were described by the annual minimum of the 7 day moving average of the mean daily discharge (7 day min). Summer-fall baseflow (BF) was calculated from the mean daily streamflow dataset following Nathan and McMahon (1990) using the EcoHydRology package (Fuka et al 2014) in R. A filter parameter of 0.925 and three filter passes were employed for this study, following the recommendations of Nathan and McMahon (1990). Summer-fall baseflow index (BFI) was calculated by dividing the summer-fall baseflow by the total streamflow for the same period. Annual and seasonal climate statistics were calculated for comparison to the streamflow statistics. Total precipitation and mean temperature were calculated on a seasonal and WY basis from the climate datasets for each watershed. In addition, annual runoff ratio was calculated for each watershed.

Change and trend detection
The Wilcoxon rank-sum test and Mann-Kendall trend test were used to test for significant differences and trends in streamflow and climate statistics prior to and during the bark beetle outbreaks. These tests are nonparametric and do not require that the data conform to the normal distribution. For this study, p values <0.05 were considered statistically significant. The Mann-Kendall test was used to test for monotonic increasing or decreasing trends in the streamflow and climate statistics over the entire time series. The Wilcoxon rank-sum test was applied to streamflow and climate statistics for the red/mixed phase and early gray phase periods to test for evidence of significant differences compared to those for the preoutbreak period. The pre-outbreak period included available climate and streamflow data from WY 1970 until the first year of the bark beetle attack (green phase). No correction for multiple comparisons was made in determining significance due to the conservative nature of these corrections for large numbers of samples.
The Durbin-Watson test found some evidence of serial correlation in the climate and streamflow statistics (see supplemental data). Serial correlation increases the probability that the test detects a significant trend (using the Mann-Kendal test) or significant difference (using the Wilcoxon rank sum test) when a one is not actually present (type I error). No correction for serial autocorrelation was made because few significant trends or differences were detected in the hydrologic times series data, as discussed further in the next section. Prewhitening the data following von Storch (1995) prior to statistical analysis was attempted and found to have little effect on the study findings. Prewhitening resulted in fewer significant differences, but with increased type II error rates (the probability of not detecting a significant trend or difference when one is actually present.

Results and discussion
Precipitation and temperature data were analyzed for trends over the 35 year period of record and for significant changes after onset of bark beetle outbreak. Wilcoxon rank sum test results (table 1) for mean annual precipitation indicate that there were no significant differences in precipitation in the red/ mixed phase and early gray phase periods compared to pre-outbreak years. This is consistent with Chapman et al (2012), who found that the MPB epidemic was triggered during a short, intense drought from 2002 to 2003 but that rainfall returned to normal levels after 2003. The Mann-Kendall test results (table 1) indicate that a significant decreasing monotonic precipitation trend over the 35 year period of record  is present in only four of the 33 watersheds. No significant increasing trends were identified. Insignificant trends in precipitation largely follow regional patterns. Total annual precipitation in watersheds located in the Rocky Mountain range decreased over the 35 year record, except for basins in Colorado, while precipitation in the Colorado, Cascadian, and western South Dakota watersheds increased. Similar to other studies (e.g., Kunkel et al 2013a, 2013b), a warming trend was observed over the study area. Mann Kendall test results indicate that 26 of the 33 watersheds show evidence of a significant increasing trend in temperature over the 35 year study period. Wilcoxon rank sum test results indicate that mean annual temperature is higher in nearly all watersheds after onset of bark beetle. The increases were significant in seven watersheds for the red/mixed phase and in 19 watersheds for the early gray phase.
The Wilcoxon rank sum test results (table 1) indicate that there are few significant changes to discharge statistics following bark beetle outbreak. Changes to daily streamflow statistics (mean, median, min, max, and total discharge) are split nearly evenly between those with higher discharge statistics after onset of bark beetle infestation and those with lower statistics. Results for runoff ratio are also split, generally following the changes to discharge and mean streamflow with few statistically significant changes. The low flow (7 day min) and baseflow (BF and BFI) streamflow statistics do not show consistent patterns of higher or lower values after onset of bark beetle infestation. While our study shows significant trends in earlier peak day in seven watersheds and earlier CT in four watersheds over the 35 year period of record, the Wilcoxon rank sum tests found few significant changes in the high flow (peak day, peak, and CT) streamflow statistics after onset of bark beetle infestation. The insignificant changes in spring melt timing do not suggest a consistent pattern of higher or lower values following bark beetle infestation, and, in some cases, do not follow the 35 year trend.
While forest cover type, percent of watershed area impacted, and percent of precipitation that falls as snow result in a distributional change in the outcome, there are few statistically significant differences in the outcome (see supplementary data). The majority of significant changes were detected in watersheds with smaller catchment areas (table 2). The two watersheds with significantly higher discharge, mean, and min statistics during the red/mixed phase concurrently experienced higher precipitation with p-values around 0.10 (see supplemental data). Higher precipitation can increase snowpack, which could be further influenced by canopy loss (Boon 2007, Pugh and Small 2012, 2013. The majority of the significantly higher daily streamflow, low flow, and baseflow statistics are associated with periods of higher precipitation and significantly lower statistics are associated with lower precipitation (table 3). Further, significantly higher peak flows are associated with periods with higher spring precipitation (table 4), significantly higher CT and peak day are associated with periods of lower spring temperature, and significantly lower CT Table 1. Number of watersheds with lower, higher, and no change in climate and streamflow statistics, as indicated by the Wilcoxon rank sum test results, for two time period comparisons: (i) red/mixed phase compared to the pre-onset period (n = 33) and (ii) early gray phase compared to the pre-onset period (n = 32). Number of watersheds with increasing, decreasing, or no trend in climate and streamflow statistics, as indicated by the Mann Kendall trend test results for the full study record (1970-2014, n = 33). Numbers of watersheds with significant change (α = 0.05) are indicated in parentheses.
and peak day area associated with periods of higher spring temperature (table 5). This is not surprising, as earlier snowmelt have been linked to increasing temperature trends in the Western US (Barnett et al 2004, Stewart et al 2005, Clow 2010, Harpold et al 2012. These results suggest that the most significant changes in streamflow statistics are likely climate related, rather than due to bark beetle impacts and that these changes are more pronounced in watersheds with smaller catchment areas. The autocorrelation in the data and the absence of a correction for multiple comparisons in determining significance, discussed in section 2.6, reinforce this result, as they serve to overestimate the number of significant results in table 1. Similar to results reported by Stednick and Jenson (2007), Somor (2010), and Biederman et al (2015), this study does not find evidence of consistent increased streamflow or baseflow following bark beetle disturbance. Instead, the results suggest that the post-outbreak streamflow statistics generally fall within the range of the pre-outbreak variability.

Conclusions
The current bark beetle infestations in the Western US were predicted to modify watershed hydrology, assuming the reduced transpiration and interception of the attacked trees would increase stream discharge and the reduction in forest canopy would lead to higher peak flows and earlier snowmelt timing (Edburg et al 2012, Mikkelson et al 2013. While bark beetles immediately and severely reduce the transpiration of the attacked tree, there is a growing body of empirical evidence indicating that the excess moisture is utilized by the ecosystem response, resulting in little net impact to discharge (Biederman et al 2014b, Brown et al 2014, Reed et al 2014. In addition, canopy reduction has been shown to have inconsistent impact on snow water equivalent in the snow pack and on snowmelt timing (Boon 2012. Our results show no significant change in daily streamflow statistics, peak flow, or snowmelt timing relative to the current bark beetle outbreaks. Climate variability may be a stronger driver of streamflow patterns and snowmelt timing than chronic forest disturbance, as post-outbreak flows generally fall within the range of the pre-outbreak variability. The impact of bark beetles on water quality, however, remains an open question (Mikkelson et al 2013). This work expands upon previous stand-and plot-scale findings to the watershed scale, providing more evidence that the current bark beetle outbreak is not significantly altering streamflow hydrology across western US watersheds.