Forest Treatment Effects on Watershed Responses Under Warming

The effects of forest treatments on watershed hydrology have often been studied in isolation from climate change. Consequently, under a warming climate, it is unclear how forest thinning will impact snowpacks, evapotranspiration, and streamflow availability. In this study, we used a distributed hydrologic model to provide insight into the effects of warming and forest treatment on the hydrologic response of the Beaver Creek watershed (∼1,100 km2) of central Arizona. Prior to the numerical experiments, confidence in the hydrologic model performance was established by comparisons to long‐term observations (2003–2018) of snow water equivalent and streamflow using station observations and through spatially distributed estimates. Results indicated that warming during the 21st century could increase mean annual streamflow by 1.5% for warming levels up to +1°C, followed by a −29% decrease for continued warming up to +6°C, due to the varying effects of warming on snow sublimation, soil evaporation, and plant transpiration. On average, forest thinning increased streamflow by +12% (or 7 mm/yr) through lower plant transpiration by −19% (or −18 mm/yr), while also increasing the change in soil water storage by +42% (or 11 mm/yr). Forest thinning delayed the detrimental effects of warming on streamflow until +4°C, as compared to +2°C without forest treatment. Furthermore, model results suggested that forest cover reductions laterally displace water availability and evapotranspiration to downstream sites. These model‐derived mechanisms provide insights on the potential for water resilience toward warming effects afforded through treatment projects in southwestern US ponderosa pine forests.


Introduction
Forested landscapes are critical for providing ecosystem services, including freshwater, carbon sequestration and storage, soil conservation, and habitat for different species (e.g., Foley et al., 2005;Trumbore et al., 2015).Often located in headwater regions, forests supply water to downstream agricultural and urbanized areas that depend on these landscapes for socio-economic activities (Mueller et al., 2013;Simonit et al., 2015).Because of this importance, the relationship between forest treatments and water supply has been the subject of many experiments (e.g., Baker, 1986;Bosch & Hewlett, 1982;Rich, 1972).Prior field efforts have investigated how forest disturbances and treatments can affect water supply (see, for instance, the reviews of Goeking and Tarboton (2020) and del Campo et al. (2022)), revealing that a set of complex feedbacks exists between the forest condition, its climatic setting, and the underlying hydrologic processes.In semiarid watersheds of the southwestern US, the hydrologic response to forest disturbances have been mixed, with negligible, positive, and negative changes in streamflow all noted across different locations (e.g., Baker, 1986;Biederman et al., 2015;Guardiola-Claramonte et al., 2011;Ren et al., 2021;Rich & Gottfried, 1976).Furthermore, the effects of forest treatments (e.g., prescribed burning, thinning, grazing) on watershed hydrology have often been studied in isolation from global change stressors (i.e., warming, precipitation shifts), though recent efforts have begun to address these compound effects (e.g., Bennett et al., 2018;Robles et al., 2017;Sexstone et al., 2018;Whitney, Vivoni, Wang, et al., 2023).It is critical to address these compounding factors since forest treatment measures are being considered in the context of climate change mitigation and warming conditions will be superimposed on treatment projects spanning years to decades.
Ponderosa pines (Pinus ponderosa) are the most expansive forest ecosystems in the southwestern US (Garrison et al., 1977) and have exhibited an increase in tree density since the late 1800s due to forest management activities such as fire suppression (Allen et al., 2002;Covington & Moore, 1994).Under current conditions, ponderosa pine landscapes are susceptible to large crown fires (e.g., Hessburg et al., 2005;Savage & Mast, 2005) that lead to high tree mortality and impact nearby and downstream socio-economic activities.For example, forests in Arizona have recently experienced several destructive crown fires (Margolis et al., 2011;McHugh & Kolb, 2003) which have become more frequent, larger, and severe over time (Singleton et al., 2019).As a result, there has been an increasing interest in forest treatment to reduce stand density and decrease the risks of wildfire and droughtinduced mortality (Fulé et al., 1997;McCauley et al., 2022;Stephens et al., 2009).While forest fires lead to complex hydrologic impacts (e.g., Biederman et al., 2022;Goeking & Tarboton, 2020;Huffman et al., 2001;Ice et al., 2004;Leonard et al., 2017), it is also important to determine how proposed forest treatments intended for fire risk reduction might affect hydrologic conditions.After forest thinning, ponderosa pine areas in the southwestern US often exhibit reduced evaporative losses (Dore et al., 2010;Svoma, 2017), increased snowpacks (Sankey et al., 2015), higher streamflow (Baker, 1986;Kaye et al., 1999), and enhanced recharge (Schenk et al., 2020), among other effects.
Local climate interacts with ponderosa pine forests and their structure to determine cold season processes of snow accumulation, sublimation, and melt on the ground surface and on tree canopies (e.g., Harpold et al., 2015;Rinehart et al., 2008;Sun et al., 2022).Climate warming, in particular, has been documented as an overriding factor for changes in snow cover, magnitude, and timing in the forested mountains of the western US (e.g., Clow, 2010;Das et al., 2011;Ficklin et al., 2013;Harpold et al., 2012;McCabe et al., 2017;Milly & Dunne, 2020;Mote et al., 2005;Musselman et al., 2017;Robles et al., 2014;Xiao et al., 2018).In forests of central Arizona, Robles et al. (2017) found that recent warming trends  of +1 to +2.5°C resulted in a shift to an earlier streamflow in the spring and argued that earlier snow melt could be an explanation.In the same region and over a similar period, Svoma (2011) found that the snow level (i.e., altitude that rainfall transitions to snowfall) increased in elevation, indicating a long-term snow cover decline related to warming.Mountainous areas in central Arizona are often characterized with shallow, ephemeral, and discontinuous snowpacks (e.g., Broxton et al., 2020;Sankey et al., 2015).In addition, snowpacks are susceptible to earlier melting, as compared to higher latitudes, due to rainon-snow events and warm periods during the winter (e.g., Mascaro et al., 2023;Robles et al., 2021).Despite the ample evidence for hydrologic changes related to warming in the southwestern US, the combined effects of forest treatments and higher air temperatures on cold-season vegetation-hydrology interactions are still unknown.
A useful approach to quantify compounding effects on the watershed response is through numerical modeling that incorporates the relevant vegetation and hydrologic processes (e.g., Moreno et al., 2016;Sexstone et al., 2018;Wang & Vivoni, 2022).Process-based, distributed hydrologic models are applicable due to their ability to: (a) resolve terrain, soil, vegetation, and channel features at high resolution and their effects on hydrologic dynamics, (b) capture the differential impact of elevation on meteorological forcing during the cold and warm seasons, and (c) account for the carryover of snow, soil moisture, and groundwater storage from seasonal to interannual periods.Furthermore, process-based models at high resolutions allow for realistic depictions of different forest treatments and their effects on vegetation parameters (e.g., Moreno et al., 2016;Sun et al., 2018;Tague & Band, 2001;Tsamir et al., 2019;Waichler et al., 2005).When applied in a distributed fashion in a watershed, models allow for the investigation of the spatial pattern of the hydrologic response within forest treatment areas as well as downstream impacts in areas not experiencing disturbances, including streamflow variations.Due to their ability to conduct scenarios that address both independent and combined effects, numerical models can explore the outcomes from interacting hydrologic processes and feedbacks.
In this study, we used a process-based, distributed hydrologic model to quantify the impacts of a large-scale forest treatment plan known as the Four Forest Restoration Initiative (4FRI, Hampton et al., 2011;Robles et al., 2014, Figure 1a shows the full extent of 4FRI in Arizona) in the Beaver Creek watershed under different warming levels that represent plausible future conditions in the 21st century, ranging from +1 to +6°C (Taylor et al., 2012).Specifically, we applied the Triangulated Irregular Network (TIN)-based Real-time Integrated Basin Simulator (tRIBS, Ivanov et al., 2004;Rinehart et al., 2008;Vivoni et al., 2011) in the Beaver Creek watershed whose upper elevations are populated with dense ponderosa pine stands.The watershed drains to the Verde River, a key source of renewable surface water for the Phoenix metropolitan area.We built upon previous efforts using tRIBS in the Beaver Creek (Hawkins et al., 2015) and other ponderosa pine areas (Mahmood & Vivoni, 2008, 2011;Moreno et al., 2016;Rinehart et al., 2008).Here, we performed modifications to improve the representation of cold season processes, to enhance grid-based parameters and meteorological forcings, and to extend simulations over a 16year period to explore the interannual variability in hydrologic responses.A set of synthetic modeling scenarios were designed to identify the independent and combined hydrologic effects of warming conditions and forest treatment under a set of simplifying assumptions.These included the full implementation and maintenance of the 4FRI plan within the Beaver Creek for the different levels of air temperature increases imposed on a historical record without considering precipitation or relative humidity changes.Our goal was to identify if the 4FRI treatment project has the potential to lessen the impacts of warming on the hydrologic processes and the overall streamflow production from the Beaver Creek watershed.

