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

Abstract 36 Mountain pine beetle (MPB) outbreaks in western United States result in widespread tree 37 mortality, transforming forest structure within watersheds. While there is evidence that these 38 changes can alter the timing and quantity of streamflow, there is substantial variation in both the 39 magnitude and direction of responses and the climatic and environmental mechanisms driving 40 this variation are not well understood. Herein, we coupled an eco-hydrologic model (RHESSys) 41 with a beetle effects model and applied it to a semiarid watershed, Trail Creek, in the Bigwood 42 River basin in central Idaho to evaluate how varying degrees of beetle-caused tree mortality 43 influence water yield. Simulation results show that water yield during the first 15 years after 44 beetle outbreak is controlled by interactions among interannual climate variability, the extent of 45 vegetation mortality, and long-term aridity. During wet years, water yield after beetle outbreak 46 increases with greater tree mortality. During dry years, water yield decreases at low to medium 47 mortality but increases at high mortality. The mortality threshold for the direction of change is 48 location-specific. The change in water yield also varies spatially along aridity gradients during 49 dry years. In relatively wetter areas of the Trail Creek basin, water yield switches from a 50 decrease to an increase when vegetation mortality is greater than 40 percent. In more water- 51 limited areas on the other hand, water yield typically decreases after beetle outbreaks, regardless 52 of mortality level. Results suggest that long-term aridity can be a useful indicator for the 53 direction of water yield changes after disturbance.

