How does water yield respond to mountain pine beetle infestation in a semiarid forest?

Mountain pine beetle (MPB) outbreaks in the western United States result in widespread tree mortality, transforming forest structure within watersheds. While there is evidence that these changes can alter the timing and quantity of streamflow, there is substantial variation in both the magnitude and direction of hydrologic responses, and the climatic and environmental mechanisms driving this variation are not well understood. Herein, we coupled an eco-hydrologic model (RHESSys) with a beetle effects model and applied it to a semiarid watershed, Trail Creek, in the Bigwood River basin in central Idaho, USA, to examine how varying degrees of beetle-caused tree mortality influence water yield. Simulation results show that water yield during the first 15 years after beetle outbreak is controlled by interactions between interannual climate variability, the extent of vegetation mortality, and long-term aridity. During wet years, water yield after a beetle outbreak increased with greater tree mortality; this was driven by mortality-caused decreases in evapotranspiration. During dry years, water yield decreased at low-to-medium mortality but increased at high mortality. The mortality threshold for the direction of change was location specific. The change in water yield also varied spatially along aridity gradients during dry years. In wetter areas of the Trail Creek basin, post-outbreak water yield decreased at low mortality (driven by an increase in ground evaporation) and increased when vegetation mortality was greater than 40 % (driven by a decrease in canopy evaporation and transpiration). In contrast, in more water-limited areas, water yield typically decreased after beetle outbreaks, regardless of mortality level (although the driving mechanisms varied). Our findings highlight the complexity and variability of hydrologic responses and suggest that long-term (i.e., multi-decadal mean) aridity can be a useful indicator for the direction of water yield changes after a disturbance.

1 Introduction

In recent decades, mountain pine beetle (MPB) outbreaks in the Western US and Canada have killed billions of coniferous trees (Bentz et al., 2010). Coniferous forests can provide essential ecosystem services, including a water supply for local communities (Anderegg et al., 2013). Therefore, it is essential to understand how ecosystems and watersheds respond to beetle outbreaks and to identify the dominant processes that drive these responses (Bennett et al., 2018). A growing number of studies have qualitatively examined hydrologic responses to beetle outbreaks and disturbance; however, these studies have produced conflicting results (Adams et al., 2012; Goeking and Tarboton, 2020). While some studies show increases in water yield following beetle outbreaks (e.g., Bethlahmy, 1974; Livneh et al., 2015; Potts, 1984), many others show no change, or even decreases (e.g., Biederman et al., 2014; Guardiola-Claramonte et al., 2011; Slinski et al., 2016). To determine which mechanisms control the change in water yield following a beetle outbreak, more quantitative approaches are needed.

Water yield is often thought to increase after vegetation is killed or removed by disturbances such as fire, thinning, and harvesting (Buma and Livneh, 2017; Chen et al., 2014; Hubbart, 2007; Robles et al., 2014; Wine et al., 2018). In the Rocky Mountain West, beetle outbreaks have increased water yield through multiple mechanisms. First, defoliation/needle loss can reduce plant transpiration, canopy evaporation, and canopy snow sublimation losses to the atmosphere (Montesi et al., 2004). Snow sublimation is an important process in snow-dominated forest systems. Beetle-caused decreases in total sublimation can increase water yield, especially since canopy sublimation is more sensitive to disturbances than ground snow sublimation (Frank et al., 2019). Increased canopy openings can also enable snow accumulation and allow more radiation to reach the ground surface, resulting in earlier and larger peak snowmelt events, which can in turn reduce soil moisture and therefore decrease summer evapotranspiration (ET).

Several studies have documented decreases in water yield following disturbances (e.g., mortality, fire, beetle outbreaks; Bart et al., 2016; Biederman et al., 2014; Goeking and Tarboton, 2020; Slinski et al., 2016). For example, in the southwestern US, beetle outbreaks have decreased streamflow by opening forest canopies and increasing radiation to the understory and at the ground surface, which leads to increases in understory vegetation transpiration (Guardiola-Claramonte et al., 2011) and soil evaporation, therefore increasing total ET (Bennett et al., 2018). Tree mortality or removal can reduce streamflow because surviving trees and/or understory vegetation compensate by using more water (Tague et al., 2019).

In a review of 78 studies, Goeking and Tarboton (2020) concluded that the decrease in water yield after tree mortality mainly happens in semiarid regions. Previous studies also provide rule-of-thumb thresholds above which water yield will increase: a loss of vegetation cover of at least 20 % and a mean precipitation of 500 mm yr⁻¹ (Adams et al., 2012). However, many watersheds in the western US experience high interannual climate variability (Fyfe et al., 2017), and local environmental gradients (e.g., long-term aridity gradients) may strongly influence vegetation and hydrologic responses to disturbances such as beetle outbreaks, making such rules of thumb difficult to apply in practice (Winkler et al., 2014). Given the possibility of either increases or decreases in water yield following beetle outbreaks, modeling approaches are crucial for identifying the specific mechanisms that control these responses.

Our overarching goal was to identify mechanisms driving the direction of change in annual water yield after beetle outbreaks in semi-arid regions (note that in the following text, “water yield” refers to the mean annual water yield). To accomplish this goal, we asked the following questions:

• Q1: What is the role of interannual climate variability in water yield response?

• Q2: What is the role of mortality level in water yield response?

• Q3: How does long-term aridity (defined as the temporally averaged potential evapotranspiration relative to precipitation) modify these responses, and how do responses vary spatially within a watershed along aridity gradients?

We hypothesized that multiple ecohydrologic processes (e.g., snow accumulation and melt, evaporation, transpiration, drainage, and a range of forest structural and functional responses to beetles) could interactively influence how water yield responds to beetle outbreaks. However, one or more processes may dominate in certain locations. In addition, the dominant ecohydrologic processes may vary over space and time due to interannual climate variability (i.e., precipitation), vegetation mortality, and long-term aridity. In Sect. 2 of this paper, we present a conceptual framework for identifying and depicting dominant hydrologic processes through which forests respond to beetle infestation. We used this framework to interpret the modeling results. In Sect. 3, we describe our mechanistic modeling approach, i.e., using the Regional Hydro-Ecological Simulation System (RHESSys), which can prescribe a range of vegetation mortality levels, capture the effects of landscape heterogeneity and the role of lateral soil moisture redistribution, and project ecosystem carbon and nitrogen dynamics, including post-disturbance plant recovery. In Sects. 4 and 5, we then present modeling results that explore how multiple mechanisms influence water yield responses. This study can help inform management in beetle-affected watersheds by providing a tool for identifying locations that should be prioritized to mitigate flooding and erosion risk under different climate conditions.

Figure 1 Mechanism of water yield responses to beetle-caused mortality during the red and gray phases (0–10 years after beetle outbreak); semicircular boxes represent understory responses and square boxes represent overstory responses.

2 Conceptual framework