Study Watershed and Data Sets
The Beaver Creek is a large sub-watershed of the Verde River with an area of 1,100 km 2 (Figure 1).The watershed has variable terrain and ecosystem conditions representative of the central Arizona highlands with a 1,625-m elevation gradient due to the Mogollon Rim, a rugged escarpment that forms the southern border of the Colorado Plateau (Hawkins et al., 2015).Two internal stream gauges were used to delineate the internal areas of the Wet Beaver Creek (WBC, 286 km 2 ) and Dry Beaver Creek (DBC, 367 km 2 ).Ecosystems vary with elevation from desert shrublands in the lowlands, through pinyon-juniper woodlands, and up to ponderosa pines at the highest elevations (Figure 2a).To complement this, Figure 2b shows the distribution of vegetation fraction, which is low in most of the watershed, except in forested regions and in developed areas.Soils are composed primarily of clay, clay loam, and loam, developed on basalts and cinders of volcanic origin (Soil Survey Staff, 2020), with some bedrock-lined channels that have little subsurface hydrologic connection (Figure 2c, Baker, 1982).As shown in Figure 2d, exposed bedrock leads to shallow soils, while other areas have deeper soil profiles.To summarize this data, Table 1 describes the areal coverage of the land cover and soil types in the study area.
Climate in the Beaver Creek watershed is classified as semiarid in a transition between hot (BSh) and cold (BSk) types according to Köppen-Geiger (Peel et al., 2007).During the study period defined as the Water Years (WYs, starting 1 October and ending 30 September) from 2003 to 2018, the areal averaged, mean annual precipitation was 524 mm, with a standard deviation of 145 mm.These estimates were obtained using the Next Generation Weather Radar (NEXRAD) Stage IV product (Du, 2011) at 4-km resolution, that was bias-corrected using five precipitation gauges operated by the Salt River Project (SRP) in and around the Beaver Creek (see Figure 1, Cederstrom, 2021).The annual cycle of precipitation (P) is bimodal (Mascaro, 2017;Vivoni et al., 2008), with a mean cold season (November to April) precipitation of ∼290 mm, primarily as snowfall, and a mean warm season (May to October) rainfall of ∼230 mm.The Beaver Creek is characterized by an areal averaged mean annual air  temperature (T air ) of 9.8°C during the study period, with the cold and warm seasons exhibiting 2.8°C and 24.9°C, respectively.These estimates were obtained from a spatial interpolation method that used bias-corrected 12-km resolution North American Land Data Assimilation (NLDAS-2, Cosgrove et al., 2003), three Remote Automatic Weather Stations (RAWS), and one of the Snow Telemetry (SNOTEL) sites at Happy Jack (Figure 1).In addition, observations of multiple variables at the SNOTEL sites at Happy Jack and Bar-M as well as daily estimates from a gridded (1-km) product referred to as Snow Water Artificial Neural Network (SWANN, Broxton, Leeuwen, & Biederman, 2019) were used to quantify climate and snow dynamics during the study period (Cederstrom, 2021).

Distributed Hydrologic Model Description
tRIBS is a continuous, spatially-explicit model of cold and warm season hydrologic processes occurring on complex terrain (Ivanov et al., 2004;Rinehart et al., 2008).To make full use of available geospatial data sets, tRIBS ingests elevation, soil, land cover, and meteorological conditions and resamples each to a domain consisting of a Triangulated Irregular Network (TIN, Vivoni et al., 2004).In tRIBS, Voronoi polygons are associated with each TIN node and serve as the finite-volume domain for water and energy balance calculations.For each Voronoi polygon, the model tracks the hydrologic response, including: (a) canopy snow and rainfall interception; (b) evapotranspiration from bare soil and vegetation; (c) snow accumulation, sublimation, and melt; (d) infiltration and soil moisture redistribution; (e) shallow subsurface flow; and (f) overland and channel flow, as detailed in Table 2.In prior studies, tRIBS has shown a good model performance with respect to hydrologic data in semiarid forested watersheds (e.g., Hawkins et al., 2015;Ko et al., 2019;Mahmood & Vivoni, 2008, 2011;Xiang et al., 2014).In addition, a range of different climate change assessments (e.g., Hawkins et al., 2015;Liuzzo et al., 2010;Piras et al., 2014;Schreiner-McGraw et al., 2020) and vegetation change evaluations (e.g., Moreno et al., 2016;Schreiner-McGraw et al., 2020) have been conducted with the model.
Due to the importance of the cold season in the Beaver Creek, we focused model calibration and validation efforts on the physical processes described by Rinehart et al. (2008) and Moreno et al. (2016), specifically the: (a) precipitation phase thresholds; (b) canopy wind profile for snow sublimation; and (c) snow albedo parameterization.Snow water equivalent (SWE) from the single layer energy and mass balance scheme on the ground surface was an important model performance target, as was the streamflow (Q) generated from snow melt at the two internal gauging stations.To improve the streamflow simulations, we added a frozen soil component that reduced the surface saturated hydraulic conductivity (K s ) in the model using information on T air as (Koren et al., 2014): where T min and T max are minimum and maximum air temperature thresholds, α is the conductivity reduction coefficient, and K′ s is the frozen soil value for K s .In the analysis of the simulations, we inspected the spatiotemporal variability of the water balance components (ΔS/Δt = P ET Q, where ΔS/Δt is the change in total soil water storage computed for the saturated and unsaturated as a residual, and ET is the evapotranspiration), and the partitioning of ET into soil evaporation (E soil ) and the canopy contribution (ET C = T + I ) from plant transpiration (T ) and evaporation from rainfall interception (I ), as well as the total sublimation (Sub = Sub C + Sub G ) from canopy (Sub C ) and ground (Sub G ) components.ET C that should be interpreted as primarily the effect of plant transpiration since T >> I (Vivoni, 2012a).Total ET is thus E soil + ET C + Sub.