climate variability is an important driver of water yield response to beetle-caused tree 23 mortality. 24  A long-term (multi-decade) aridity index is a reliable indicator of water yield response to 25 MPBs: in a dry year, decreases occur mainly in "water-limited" areas and vegetation 26 mortality levels have only minor effects; in wetter areas, decreases only occur at low 27 mortality levels. 28  Generally, in a dry year, low to medium MPB-caused vegetation mortality decreases 29 water yield, and high mortality increases water yield; this response to mortality level is In a review of 78 studies, Goeking and Tarboton (2020) concluded that the decrease in water 84 yield after tree-mortality mainly happens in semiarid regions. Previous studies also provide rule-85 of-thumb thresholds above which water yield will increase: at least 20 percent loss of vegetation 86 cover and mean precipitation of 500 mm/year (Adams et al. 2012). However, many watersheds 87 in the western U.S. experience high interannual climate variability (Fyfe et al. 2017), and local 88 environmental gradients (e.g., long-term aridity gradients) may strongly influence vegetation and 89 hydrologic responses to disturbances, including beetle outbreaks, making predictions difficult 90 (Winkler et al. 2014). Given the possibility of either increases or decreases in water yield 91 following beetle outbreaks, modeling approaches are crucial for identifying the specific 92 mechanisms that control these responses. 93 The overarching goal of this study is to identify mechanisms driving the direction of change in 94 annual water yield after beetle outbreaks in semi-arid regions (note that in the following text, 95 "water yield" refers to means annual water yield). The following specific questions address this 96 goal: 97  Q1: What is the role of interannual climate variability in water yield response? 98  Q2: What is the role of mortality level in water yield response? 99  Q3: How does long-term aridity (defined as temporally averaged potential 100 evapotranspiration relative to precipitation for a period of 38 years) modify these 101 responses, and how do responses vary spatially within a watershed along aridity 102 gradients? 103 We hypothesize that multiple ecohydrologic processes (e.g., snow accumulation and melt, 104 evaporation, transpiration, drainage, and a range of forest structural and functional responses to 105 beetles) could interactively influence how water yield responds to beetle outbreaks-however, in 106 certain locations one or more processes may dominate. In addition, the dominant ecohydrologic 107 processes may vary over space and time due to interannual climate variability (i.e., 108 precipitation), vegetation mortality, and long-term aridity. In Sect 2, we present a conceptual 109 framework for identifying and depicting dominant hydrological processes through which forests 110 respond to beetle infestation. We use this framework to interpret the modeling results. In Sect 3, 111 we describe our mechanistic modeling approach, i.e., using the Regional Hydro-Ecological 112 Simulation System (RHESSys), which can prescribe a range of vegetation mortality levels, 113 capture the effects of landscape heterogeneity and the role of lateral soil moisture redistribution, 114 and project ecosystem carbon and nitrogen dynamics, including post-disturbance plant recovery. 115 In Sects 4 and 5, we then present modeling results that explore how multiple mechanisms 116 influence water yield responses. 117 2 Conceptual framework 118 2.1 Vegetation response to beetle outbreaks 119 Mountain pine beetles (MPB) introduce blue stain fungi into the xylem of attacked trees, which 120 reduces water transport in plants and eventually shuts it off (Paine et al. 1997 Model (Bennett et al. 2018) also have some limitations. For example, 1) they may assume 167 constant LAI after disturbances and static vegetation growth (e.g., VIC and DHSVM), 2) they 168 may not include lateral flow to redistribute soil moisture (VIC), and 3) in some cases, the 169 approach to represent the effects of beetle outbreaks may be too simplified (e.g., changing only 170 LAI and conductance without considering two-way beetle-vegetation interactions in post-171 disturbance biogeochemical and water cycling e.g., as in CLM-ParFlow). Thus, improving 172 current fully distributed process-based models to capture the coupled dynamics between 173 hydrology and vegetation at multiple scales is a critical step for projecting how beetle outbreaks 174 will affect water yield in semiarid systems (Goeking and Tarboton 2020). Here we use 175 RHESSys7.1, which captures these processes. Trail Creek has cold, wet winters and warm, dry summers; mean annual precipitation is 185 approximately 978 mm with 60% snow (Frenzel 1989). The soil is mostly permeable coarse 186 alluvium (Smith 1960). Vegetation is clustered into two major groups along the elevation which 187 ranges from 1760 to 3478 m: sagebrush, riparian species, and grasslands in lower to middle 188 elevation areas and Douglas-fir (Pseudotsuga menziesii), lodgepole pine (Pinus contorta var. 189 latifolia), subalpine fir (Abies lasiocarpa), and mixed shrub and herbaceous vegetation in 190 middle-to-higher elevations (Buhidar 2002). 191 A strong upper to lower vegetation and long-term aridity gradient exists for Trail Creek (Fig. 3). 192 The northern (higher elevation) portion of the basin is mesic and covered principally by 193 evergreen forest; the southern (lower elevation) portion is xeric and covered by shrubs, grasses, 194 and mixed herbaceous species. In total, Trail creek contains 72 sub-basins and two of them (e.g., 195 and high elevation area is balanced (i.e., PET/P between 0.8 and 2) and evergreen tree coverage 199 is more than 50%; the southern part is water-limited (i.e., PET/P > 2) and evergreen tree 200 coverage is less than 30% (Figs. 2 and 3).  Fig. 4). Here we integrated this beetle effects model into RHESSys (Fig. 4). Beetles attack 257 trees mainly during late summer, and needles will turn from green to red at the beginning of the 258 following summer. We simplify this process with prescribed tree mortality on September 1 to 259 represent a beetle outbreak of the current year. The advantage of this integration is that RHESSys 260 accounts for the lateral connectivity in water and nitrogen fluxes among patches which is not Climate inputs for this study, including maximum and minimum temperatures, precipitation, 301 relative humidity, radiation, and wind speed, were acquired from gridMET for years from 1980 302 to 2018. GridMET provides daily high-resolution (1/24 degree or ~4 km) gridded meteorological 303 data (Abatzoglou 2013). It is a blended climate dataset that combines the temporal attributes of 304 gauge-based precipitation data from NLDAS-2 (Mitchell et al. 2004) with the spatial attributes of 305 gridded climate data from PRISM (Daly et al. 1994). 306