2.1 Vegetation response to beetle outbreaks

Mountain pine beetles (MPB) introduce blue stain fungi into the xylem of attacked trees, which reduces water transport in plants and eventually shuts it off (Paine et al., 1997). During outbreaks, MPBs prefer to attack and kill larger host trees that have greater resources (e.g., carbon), while smaller-diameter host trees and non-host vegetation (including the understory) remain unaffected (Edburg et al., 2012). After an MBP outbreak, trees experience three phases (i.e., red, gray, and old) over time (Hicke et al., 2012). During the red phase, the trees' needles turn red. During the gray phase, there are no needles in the canopy. During the old phase, killed trees have fallen and understory vegetation and new seedlings experience rapid growth (Hicke et al., 2012; Mikkelson et al., 2013).

2.2 Hydrologic response to beetle outbreaks

Figure 1 describes the main processes that alter evapotranspiration and either decrease or increase water yield, depending on which processes dominate (Adams et al., 2012; Goeking and Tarboton, 2020). During the red and gray phases, needles fall to the ground, the leaf area index (LAI) is lower, and the canopy is more open (Hicke et al., 2012). During these phases, changes in canopy cover can interact with hydroclimatic conditions across a watershed to alter transpiration and evaporation in a variety of ways.

Canopy mortality can reduce transpiration rates in infected trees, though in water-limited environments, remaining trees may compensate to some extent by increasing transpiration (Adams et al., 2012; Tague et al., 2019). More open canopies (i.e., following tree mortality) intercept less precipitation than closed canopies, reducing canopy evaporation but potentially increasing it from soil and litter layers (Montesi et al., 2004; Sexstone et al., 2018). Meanwhile, opening canopies can also increase the proportion of snow falling to the ground and, therefore, increase snowpack accumulation. With more solar radiation reaching the ground, earlier and larger peak snowmelt can also occur (Bennett et al., 2018). Generally, earlier snowmelt increases water for spring streamflow and decreases water for summertime ET (Pomeroy et al., 2012). However, once snags fall, reductions in longwave radiation can actually lead to later snowmelt (Lundquist et al., 2013). An open canopy, combined with less competition for resources such as solar radiation and nutrients, can also promote understory vegetation growth, which may increase understory transpiration (Biederman et al., 2014; Tague et al., 2019). In some riparian corridors, the regreening of surviving vegetation and the compensatory response of remaining tissues could diminish the reduction in ET caused by foliage fall, leading to no significant water yield response to beetle-caused mortality (Nagler et al., 2018; Snyder et al., 2012). Therefore, whether water yield increases or decreases following beetle outbreak ultimately depends on how these processes interact.

Figure 2 Land cover, elevation, and tree mortality for Trail Creek: (a) is a land cover map showing the main vegetation types, (b) shows the elevation gradient, and (c) shows the severity of beetle-caused tree mortality (during the period 2003–2012; Meddens et al., 2012). Note that, for our modeling experiments, we prescribed that the beetle outbreak occurred uniformly across the evergreen patches instead of using historical beetle outbreak data.

Interannual variability in climate (e.g., dry versus wet years) can affect how hydrologic processes interact in forested watersheds (Goeking and Tarboton, 2020; Winkler et al., 2014). For instance, during wet years, remaining plants are not water limited, and reductions in plant transpiration due to beetle-caused mortality dominate over increases in soil evaporation or remaining plant transpiration, resulting in a higher water yield. In contrast, during dry years, plants are already under water stress and decreases in plant transpiration caused by tree mortality may be compensated for by increased soil evaporation and transpiration by remaining trees or understory vegetation, leading to declines in water yield. Moreover, these responses are also affected by land cover types (e.g., young versus old trees, different tree species), which is not currently well documented (Morillas et al., 2017; Perry and Jones, 2017).

2.3 Review of modeling approaches

Many models, ranging from empirical and lumped to physically based and fully distributed, have been used to study hydrologic responses to disturbances. Goeking and Tarboton (2020) argue that only physically based and fully distributed models can capture how disturbances alter water yield because they represent fine-scale spatial heterogeneity and physical process that vary over space and time. Despite their advantages, many process-based models such as the coupled CLM-ParFlow model (Mikkelson et al., 2013; Penn et al., 2016), the Distributed Hydrology Soil Vegetation Model (Livneh et al., 2015; Sun et al., 2018), and the Variable Infiltration Capacity Model (Bennett et al., 2018) also have some limitations. For example, (1) they may assume a constant LAI after disturbances and static vegetation growth (e.g., VIC and DHSVM), (2) they may omit lateral flow that redistributes soil moisture (VIC), and (3) in some cases, the approach used to represent the effects of beetle outbreaks may be too simplified (e.g., changing only LAI and conductance without considering two-way beetle–vegetation interactions in post-disturbance biogeochemical and water cycling, as in CLM-ParFlow). Thus, improving current fully distributed process-based models to capture the coupled dynamics between hydrology and vegetation at multiple scales is critical if we are to project how beetle outbreaks will affect water yield in semiarid systems (Goeking and Tarboton, 2020). Here we use the RHESSys–Beetle model, which captures these processes.

3 Model, data, and simulation experiment design

3.1 Study area

Our study focused on the Trail Creek watershed, which is located in Blaine County between the Sawtooth National Forest and the Salmon-Challis National Forest (43.44^∘ N, 114.19^∘ W; Fig. 2). This is a 167 km² sub-catchment in the south part of Big Wood River basin, and is within the wildland–urban interface, where residents are vulnerable to the flood and debris flows caused by forest disturbances (Skinner, 2013). Trail Creek has frequently experienced beetle outbreaks, notably in 2004 and 2009, when beetles killed 7 and 19 km² of trees, respectively (Berner et al., 2017).

Figure 3 Evergreen forest cover percentage and ID of each sub-basin in Trail Creek, and long-term aridity index for this watershed. Aridity index is defined as annual mean potential evapotranspiration (PET)/precipitation (P) from 38 years of data (see Sect. 3.4); PET/P>2 is water limited, PET/P<0.8 is energy limited, and PET/P between 0.8 and 2 is balanced. Recall that only evergreen forest trees are attacked during beetle outbreaks.

Trail Creek has cold, wet winters and warm, dry summers; mean annual precipitation is approximately 978 mm, 60 % of which falls as snow (Frenzel, 1989). The soil is mostly permeable coarse alluvium (Smith, 1960). Elevations range from 1760 to 3478 m. Along this elevation gradient, there are also strong vegetation and aridity gradients (Fig. 3). The northern (higher elevation) portion of the basin is mesic and covered principally by evergreen forest containing Douglas fir (Pseudotsuga menziesii), lodgepole pine (Pinus contorta var. latifolia), subalpine fir (Abies lasiocarpa), and mixed shrub and herbaceous vegetation. The southern (lower elevation) portion is xeric and covered by shrubs, grasses, and mixed herbaceous species, including sagebrush, riparian species, and grasslands (Buhidar, 2002).