Model Domain, Parameterization, and Initialization
Model simulations in the Beaver Creek were informed by the efforts of Hawkins et al. (2015).The watershed domain was converted from a 30-m resolution United States Geological Survey (USGS) digital elevation model (Figure 1) into a TIN using the hydrographic procedure of Vivoni et al. (2004), resulting in 76,624 Voronoi polygons (∼120-m resolution) that contained a stream network that mimicked the available hydrography.Elevation values were used to compute precipitation and air temperature lapse rates when preparing the meteorological forcing.The topographic index (TI) of each Voronoi polygon was computed as (Beven & Kirkby, 1979): where A c is the upslope contributing area and tan β is the terrain slope in the direction of downstream flow.TI is a useful hydrologic index to classify watershed locations, with low values often in areas near the watershed boundary and higher values found along the stream network.Here, the static TI is used to determine the locations within the watershed where changes occurred in hydrologic responses due to the warming and forest treatment scenarios.
Model parameterization in terms of soil and vegetation conditions followed previous tRIBS applications where initial values were obtained from literature sources (e.g., Hawkins et al., 2015;Ivanov et al., 2004;Mahmood & Vivoni, 2011;Moreno et al., 2016) and assumed to be spatially uniform within each class (Figures 2a and 2c).
Table 3 presents the model parameter values for the major soil and land cover types.Soil hydraulic and thermal parameters were related to the surface soil texture, obtained from the Soil Survey Geographic Database (SSURGO, Soil Survey Staff, 2020), and assumed invariant in time.Similarly, no temporal variability was included in the vegetation parameters due to our focus on evergreen forests, though the model can represent seasonal changes for deciduous ecosystems (e.g., Vivoni, 2012a;Xiang et al., 2014).We related a subset of the vegetation parameters (v f and LAI) to ponderosa pine basal area (BA) estimates for pre-and post-treatment conditions, following Moreno et al. (2016).Most vegetation parameters were aggregated to the polygons of the different land cover types, but we opted to map a set of vegetation parameters (v f and h) at 30-m resolution to the coarser model domain (i.e., ∼120-m resolution) to better capture the variations of the forest treatment in the ponderosa pine areas.
Model calibration and validation were conducted using a sequential approach at different spatial scales based on: (a) hourly SWE observations at two SNOTEL stations, (b) daily SWE estimates from SWANN, and (c) hourly observations at the two internal stream gauges.The Beaver Creek outlet gauge (Figure 1) did not have an updated rating curve after 2010 and was excluded due to a lack of high-quality data.For the SNOTEL site at  (2014).In these comparisons (Cederstrom, 2021), the simulated monthly average surface soil moisture and temperature (averaged over the top 10 cm) was compared to observed values from the near surface sensor at 5 cm depth.Model calibration at the Happy Jack SNOTEL station showed that the most sensitive cold season processes were the precipitation partitioning into liquid and solid fractions, the aging of snow albedo, and liquid water retention in the snowpack.This calibration effort resulted in cold season parameter values that varied slightly from Rinehart et al. (2008), as shown in Table 4.
After building confidence in the cold season processes, we performed model calibration and validation with respect to the internal stream gauges by focusing on the main soil and land cover types (Table 3).Manual calibration involved varying the most sensitive parameters within acceptable ranges, specifically: saturated hydraulic conductivity (K s ), hydraulic conductivity decay parameter with depth ( f ), air entry bubbling pressure (Ψ b ), pore size distribution index (m), leaf area index (LAI), optical transmission coefficient (k t ), stomatal resistance (r s ), and the stress thresholds for soil evaporation (θ * e ) and plant transpiration (θ * t ).Initial parameter values from Hawkins et al. (2015) were updated to account for: (a) a revised spatial distribution of soil depth from Soil Survey Staff (2020) that resulted in shallower soils in the Beaver Creek, (b) higher spatial resolution and bias-corrected meteorological forcing, (c) a revised spatial map of land cover types, and (d) the inclusion of a frozen soil effect deemed to be present in the forested clay soils of the Beaver Creek (NRCS, 2010).The manual calibration was not designed to reproduce low baseflows in the WBC, nor channel losses in the DBC.
Note.K s is the surface saturated hydraulic conductivity, θ s and θ s are the soil moisture contents at saturation and residual values, m is the pore size distribution index, Ψ b is the air entry bubbling pressure, f is the hydraulic conductivity decay parameter with depth, A s and A u are the saturated and unsaturated anisotropy ratios, n is the soil porosity, k s and C s are the soil heat conductivity and heat capacity, p is the free throughfall coefficient, S c is the canopy storage capacity, K and g are the drainage coefficient and exponential parameters, a is albedo, k t is the optical transmission coefficient, r s is the canopy minimum stomatal resistance, LAI is the leaf area index, and θ * e and θ * t are stress thresholds for soil evaporation and plant transpiration, respectively.Two vegetation parameters (v f , vegetation fraction and h, vegetation height) were obtained from Picotte et al. (2019).Model initialization consisted of specifying a spatially-distributed depth to the water table which sets the initial soil moisture profile (Ivanov et al., 2004).In the absence of field data, the initial groundwater depth was obtained through a spin-up simulation over the period of 1 October 1997, to 30 September 2002 (4 years).This simulation used the observed meteorological forcing to allow the watershed hydrologic processes to occur under the influence of the specified terrain, soil, and vegetation properties.Since the soil depths were shallow, the groundwater depth stabilized after 4 years, as demonstrated by longer trials of up to 8 years in duration.We then used the final model state of the spin-up period as the initial condition for the WYs 2003 to 2018 simulation (16 years) and the warming and forest treatment scenarios.This was performed using a model restarting capability from previously outputted model states (Vivoni et al., 2011).

Model Meteorological Forcing
Meteorological forcings consisted of gridded, hourly variables at 2-m height above the canopy for precipitation, incoming solar radiation, air temperature, relative humidity, atmospheric pressure, and wind speed.The process to generate forcings involved bias-correction and spatial interpolation, where applicable, using ground stations and gridded products.For P, the NEXRAD Stage IV product had a spatial resolution of 4-km deemed sufficient for the model simulations (about 90 pixels in the Beaver Creek).As such, only a bias correction was applied using two multiplicative factors: (a) a monthly lapse rate, γ m (z), derived using the observations and the elevations (z) of five precipitation gauges and NEXRAD pixels (Mascaro et al., 2023), and (b) a monthly mean-field bias correction factor (ω m , Robles-Morua et al., 2012).Both γ m (z) and ω m had larger values in winter months as compared to the summer season, indicating larger elevation controls and more substantial NEXRAD errors during snowfall.To illustrate the outcomes, Figure 3 shows the areal averaged daily P in the Beaver Creek over the WYs 2003 to 2018, ranging from 0 to 68 mm/d, as well as mean seasonal maps of winter and summer P at 4-km resolution during the same period of WYs 2003 to 2018.Note the higher winter P, ranging from 100 to 500 mm per season, and the elevation variations for both seasons.
Other meteorological variables were spatially interpolated from three RAWS stations and the SNOTEL site at Happy Jack, and gap-filled with the 12-km, hourly resolution NLDAS-2 gridded products.The spatial interpolation methods of Liston and Elder (2006) were applied to the station data.Quality control at the ground stations resulted in gaps that were subsequently filled in with bias-corrected NLDAS-2 time series at the co-located pixels.
Bias-correction of NLDAS-2 was conducted using quantile mapping at the hourly resolution for the co-located pixels to the weather stations, applied separately for each month (Cannon et al., 2015).Due to the sparse station coverage and coarse resolution of NLDAS-2, local air temperature lapse rates in the Beaver Creek were calculated for each month using daily, 4-km maps from the Parameter-Elevation Regression on Independent Slope Model (PRISM, Daly et al., 1994;PRISM Climate Group, 2020).Estimated lapse rates varied from 6.78°C/km in December to 8.47°C/km in May (Cederstrom, 2021).Other meteorological variables were extracted from stations and gap-filled with bias-corrected NLDAS-2 products, where necessary, except atmospheric pressure that was unavailable and was obtained directly from NLDAS-2 and adjusted for elevation.Subsequently, the interpolation method of Liston and Elder (2006) was applied to derive 1-km resolution meteorological variables using: (a) the Barnes objective analysis scheme (Barnes, 1994), (b) adjustments with elevation, slope, and aspect, and (c) incorporation of local air temperature lapse rates.Figure 3 also shows the spatially averaged daily T air , ranging from 14 to 27°C, and the mean winter and summer seasonal averaged T air from WYs 2003 to 2018.Elevation, slope, and aspect impact T air , with lower values at higher elevations, for instance.