Simulation experiments 307
To quantify how water yield responds to beetle-caused mortality, we designed the following 308 simulation experiment. We prescribed a beetle outbreak in September 1989, the mortality level 309 (%) is applied to all evergreen patches for each sub-basin. After beetle outbreak, red needles stay 310 on the trees for one year before they start to fall (transferred to the litter pool) at an exponential 311 rate with a half-life of two years. The snag pools stay in the standing trees for five years and then 312 start to fall and are added to the CWD pool which decays at an exponential rate with a half-life 313 of ten years. To address Q2 (i.e., the role of vegetation mortality), we prescribe a range of infestation-caused 333 mortality levels (i.e., from 10% to 60% by a step of 10% in terms of carbon, uniformly applied to 334 all evergreen patches for each sub-basins ) and a control run (no mortality) to quantify the 335 response of forests in water yield to vegetation mortality level (for each sub-basin vegetation 336 mortality is evergreen mortality multiplied by evergreen coverage of that basin). The 337 differences in water yield between each mortality level and the control run represent the effects 338 of beetle kill: a positive value means that mortality increased water yield, and vice versa. 339 We quantified the water budget for each sub-basin to examine which hydrological process dropped immediately after beetle outbreak, then gradually recovered to pre-outbreak levels 366 during following years (Fig. 5a). Total LAI (i.e., including dead foliage) showed a slight increase 367 during the first ten years after beetle outbreak (1990 -2000), which is due to the retention of 368 dead leaves in the canopy and the simultaneous growth of residual (unaffected) overstory and 369 understory vegetation (Fig. 5b). The dead foliage pool (Fig. 5c)  outbreak, scenarios where the evergreen mortality level was larger than zero had higher basin-379 scale water yield than the control scenario (where the evergreen mortality level was zero). This 380 was especially true during wet years; however, there was no significant increase during dry years 381 (i. e., 1992, 1994, 2001, and 2004; Fig. 6a). The year-to-year soil storage fluxes responded 382 strongly in the first two years after beetle outbreak, then stabilized to the pre-outbreak condition 383 (Fig. 6b). Note that year-to-year soil storage change is not the same as soil water storage. After 384 beetle outbreak, the soil can hold some portion of water that not being up taken by the plants, but 385 it was confined by the soil water holding capacity. This phenomenon indicates that the soil has 386 some resilience to vegetation change. 387 Beetle outbreaks reduced transpiration during wet years but did not have significant effects in 388 dry years (Fig. 6c). This is because transpiration in dry years was water-limited and so was much 389 lower than the potential rate (more water is partitioned to evaporation; Biederman et al. 2014). 390 Thus, killing more trees had little effect on stand scale transpiration because remaining trees 391 utilized any water released by the dead trees in dry years. On the other hand, plant transpiration 392 in wet years was close to the potential rate; therefore, decreases in canopy cover reduced 393 transpiration. The simulation results did not show any apparent effect on snowmelt after beetle 394 outbreak. 395 The evaporation response was opposite in dry and wet years: evaporation increased in dry years, 396 while it decreased in wet years (Fig. 6d). This phenomenon is caused by tradeoffs and 397 interactions among multiple processes, as will be explained in more detail in the next section. Beetle outbreak had opposite effects on evaporation between a dry year and a wet year (Fig. 7). 402 In the dry year, most sub-basins experienced higher evaporation for beetle outbreak scenarios 403 than in the control scenario (Fig. 7a). This was the cumulative consequence of decreased canopy 404 evaporation and increased ground (soil, litter, pond) evaporation due to decreases in LAI (caused 405 by mortality). In the dry year, the latter effect (i.e., increased ground evaporation) dominated 406 Beetle outbreak decreased transpiration in both dry and wet years, and with higher mortality 419 levels the decrease became larger (Fig. 8). However, during the dry year, the water-limited area 420 showed less change than the balanced area; some sub-basins even showed slight increases. This  In the dry year (1994), beetle-caused vegetation mortality affected water yield (Fig. 10), but the 440 responses differed between the balanced and water-limited areas. For the balanced area, most 441 sub-basins showed slight decreases in water yield after beetle outbreak and no significant 442 differences among low vegetation mortality level (<=40%, Fig. 10a). However, with increased 443 mortality levels, more sub-basins showed increases in water yield, particularly with vegetation 444 mortality higher than 40% (Fig. 10a). Moreover, the vegetation mortality threshold that changed 445 the direction of water yield response was altered by long-term aridity, e.g., it was 40% for aridity 446 2.0 but 20% for aridity 1.0. For the water-limited area, water yield decreased and was 447 independent from mortality level (Fig. 10a). In the wet year (1995), the water yield in most sub-448 basins increased after beetle outbreak, and the balanced area increased more significantly than 449 the water-limited area. Furthermore, for the balanced area, higher mortality levels caused larger 450 increases in water yield which responded more linearly (Fig. 10b). In summary, for a wet year, 451 increases in water yield occurred for most sub-basins, driven by a decrease in ET. However, 452 during dry years, the water yield and ET responses were spatially heterogeneous, and the 453 competing changes in evaporation and transpiration changed the direction and magnitude of ET 454 and thus water yield response. The competing effect among different hydrologic fluxes for a dry 455 year is explored in more detail in the next section. with low evergreen mortality (<=30%), the major response types were D1 and D2, in which the 471 increase in ground evaporation dominated over the decrease in transpiration and canopy 472 evaporation (Fig. 11a, b, and c). However, with higher evergreen mortality (>30%), the major 473 response type became W2, where the increase in ground evaporation did not exceed the decrease 474 in canopy evaporation and transpiration (Fig. 11e, f, and g). This indicates that, in a dry year, 475 when more evergreen stands are killed, the increase in ground evaporation reaches a limit while 476 transpiration and canopy evaporation continue to decrease with decreasing LAI. The increase in 477 ground evaporation was triggered either by decreased Total LAI and open canopy, which 478 allowed more solar radiation penetration to the ground for evaporation (Fig. S5c), or less 479 transpiration from plants, which left more water available to evaporate (Fig. 8a). The decrease in 480 plant transpiration and canopy evaporation was driven by a lower Live LAI and a lower Total 481 LAI, respectively (Fig. S5 a&c and Fig. 8a). 482 The decrease in water yield in the water-limited area (lower part of the basin) was driven by 483 different hydrologic flux competitions in different mortality levels. When evergreen stand 484 mortality level was low (<=30%), the response types were D5 and D7, in which the increase in 485 ground and canopy evaporation dominated over the decrease of transpiration (Fig. 11a, b, and c). 486 However, with high evergreen stand mortality (>30%), the response types became D1 and D2 487 ( Fig. 11e, f, and g), in which the canopy evaporation changed from an increase to a decrease that 488 was driven by a decrease in Total LAI (Fig. S5c). When mortality was low, the increases in 489 growth from residual plants and understory outstripped the litter fall of dead foliage; thus, Total 490 LAI increased, and vice versa when mortality was high. 491 Goeking and Tarboton 2020). Our results show that mainly decreases in water yield occurred in 500 dry years, while increases occurred in wet years. During a wet year, plant ET can reach its 501 potential so that any reductions in actual plant ET will dominate over any increases in ground 502 evaporation, resulting in a net increase in water yield. During a dry year, the relative dominance 503 of these competing effects had greater spatial heterogeneity because the water stress status of the 504 plants varied across the basin (as explained in Sect 4.2.2; Fig. 11). 505 However, the responses we observed in the dry year (1994) and in the wet year (1995) were also 506 affected by the previous year's climate (mainly precipitation) and its effects on hydrologic and 507 biogeochemical processes, which set the initial conditions for the dry and wet year (e.g., soil 508 moisture, nitrogen availability, etc.). Therefore, we also analyzed other water years during the 509 first ten years after beetle outbreak to examine whether our findings for dry and wet years follow 510 a general pattern and to what extent they are influenced by antecedent conditions. Results 511 indicate that our findings are robust through the study time period. For example, water yield 512 generally decreased during dry years (1992,1994 results corroborate these earlier studies by revealing that there are precipitation thresholds above 520 which tree removal increases water yield (Figs. 10, S1 and S2). 521