In total, Trail Creek contains 72 sub-basins, and two of them (sub-basins 412 and 416 in Fig. 3) are urban areas. When we classified this basin into different zones according to an aridity index, i.e., the ratio of the 38 year average annual potential evapotranspiration (PET) to the precipitation (P) (Sect. 3.4), two distinct hydrologic–vegetation cover types emerged: the high-elevation northern area is balanced (i.e., PET/P is between 0.8 and 2) and evergreen tree coverage is more than 50 %; the southern part is water limited (i.e., PET/P>2) and evergreen tree coverage is less than 30 % (Figs. 2 and 3).

3.2 Model descriptions

3.2.1 Ecohydrologic model

The Regional Hydro-Ecological Simulation System (RHESSys; Tague and Band, 2004) is a mechanistic model designed to simulate the effects of climate and land-use changes on ecosystem carbon and nitrogen cycling and hydrology. RHESSys fully couples hydrologic processes (including streamflow, lateral flow, ET, and soil moisture, etc.), plant growth and vegetation dynamics (including photosynthesis, maintenance respiration, and mortality, etc.), and soil biogeochemical cycling (including soil organic matter decomposition, mineralization, nitrification, denitrification, and leaching, etc.). It has been widely tested and applied in several mountainous watersheds in western North America, including many in the Pacific and Inland Northwest (e.g., Tague and Band, 2004; Garcia and Tague, 2015; Hanan et al., 2017, 2018, 2021; Lin et al., 2019; Son and Tague, 2019).

RHESSys represents a watershed using a hierarchical set of spatial units, including patches, zones, sub-basins, and the full basin, to simulate various hydrologic and biogeochemical processes occurring at these scales (Tague and Band, 2004). The patch is the finest spatial scale at which vertical soil moisture and soil biogeochemistry are simulated. In every patch, there are multiple canopy strata layers to simulate the biogeochemical processes related to plant growth and nutrient uptake. Meteorological forcing inputs (e.g., temperature, precipitation, humidity, wind speed, and solar radiation) are handled at the zone level and spatially interpolated and downscaled for each patch based on elevation, slope, and aspect. Sub-basins are closed drainage areas entering both sides of a single stream reach (the water budget is closed in sub-basins). The largest spatial unit is the basin, which aggregates the streamflow from sub-basins (Hanan et al., 2018; Tague and Band, 2004). In RHESSys, streamflow is the sum of overland flow and baseflow, and we consider the streamflow to be the water yield of each sub-basin.

RHESSys models vertical and lateral hydrologic fluxes, including canopy interception, plant transpiration, canopy evaporation/sublimation, snow accumulation, snowmelt and sublimation, soil evaporation, soil infiltration, and subsurface drainage. Canopy interception is based on the water-holding capacity of vegetation, which is also a function of plant area index (PAI). Both the canopy evaporation and transpiration are modeled using the standard Penman–Monteith equation (Monteith, 1965). Snow accumulation is calculated from incoming precipitation and is assumed to fall evenly across each zone. Snowmelt is based on a quasi-energy budget approach accounting for radiation input, sensible and latent heat fluxes, and advection. Soil evaporation is constrained by both energy and atmospheric drivers, as well as a maximum exfiltration rate controlled by soil moisture (Tague and Band, 2004). Vertical drainage and lateral flow are a function of topography and soil hydraulic conductivity, which decays exponentially with depth (Hanan et al., 2018; Tague and Band, 2004). Section S1 in the Supplement contains a more detailed synopsis of the soil hydrologic model.

Vegetation carbon and nitrogen dynamics are calculated separately for each canopy layer within each patch, while soil and litter carbon and nitrogen cycling are simulated at the patch level. Photosynthesis is calculated based on the Farquhar model, which considers the limitations of nitrogen, light, stomatal conductance (influenced by soil water availability), vapor pressure deficit, atmospheric CO₂ concentration, radiation, and air temperature (Farquhar and von Caemmerer, 1982; Tague and Band, 2004). Maintenance respiration is based on Ryan (1991), in which respiration is computed as a function of nitrogen concentration and air temperature. Growth respiration is calculated as a fixed ratio of new carbon allocation for each vegetation component (Ryan, 1991; Tague and Band, 2004). Net photosynthesis is allocated to leaves, stems, and roots at daily steps based on the Dickinson partitioning method, which varies with each plant development stage (Dickinson et al., 1998). LAI is estimated from leaf carbon and specific leaf area for each vegetation type. The soil and litter carbon and nitrogen cycling (heterotrophic respiration, mineralization, nitrification, and denitrification, etc.) are modified from the BIOME_BGC and CENTURY-NGAS models (Parton, 1996; Tague and Band, 2004; White and Running, 1994). A detailed description of RHESSys model algorithms can be found in Tague and Band (2004).

3.2.2 Beetle effects model

Edburg et al. (2011) designed and developed a model of MPB effects on carbon and nitrogen dynamics for integration with the Community Land Model Version 4 (CLM4) (Lawrence et al., 2011, Fig. 4). We integrated this beetle effects model into RHESSys (Fig. 4). Beetles attack trees mainly during late summer, with needles turning from green to red at the beginning of the following summer. We simplified this process by prescribing that tree mortality occurred on 1 September when representing a beetle outbreak in a given year. The advantage of this integration is that RHESSys accounts for the lateral connectivity in water and nitrogen fluxes among patches, which is not represented in CLM4 (Fan et al., 2019). Our approach differs from other hydrologic models of beetle effects (e.g., VIC, CLM-ParFlow, and DHSVM) because it includes dynamic changes in plant carbon and nitrogen cycling caused by beetle attack, plant recovery, and their effects on hydrologic responses. Previous studies of hydrologic effects of beetle outbreaks have mainly focused on consequences of changes in LAI and stomatal resistance during each phase of mortality while missing feedbacks between carbon and nitrogen dynamics, vegetation recovery, and hydrology (Bennett et al., 2018; Livneh et al., 2015; Mikkelson et al., 2013; Penn et al., 2016; Sun et al., 2018).

Figure 4 Conceptual framework of the beetle effects model. (a) Normal background mortality routine in RHESSys before the beetle outbreak. (b) Mortality from bark beetles. We added a snag (standing dead tree) and a dead foliage (needles still on dead tree) pool, as shown in the dashed and dotted circle, respectively. After a beetle outbreak, carbon (C) and nitrogen (N) move from the stem to the snag pool (dashed arrow). After staying in the snag pool for m years, C and N move on to the coarse wood debris pool (CWD) with an exponential decay rate to represent the snag fall (dashed arrow). There is a similar process for leaf C and N, which move from leaf to dead foliage to litter pool (dotted arrow). Furthermore, C and N in the CWD and fine root pools move to the litter pool immediately after the outbreak. Figure modified from Edburg et al. (2011).