Numerical Experiments and Scenarios
We used the computational savings afforded by the domain representation through a TIN (Vivoni et al., 2004) and the parallel model capabilities (Vivoni et al., 2011) to conduct 16-year long simulations (10 runs in total) in the Beaver Creek watershed.The base case (BC) simulation consisted of the calibrated model parameters with the original meteorological forcing and forest characteristics, whereas the nine change scenarios evaluated the independent or combined effects of warming and forest treatment over the same period.Each simulation was conducted on 52 processors in a high-performance computing cluster at Arizona State University (ASU) using a sub-basin partitioning method (Vivoni et al., 2011), with an approximate total run time of 30 hr for each 16-year period.The initial model state before each simulation was identical and based on the restart file saved after the 4year spin-up period.A high interannual variability was noted in the meteorological conditions between WYs 2003 Water Resources Research 10.1029/2023WR035627 CEDERSTROM ET AL. and 2018, with several below-average years having P < 450 mm/yr (2004,2006,2007,2009,2014,2018) and a few years with P > 600 mm/yr (2003,2005,2010), resulting in the sampling of a wide range of hydrologic conditions.
Table 5 describes the warming and post-treatment (PT) scenarios and introduces the acronyms subsequently used.Warming refers to increases in air temperature (T air ) applied uniformly in space and in time at four different levels (i.e., +1, +2, +4, and +6°C) to the meteorological forcing derived from WYs 2003 to 2018.These values were selected to represent the range of robust temperature increases simulated in the 21st century within the Coupled Model Intercomparison Project Phase 5 (CMIP5).For the central Arizona highlands, Whitney, Vivoni, Bohn, et al. (2023) presented CMIP5 scenarios with T air increases from 1 to 6°C from 2025 to 2100.This wide range is due to variations among the climate models and emissions scenarios over several decades.We tested the  In practice, T air was increased in each hourly, 1-km grid derived from the spatial interpolation method of Liston and Elder (2006) resulting in uniform upward shifts of the diurnal and seasonal cycles.Meteorological forcings dependent on T air , for instance the incoming longwave radiation, were changed accordingly.Relative humidity, however, was not altered in order to focus on the warming effects alone.This is consistent with observational and modeling studies that indicate a near constant relative humidity under higher global temperatures (e.g., Held & Soden, 2006;Willett et al., 2007).For the range of T air and RH experienced in the region, increases of 1°C in T air lead to percent changes of 6%-9% in absolute humidity, while at +6°C, the absolute humidity varies from +40% to 62%.Other meteorological variables were not modified.While simplified, this approach has been used previously to assess warming impacts on hydrologic conditions (e.g., Nemec & Schaake, 1982;Xu & Halldin, 1997; J. Y. Zhang et al., 2013) and can isolate the effects of warming on cold and warm season processes independent of precipitation changes which can introduce further variations.
In Table 5, PT refers to the cases where ponderosa pine forests have been thinned according to the specifications of 4FRI (Hampton et al., 2011;Robles et al., 2014).Note that the treatment area covers the entirety of the ponderosa pine forests in the Beaver Creek (Figure 2a) and was implemented instantaneously and maintained at that level over time (i.e., no vegetation regrowth is allowed).The pre-treatment conditions were based on estimated basal areas from 2011 at 250-m resolution over the entire forested area (Wilson et al., 2013).However, the 4FRI treatment case was only available for the DBC, as derived from the Forest Vegetation Simulator (Crookston & Dixon, 2005) applied by the US Forest Service using anticipated forest thinning practices and constraints.As a result, we needed to expand the posttreatment case BA values from the large polygon areas of variable sizes (e.g., 1-10 3 ha) in the DBC to the remaining forested areas.We resampled the posttreatment BA to a 250-m grid that matched the data of Wilson et al. (2013).
Then, we fit a linear regression to the pre-and post-treatment BA values from co-located pixels, resulting in a Pearson correlation coefficient (CC) of 0.50.As shown in Figure 4a, the pre-treatment BA distribution (black line) had larger values (mean of 17.4 m 2 /ha, Std of 10.4 m 2 /ha), while the posttreatment BA distribution (blue line) had significant reductions in BA (12.5 ± 4.5 m 2 /ha).The post-treatment case based on the regression (labeled "regression" in Figure 4a) preserved the mean of the post-treatment BA distribution from 4FRI (labeled "original" in Figure 4a) but reduced its variance.This implies that the post-treatment scenario in the areas outside of the DBC does not capture the range of BA values that would result from the application of the US Forest Service method.The spatial maps shown in   2011), with an R 2 = 0.99.This approach allows for the effects of forest treatment to be captured across all relevant vegetation-related hydrologic processes.

Testing of Snowpack and Streamflow Simulations
Figure 5 illustrates the hydrologic model performance with respect to the basin-averaged SWE time series from SWANN and the distributed SWE estimates from the same product.Prior to performing this test, the cold season processes in tRIBS were calibrated at the Happy Jack SNOTEL and validated at the Bar-M SNOTEL sites against daily SWE observations.Both SNOTEL sites are in a ponderosa pine forest clearing on flat terrain in the high elevations of the Beaver Creek.Model performance was deemed good due to a high CC, a modest to low root mean squared error (RMSE), and a Bias near 0.70, as shown in Table 6.Underestimation of SWE by the model primarily occurred during years with lower snow accumulation (i.e., <150 mm of peak annual SWE).Parameter transfers from the Happy Jack SNOTEL station to the Bar-M site showed an improvement in performance (Table 6), thus providing confidence in the calibrated snow parameters (Table 4).As at the SNOTEL stations, the basin-averaged SWE was simulated well as compared to SWANN over the calibration and validation periods which exhibited different precipitation amounts, indicating robustness to meteorological forcing (Figure 5a).Overall snow model performance was as good as in prior studies (e.g., Mahmood & Vivoni, 2011;Moreno et al., 2016).While the simulated peak annual SWE was generally lower than SWANN estimates, the model performed well in capturing the annual timing of snow accumulation and melt, leading to a CC greater than 0.90 (Table 6).Large interannual changes in SWE, as driven by variations in P and T air , were well represented.Underestimations in simulated SWE were attributed to forests from 2,000 to 2,400 m (Figures 5b and 5c), where precipitation partitioning and snowpack dynamics are sensitive to small changes in T air .
Capturing snow processes adequately was important prior to evaluating streamflow in the Beaver Creek where cold season (November-April) and warm season (May to October) runoff mechanisms differ significantly (Hawkins et al., 2015;Mascaro et al., 2023;Moreno et al., 2016;Robles et al., 2017).Figure 6 compares modeled and observed hourly streamflow at both gauges, while the insets compare monthly average streamflow.Streamflow performance metrics are shown in Table 6 for the Wet and DBC gauges, with mean flows of 0.70 and 0.85 m 3 /s, respectively.Larger peak flows and longer periods of low flows occur in the DBC, whereas the WBC has sustained baseflow during the year.Overall, the model captures the magnitude and timing of the streamflow events at both sites.While RMSE values are larger than mean flows, these are highly affected by a few large peaks (Clark et al., 2021).Therefore, the streamflow metrics indicate the model can capture well the complex runoff generation in the region (Mascaro et al., 2023).Under-estimations (Bias < 1) were common in the calibration period containing larger events, while over-estimations (Bias > 1) were more typical in the validation period at all time scales.While RMSE values were low overall, the DBC exhibited larger RMSEs due to its higher peak flows.
The model simulates monthly streamflow peaks to occur slightly earlier than observed, especially in the WBC, where CC values were often lower.Since snow melt timing was captured well (Figure 5), we attribute this shift in timing to deep subsurface processes, likely flow through fractured bedrock not captured in the model, which delay streamflow to the late spring.Similarly, the model does not reproduce low flows in the WBC arising from springs (Schenk et al., 2018).