Role of vegetation mortality 522
Vegetation mortality is another important factor that influences water yield response. We found 523 that during the wet year, beetle outbreak increased water yield across the basin and the 524 magnitude of these increases grew linearly with the level of vegetation mortality (Fig. 10b). In 525 the dry year, however, the response of water yield to the level of vegetation mortality was more 526 complicated because mortality influenced not only the magnitude of change but also the 527 direction (Fig. 10a). These opposing results (due to mortality level) mainly occurred in the 528 "balanced" northern part of the basin, where the competing effects of mortality (i.e., increases in 529 ground evaporation versus decreases in transpiration) are more balanced (Fig. 11). The level of 530 vegetation mortality played a less significant role in changing water yield in the southern "water-531 limited" area. Vegetation mortality level determined the magnitudes of Live LAI, Total LAI, 532 transpiration, canopy evaporation, and ground evaporation in such a way that it governed the 533 direction of change in both ET and water yield. Thus, when vegetation mortality level was higher 534 than 40%, its effect of decreasing transpiration became the dominant process and its effect of 535 increasing soil evaporation became minor (Fig. 11  that when at least 20% of vegetation cover is removed, water yield can increase. According to 538 previous analysis (Sect 4.1), for a dry year, water yield increases when more than 40% of 539 vegetation is removed (Fig. 10a). Our model simulations indicate similar mortality thresholds 540 exist for driving water yield increases during the dry year, however, we did not find evidence 541 that such a threshold exists during wet years. These differences between dry and wet years 542 suggest that the effects of mortality on water yield depend on climate variability. Other studies 543 corroborate this finding by demonstrating that the relationship between mortality level and water 544 yield response is complicated and nonlinear (Moore and Wondzell 2005). 545

Role of long-term aridity index (PET/P) 546
Long-term aridity indices can be used to predict where water yield will decrease after 547 disturbance. We found that water yield always increased in a wet year, irrespective of the 548 climatic aridity index (Fig. 10a). For dry years, long-term aridity index became important in 549 driving the direction of water yield responses to beetle outbreak. In areas that are less water-550 limited (balanced areas), the direction of water-yield responses to beetle outbreak in a dry year 551 was mixed and depended on mortality level. For water-limited areas, in a dry year, water yield 552 showed a more consistent decrease trend, and it was also less affected by mortality level. These The decrease in water yield for water-limited area can be driven by increases in canopy 556 evaporation or transpiration, which were different in the hydrologically-balanced area (driven by 557 increase of ground evaporation). There, the increase in canopy evaporation was due to an 558 increase in total LAI which is a combined effect of delayed decay of dead foliage and fast 559 growth of residual and understory plants (Fig. 11d type D5 found an increase in photosynthesis and transpiration after thinning in a semi-arid forest. These 563 findings illustrate that in addition to top-down climate variability, the long-term aridity index 564 (which also varies with bottom-up drivers such as vegetation and local topography) can be 565 another useful indicator of how water yield will respond to disturbances. 566

Uncertainties 567
While our findings revealed how topoclimatic gradients influenced water yield responses to 568 beetle infestation, some uncertainties remain. For one, we used uniform mortality levels for all 569 patches across the watershed rather than location and vegetation-specific mortality levels. 570 However, in reality beetles usually attack older trees first (Edburg et al. 2011). Thus, 571 incorporating a more mechanistic understanding of beetle attack patterns with our beetle effects 572 model could enable us to simulate more realistic outbreak scenarios moving forward. Another 573 source of uncertainty came from the model treatment of litter pools. In the current 574 implementation, we ignored the effects of litter on ground albedo and snowmelt (Lundquist et al. 575 2013), which could have an effect on rates of AET and PET and therefore our calculated long-576 term aridity index. Also, because we focused on water yield responses during the first 15 years 577 after beetle outbreak, we may have missed some of the long-term effects (e.g., after the 578 ecosystem has begun to recover) on forest hydrology. Future research should integrate the short-579 term and long-term effects and interactions among beetle outbreak, vegetation dynamics, and 580 hydrology. Since Trail Creek is either "balanced" or "water-limited" in terms of aridity, other 581 "energy-limited" regions could also be investigated. 582

Conclusion 583
We tested a coupled ecohydrologic and beetle effects model in a semi-arid basin in southern 584 Idaho to examine how watershed hydrology responds to beetle outbreak and how interannual 585 climatic variability, vegetation mortality, and long-term aridity influence these responses. 586 Simulation results indicate that each factor can play a discrete role in driving hydrological 587 processes (e.g., the direction and magnitude of changes in plant transpiration, canopy and soil 588 evaporation, soil and litter moisture, snow sublimation, etc.). These combined effects determine 589 the overall water budget and water yield of the basin. While interannual climate variability is the 590 key factor driving the direction of change in water yield, vegetation mortality levels and long-591 term aridity modify water yield responses. 592 In dry years, the water yield of most sub-basins slightly decreased after beetle outbreak when 593 vegetation mortality level was lower than 40%; while during wet years in most sub-basins it 594 increased. Our results show that long-term aridity index is a reliable indicator of the water yield 595 decreases that occur during dry years due to the fact that there is a consistent decrease in water 596 yield in the most water-limited portion of the basin. Generally, the effects of vegetation mortality 597 on water yield during dry years is less uniform and depends on local, long-term aridity 598 conditions. During wet years, on the other hand, mortality typically causes increases in water 599 yield. This illustrates that together interannual climate variability and mortality can have a 600 stronger effect on the direction of water yield response in water-limited regions than interannual 601 climate variability alone. Future studies to predict water yield response to disturbance should 602 consider the interactions of these factors and capture the fluctuations of competing water fluxes 603 and storage change that control overall water budget and water yield. 604 Using our novel RHESSys-beetle effects modeling framework, we demonstrate that the direction 605 of hydrologic response is a function of multiple factors (e.g., interannual climate variability, 606 vegetation mortality level, and long-term aridity) and that these results do not necessarily conflict 607 with each other but are representative of different conditions. The mechanisms behind these 608 changes compete with each other resulting in a water yield increases or decreases (Fig. 1). 609 Contradictory findings in previous studies may result from differing mortality levels (disturbance 610 severity), or differences in aridity, and consequently, the emergent drivers that dominate water 611 yield responses differ. Disentangling these drivers is difficult or impossible using a purely 612 empirical approach where it can be challenging or cost-prohibitive to experiment under a broad 613 range of controlled conditions. Distributed process-based models on the other hand, provide a 614 useful tool for examining these dynamics. 615 Findings from this study can assist water supply stakeholders in risk management in beetle 616 outbreak locations. For example, during wet years, more attention might be focused on 617 "balanced" areas, i.e., wet regions, for flooding and erosion risks after beetle outbreaks since 618 these regions may experience large increase in runoff due to decreases in plant transpiration and 619 increases in soil moisture. During the dry years, attention might need to shift to "water-limited" 620 areas for managing wildfire risk since these regions will experience elevated ET and lower soil 621 and litter moisture. Because multiple factors interact to influence hydrological processes after 622 beetle outbreak, water and forests management must respond to spatial and temporal variations 623 in climate, aridity, and vegetation mortality levels. 624     The green background color is the period before beetle outbreak, and the red background color 960 is after the beetle outbreak. 961 962 Long-term aridity is defined as temporally averaged (38 years) potential evapotranspiration 978 relative to precipitation. 979 980 water yield for a dry year (1994, a) and wet year (1995, b