To better represent the effects of beetle-caused tree mortality, we added a snag pool (standing dead tree stems) and a dead foliage pool (representing the red needle phase) in RHESSys (Fig. 4). All leaf biomass (including carbon and nitrogen) become part of the dead foliage pool. After 1 year, the dead foliage is transferred to the litter pool at an exponential rate with a half-life of 2 years (Edburg et al., 2011, 2012; Hicke et al., 2012). Similarly, stem carbon and nitrogen are moved to the snag pool immediately after the outbreak. After 5 years, carbon and nitrogen in snags begin to move into the coarse woody debris (CWD) pool at an exponential decay rate with a half-life of 10 years (Edburg et al., 2011, 2012). After the outbreak, the coarse roots that are killed move to the CWD and fine roots move to the litter pool. To simplify, we assume a uniform mortality level for all evergreen patches across the landscape. Due to the limitations of land cover data, we cannot separate the pine and fir in these evergreen patches. However, this does not affect the interpretation of our results because we analyze them based on mortality level and evergreen vegetation coverage rather than based on species.

In the integrated model, the reductions in leaf carbon and nitrogen after a beetle outbreak can directly decrease LAI and canopy height, which consequently affects energy (i.e., longwave radiation and the interception of shortwave radiation) and hydrologic (i.e., transpiration and canopy interception) fluxes. The model calculates two types of LAI: live LAI (i.e., only live leaves are included), and total LAI (i.e., both live and dead leaves are included). Plant transpiration is a function of the live LAI, while other canopy properties, including interception and canopy evaporation, are a function of the total LAI. The calculation of canopy height includes living stems and the snag pool.

3.3 Input data

We used the US Geologic Survey (USGS) National Elevation Dataset (NED) at 10 m resolution to calculate the topographic properties of Trail Creek, including elevation, slope, aspect, basin boundaries, sub-basins, and patches. Using NED, we delineated 16705 patches with 100 m resolution within 72 sub-basins. We used the National Land Cover Database (NLCD) to identify five vegetation and land cover types, i.e., evergreen, grass/herbaceous, shrub, deciduous, and urban (Homer et al., 2015). We determined soil properties for each patch using the POLARIS database (probabilistic remapping of SSURGO; Chaney et al., 2016). Parameters for soil and vegetation were based on previous research and literature (Ackerly, 2004; Berner and Law, 2016; Hanan et al., 2016, 2021; Law et al., 2003; White et al., 2000).

Climate inputs for this study, including maximum and minimum temperatures, precipitation, relative humidity, radiation, and wind speed, were acquired from gridMET for the years from 1980 to 2018. GridMET provides daily high-resolution ($\mathrm{1}/\mathrm{24}$^∘ or ∼4 km) gridded meteorological data (Abatzoglou, 2013). This is a blended climate dataset that combines the temporal attributes of gauge-based precipitation data from NLDAS-2 (Mitchell et al., 2004) with the spatial attributes of gridded climate data from PRISM (Daly et al., 1994).

3.4 Simulation experiments

To quantify how water yield responds to beetle-caused mortality, we prescribed a beetle outbreak in September 1989, with the same mortality level (%) applied to all evergreen patches for each sub-basin. After the beetle outbreak, red needles remained on the trees for 1 year before they started to fall (i.e., were transferred to the litter pool) at an exponential rate with a half-life of 2 years. The snag pool remained as standing trees for 5 years and then began to fall and were added to the CWD pool, which decayed at an exponential rate with a half-life of 10 years.

To address Q1 (i.e., the role of interannual variability), we compared water yield responses during a dry water year, 1994 (i.e., 5 years after the beetle outbreak, with an annual precipitation of 611 mm), to responses during a wet year, 1995 (i.e., 6 years after the beetle outbreak, with an annual precipitation of 1394 mm). This enabled us to estimate the role of interannual climate variability in driving changes in water yield following a beetle attack. The dry year was selected from years that had precipitation below the 15th percentile of the annual precipitation data (1980–2018; Searcy, 1959; see Fig. S1 in the Supplement). During the early period after the beetle outbreak (e.g., 1994 and 1995), the forest experienced large changes in vegetation canopy cover, plant transpiration, and soil moisture. We chose these two successive years because their canopy and vegetation statuses were similar in terms of fallen dead foliage and residual vegetation regrowth, which makes this comparison reasonable. However, it is possible that antecedent climate conditions may affect the following year's response. For example, soil moisture can be depleted during a drought year, affecting initial conditions the following year. Moreover, under drought conditions, less reactive nitrogen is taken up by the plants or leaching is reduced, so more nitrogen will be left for the following year. Therefore, the difference in water yield response between 1994 and 1995 might be affected by not only climate variations but also the initial hydrologic and biogeochemical conditions. To consider the time lag effect (antecedent conditions affecting the current year's response), we also analyzed other dry and wet years.

To address Q2 (i.e., the role of vegetation mortality), we prescribed a range of beetle-caused mortality levels (i.e., from 10 % to 60 % in steps of 10 % in terms of the reduction in carbon, which was uniformly applied to all evergreen patches for each sub-basin) and a control run (no mortality). This enabled us to quantify how forest water yield responded to the level of vegetation mortality (for each sub-basin, vegetation mortality was the evergreen mortality multiplied by the evergreen coverage of that basin). The differences in water yield between each mortality level and the control run represent the effects of beetle kill: a positive value means that mortality increased water yield, and vice versa.

We quantified the water budget for each sub-basin to examine which hydrologic processes contribute to the water yield responses: water yield (Q), precipitation (P), canopy evaporation (E_canopy: canopy evaporation and snow sublimation), transpiration (T), ground evaporation (E_ground, which includes bare soil evaporation, pond evaporation, and litter evaporation), snow sublimation (Sublim, ground), soil storage change ($\mathrm{d}{S}_{\mathrm{soil}}/\mathrm{d}t$), litter storage change ($\mathrm{d}{S}_{\mathrm{litter}}/\mathrm{d}t$), snowpack storage change ($\mathrm{d}{S}_{\mathrm{snowpack}}/\mathrm{d}t$), and canopy storage change ($\mathrm{d}{S}_{\mathrm{canopy}}/\mathrm{d}t$). We summarized these rate variables at annual time steps.

The storage components include soil, litter, and canopy. According to Eq. (1), if the storage increases, water yield decreases.

$$\begin{array}{}\text{(1)}& \begin{array}{rl}Q& =P-E_{\mathrm{canopy}}-E_{\mathrm{ground}}-\mathrm{Sublim}-T\\ & -\mathrm{d}\left(S_{\mathrm{soil}}+S_{\mathrm{litter}}+S_{\mathrm{canopy}}+S_{\mathrm{snowpack}}\right)/\mathrm{d}t.\end{array}\end{array}$$