Interannual Variability of Water Balance
Given the agreement between modeled and observed hydrologic processes, we analyzed the interannual variability in the meteorological forcing and water balance components.Table 7 shows the basin-averaged annual hydrologic fluxes (P, ET, Q), change in total soil water storage (ΔS/Δt), and water balance ratios (ET/P, Q/P) and their interannual metrics.As quantified by the coefficient of variation (CV), the simulated outlet streamflow has an amplified interannual variability relative to P (CV = 0.63 and 0.28, respectively), attributed to the watershed response.Among these are the large yearto-year changes in snowpack conditions (Figure 5), conditioned by P and T air , which lead to variability in sublimation (Sub, CV = 0.42).In contrast, components of ET C + E soil (CV = 0.10) have a moderating effect due to biophysical limits to evaporative losses (e.g., ET is often capped during wet years, Perez-Ruiz et al., 2022), leading to a wide variation of ET/P from 0.54 to 1.41 with a mean of 0.88.This range of ET/P is also explained by the carryover of subsurface water from wet to dry years, as quantified by ΔS/Δt computed as a residual.Positive ΔS/Δt indicated a net increase of subsurface water for that WY, while negative values resulted from sustained evaporative losses during subsequent drier periods.In addition, WYs with a high P typically yield a high Q, with runoff ratio (Q/P) values varying from 0.03 to 0.19 with an interannual mean of 0.11, consistent with regional estimates (Grimm et al., 1997).

Impacts of Warming and Forest Treatment
The independent and combined modeling scenarios were first analyzed for their impact on the snowpack conditions.Figure 7 presents the mean annual snowfall fractions (a-c), snow covered area (d), and snow melt (e) for a subset of the scenarios.As warming is applied, the mean snowfall fractions showed large decreases for higher elevations in the Beaver Creek (e.g., above 2,000 m).For BC6 (+6°C), most of the cold season precipitation is in the form of rainfall.A trend toward less snowfall has been identified in both the historical record and in climate change projections (e.g., Klos et al., 2014;Knowles et al., 2006;Whitney, Vivoni, Wang, et al., 2023).The snowto-rain transition is controlled by the precipitation partitioning algorithm and its calibrated parameters (Table 4) and is thus independent of forest treatment.With warmer temperatures and less snowfall, the snow covered area was reduced in magnitude and duration.For BC0, monthly mean peak values occurred in January at 28% of the area, whereas warming led to a lower snow covered areas (18% at BC2, 6% at BC6) and an earlier snow free   period in the spring.Peak snow melt also occurred earlier, shifting from March to December, and snow melt ended earlier in the spring.Lower snow melt rates under warming were due to less available snowpacks.Overall, the forest treatment showed a limited impact on snowpack conditions, except for increasing snow melt rates at +0 and +2°C by 1-3 mm/month.This was attributed to reductions in canopy and ground sublimation which led to higher ground snowpacks, as detailed in the following.
Figure 8 shows mean monthly variations of soil moisture and ET components for a subset of the scenarios to illustrate the impact of changes in cold season processes.Surface (θ surf , top 10 cm) and root zone (θ root , top 1 m) soil moisture both increased due to warming, with a more prolonged effect into the warm season in the deeper soil.
In addition, a 1 month shift was noted in the timing of the peak soil moisture due to warming (e.g., from January to December for θ surf and from March to February for θ root for +6°C).Forest treatments further increased soil moisture, with a larger change for θ root that persisted year-round.As a result, both effects altered soil moisture conditions which feedback onto the evaporative losses.E soil and ET C increased with warming throughout the year, but larger changes were noted in the cold season in response to available θ surf and θ root from increased snow melt.
Note that E soil is greater than ET C due to the low vegetation fraction in the watershed (Figure 2b).Post-treatment cases showed a negligible change E soil for all values of T air .In contrast, ET C was reduced by forest thinning, in particular at higher warming, with decreases that tended to be larger during the warm season.This suggests that forest treatments reduced plant water uptake from the root zone more effectively when the atmospheric demand for water is higher, thus retaining soil moisture at greater soil depths.
Figure 9 summarizes the independent effects of warming and its combination with forest treatment on the mean annual water balance components averaged over the Beaver Creek, relative to the BC0.Positive percent changes (%), defined as 100 × [(Scenario BC0)/BC0], indicate higher amounts of a water balance component.This is complemented with Table 8 that presents absolute mean annual values for each water balance component in each scenario.With warming alone, an increase in total evapotranspiration (ET BC ) occurred, reaching up to 11% for BC6, within which ET C exhibited a larger percent change than E soil (whose absolute values are larger).This was due to the capacity of ponderosa pines to access deep soil moisture (θ root ) that carried over in greater amounts to the warm season to sustain transpiration, whereas this effect was more limited for shallow soils (θ surf ).In response, the change in storage (ΔS BC ) and streamflow (Q BC ) experienced large ( 46% or 12 mm/yr) and modest ( 8% or 5 mm/yr) changes, respectively, with higher warming alone.Note that Q/P across the scenarios was generally low, ranging from 0.08 to 0.13.Large percent changes in ΔS BC are attributed to the low mean annual value in BC0 (Table 7), whereas its larger absolute values as compared to Q BC indicate that soil water storage is impacted significantly.Q BC first exhibited an increase with warming of +1.5% or 1 mm/yr for BC1, followed by a decrease with warming, reaching 29% (or 18 mm/yr) for BC6 (Table 8).This behavior was due to counteracting effects of warming on the different ET components.Increases in ET C,BC and E soil,BC with T air occurred at the same time that snow sublimation (Sub BC ) and its components decreased significantly due to lower amounts of snowfall and less available snowpacks (note that absolute values of Sub are less than both E soil and ET C , Table 8).This is consistent with a significant reduction of the mean annual value of the maximum SWE (SWE max in mm) experienced each year as a basin average (Table 8).After a +1°C increase in T air , the positive warming effects on ET C,BC and E soil,BC outpaced the reductions in Sub BC , such that total ET was significantly increased and Q BC decreased overall.
Forest treatment led to shifts in water balance components that were nearly constant for the warming levels.The removal of ponderosa pines reduced ET relative to the BC (ET PT < ET BC ) due to a reduction in canopy contributions (ET C , PT < ET C , BC , 19%) and in sublimation (Sub PT < Sub BC , 4%), with a small overall impact on soil evaporation (+0.5%), which had an overall larger magnitude.Forest treatment effects on Sub were attributed to: (a) less canopy snow interception reduced sublimation from trees and increased the ground snowpack, and (b) more radiation reaching the ground promoted faster snowpack melting that reduced the time available for sublimation.SWE max exhibited small increases (∼1 mm) due to forest thinning, relative to no treatment, up to +4°C  (Table 8).In response to reductions in ET components, ΔS PT and Q PT experienced large to modest percent increases, respectively, with forest thinning at all warming levels (Table 8 presents absolute values).Across all warming levels, Q PT was 12% (or 7 mm/yr) higher than Q BC while ΔS PT was 42% (or 11 mm/yr) greater than ΔS BC , on average.An important warming threshold was identified at +4°C where the thinning effects leading to higher streamflow no longer counteracted the warming effects leading to lower streamflow.As a result, the scenarios suggest forest thinning increases streamflow relative to the BC only up to +4°C.Note, however, that storage decreases (ΔS PT = 32% for PT4) due to warming still occurred at this threshold after forest treatment.This implies that forest thinning led to increased subsurface water storage, relative to no treatment, but that the watershed still experienced an overall loss due to the higher evaporative demand (E soil,PT and ET C,PT ) at +4°C.
The sensitivity of streamflow to warming and forest treatment was further explored using double-mass plots (Figure 10a) of the outlet streamflow.These plots are useful for identifying hydrologic changes due to forest disturbances (e.g., Biederman et al., 2015;M. Zhang & Wei, 2012).Annual cumulative Q (mm) from each scenario were compared to the base case (BC0), with linear regressions shown and metrics provided in Table 9. Deviations from the line of perfect agreement (labeled 1:1 line) indicated if a