Q is the water yield (mm yr⁻¹); P is the precipitation (mm yr⁻¹); E_canopy is the canopy evaporation (including the canopy snow sublimation; mm yr⁻¹); E_ground is the ground evaporation (including bare soil evaporation, pond evaporation, and litter evaporation; mm yr⁻¹); T is the transpiration (mm yr⁻¹); Sublim is the ground snow sublimation (mm yr⁻¹); $\mathrm{d}{S}_{\mathrm{soil}}/\mathrm{d}t$ is the change in soil water storage calculated at yearly intervals (mm yr⁻¹); $\mathrm{d}{S}_{\mathrm{litter}}/\mathrm{d}t$ is the change in litter water storage calculated at yearly intervals (mm yr⁻¹); ${\mathrm{S}}_{\mathrm{canopy}}/\mathrm{d}t$ is the change in canopy water storage calculated at yearly intervals (mm yr⁻¹), and $\mathrm{d}{S}_{\mathrm{snowpack}}/\mathrm{d}t$ is the change in snowpack water storage calculated at yearly intervals (mm yr⁻¹).

Water balance differences between different mortality scenarios and control scenario results were calculated using Eq. (2). (Note that precipitation is a model input and is unaffected by mortality and so ΔP=0.)

$$\begin{array}{}\text{(2)}& \begin{array}{rl}\mathrm{\Delta }Q& =\mathrm{\Delta }{E}_{\mathrm{canopy}}+\mathrm{\Delta }{E}_{\mathrm{ground}}+\mathrm{\Delta }\mathrm{Sublim}+\mathrm{\Delta }T\\ & +\mathrm{\Delta }\left(\mathrm{d}\left(S_{\mathrm{soil}}+S_{\mathrm{litter}}+S_{\mathrm{canopy}}+S_{\mathrm{snowpack}}\right)/\mathrm{d}t\right)\end{array}\end{array}$$

To address Q3 (i.e., the role of long-term aridity), we calculated a long-term aridity index (PET/P, Fig. 3) across the basin and analyzed the relationship between long-term aridity index and hydrologic response. As mentioned earlier, the long-term aridity index is defined as the ratio of the mean annual potential ET (PET) to the annual precipitation (P) averaged over 38 years (water year 1980–2018) of historical meteorological data. Based on the long-term aridity index, we classified our sub-basins into three types (i.e., water limited, balanced, and energy limited; McVicar et al., 2012; Table 1).

Table 1 Classification of the aridity index.

4 Results

4.1 Simulated vegetation response to a beetle outbreak at the basin scale

4.1.1 Vegetation response to a beetle outbreak

Figure 5 shows the basin-scale vegetation response after a beetle outbreak in 1989. The live LAI dropped immediately after the beetle outbreak, then gradually recovered to pre-outbreak levels during the following years (Fig. 5a). The total LAI (including dead foliage) slightly increased during the first 10 years after the beetle outbreak (1990–2000), which was due to the retention of dead leaves in the canopy and the simultaneous growth of the residual (unaffected) overstory and understory vegetation (Fig. 5b). The dead foliage pool (Fig. 5c) remained in place for 1 year and then began to fall to the ground (was converted to litter) exponentially with a half-life of 2 years, and the snag pool (Fig. 5d) remained in place for 5 years and then began to fall to the ground (was converted to CWD) exponentially with a half-life of 10 years. These behaviors of the dead foliage and snag pools are similar to those noted in Edburg et al. (2011), which demonstrates that the integrated model is simulating expected vegetation dynamics following a beetle outbreak.

Figure 5 Basin-scale vegetation responses after a beetle outbreak for different evergreen mortality levels. (a) Annual live leaf area index (live LAI), (b) annual total LAI (calculated including the dead foliage pool), (c) daily dead foliage pool, and (d) daily snag pool after the outbreak. The green and red background colors denote the periods before and after the beetle outbreak, respectively.

4.1.2 Time series of the hydrologic response to the beetle outbreak

Figure 6 shows the changes in simulated water fluxes and soil moisture over the basin after the beetle outbreak with various evergreen mortality levels. During the first 15 years after the beetle outbreak, scenarios where the evergreen mortality level was larger than zero had higher basin-scale water yields than the control scenario (where the evergreen mortality level was zero). This was especially true during wet years; however, there was no significant increase during dry years (i.e., 1992, 1994, 2001, and 2004; Fig. 6a). The year-to-year soil storage fluxes responded strongly in the first two years after the beetle outbreak and then stabilized to the pre-outbreak condition (Fig. 6b). Note that year-to-year soil storage change is not the same as soil water storage. After the beetle outbreak, the soil held some portion of the water that was not taken up by plants, but this was constrained by the soil water holding capacity. This phenomenon indicates that the soil has some resilience to vegetation change.

Figure 6 Basin-scale annual sum of hydrologic flux responses after the beetle outbreak (1989) for different evergreen mortality levels. (a) Annual water yield calculated as the annual sum of the basin streamflow and (b) annual soil water storage change calculated as the soil water storage at the end of the water year minus the soil water storage at the beginning of the water year. (c) Transpiration is the annual sum of the transpiration for both the overstory and the understory. (d) Evaporation is calculated as the annual sum of canopy evaporation, ground evaporation, and snow sublimation.

The beetle outbreak reduced transpiration during wet years but did not have significant effects in dry years (Fig. 6c). This occurred because transpiration in dry years was water limited and therefore much lower than the potential rate (more water was partitioned to evaporation; similar to Biederman et al., 2014). Thus, killing more trees had little effect on stand-scale transpiration because the remaining trees used any water released by the dead trees in dry years. On the other hand, plant transpiration in wet years was close to the potential rate; therefore, decreases in canopy cover reduced transpiration. There was no apparent effect of the beetle outbreak on snowmelt.

Snow sublimation played an essential role in driving the evaporation responses we observed. In the Trail Creek watershed, snow sublimation accounted for around 50 % of the total evaporation (not shown in the figure), and around 60 % came from the canopy. Canopy sublimation accounted for an even larger proportion of the total sublimation during high snow years (Figs. S7d and S1). These results are similar to other western US forests, where 50–60 % of the total sublimation has been found to come from canopy sublimation, which is more sensitive to beetle kill than ground snow sublimation (Frank et al., 2019; Molotch et al., 2007). We also found that during the first three years after beetle outbreak, when dead foliage was still on the canopy, canopy sublimation increased by approximately 6 % due to an increase in the total LAI as new needles grew and dead foliage remained on the canopy. This increased canopy snow interception and subsequent sublimation (Fig. 5). However, when the dead foliage fell to the ground and snags began to fall, the canopy sublimation decreased by approximately 10 % for the most severe mortality scenario (60 % evergreen mortality) compared to the no-outbreak scenario. This occurred because the canopy total LAI decreased and there was less canopy interception (Fig. 5).

Ground snow sublimation was less sensitive to beetle kill (Fig. S7b). In the first 3 years after beetle kill (at 60 % mortality), ground snow sublimation increased by approximately 7.5 % due to an increase in aerodynamic conductance caused by the higher understory canopy height. However, from 1993 to 2002, there were no obvious changes in ground snow sublimation after beetle outbreak. When all dead foliage and more than 50 % of snags had fallen to the ground, ground snow sublimation decreased because snowmelt increased as the canopy opened (Figs. 5 and S7b). In general, for the 60 % mortality scenario, the ground snow sublimation first increased by approximately 5 % when dead foliage was still on the trees, and then decreased by approximately 6 % when the canopy was open.

The evaporation response was opposite in dry and wet years: evaporation increased in dry years, while it decreased in wet years (Fig. 6d). This phenomenon is caused by tradeoffs and interactions among multiple processes, as will be explained in more detail in the next section.

Figure 7 Relationships between long-term aridity, vegetation mortality level, and differences in evaporation for a dry year (1994, a) and a wet year (1995, b). Differences are calculated as the normalized differences (%) in evaporation between each evergreen mortality scenario and the control run for no beetle outbreak. Vegetation mortality for each sub-basin is calculated as the percentage of evergreen patches multiplied by the level of evergreen mortality caused by beetles. Long-term aridity is defined as the temporally averaged (38 years) potential evapotranspiration relative to the precipitation.

4.2 The role of spatial heterogeneity in the water yield response

4.2.1 Spatial patterns of hydrologic response along the long-term aridity gradient

Evaporation

The beetle outbreak had opposite effects on evaporation in a dry year and in a wet year (Fig. 7). In the dry year, most sub-basins experienced higher evaporation for beetle outbreak scenarios than in the control scenario (Fig. 7a). This was the cumulative consequence of decreased canopy evaporation and increased ground (soil, litter, pond) evaporation due to decreases in the LAI (caused by mortality). In the dry year, the latter effect (i.e., increased ground evaporation) dominated over the former, leading to an overall increase in evaporation. When the vegetation mortality level (calculated as the percentage of evergreen patches in a sub-basin multiplied by the level of evergreen mortality caused by beetles) was higher than 20 %, a few sub-basins in the balanced (more mesic) area experienced a decrease in evaporation, indicating that the effects of decreasing canopy evaporation outstripped the effects of increasing ground evaporation. In the wet year, most of the sub-basins located in the balanced area (where canopy evaporation decreases dominated) experienced decreases in evaporation. The decrease in evaporation responded linearly to the level of vegetation mortality (Fig. 7b). However, sub-basins located in much drier regions (aridity > 3.5) had relatively minimal responses to the vegetation mortality level, and some even had slight increases in evaporation (where ground evaporation increases were dominant due to the drier long-term climate and lower canopy mortality resulted from less evergreen coverage).

Transpiration

The beetle outbreak decreased transpiration in both dry and wet years, and the decrease became larger with higher mortality levels (Fig. 8). However, during the dry year, the water-limited area experienced less change than the balanced area; some sub-basins even showed slight increases. This increase in the water-limited part of the basin occurred because, after the beetles have killed some overstory evergreen trees, the living trees and understory plants can together exhibit higher transpiration rates in dry years (Tsamir et al., 2019). In the wet year, when most canopies reach their potential transpiration rates (less competition for water), beetle outbreaks can reduce transpiration rates by decreasing the live LAI.

Figure 8 Relationships between long-term aridity, vegetation mortality, and difference in transpiration for a dry year (1994, a) and a wet year (1995, b).

Total ET

In a dry year, the balanced and water-limited areas had opposite responses to mortality: the balanced area experienced a decrease in ET and the water-limited area experienced a slight increase (Fig. 9). In the balanced area, larger ET decreases occurred with higher mortality levels. However, increases in ET in water-limited regions were less sensitive to vegetation mortality level; even when mortality was high (>40 %), ET still increased (Fig. 9a). During the wet year, most sub-basins experienced decreasing ET after the beetle outbreak, and the magnitude was larger with higher vegetation mortality. The different responses of ET were driven by different hydrologic responses (transpiration, ground evaporation, and canopy evaporation) competing with each other; this competition was influenced by climate conditions, mortality level, and spatial heterogeneity in long-term aridity.

Figure 9 Relationships between long-term aridity, vegetation mortality level, and difference in ET for a dry year (1994, a) and a wet year (1995, b).

Water yield

In the dry year (1994), beetle-caused vegetation mortality affected water yield (Fig. 10), but the responses differed between the balanced and water-limited areas. For the balanced area, most sub-basins showed slight decreases in water yield after the beetle outbreak (Fig. 10a). The water yield in sub-basins with low vegetation mortality levels (≤40 %) did not change significantly; however, with increased mortality, more sub-basins showed increases in water yield, particularly the sub-basins with vegetation mortality levels higher than 40 % (Fig. 10a). Moreover, the vegetation mortality threshold that changed the direction of water yield response was affected by long-term aridity, e.g., it was 40 % for an aridity of 2.0 but 20 % for an aridity of 1.0. For the water-limited area, the water yield decreased and was independent of the mortality level (Fig. 10a). In the wet year (1995), the water yield in most sub-basins increased after the beetle outbreak, and the water yield in the balanced area increased more significantly than that in the water-limited area. Furthermore, in the balanced area, higher mortality levels caused larger increases in water yield, which responded more linearly (Fig. 10b). In summary, for a wet year, increases in water yield occurred for most sub-basins, driven by a decrease in ET. However, during dry years, the water yield and ET responses were spatially heterogeneous, and the competing changes in evaporation and transpiration changed the direction and magnitude of ET and thus the water yield response. The competing effects of different hydrologic fluxes for a dry year are explored in more detail in the next section.

Figure 10 Relationships between long-term aridity, vegetation mortality level, and difference in water yield for a dry year (1994, a) and a wet year (1995, b).

4.2.2 Analysis of water budgets to understand the decreasing water yield in the dry year

We analyzed the fluxes in a dry year (1994) in greater detail to understand the responses of the hydrologic fluxes and the resulting water yield. Based on Eq. (2), we identified four hydrologic fluxes that can potentially affect the water yield: canopy evaporation (canopy evaporation and canopy snow sublimation), ground evaporation (bare soil evaporation, ground snow sublimation, litter evaporation, pond evaporation), plant transpiration, and year-to-year storage change (soil, canopy, litter, snowpack). These three storage terms (canopy, litter, snowpack) were considered together with soil storage because their contribution was minor in comparison with other fluxes. Figure 11 summarizes different combinations of the four dominant processes during the dry year (1994) based on their directions (increase or decrease in water yield) after the beetle outbreak. In total, 14 combinations of changes in these fluxes (referred to as “response types”) were found. Five of those resulted in an increase in water yield; the others resulted in a decrease.