Water Resources Research
10.1029/2023WR035627 scenario had higher (BC1, PT0, PT1) or lower (BC4, BC6, PT6) Q as compared to BC0.Progressive changes were noted in annual Q as warming was applied, with slight increases from BC0 to BC2, followed by large decreases at BC4 and BC6.Also evident in the double-mass plots is increased variability with low levels of warming (e.g., BC1 and BC2) with the drier years falling below the 1:1 line and wetter years above the 1:1 line (Figure 10a and Table 9).The largest snow inputs tend to occur in wet years, therefore, reduced Sub from warming overcomes increased ET only in years with substantial snowfall.Forest treatment led to maximum values of annual Q for PT0, followed by progressive decreases up to PT6, as shown by the regression slope (m) indicating the percent changes in Q relative to BC0.When compared across scenarios, we found that at a warming level of +6°C, forest thinning (PT6) leads to a +16% increase in Q relative to no treatment (BC6).
Outlet streamflow changes were further diagnosed through monthly mean variations (Figure 10b).A comparison of monthly Q from a subset of the scenarios revealed the higher sensitivity to warming in seasonal streamflow as compared to forest treatment.Note that warming initially increased Q in DJF and decreased Q in the transition periods of ON and AM between BC0 and BC2.This is attributed to: (a) higher rainfall fraction and larger snow melt leading to higher Q in DJF, and (b) lower snow melt and higher ET in the transition periods both reducing Q.At a warming of +6°C, monthly Q was reduced in months with snow relative to BC0, except in January.Overall, warming decreased the duration of the winter streamflow period and concentrated it within January.In contrast, forest treatment led to a small increase in monthly Q at each level of warming, with the largest impacts during DJF for +6°C when a reduction of trees decreasing the plant transpiration which promoted higher soil moisture values in the surface and root zones, with a greater effect for deeper soils.Streamflow during the months of JJAS had a small sensitivity to warming and showed limited increases after forest treatment.Importantly, the independent and combined scenarios had higher impacts on streamflow during the cold season (November to April), as compared to the warm season, in response to larger effects on snow dynamics, soil moisture, and evapotranspiration components in that season (Figure 8).

Effects of Forest Treatment on Downstream Locations
When averaged over the Beaver Creek, forest treatment decreased ET, mostly through reductions in plant transpiration, which led to increases in Q and ΔS/Δt across all warming levels.To illustrate these spatial dynamics, Figure 11 shows the mean annual evapotranspiration (ET = E soil + ET C + Sub) as difference maps in the ponderosa pine areas for a subset of the modeling scenarios.The forest treatment effect alone is shown by PT0-BC0, where large regions of negative (red) values indicated that ET is reduced by upwards of 70 mm/yr.Nevertheless, small positive (blue) areas were also present in PT0-BC0 in or near low-lying areas and channels downstream of the forest treatment.Warming effects at +1°C are shown in the difference maps of BC0-BC1 and PT0-PT1, with the latter having forest treatments superimposed.In both cases, opposite effects occurred to ET at +1°C for elevations higher than or lower than 2200 m.Above this transition zone between snowfall and rainfall (Figure 7), modest decreases in ET ( 10 to 20 mm/yr) are due to the presence of less snow which reduced Sub.Below this elevation, small increases in ET (+10 to +20 mm/yr) are due to enhanced E soil and ET C at +1°C.At higher levels of warming (BC6-BC4 and PT6-PT4), increases in E soil and ET C overwhelmed the effect of reduced Sub and ponderosa pine area experienced large increases in ET (+10 to +70 mm/yr).Forest treatments moderated ET enhancement due to warming and reconfigured the spatial distribution.For instance, forest thinning at high elevations displaced ET toward lower-lying channel areas.Figure 12 illustrates this displacement by comparing the mean annual ET as function of the topographic index (TI).Post-treatment ET for PT0 was lower than BC0 at all values of TI.However, the difference between the two cases was amplified for low TI (values of 4-10) depicting high slope regions with low contributing area that are associated with forested regions.In contrast, ET differences were significantly reduced for high TI (values of 20-27) found in flatter areas with high contributing area, often near the channel network.This suggests that forest thinning in upslope areas with low TI led to large ET reductions, which in turn increased soil moisture in the root zone (Figure 8b).Since downstream areas with high TI had smaller ET reductions after thinning, we hypothesize that lateral connectivity through the subsurface is responsible for the displacement of soil water from low to high TI areas as noted in the ET patterns (Figure 11).In essence, the lateral hydrologic fluxes captured in the model allow for the displacement or translocation of forest thinning effects toward downstream areas.