Figure 11 Water yield response types after the beetle outbreak for different evergreen mortality scenarios compared with the control scenario. D1–D9 are water yield decrease types and W1–W5 are water yield increase types. In panels (d) and (h), the left side of each type are increasing fluxes that cause water yield decreases, and the right side of each type are decreasing fluxes that cause water yield increases. If the left side is larger than the right side, the water yield increases, and vice versa. (Note: this mortality is evergreen mortality, which is different from vegetation mortality.)

Water yield responses caused by competition between hydrologic fluxes showed different patterns across the aridity gradient (Figs. 3 and 10). For the balanced area (upper part of the basin) with low evergreen mortality (≤30 %), the major response types were D1 and D2, in which the increase in ground evaporation dominated over the decreases in transpiration and canopy evaporation (Fig. 11a–c). However, with higher evergreen mortality (>30 %), the major response type became W2, where the increase in ground evaporation did not exceed the decrease in canopy evaporation and transpiration (Fig. 11e–g). This indicates that, in a dry year, when more evergreen stands are killed, the increase in ground evaporation reaches a limit while transpiration and canopy evaporation continue to decrease with decreasing LAI. The increase in ground evaporation was triggered either by decreased total LAI and open canopy, which allowed more solar radiation to penetrate to the ground for evaporation (Fig. S5c), or by less transpiration from plants, which left more water available to evaporate (Fig. 8a). The decreases in plant transpiration and canopy evaporation were driven by a lower live LAI and a lower total LAI, respectively (Figs. S5a, c and 8a).

The decrease in water yield in the water-limited area (i.e., the lower part of the basin) was driven by differences in how competing hydrologic responses interacted at different mortality levels. When evergreen stand mortality level was low (30 %), the response types were D5 and D7, in which the increase in ground and canopy evaporation dominated over the decrease in transpiration (Fig. 11a–c). However, with high evergreen stand mortality (>30 %), the response types became D1 and D2 (Fig. 11e–g), in which the canopy evaporation changed from an increase to a decrease that was driven by a decrease in the total LAI (Fig. S5c). When mortality was low, the increased growth of residual plants and the understory outstripped the litter fall of dead foliage, thus increasing the total LAI, and vice versa when mortality was high.

5 Discussion

5.1 Role of interannual climate variability

During the first 15 years after a beetle attack, various hydrologic processes oppose and/or reinforce one another to either increase or decrease the water yield: a decrease in live LAI can reduce transpiration, while a decrease in total LAI can enhance ground evaporation but diminish canopy evaporation (Montesi et al., 2004; Tsamir et al., 2019). Interannual climate variability plays an important role in determining which of these competing effects dominate, and therefore drives the direction of water yield response to beetle outbreak (Goeking and Tarboton, 2020; Winkler et al., 2014). Our results show that mainly decreases in water yield occurred in dry years, while the water yield increased in wet years. During a wet year, we found that plant ET reached its potential, so reductions in actual plant ET dominated over increases in ground evaporation, resulting in a net increase in water yield. During a dry year, the relative dominance of these competing effects had greater spatial heterogeneity because the water stress status of the plants varied across the basin (as explained in Sect. 4.2.2; Fig. 11).

However, the responses we observed in the dry year (1994) and in the wet year (1995) were also affected by the previous year's climate (mainly precipitation) and its effects on hydrologic and biogeochemical processes, which set the initial conditions for the dry and wet years (e.g., soil moisture and nitrogen availability). Therefore, we also analyzed other water years during the first 10 years after the beetle outbreak to examine whether our findings for dry and wet years followed a general pattern and to what extent they were influenced by antecedent conditions. Results indicate that our findings are robust throughout the study time period. For example, water yield generally decreased during dry years (1992, 1994, and 2001, see Figs. S1 and S2) and always increased during wet years (1993 and 1995–2000, see Figs. S1 and S2).

Adams et al. (2012) provide a threshold of precipitation under which the water yield increases after disturbances: at least 500 mm yr⁻¹. The average annual precipitation over this study basin was 600–900 mm in dry years and higher than 900 mm in wet years. Recent field observations also suggest that annual climate variability can affect the magnitudes of evapotranspiration fluxes, which have the potential to change the water yield direction (Biederman et al., 2014). Our results corroborate those earlier studies by revealing that there are precipitation thresholds above which tree removal increases water yield (Figs. 10, S1, and S2).

5.2 Role of vegetation mortality

Vegetation mortality is another important factor that influences water yield response. We found that during the wet year, beetle outbreak increased water yield across the basin, and the magnitude of these increases grew linearly with the level of vegetation mortality (Fig. 10b). In the dry year, however, the response of water yield to the level of vegetation mortality was more complicated because mortality level influenced not only the magnitude of change but also the direction (Fig. 10a). These opposing results (due to mortality level) mainly occurred in the balanced northern part of the basin, where the competing effects of mortality (i.e., increases in ground evaporation versus decreases in transpiration) were more balanced (Fig. 11). The level of vegetation mortality played a less significant role in water yield changes in the southern water-limited area. Vegetation mortality level determined the magnitudes of the live LAI, total LAI, transpiration, canopy evaporation, and ground evaporation in such a way that it governed the directions of change in both ET and water yield. Thus, when the vegetation mortality level was higher than 40 %, its capacity to reduce transpiration became the dominant effect, whereas its capacity to increase soil evaporation became a minor effect (Fig. 11f and g; Guardiola-Claramonte et al., 2011).

Besides the precipitation threshold of at least 500 mm yr⁻¹, Adams et al. (2012) also estimate that when at least 20 % of the vegetation cover is removed, the water yield can increase. According to our previous analysis (Sect. 4.1), for a dry year, water yield increases when more than 40 % of the vegetation is removed (Fig. 10a). Our model simulations indicate that there are similar mortality thresholds that drive water yield increases during the dry year, but we did not find evidence of such a threshold during wet years. These differences between dry and wet years suggest that the effects of mortality on water yield depend on climate variability. Similarly, other studies demonstrate that the relationship between mortality level and water yield response is complicated and nonlinear (Moore and Wondzell, 2005).

5.3 Role of the long-term aridity index (PET/P)

Long-term aridity indices can be used to predict where the water yield will decrease after a disturbance. We found that water yield always increased in a wet year, irrespective of the climatic aridity index (Fig. 10a). For dry years, the long-term aridity index was important in driving the direction of water yield responses to a beetle outbreak. In areas that were less water limited (balanced areas), the direction of water yield responses to a beetle outbreak in a dry year was mixed and depended on mortality level. For water-limited areas, in a dry year, water yield showed a more consistent decrease and was also less affected by mortality level. These results agree with previous studies which found that water yield decreases largely happen in semi-arid areas (Biederman et al., 2014; Guardiola-Claramonte et al., 2011).

The decrease in water yield for a water-limited area can be driven by increases in canopy evaporation and transpiration, in contrast to a hydrologically balanced area, where it was driven by increased ground evaporation. In the water-limited area, the increase in canopy evaporation was due to an increase in total LAI, which is a combined effect of delayed decay of dead foliage and fast growth of residual and understory plants (Fig. 11d, types D5, D7, D8, and D9; Fig. S5). The surviving and understory plants in the water-limited area also had higher transpiration rates after mortality (Fig. 11d, type D6; Fig. 8). Similarly, in field studies, Tsamir et al. (2019) found that photosynthesis and transpiration increased after thinning in a semi-arid forest. These findings illustrate that, in addition to top-down climate variability, the long-term aridity index (which also varies with bottom-up drivers such as vegetation and local topography) can be a useful indicator of how water yield will respond to disturbances.

In addition to evaporation and transpiration, snow sublimation can also influence the direction of hydrologic response. Similar to other process-based snow models, we found that once dead foliage fell to the ground, canopy sublimation decreased (e.g., Koeniger et al., 2008; Sexstone et al., 2018), which in turn increased the water yield relative to the period when dead needles remained on the trees (Figs. 5 and S7). In water-limited regions, the decrease in canopy sublimation was much smaller than in the balanced regions because there were smaller changes in total LAI (Fig. S5c and d). However, immediately after beetle outbreak (e.g., 1990–1992), we found that canopy sublimation increased in both regions due to an increase in total LAI (Fig. S7). This finding is supported by observational studies showing that canopy sublimation can increase with increasing leaf area (Koeniger et al., 2008).

We also found that ground/snowpack sublimation decreased when all dead foliage fell to the ground because snowmelt increased with the opening of the canopy. However, this finding differs from that of other studies suggesting that snowpack sublimation can increase with a more open canopy (Biederman et al., 2014; Harpold et al., 2014). The latter can occur because open canopies may allow more snow to reach the ground, which can increase sublimation. However, in our study, faster snowmelt appeared to dominate over increases in ground sublimation. These contrasts between our research and previous studies illustrate a sophisticated balance between canopy–atmosphere–environmental processes that must be accounted for when studying the sublimation response to disturbances (Edburg et al., 2012; Frank et al., 2019). Although RHESSys is a powerful tool for representing these complex interactions, some process representations warrant further analysis. For example, RHESSys currently ignores the effects of litter on ground albedo and snowmelt (Lundquist et al., 2013), which could affect actual ET and PET rates and therefore the long-term aridity index.

5.4 Uncertainties and recommendations for future research

We found the long-term (38 year) aridity index for our study region to be a key driver influencing hydrologic responses to beetle outbreaks. While this trend is likely to continue in the future as climate change intensifies aridity in the western US (Livneh and Badger, 2020), the classification of water-limited/balanced regions based on the 38 year aridity index may change. Thus, projecting how responses will change under future aridity scenarios requires further modeling research. We used historical 38 year (1980–2018) data to calculate the aridity index (PET/P). This method can be extended to project future responses to beetle outbreaks by using future climate data from generalized circulation models (GCMs) to drive the process-based ecohydrologic beetle effects model. Another consideration, however, is that as aridity continues to increase, vegetation may shift from evergreen to more drought-tolerant shrub or grass species. This would in turn alter beetle outbreak patterns and the corresponding water yield responses (Abatzoglou and Kolden, 2011; Bart et al., 2016). However, this process is not well understood and is not currently represented in our modeling framework. A key uncertainty in predicting future beetle effects is how vegetation will respond to climate change.

Another key uncertainty is how beetle attacks will change in the future. We used uniform mortality levels for all patches across the watershed and focused our analyses on potential beetle effects. However, in reality, beetles usually attack older trees first (Edburg et al., 2011). Thus, incorporating a more mechanistic understanding of beetle attack patterns with our beetle effects model could enable us to simulate more realistic outbreak scenarios moving forward. We also focused on water yield responses during the first 15 years after a beetle outbreak in a watershed that contained balanced or water-limited sub-basins. Future research should analyze long-term effects (e.g., after the ecosystem begins to recover) on forest hydrology and investigate wetter, energy-limited regions.

6 Conclusion

We tested a coupled ecohydrologic and beetle effects model of a semi-arid basin in southern Idaho to examine how the watershed hydrology responds to a beetle outbreak and how interannual climatic variability, vegetation mortality, and long-term aridity influence these responses. Simulation results indicate that each factor can play a discrete role in driving hydrologic processes (e.g., the directions and magnitudes of changes in plant transpiration, canopy and soil evaporation, soil and litter moisture, and snow sublimation). These combined effects determine the overall water budget and water yield of the basin. While interannual climate variability is the key factor driving the direction of change in water yield, vegetation mortality levels and long-term aridity modify water yield responses.

In dry years, the water yields of most sub-basins slightly decreased after the beetle outbreak when vegetation mortality was lower than 40 %, while they increased during wet years in most sub-basins. Our results show that long-term aridity index is a reliable indicator of the water yield decreases that occur during dry years due to the fact that there is a consistent decrease in water yield in the most water-limited portion of the basin. Generally, the effect of vegetation mortality on water yield during dry years is less uniform and depends on local aridity. During wet years, on the other hand, mortality typically causes increases in water yield. This illustrates that, together, interannual climate variability and mortality can have a stronger effect on the direction of water yield response in water-limited regions than interannual climate variability alone. Future studies to predict water yield response to disturbance should consider the interactions of these factors and capture the fluctuations of competing water fluxes and storage changes that control the overall water budget and water yield.

Using our novel RHESSys–beetle effects modeling framework, we have demonstrated that the direction of hydrologic response is a function of multiple factors (e.g., interannual climate variability, vegetation mortality level, and long-term aridity) and that these results do not necessarily conflict with each other but are representative of different conditions. The mechanisms behind these changes compete with each other, resulting in water yield increases or decreases (Fig. 1). Contradictory findings in previous studies may result from differing mortality levels or differences in aridity because the emergent drivers that dominate water yield responses differ. Disentangling these drivers is difficult or impossible using a purely empirical approach, where it can be challenging or cost prohibitive to experiment under a broad range of controlled conditions. Distributed process-based models, on the other hand, provide a useful tool for examining these dynamics.

This study can assist water supply stakeholders in risk management in beetle outbreak locations. For example, during wet years, more attention might be focused on balanced areas, i.e., wet regions, regarding flooding and erosion risks after beetle outbreaks, since these regions may experience large increases in runoff due to decreases in plant transpiration and increases in soil moisture. During dry years, attention might need to shift to water-limited areas for managing wildfire risk, since these regions will experience elevated ET and lower soil and litter moisture. Because multiple factors interact to influence hydrologic processes after a beetle outbreak, water and forest management must respond to spatial and temporal variations in climate, aridity, and vegetation mortality levels.