Discussion and Conclusions
Ponderosa pine forests in headwater regions of the southwestern US are important for provisioning freshwater to downstream areas.In addition to their susceptibility to warming (Adams & Kolb, 2005;Allen et al., 2002), forested regions are also prone to destructive crown fires which modify large portions of the landscape (Dore et al., 2010;Huffman et al., 2001;Neary et al., 2012).In response, various types of forest treatments (e.g., thinning, grazing, prescribed burning) are underway in the southwestern US to reduce severe wildfires (e.g., Johnston et al., 2021;Sankey & Tatum, 2022).These ongoing or planned landscape-scale forest treatments require quantifying the induced watershed responses as part of initial treatment planning and in their posttreatment verification (e.g., O'Donnell et al., 2021;Simonit et al., 2015).This is important given that prior studies have found mixed outcomes regarding the hydrologic responses to forest thinning, as reviewed by Goeking andTarboton (2020) anddel Campo et al. (2022).In semiarid ponderosa pine forests, streamflow increases observed after different types of forest treatments are often cited to be due to reductions in evapotranspiration and interception and their effects on increased snowpacks and soil moisture (e.g., Baker, 1986).However, tracking the hydrologic mechanisms responsible for streamflow changes after forest thinning is difficult from field observations alone due to the: (a) large year to year variations in meteorological conditions that create differing hydrologic responses, (b) recovery of understory vegetation that often occurs over a several year period after treatments, and (c) limitations in the measurement of the hydrologic states and fluxes within a watershed, among other factors.
Here, we developed a set of spatially-explicit modeling scenarios that combined a prescribed forest treatment for ponderosa pine areas in the Beaver Creek watershed of central Arizona (Hampton et al., 2011) with warming experiments.We adapted treatment plans from the regional 4FRI effort, but only considered its effects within the Beaver Creek and its local terrain, soil, and vegetation conditions.Independent and combined scenarios were conducted after establishing the model performance with respect to available observations and estimates over a 16-year period to capture interannual variability.This included calibration and verification efforts at two SNOTEL stations and at two stream gauges and with respect to a spatiotemporal SWE product.In addition, modeling efforts improved the parameterization, meteorological forcing, process representation, and performance evaluation for the distributed hydrologic model in relation to prior studies in the region (Hawkins et al., 2015;Moreno et al., 2016).As noted by Mascaro et al. (2023), hydrologic simulations in forested regions of central Arizona are challenging due to the variable nature of the meteorological forcing during the cold and warm seasons.The long-term evaluation performed here over a set of highly variable meteorological conditions indicated that the modeled snowpack dynamics are captured well in terms of timing, but often have lower magnitudes which translate to underestimations in cold season streamflow and shifts in streamflow timing.These limitations, however, do not preclude an evaluation of the model sensitivity to forest treatments under different levels of warming.Indeed, we would anticipate a larger sensitivity to the hydrologic processes that are summarized next, especially during the cold season: 1. Increases in air temperature have a significant impact on snowpack conditions, including the snow covered area and snow melt rates, due to a reduction in snowfall relative to rainfall.At +2 and +6°C increases, mean annual SWE was reduced by 11% and 99%, respectively, with an accelerating decrease at higher warming levels.Furthermore, snowpack duration was reduced by 1-2 months and peak snow timing shifted earlier toward December by +6°C.In response, snow sublimation rates from the ground and tree canopies were reduced significantly.Forest treatment had a small impact on snowpack conditions at +0 and +1°C by reducing sublimation rates, but the benefits were negligible for greater warming.This implies that the oftencited effects of forest thinning on snowpack dynamics are relatively small.2. Earlier snow melt and reduced sublimation induced by warming promoted increases in soil moisture, with a larger effect over the root zone.Higher soil moisture, in turn, increased soil evaporation and plant transpiration by +9% and +16%, respectively, at +2°C, counteracting sublimation reductions to cancel out changes in evapotranspiration.However, at higher levels of warming, evapotranspiration increased by +11% at +6°C, driven by plant transpiration.Spatial patterns of evapotranspiration changes showed an elevation gradient due to the differential warming effects on sublimation and other components.Forest treatment impacted the spatial distribution by reducing plant transpiration by 17% at +6°C relative to no treatment which led to more evapotranspiration concentrated downstream in channels.Thus, plant transpiration was responsible for most of the thinning impacts and exhibited both local and downstream effects.

Water Resources Research
10.1029/2023WR035627 3. Counteracting effects of warming on different evapotranspiration components also explained how streamflow was increased by +1.5% up to +1°C followed by strong reductions of 29% up to +6°C.Low levels of warming (1-2°C) also increased streamflow variability, as only years with substantial snowfall showed reductions in sublimation that overcame increases in evapotranspiration.Warming shortened the cold season streamflow period from 6 to 4 months and shifted the peak streamflow to occur earlier in January.Forest treatment resulted in an average increase in streamflow of 12% and in the change in soil water storage of 42% relative to no treatment when averaged across all warming levels, with increases occurring exclusively in the cold season for streamflow and throughout the year for soil water.A threshold was identified at +4°C at which the forest treatment effects leading to higher streamflow no longer counteracted the streamflow reductions induced by warming.
Here, we combined the capabilities of a process-based, distributed hydrologic model operating at high resolution with a set of scenarios to provide insights into the compounding hydrologic processes that are critical for understanding forest treatment effects under a warming climate (e.g., Stephens et al., 2013).To isolate the warming effect from other climate change factors, such as precipitation and relative humidity changes, we constructed scenarios that increased air temperature uniformly in space and time by fixed amounts to represent plausible ranges during the 21st century in the region.In addition to ignoring precipitation changes, this approach kept the same relative humidity across all warming levels, which resulted in small increases in absolute humidity (e.g., 6%-9% at +1°C, and 40%-62% at +6°C).We note that Simpson et al. (2024) identified discrepancies between modeled and observed trends in absolute humidity for the southwestern US (e.g., increasing climate model trends that are not present in data).However, the largest discrepancies are present during the warm season, while streamflow in the Beaver Creek is generated primarily in the cold season.In future work, improvements to the meteorological forcing scenarios can be addressed using statistically downscaled climate model outputs, as performed by Whitney, Vivoni, Bohn, et al. (2023) and Whitney, Vivoni, Wang, et al. (2023) for the Colorado River.For these scenarios, changes in relative humidity are small (±5%) for the low mean annual values (∼40%) in the region.Future work can also assess the model sensitivities to changes in relative humidity.
Other limitations are noteworthy with respect to the forest treatment scenarios.First, the model resolution (∼120m) is limited in its ability to capture fine-scale dynamics in snow and soil moisture that can impact responses to forest disturbance and warming (e.g., Moeser et al., 2020;Sexstone et al., 2018;Svoma, 2017).Second, we applied a simple transformation of BA changes to modify the vegetation parameters in ponderosa pine areas.
While the parameter transformations were based on previous studies in ponderosa pine regions, these could be improved through site-specific regressions using, for instance, the products of Picotte et al. (2019) and Hampton et al. (2011).Furthermore, post-treatment BA estimates did not capture the full range of observed variability outside of the DBC.Third, we assumed that thinning occurred instantaneously and was maintained throughout the study period and for the different levels of warming.As a result, the forest treatment scenarios did not account for vegetation regrowth (e.g., Minott & Kolb, 2020) or the interaction of warming with forest structure (e.g., Shriver et al., 2021).For future studies at specific treatment sites, the distributed model could be applied at higher resolutions (1-2 m) using elevation and canopy height data from Light Detection and Ranging (LiDAR) to capture fine-scale snow and soil moisture dynamics (e.g., Mahmood & Vivoni, 2011).Using LiDAR data sets, model vegetation parameters could be more faithfully tied to the properties of individual trees and hydrologic simulations could capture the wide range of natural variability within a forested landscape.Furthermore, modeling scenarios could be developed to: (a) estimate the watershed responses of phased treatments occurring at different times, (b) implement changes in model parameters to account for post-treatment regrowth trajectories, and (c) guide the selection of treatment areas and types that maximize hydrologic benefits for a targeted level of warming.
As noted in this model application, the fidelity of process-based, distributed hydrologic model is enhanced with high-resolution meteorological forcings and descriptors of terrain, soil, and vegetation properties, as well as with the ability to confront the model with a coordinated set of long-term observations of different cold and warm season processes (e.g., Fatichi et al., 2016;Schreiner-McGraw & Vivoni, 2018;Vivoni, 2012b).When confidence is built iteratively in the process-based simulations, modeling scenarios can then help explain the hydrologic mechanisms operating under confounding factors and provide novel insights that are directly applicable to forest treatment decisions.For instance, the modeling scenarios indicated that forest thinning primarily acted to reduce ponderosa pine transpiration (average of 19% relative to no treatment) which propagated to higher soil water storage year-round (+42%), especially at larger soil depths, and increased streamflow (+12%) in the cold season.

10.1029/2023WR035627
Other hydrologic mechanisms that are often associated with forest treatment, such as effects on snowpack sublimation ( 4%) or soil evaporation (+0.5%), played minor or negligible roles, respectively.
Modeling scenarios also revealed that local forest management for wildfire risk reduction might compensate for the warming effects on hydrologic processes up to a certain temperature threshold beyond which their effectiveness is negligible.As a result, the costs and benefits of treatment projects need to be considered under the current climate and for the warmer conditions imposed onto forested ecosystems for decades into the future.In addition, projects need to be cognizant of the varying time scales between the slow, long-term warming effect and the staggered pace of treatments, their maintenance, and the reduced hydrologic benefits that are often observed after plant regrowth in large forests.As demonstrated in this study, extensive forest treatments that are maintained over time can provide short-term resilience toward anticipated warming effects in ponderosa pine forests of central Arizona, with implications toward other forested regions in the southwestern US.

Figure 1 .
Figure 1.(a) Location of Beaver Creek watershed in the Verde River of central Arizona (AZ), shown with black polygon and red box.(b) Elevation distribution from a 30-m resolution digital elevation model (DEM, USGS, 2017), including the boundaries of Wet and Dry Beaver Creek (DBC), and locations of United States Geological Survey (USGS) stream gauges, Snow Telemetry (SNOTEL) stations at Bar-M and Happy Jack, precipitation gauges, and weather stations.The stream gauges are Beaver Creek at Camp Verde, AZ (USGS 09505500), Wet Beaver Creek near Rimrock, AZ (09505200), and DBC near Rimrock, AZ (USGS 09505350).

Figure 2 .
Figure 2. Beaver Creek watershed characteristics.(a and b) Land cover types and vegetation fraction at 30-m resolution (Picotte et al., 2019).(c and d) Soil type polygons and soil depth at 30-m resolution (Soil Survey Staff, 2020).
SWANN, and the internal gauges, we divided the study period into a calibration(WY 2003 to WY  2010)  and a validation (WY 2011 to WY 2018) set of equal length.The entire period was used as independent validation for SWE comparisons at the Bar-M site.Model setups at the SNOTEL stations used a single Voronoi polygon, followingHawkins et al. (2015) and Méndez-Barroso et al.

Figure 3 .
Figure 3. Meteorological forcing for Water Years 2003 to 2018.(a) Time series of daily values of basin-averaged precipitation (P, mm/d) and air temperature (T air , °C).Maps of mean seasonal total or average P and T air , respectively, for (band c) winter and (d and e) summer.Here, winter was defined as the cold season (November-April) and summer as the warm season (May-October).P is from bias-corrected NEXRAD data and T air is from interpolated stations and NLDAS-2 after applying lapse rates.
aforementioned four warming levels to capture the full range of variability, with lower values projected earlier in the 21st century and larger values near 2100.

Figure 4 .
Figure 4. Basal area (BA) (m 2 /ha) for pre-and post-treatment ponderosa pine regions in the Beaver Creek.(a) Relative frequency of BA for pre-and posttreatment cases (original and regression).(b and c) Spatial maps of pretreatment and post-treatment (regression) BA which show a mean reduction of 4.9 m 2 /ha.

Figure 5 .
Figure 5. Modeled and SWANN-estimated snow water equivalent (SWE) from the calibration (Water Years (WYs) 2003 to 2010) and validation (WYs 2011 to 2018) periods.(a) Daily comparison of basin-averaged SWE.(b and c) Mean winter SWE along elevation bands aggregated to 1 km resolution.

Figure 6 .
Figure 6.Modeled and observed hourly and monthly mean Q for (a and b) Wet Beaver Creek and (c and d) Dry Beaver Creek during the entire study period.Shaded areas represent the interquartile range of monthly streamflow.Monthly precipitation shown as an average over the entire Beaver Creek in (b) and (d) for reference.

Figure 10 .
Figure 10.(a) Double-mass plots of cumulative outlet streamflow for each warming and post-treatment scenario with respect to the base case (BC0) during the study period.(b) Monthly mean outlet streamflow for the base case (BC0), warming only (BC2, BC6), post-treatment only (PT0) and combined warming and post-treatment (PT2, PT6) scenarios organized by Water Year from October to September.Interannual variations in (b) are omitted for clarity.

Figure 11 .
Figure 11.Spatial distribution of the differences in mean annual evapotranspiration (ET = E soil + ET C + Sub) across a set of scenarios for the ponderosa pine areas of the Beaver Creek.Positive values (blue) indicate ET is larger in the first scenario, whereas negative values (red) indicate the opposite.PT0-BC0 shows effects of forest treatment with no warming.BC1-BC0 shows the effect of +1°C under no forest treatment.Other cases combine warming and forest treatment.

Figure 12 .
Figure 12.Topographic index (TI) control on mean annual evapotranspiration (ET = E soil + ET C + Sub) in the base case (BC0) and post-treatment (PT0) simulations under no warming.

Table 1
Areal Coverage (%) of Soil and Land Cover Types in the Beaver Creek CEDERSTROM ET AL.

Table 2
Hydrologic Components of the tRIBS Model With Asterisk Indicating New Development

Table 3
Calibrated Model Parameters for the Major Soil and Land Cover Types

Table 4
Calibrated Snow Model Parameters Obtained at the Happy Jack SNOTEL Station min and T max are the minimum and maximum air temperature thresholds for frozen soil effects on K s , α is the conductivity reduction coefficient, T liq and T sol are the air temperature thresholds for precipitation partitioning into liquid and solid phases, ϕ liq is the liquid water holding capacity of the snowpack, α a and α m are snow albedo age values for accumulating and melting snow, and λ a and λ m are snow albedo age decay coefficients for accumulating and melting snow.

Table 5
Base Case Simulation and Scenarios With Independent and Combined Warming and Forest Treatment, Each Over a 16-Year Duration (Water Years 2003 to 2018)Precipitation forcings are from bias-corrected NEXRAD data over the same period for all cases.

Table 6
Performance Metrics Between Modeled and Estimated or ObservedVariables at Different Sites During Overlapping Periods Bias is the ratio of the mean of the simulations divided by the mean of the observations.RMSE is the root mean squared error.CC is the correlation coefficient.

Table 7
Simulated Annual Water Balance Components Averaged Over the Beaver Creek for the Base Case (BC0), Including Precipitation (P), Canopy Contribution and Soil Evaporation (ET C + E soil ), Total Sublimation From Ground and Canopy (Sub), Outlet Streamflow (Q), and Change in Total Soil Water Storage (ΔS/Δt) Water year P (mm/yr) ET C + E soil (mm/yr) Sub (mm/yr) Q (mm/yr) ΔS/Δt (mm/yr) ET/P ( ) Q/P ( ) Note.The annual total evapotranspiration (ET = ET C + E soil + Sub) ratio (ET/P) and the runoff ratio (Q/P) are also shown.Interannual mean, standard deviation (Std), and coefficient of variation (CV) values for all WYs (2003-2018) are shown at the bottom.

Table 9
Linear Regressions of Double-Mass Plots of Cumulative Outlet Streamflow Relative to the Base Case (BC0) for All ScenariosNote.ΔQ BC0 is the difference in annual streamflow relative to BC0.