Modeling the distributed effects of forest thinning on the long-term water balance and streamflow extremes for a semi-arid basin in the southwestern US

To achieve water resource sustainability in the water-limited southwestern US, it is critical to understand the potential effects of proposed forest thinning on the hydrology of semi-arid basins, where disturbances to headwater catchments can cause significant changes in the local water balance components and basinwise streamflows. In Arizona, the Four Forest Restoration Initiative (4FRI) is being developed with the goal of restoring 2.4 million acres of ponderosa pine along the Mogollon Rim. Using the physically based, spatially distributed triangulated irregular network (TIN)-based Real-time Integrated Basin Simulator (tRIBS) model, we examine the potential impacts of the 4FRI on the hydrology of Tonto Creek, a basin in the Verde–Tonto–Salt (VTS) system, which provides much of the water supply for the Phoenix metropolitan area. Long-term (20-year) simulations indicate that forest removal can trigger significant shifts in the spatiotemporal patterns of various hydrological components, causing increases in net radiation, surface temperature, wind speed, soil evaporation, groundwater recharge and runoff, at the expense of reductions in interception and shading, transpiration, vadose zone moisture and snow water equivalent, with south-facing slopes being more susceptible to enhanced atmospheric losses. The net effect will likely be increases in mean and maximum streamflow, particularly during El Niño events and the winter months, and chiefly for those scenarios in which soil hydraulic conductivity has been significantly reduced due to thinning operations. In this particular climate, forest thinning can lead to net loss of surface water storage by vegetation and snowpack, increasing the vulnerability of ecosystems and populations to larger and more frequent hydrologic extreme conditions on these semi-arid systems.


Introduction
Quantifying the hydrological effects of extensive, humandriven forest thinning is of primary importance for sustainable water resource management in semi-arid basins, where disturbances in the upland vegetation density and architecture can trigger zonal alterations to the components of the water balance (Biederman et al., 2014), resulting, sometimes, in streamflow shifts along an entire basin (MacDonald, 2000;Reid, 1993;Webb and Kathuria, 2012). Because precipitation is cycled through forests and soil, upland modifications in vegetation cover are expected to affect the dynamics of the entire basin in terms of water yield quantity and quality, and peak and low flows (Jones, 2000;Moore and Wondzell, 2005;Alila, 2004, 2013).
In north-central Arizona, the US Forest Service is leading a collaborative effort known as the Four Forest Restoration Initiative (4FRI), a large-scale restoration of ponderosa pine (Pinus ponderosa) along the Mogollon Rim, with the primary goal to mitigate fire risk through forest thinning (Hampton et al., 2011;Stephens et al., 2013). In addition to the Phoenix H. A. Moreno et al.: Hydrologic effects of forest thinning metropolitan area (PMA), as well as other towns and cities in the region, a number of ecological communities depend upon the freshwater derived from basins, whose headwaters extend along the restoration areas (Arizona Department of Water Resources, 2010;Baker, 1986). Besides changes in mean water yields, projected forest removal could potentially affect base flows during dry periods (Dung et al., 2012;Lin et al., 2007), while increasing the risks of downstream flooding in the rapidly responsive, steep-slope mountain basins (Eisenbies et al., 2007;Jones, 2000;Jones and Grant, 1996;Jones and Post, 2004). It is, therefore, critical to understand the hydrologic effects of forest thinning, in conjunction with the cumulative effects of climate change and other stressors (e.g., population increase, urbanization) that can be expected to exacerbate the impacts of human interventions in current basin land cover (Barnett et al., 2005;Dale et al., 2001;National Research Council, 2008).
Traditionally, evidence of the connections between forest thinning and water yield responses has been based on paired watershed studies. Most of these studies have identified immediate increases in runoff and sediment production (Bosch and Hewlett, 1982;Brown et al., 2005;Hibbert, 1983;Hornbeck et al., 1993;Sahin and Hall, 1996). However, in basins where water yield depends mainly on snow accumulation and melt, researchers have reported high variability and uncertainty tied to site-specific topography, forest structure and microclimatic conditions (Cline et al., 1977;Lundquist et al., 2013;Schelker et al., 2013;Stottlemyer and Troendle, 2001;Troendle and Reuss, 1997;Venkatarama, 2014;Woods et al., 2006). Multiple authors have found a direct relationship between thinning, snow interception reduction and ablation increase (Link and Marks, 1999;Lundquist et al., 2013;Varhola et al., 2010;Venkatarama, 2014). In Arizona, most of the data and knowledge regarding hydrologic response to treatments in piñon-juniper and ponderosa pine forests have been obtained from the Beaver Creek research watershed, located within the Verde River basin (Baker, 1984(Baker, , 1986Brown et al., 1974). Results indicate that the thinning of poderosa pine leads to statistically significant short-term increases in runoff, particularly in steep north-facing slopes. In addition, the duration of snow on south-facing slopes is affected by thinning intensity, overstory removal and higher exposure to wind and solar radiation (Baker, 1986).
More recently, physically based, spatially distributed hydrological models have complemented the experimental approach to provide new insights into the processes undergoing change, both prior and post-forest removal (Bathurst et al., 2004;Legesse et al., 2003;Li et al., 2007). Such work indicates that, due to shifts in evapotranspiration and soil hydraulic properties and moisture, increases in water yield can be expected after forest thinning (Hundecha and Bardossy, 2004;Li et al., 2007;Serengil et al., 2007;Webb and Kathuria, 2012)

Goals, organization and scope of this paper
While much has been learned from the Beaver Creek experiments, greater understanding is still necessary to provide the long-term estimates of water yield needed by water managers and land and water decision makers for semi-arid basins in Arizona. In this regard, the application of highly realistic, physically based, spatially distributed models that appropriately simulate the detailed behavior of catchment dynamics at relevant spatial and temporal scales can provide valuable insights.
Here we examine the potential impacts of extensive forest thinning on the hydrology of Tonto Creek, selected as a prototypical semi-arid watershed suitable for the inference of long-term impacts on water yield and extreme conditions on neighboring basins. Additionally, we explore the mechanisms responsible for change due to forest removal from local to basin scale. Specifically, we examine the following three questions related to the sustainability of water resources in this region: 1. Is the 4FRI likely to produce significant alterations in streamflow and the components of the water balance at the basin scale?
2. If so, what are the expected magnitudes of annual and seasonal water changes?
3. What are the physical mechanisms likely to be responsible for observed hydrologic shifts at the element (smallest computational unit) scale and how do they alter the soil column water balance in hillslopes having contrasting aspects?
We address these questions using a calibrated, highresolution, catchment-scale hydrological model (see Sect. 3) as a tool to reproduce the spatiotemporal dynamics of the Tonto Creek basin, both prior and post-forest treatment, under long-term (20-year) historic climate forcing. Using 20 consecutive years provides an ample range of climate variability (including El Niño-Southern Oscillation (ENSO) phases), while the study of feasible forest thinning scenarios within a distributed model provides management and policy relevance to the research questions in this study. In particular, we analyze the shifts in the probability distribution functions of mean and extreme (low and peak) streamflow values, and the implications for water security and flood risk of downstream communities. Further, we investigate the inter-annual and seasonal mechanisms that explain effects of forest thinning on river flows, snow water equivalent, basin evaporation and transpiration and soil water storage in the vadose and saturated zones. Subsequently, a closer look to the spatially distributed hydrological fields evidence their relation to the areas where restoration occurred and the physical mechanisms responsible for such responses. Finally, a more detailed analysis of the changes triggered at the element scale is Figure 1. Elements of the water balance in a forested hillslope: A fraction of the gross precipitation or snow (P) is intercepted by vegetation (Int) and the remaining volume reaches the ground as net precipitation or snow (P net ). Intercepted water (Int) is either unloaded from leaves and branches (P unl ) or temporarily stored for evaporation (E int ) back to the atmosphere. If snow occurs, P net builds up snowpack. When thermodynamic conditions allow, retained water in the snow can sublimate (S snow ), or after melting, it can infiltrate (Inf melt ), runoff (R) or be transpired by plants (T), evaporated from soil (E soil ), serve as groundwater recharge (G R ) or remain as soil moisture in the vadose zone. Analogously, if only liquid precipitation occurs, P net redistributes according to the processes mentioned above, except for the snow-related mechanisms. Runoff (R) can be produced through infiltration excess, saturation excess, perched return flow and/or groundwater contribution (Exf). Subsurface flow can occur through lateral vadose zone flow (S R ) and/or groundwater flow (GW flow ). An aerodynamic component has been added to this plot to mark the importance of the surface roughness by trees on evaporation and sublimation water fluxes.
performed at sites having contrasting (north and south) hillslope aspect.

Effects of forest thinning on hydrology
Forest disturbance and management activities have been shown to influence nearly all components of the water budget from the plot to the entire basin scale (Ice and Stednick, 2004;Waring and Schlesinger, 1985). Figure 1 illustrates the components of the water balance in a typical forested hillslope in the semi-arid southwestern US (with snow presence during the winter months). Liquid and solid precipitation (P) are the principal controls on spatial distribution, timing and magnitude of runoff, evapotranspiration, snow accumulation, soil water fluxes and storage. Forest reduction will impact mostly surface water storage and flow, and sub-surface flow within the vadose zone. Removal of trees reduces leaf area and, thus, plant interception (Int) allowing for more net pre-cipitation (P net ) to reach the ground surface (National Research Council, 2008;Verry et al., 1983). During the winter, reductions in Int lead to increases in snowpack depth and cover (Woods et al., 2006). Increases in P net result in increases in soil moisture, plant water availability and rapid runoff production, particularly during intense rainfall events (Helvey and Patric, 1965). In contrast, reduced biomass consumes less water volume through plant transpiration (T) but enhances evaporation from the soil, melted water and/or sublimation from frozen surfaces (E soil , S snow ) due to reduced shading of clear-sky shortwave solar radiation and sheltering for turbulent moment transfer by wind gusts (Biederman et al., 2012;Gustafson et al., 2010;Harpold et al., 2012a, b;Musselman et al., 2008;Veatch et al., 2009). Thus, water yield increases are expected earlier in the year due to a premature snowmelt season caused by increased wind and shortwave radiation exposure in this semi-arid, high elevation forest (Helvey, 1980;Hornbeck and Smith, 1997;Jones and Post, 2004;Link and Marks, 1999;Mahmood and Vivoni, 2013;Megahan, 1983). It has also been shown that silvicultural manipulations in forests, via prescribed fires, have produced changes in the hydraulic properties of the underlying soil that can persist for several years, depending on the fire intensity and soil composition (Benavides-Solorio and MacDonald, 2005;DeBano, 2000;Moody et al., 2005;Neary et al., 1999;Robichaud, 2000;Shakesby and Doerr, 2006;Lear and Danielovich, 1988;Woods et al., 2007). Previous studies report reductions of between 10 and 40 % in soil hydraulic conductivity during post-fire conditions (Leighton-Boyce et al., 2007;Robichaud, 2000;Shakesby and Doerr, 2006). Additionally, effects of forest operations for mechanical thinning, such as logging and carrying of heavy material on roads, trails and hillslope contours, favor the occurrence of faster and larger volumes of overland flow due to soil compaction (Bowling and Lettenmaier, 2001;Cline et al., 2010;Cuo et al., 2006;Fatichi et al., 2014;Harr et al., 1975;Jones and Grant, 1996;Marche and Lettenmaier, 2001;Wemple and Jones, 2003). Field studies conducted during pre-and post-treatment conditions reveal reductions of up to 67 % in soil hydraulic conductivity for randomly distributed locations within an area mechanically restored with heavy equipment (Grace et al., 2006;Grace III et al., 2007). The duration of this disturbance to soil conditions has received very little attention in the literature; however, a few authors consider it to be highly variable (from months to years) and dependent on both climate conditions and whether recurrent operations are maintained (Cline et al., 2010;Robichaud, 2000). The overall effects of human-driven forest modifications include induced changes in the basin hydrology through direct forest effects and soil collateral effects, which then determine the total hydrological response during storm and inter-storm periods.

The Four Forest Restoration Initiative as an agent of hydrologic change for the Verde-Tonto-Salt system
The 4FRI, led by the US Forest Service, is targeting the restoration of up to 9712 km 2 of contiguous ponderosa pine of the Kaibab, Coconino, Apache-Sitgreaves, and Tonto National Forests across the Mogollon Rim in Arizona. The primary goal of 4FRI is to improve forest resilience and function by reducing forest cover, through the use of prescribed burns and mechanical thinning to historical conditions similar to that of the early 20th century (Hampton et al., 2011;Schoennagel et al., 2004). The projected treatment areas overlap with the headwaters of important water supply basins, including the Little Colorado, an important tributary of the Colorado River, and the Verde-Tonto-Salt (VTS) system, whose surface waters serve important cities and villages in north-central Arizona, including the PMA (see Fig. 2). Agency representatives and stakeholder groups recently agreed on future reductions in the current basal area conditions of the ponderosa pine from an average of 2755 to 1332 m 2 km −2 , by focusing in the removal of small-diameter trees (Hampton et al., 2011;Sisk et al., 2006) to reduce the threat of intense fire events to human communities, wildlife habitat and key ecosystem components (Allen et al.,  2002; Chambers and Germaine, 2003). Figure 3a and b illustrates the current pre-treatment and projected post-treatment change basal area of ponderosa pine for Tonto Creek. The post-treatment scenario was obtained from the Four Forest Restoration Initiative implementation plan (http://www. fs.usda.gov/4fri). The reader is referred to (Hampton et al., 2011) for more details on the density criteria and projections. Restoration of sensitive areas is discouraged, including those with steep slopes or sensitive soils, in proximity to streams, having wildlife regulations, and areas of recent tree harvesting. However, the vast majority of the ponderosa pine covered area, classified as community protection management areas (CMPA), aquatic and municipal watersheds, Mexican spotted owl (MSO) restricted and wildlife habitat, have been declared suitable for restoration.

Study region and watershed characteristics
The VTS system is located in the central Arizona highlands, characterized by rugged mountains with steep slopes separated by narrow valleys. The headwater catchments of the VTS system lie on the Mogollon Rim, a large escarpment that holds a wide diversity of vegetation types and ecosystems (Arizona Department of Water Resources, 2010). Because of the high elevations and associated higher amounts of rainfall and snowfall, the Mogollon Rim area contains the state's most important water-producing watersheds and the greatest concentration of perennial streams, which, in turn, support riparian habitat (Arizona Department of Water Resources, 2010). Precipitation is bimodal at a mean annual value of 481 mm yr −1 , with the largest amounts during the winter months due to frontal storm systems and a secondary rainy period during summer, coincident with the highest evapotranspiration rates, via monsoon-driven precipitation (Arizona Department of Water Resources, 2010). The mean annual temperature and runoff in the region have been estimated as 17.9 • C and 79.8 mm yr −1 (Arizona Department of Water Resources, 2010). The VTS system provides groundwater to small communities and individual farmers, mostly based on the Tonto and Verde rivers, and, along with the water allocation from the lower Colorado River through the CAP (Central Arizona Project) canal, groundwater and treated effluent, supplies water for the two million inhabitants of the PMA in the Salt River Valley Water Users' Association. We use the Tonto Creek basin as a case study to ex- plore the potential impacts of the 4FRI during 20-year long simulations by imposing historic climate forcing. Although Tonto has the smallest catchment area in the VTS system, the areal fraction covered by ponderosa pine is one of the largest, and so it provides a good indication of the processes triggered by forest removal across the whole VTS system. Table 1 summarizes the major characteristics of this basin. Slopes vary around a mean of 28 % with a standard deviation of 21 % induced by drastic changes in elevation over short distances. The contrasting relief and the steep slopes lead to rapid runoff responses and short concentration (or response) times. Figure 4 shows the spatial distribution of elevation, hydrography, vegetation, soils and depth to bedrock for the study basin. Overall, the area is characterized by a dominance of sandy loam soils, forest vegetation and deep impervious rock. The projected restoration area lies between the lines of 1800 to 2400 m elevation.

Observed hydrologic data and climate forcing
We compiled regional weather and rain gauge, snow, and streamflow station data at a daily timescale from the NOAA, National Climatic Data Center (http://www.ncdc.noaa.gov/ cdoweb/search), Natural Resources Conservation Service (http://www.wcc.nrcs.usda.gov/snow/) and USGS National Water Information System (http://waterdata.usgs.gov/nwis), respectively (see Fig. 2). This set of stations was selected because of the continuous data availability from 1 January 1990 to 30 December 2010, the prevalence of stations within the VTS basin and few information gaps (< 0.5 % gaps). For regional climate forcing, we used the NASA Land Data Assimilation Systems data set (NLDAS; Mitchell et al., 2004), which includes net radiation, atmospheric pressure, air temperature, wind speed, precipitation and vapor pressure (http://ldas.gsfc.nasa.gov/nldas/). NLDAS is released on a one-eighth-degree grid over central North America on an hourly basis, constituting a superb climate forcing for continuous, distributed modeling purposes. For precipitation, NLDAS constructs its forcing data set from CPC (Climate Prediction Center) PRISM (Parameter-Elevation Regressions on Independent Slopes Model)-adjusted one-eighth-degree daily gauge analyses, temporally disaggregated using Stage II radar fields (Mitchell et al., 2004). Since the quality of distributed hydrologic simulations highly depends on the accuracy of Quantitative Precipitation Estimates (Carpenter and Georgakakos, 2004;Collier, 2007;Moreno et al., 2013Moreno et al., , 2014, we first evaluated and bias corrected NLDAS rainfall forcing to minimize model error propagation from the precipitation input (see Appendix A1). Using NLDAS, it can be seen that Tonto Creek presents a bimodal precipitation distribution with above-average values during DJFM and JAS and a unimodal temperature pattern, whose peak occurs during JJAS (Fig. 5). Further, a map with the spatial distribution of mean annual precipitation and surface air temperature is presented in Fig. 6. Compared with Fig. 3, it can be determined that projected areas for forest thinning coincide with the higher annual basin precipitation (P > 500 mm yr −1 ) and lower mean temperatures (Temp < 16 • C, see Fig. 6a, b).

Distributed hydrologic model
The triangulated irregular network (TIN)-based Real-time Integrated Basin Simulator (tRIBS) (Ivanov et al., 2004a;Vivoni et al., 2007b) is a continuous, physically based simu- lator of watershed dynamics. The model uses spatially varying topography, soil and vegetation characteristics and timeevolving distributed climate forcing to represent the processes governing movements of surface and subsurface water in a basin. tRIBS uses a TIN scheme to reduce computational workload and accurately represent topography, water flow paths and river networks (Vivoni et al., 2004). This TIN geometry determines a network of sloped Voronoi polygons that communicate through their edges by mass continuity and flux equations. Underground dynamics are constrained by spatially varying depth to bedrock, which acts as an impermeable surface that determines the lower aquifer boundary. tRIBS can be run on a multi-processor computer by taking advantage of parallelization via domain decomposition (Vivoni et al., 2011). tRIBS computes short-and long-wave radiation fluxes using geographic location, time of the year, cloudiness, aspect, emissivity, slope and albedo at each computational element. Incoming solar radiation is reduced by vegetative shading according to the Beer-Lambert law (Brantley and Young, 2007;Marshall and Waring, 1986) (see Appendix B). Effects of distant landscape on the amount of incoming radiation are accounted through radiation scattering and sheltering functions that are controlled by landview factors and hillslope albedo (Rinehart et al., 2008). Surface latent (i.e. evaporation and transpiration), sensible and ground heat fluxes are computed using meteorological conditions and soil moisture (Ivanov et al., 2004b). Snow processes are accounted for through a single-layer snow module with a coupled energy and mass balance approach that accounts for direct and diffuse solar (shortwave) and longwave radiation, snow interception and unloading, sublimation of intercepted and on-the-ground snow, accumulation and ablation of snow and infiltration of meltwater (Mahmood and Vivoni, 2013;Rinehart et al., 2008). Vegetation intercepts snow falling in solid form, based on its leaf area index, and unloads snow in relation to air temperature. Remaining on-the-ground and canopy snow can be sublimated depending on absorbed shortwave and long-wave radiation and aerodynamic conditions (Liston and Elder, 2006;Pomeroy et al., 1998;Wigmosta, 1994). Meltwater can either infiltrate or run off and eventually is routed downslope to the channel as surface or subsurface runoff. Rainfall interception follows the canopy water balance scheme (Rutter et al., 1971(Rutter et al., , 1975 including throughflow, drainage, storage and evaporation, values that are determined by plant architecture properties and vegetation fraction. Evapotranspiration processes account for (1) evaporation from wet canopy (E int ), (2) evaporation from bare soil (E soil ), and T. Total evapotranspiration (ET) is estimated using the Penman-Monteith equation that depends on the surface energy balance and aerodynamic conditions for surface and plants. The below-canopy distribution of the vertical wind speed follows a decay-exponential function depending on the biometric features of the forest determined by projected leaf area index (LAI) and vegetation height (see Appendix B) (Sypka and Starzak, 2013;Yi, 2008). Evapotranspiration partitioning depends on the ability of E soil and T to extract soil water from the surface and root zones and is determined by constant model stress factors (Ivanov et al., 2004a;Mendez-Barroso et al., 2013). A kinematic approximation for unsaturated flow is used to compute infiltration and propagate soil moisture fronts in an anisotropic soil column according to an exponentially decaying hydraulic conductivity condition (Cabral et al., 1992;Garrote and Bras, 1995;Ivanov et al., 2004a). The coupled framework of the unsaturated and saturated processes results in a set of runoff mechanisms, namely, infiltration-excess runoff (Horton, 1933), saturation excess runoff (Dunne and Black, 1970), groundwater exfiltration (Hursh and Brater, 1941), and perched return flow (Weyman, 1970). Routing of surface flow is achieved via hydrologic overland flow and hydraulic channel routing that uses a kinematic wave approximation (Vivoni et al., 2007a).

Computational domain, model parameters and initialization
We obtained a 30 m digital elevation model (DEM) from the National Elevation Dataset (Gesch et al., 2002) for the central Arizona region. A grid sensitivity analysis was performed, leading to a convenient mesh simplification through selection of a coarser grid resolution that guaranteed: (1) preservation of the spatial distributions of elevation, slope, curvature and hillslope aspect, and (2) scheduling of a multipleyear parallelized model calibration procedure in a feasible period of time. A TIN geometry was then constructed following a modified VIP (very important point) method that minimized the number of computational nodes and the Kullback-Leibler divergence between topographic density functions. This resulted in an optimum horizontal point density of d = 0.86 and n t = 1970 (d = n t / n g , where n t is the number of TIN nodes and n g is the number of DEM cells) with an equivalent cell size of r e = 964 m. The final TIN represents the basin topography with high accuracy and preserves the finest level structures of stream network, river flood plains and watershed divide through a double buffer node strategy. tRIBS requires specification of the spatially varying parameters associated with individual soil and vegetation classes, and of those that describe the properties of the hillslope and channel network routing, and the underground aquifer (Ivanov et al., 2004b;Moreno et al., 2012). Soil and vegetation parameters are assigned to the different classes represented in Fig. 4. Soil texture maps were derived from the State Soil Geographic (STATSGO) database at 1 : 250 000 scale providing full regional coverage. Similarly, vegetation type and fraction maps were obtained from the USGS National Land Cover Dataset (Homer et al., 2004) at 30 m resolution for the year 2006. Distributed land cover properties were determined by vegetation parameters extracted from ancillary 2006 Landfire products (http://www.landfire.gov/) and mathematical expressions, from the literature, depending on the pre-treatment and post-treatment forest basal area maps (see Fig. 3). Associated parameters include vegetation fraction, LAI, vegetation throughfall and canopy storage (see Appendix B). We consider only two vegetation fraction cases (pre-treatment and post-treatment) ignoring any intermediate vegetation phenology, re-growth or recurrent thinning operations (see Sect. 4.5). A spatially distributed bedrock depth map, at 1500 m spatial resolution, was obtained from the Northern Arizona Regional Groundwater-Flow Model (Pool et al., 2011) and used to set a lower impermeable aquifer boundary. Finally, a geomorphic relation between channel width (w in m) and contributing area (A in km 2 ) was derived from 21 field measurements taken during a field cam-paign along the basin main channel, resulting in the expression w = 9.303A 0.243 with R 2 = 0.76. tRIBS also requires a spatially distributed initial condition, provided by the depth to groundwater surface, to set soil moisture profiles following a hydrostatic equilibrium assumption. A 1500 m spatial resolution hydraulic head map, issued for spring 1990, from the Northern Arizona Regional Groundwater-Flow Model (Pool et al., 2011) was adopted as the distributed initial condition. The depth to groundwater then had a mean value of 248 m with a standard deviation of 183 m. The model was spun-up for 1 year (January to December, 1990) when dynamic steady-state conditions were reached in streamflow, groundwater and vadose zone moisture profiles.

Calibration and evaluation strategy
Our results are supported by calibration and evaluation tests with continually available hydrological information on the ground. First, a one-at-a-time (OAT) sensitivity analysis facilitated determination of the relative importance of model parameters as evaluated by performance criteria (Gupta et al., 2009;Gupta and Kling, 2011), revealing that watershed responses are mainly controlled by the set of soil and vegetation parameters shown in Table 2. For the case of soil parameters, those properties are the saturated hydraulic conductivity (K s ) and its decay exponent with depth (f), the air entry bubbling pressure (ψ b ) and the pore size distribution index (λ 0 ). These parameters control the infiltration, percolation, throughflow and runoff production rates, water retention and vadose zone wet front evolution. Complementary, for the vegetation classes, three parameters were found to dominate the runoff production through controls on interception, soil moisture, evapotranspiration and snowmelt rates. Those parameters are albedo (a), vegetation height (H v ) and optical transmission coefficient (K t ). Parameters, other than those listed in Table 2, were assigned reference values from the literature within feasible ranges of variation (Ivanov et al., 2004a, b;Moreno et al., 2012;Rutter et al., 1971). Subsequently, daily time series of streamflow (Q) and snow water equivalent (SW) were used as targets for model calibration, during the 10-year period N = [01/01/1991,12/31/2000], selected to include important drivers of seasonal and interannual climate variability including winter frontal, monsoonal systems, Pacific decadal oscillation (PDO) and ENSO events (Dominguez et al., 2010). For calibration, we implemented a model pre-emption framework (Razavi et al., 2010) to improve computational efficiency by terminating model runs in poorly performing parts of the parameter space. The Shuffled Complex Evolution (SCE) algorithm (Duan et al., 1993) was used to find optimum values within feasible ranges of variation that minimize the normalized residuals of simulated and observed time series of Q and SW, as dictated by the normalized objective function M(t) evaluated at each pre- Table 2. Model calibrated parameters for the period 1 January 1991 to 31 December 2010 at the Tonto Creek basin. Parameters for soil are saturated hydraulic conductivity (K s ) and its decay exponent with depth (f), pore-size distribution index (λ 0 ) and air entry bubbling pressure (ψ b ) and for vegetation, albedo (a), vegetation height (H v ) and optical transmission coefficient (K t ).

Soil type
with where w 1 and w 2 are optimization weights, σ ox is the standard deviation of observed values and x obs j , x sim j are the observed and simulated values during simultaneous time steps j.
Calibrated values, illustrated by Table 2, were then used to evaluate model robustness during the period 1 January 2001 to 31 December 2010. Figure 7 shows the daily observed and simulated time series of Q and SW for the calibration and evaluation periods with complementary information about model skill at the daily timescale, in terms of the meansquared error (MSE), Nash-Sutcliffe efficiency (NSE) and Pearson correlation coefficient ρ so . Together, these scores provide a complementary view of the model simulations in terms of the mean, variability and overall correlation. Figure 7 and skill scores suggest that despite certain discrepancies in the simulation of long recessions, and certain peak streamflows and snow water equivalent maximums, the model is able to reproduce the distinct hydrologic patterns that determine the presence of on-the-ground snow and mean and variability of stream discharge. As indicated before, the overall quality of hydrologic simulations is largely tied to the quality of hourly precipitation inputs, whose uncertainties propagate basinwise (Bardossy and Das, 2008;Borga et al., 2006;Michaud and Sorooshian, 1994). Model robustness is indicated by the evaluation scores, which summarize predictive capability during the entire 20-year period.

Design of numerical experiments
Our goal is to understand the individual and collateral effects of forest thinning and related soil disturbances due to forest removal operations on the total hydrologic response, using historic climate forcing. Modeling experiments were therefore conducted during the period 1 January 1991 to 31 December 2010 with adoption of a base-case scenario determined by 2006 soil and vegetation cover maps (Fig. 4). Forest thinning induces model changes in vegetation fraction, LAI, vegetation throughfall and canopy storage (see Appendix B). In all cases we assume that litter is also removed from the thinned areas and vegetation condition does not dynamically evolve during post-treatment conditions (see Sect. 4.5). Soil changes are fundamentally represented by modifications in the saturated hydraulic conductivity, which are triggered by compaction and water-repellency processes after mechanical thinning and prescribed burning. Given the high uncertainty in such values, as reported in the literature, we assume three additional cases of post-treatment steady reductions in original soil hydraulic conductivity (imperviousness; from Table 2) between 10 and 40 % (10, 20, 40 %) of the current values and only in the areas covered by ponderosa pine. Table 3 summarizes the simulation of scenarios and the corresponding adopted naming convention (case). While representing post-fire conditions as constant over time might be considered unrealistic, the results are indicative of the immediate sensitivity of basin response to a drastic (as planned) land cover change. Spatially distributed water footprints due to forest thinning can be understood through an elementscale view of the long-term shifts on water fluxes and stocks. This analysis is performed through the selection of multiple domain elements located within forest treated areas of different thinning intensity values; elements with contrasting solar aspects are paired according to similar elevation, slope, air temperature, wind speed, net radiation, evapotranspiration and soil moisture to compare their hydrologic evolution from pre-to post-treatment conditions. A total of eight element pairs were found to fulfill these requirements. For each element, the components of the water balance can be estimated as in the soil column schematic in Fig. 8, where surface and subsurface reservoirs and input/output fluxes have been included in annual ( t = 1 year) mass continuity equa-  Fig. 4), along with NSE, MSE and ρ so skill scores. To improve the visualization of low streamflow values, the time series of discharges have been elevated to a 0.5 exponent. Mean Areal and pixel Precipitation (MAP and P , respectively) are derived from the corrected NLDAS product. Post-treatment basal area Calibrated K s VS10 Post-treatment basal area 10 % reduction in K s across soil types in ponderosa pine areas VS20 Post-treatment basal area 20 % reduction in K s across soil types in ponderosa pine areas VS30 Post-treatment basal area 30 % reduction in K s across soil types in ponderosa pine areas tions (Eqs. 6-8). The different pre-and post-forest-thinning components of the water balance in the soil column are appraised to only evaluate the effect of thinning in contrasting hillslope aspects.
4 Results and discussion

Streamflow shifts and extreme event probability
Forest removal affects the distribution and magnitude of streamflows in a different manner depending on the seasonal magnitude of runoff generation, the shifts in INT, SW, θ and GW and the soil hydraulic conditions imposed by thinning operations. Field observations have shown an immediate decrease in soil hydraulic conductivity, but recovering to historic soil conditions with time, after forest treatment. This section addresses the concerns for increased flood risks during heavy rainstorms and sustained river water supply for urban populations and ecological processes during low discharge conditions, as a result of a vegetation-reduced system. According to the annual patterns of precipitation, temperature (Fig. 5) and streamflow (Fig. 11a), there are three differentiable conditions in Tonto Creek characterized by the (1) wetter, higher flows during winter (e.g. January) season and (2) the summer monsoon (e.g. August), and (3) drier, low flow circumstance during the pre-monsoon period Figure 8. Soil column water balance storages and fluxes of a typical hillslope computational element. The computational element's Voronoi geometry has been represented by a rectangular shape in the interest of simplification. Water is mostly stored through vegetation interception (Int), snow accumulation (SW), vadose zone soil moisture (θ ) and groundwater in the saturated zone (GW). Surface and subsurface water (in and out) fluxes include above canopy gross precipitation (P ), vegetation transpiration (T ), evaporation from intercepted water (E int ), evaporation from soil (E soil ), sublimation from intercepted (S int ) and on-the-ground snow (S snow ), net surface (R = R in − R out ) and subsurface runoff (θ f = θ in − θ out ) and net groundwater flow (GWf = GW in − GW out ). The column is constrained by an impervious bedrock layer, whose depth varies from element to element.
(e.g. June). Hourly time series from the reference and simulated cases are classified by hydrologic period (winter, premonsoon, monsoon and all months included) to understand the probability distribution shifts that forest thinning produces on quartiles, Q 1 through Q 4 (where Q 1 and Q 4 correspond to low and high flows, respectively) and low-order statistical moments (µ, σ ) of long-term (20-year) simulations (Fig. 9). Results are expressed in terms of ratios relative to distributional values obtained by the reference case for each type of hydrologic condition.
Model results indicate that Q 1 , µ, σ and Q 4 are larger across cases, confirming not only the higher runoff efficiency but also the increased flood risk for riverine communities during the winter season under post-treatment conditions. In contrast, during the monsoon season, differences in the soil hydraulic conductivity play a major role in the distribution of streamflow values. For instance, V and VS10 produce net reductions in µ, σ and Q 4 ; increases in the same statistics are observed for the most impervious cases (VS20, VS40). In the long term, if hydraulic conductivities return to normal, it might mean reductions in the mean and extreme runoff production during monsoon showers. On the other hand, during pre-monsoon conditions, forest thinning seems to be increasing the lowest streamflows, but has a mixed effect on µ, σ , Q 3 and Q 4 . In these cases, the less impermeable scenarios achieve reductions in distribution values, indicating drier hydrologic conditions, while the most permeable scenarios (VS20, VS40) evidence increases in the same distributional parameters.
Results for all months together suggest net increases in Q 1 , µ, σ and Q 4 , indicating a net distributional shift to the right, relative to the reference case. These changes in distributional values of streamflow triggered by land cover changes may support the need for decision making oriented towards water preservation during dry conditions and mitigation or adaptation of the negative effects of floods on urban settings and ecological communities.

Effects of forest thinning on mean and variability of basin-scale water balance components
Hydrologic effects of headwater forest thinning are reflected through both local changes in the mean and variability of water fluxes and stocks and basin-scale shifts in discharge yield.
The following analysis supports this statement by quantifying the magnitude and direction of the water changes that are statistically significant at the basin scale. First, an interannual examination is conducted to understand shifts in key water variables and their patterns, both in the long-term and during warm and cold phases of ENSO. Like the entire southwestern US, the Tonto Creek basin experiences increases in total annual P during El Niño (by 20 %) and reductions during La Niña (by 11 %), with both phases presenting slight reductions in mean air temperatures (Temp), as estimated from NLDAS corrected 20-year records ( Fig. 10a and b). Since water balance is affected by ENSO, alterations in the basin's response to forest thinning are also expected. In addition to inter-annual variations, seasonal shifts are expected as modifications in the below-canopy energy balance, wind speed and net precipitation impose differential effects according to each month's climate regime. Results are presented in terms of each simulated case relative to the corresponding reference scenario, and for each ENSO phase type, as 20-year mean and standard deviation ratios and monthly absolute differences.

Inter-annual trends
In the long-term, forest thinning leads to changes in water distribution that are exacerbated during an ENSO event. Results suggest increased annual average streamflows (Q) of up to 7 %, but reductions of SW of 16 % and snow covered area (SA) of about 5 % (Fig. 10c), while only slight reductions Figure 9. Long-term ratios (Q case / Q ref ) between streamflow probability distribution properties for the forest thinning scenarios and the reference case, computed from hourly simulated time series for typical winter (January), pre-monsoon (June) and monsoon (August) months and all months. Statistical properties include first, second, third and fourth quartiles (Q 1 , Q 2 , Q 3 and Q 4 ), mean (µ) and standard deviation (σ ). In all plots, the dashed line represents the reference case.
(less than 2 %) in vadose zone soil moisture and evapotranspiration (θ and ET) are observed. Similarly, 10 cm and root zone soil moisture (θ 10 and θ root ) and depth to groundwater (DG) do not show significant changes, relative to the reference case. Comparatively, thinning simulation cases differentially impact the mean Q, with VS40 being the most efficient in increasing runoff through a decrease in soil infiltration capacity. In addition, temporal hydrologic variability appears to be dampened by forest thinning, with the exception of streamflow, as illustrated by the lower time series standard deviations of Fig. 10d. Interestingly, ENSO appears to modulate these shifts by exacerbating or moderating forest thinning impacts. For instance, El Niño enhances direct surface responses in Q and θ 10 and compensates for the losses in SW and SA. In contrast, La Niña dramatically reduces Q, SW and SA (See Fig. 10c). In terms of time series variability, ENSO appears to intensify reductions in inter-annual variability due to forest thinning across the tested variables, with the exception of ET during La Niña and SA during EL Niño, as illustrated by Fig. 10d. A seasonal analysis (next subsection) facilitates identification of the emerging monthly patterns responsible for these inter-annual trends.

Seasonal shifts and emerging hydrological patterns
At the monthly scale, forest thinning increases streamflows and groundwater recharge, at the expense of reduced interception and snowpack, and a pattern emerges of a less regulated runoff system that exacerbates both higher and lower levels of river flow. At Tonto Creek, the high precipitation and low air temperatures during winter months drive the unimodal annual cycle of Q and SA with maximum values in January of each year (Fig. 11a, f). Cumulative effects of this wetter period are also observed through delayed responses of θ , DG and SW, with maximum peaks appearing in March (Fig. 11b, c, e). Comparatively, the second rainfall peak only produces Q values below the annual mean, as most water leaves the basin through higher ET rates, a typical behavior of water-limited basins (Fig. 11a, d). For the most part, forest thinning tends to increase Q, in particular for those months with already high runoff production and for those cases with the most impervious soils (i.e., DJF and VS40; see Fig. 11g). Nonetheless, during the monsoon season (JAS), changes in Q are less clear across thinning cases with the less impervious scenarios (V, VS10) instead showing net reductions in Q, even when ET values have been simultaneously reduced (Fig. 11g, j). The emerging shift in patterns of SW and SA reveal reductions that are more marked during their peak values (i.e. during MA; Fig. 11k, l). Aside from SW, vadose zone water availability (θ ) does not show significant changes during the year due to thinning (Fig. 11h). In contrast, the depth to groundwater decreases almost uniformly year round, with the least impervious scenario having the largest aquifer recharge values (Fig. 11i). On balance, reductions in snow water equivalent and, less likely, in evapotranspiration linked to vegetation removal, compensate for the increased (especially winter) runoff and groundwater recharge, resulting in an emerging pattern shifting from surface snow to groundwater storage, an issue in semi-arid basins, whose deep aquifers may remain disconnected from the channel base flows for many months of the year. A detailed spatial analysis (see the next sub-section) provides information about important local water trends for mountain ecosystems settled directly on thinned areas of the forest.

Distributed hydrologic effects of forest removal
As forest reduction will only occur in the headwaters of Tonto Creek basin (see Fig. 3), direct hydrologic effects are likely to be particularly marked in such areas, which are subject to higher annual basin precipitation and lower mean temperatures. Figure 12 presents the spatial hydrologic patterns for the reference case (first column) and projected changes for three representative cases (V, VS10, VS40; columns 2 through 4) relative to the reference. Results shown in Fig. 12 indicate that averaging over time, spatial differences due to changes in soil hydraulic conductivity (i.e. V vs. VS10 or VS40) are not salient among cases but rather that any level of forest removal imposes major changes in local water.
Runoff and soil moisture: in terms of runoff (R ref ), current rates attain the highest values in shrubland and low basal area ponderosa pine cover, as most of the water in forested areas is intercepted or bound up by snowpack for slower release to the channel network in the spring. Consistently across scales, thinning promotes increases in local runoff production, of up to 40 % in heavily thinned areas and for the most impervious case (VS40). On the other hand, storage of water in the vadose zone (θ ref ) is characterized by higher values in proximity to the channel network and riparian areas, particularly in high elevation areas, dominated by forest cover and higher rainfall values. Forest removal induces mixed shifts in θ, but a dominant reduction trend is observed in heavily thinned areas, with VS40 producing the largest reduction rates (of up to 15 %) in θ .
Evapotranspiration: coupled to soil moisture, atmospheric losses through evapotranspiration are evidently larger along the river network and riparian areas where ET consumes available surface and subsurface water through rates that equal annual precipitation (ET ∼ P ∼ 700 mm yr −1 ) in some riparian corridors. Except in heavily thinned transects with slightly higher temperature (Temp), where increases of up to 30 mm yr −1 in ET are seen, the vast majority of the thinned area indicates decreases in ET, of up to 40 mm yr −1 , presumably associated to reductions in transpiration rates (T). As impervious cases (i.e., VS40) produce increases in surface runoff production to river network, both θ and ET decrease more drastically.
Snow: in terms of snow processes, current conditions allow for the formation, accumulation and melt of on-theground snow differentially across the Mogollon Rim during the winter and spring months. In the case of the Tonto Creek Figure 13. Long-term element-scale shifts in mean water fluxes and stocks relative to the reference case during 20-year model simulations. Results are presented for (a) 7N-6S, and (b) 6N-7S, as representative element pairs with different thinning degrees and contrasting hillslope aspects. Tested cases (V, VS10, VS20, VS40) are differentiated by the geometric symbols aligned vertically for each variable with north represented by solid and south represented by hollow symbols. Water fluxes include runoff (R), groundwater flow (GWf), sublimation from on-the-ground snow (S snow ) and intercepted (S int ) snow, evaporation from soil (E soil ) and intercepted water (E int ), vegetation transpiration (T ) and total evapotranspiration (ET). Water stocks include vegetation interception (Int), on-the-ground snow water (SW), vadose zone soil moisture (θ ) and groundwater storage (GW). Auxiliary variables, including 2 m surface temperature (T s ), wind speed (WS), net radiation (NR) and soil moisture at 10 cm and root zone depths (θ 10 and θ root ), have been added to the plot to aid interpreting budget shifts. basin, exceptional wet, cold winter seasons result in a local maximum of 1000 mm snow water equivalent (SWmax ref ), with snowpack persisting (NDS ref ) for up to 170 consecutive days. Forest thinning consistently reduces NDS for as long as 60 days and SW max by as much as 350 mm, in the most intensively thinned areas by an increased forcing of shortwave energy on thinned patches. In summary, model simulations reveal that vegetation removal is the most important factor determining distributed changes in fluxes and storages of water, more so than hydraulic changes in soil. The Tonto Creek basin presents spatially distinct responses to forest thinning characterized by punctuated increases in runoff and generalized decreases in soil moisture, evapotranspiration and snow persistence and volume, compared to historically simulated levels. In the next sub-section, the physical mechanisms inducing change at the element level are explored in higher detail, through soil col-umn analysis of multiple computational elements with contrasting annual radiation differences.

Soil column water balance in hillslopes with contrasting solar aspect
This section aims to identify the effect of forest thinning in contrasting solar aspects. Figure 13 (a and b), summarizes two examples of the typical shifts in the soil column water balance terms as a proportion of the reference case. Although only two element pairs (7N, 6S and 6N, 7S) are shown, the balance of evidence indicates that forest thinning induces local increases in below-canopy P net (P -Int), NR, T s and WS, which trigger increases in R (exacerbated by soil imperviousness), S snow , and E soil , at the expense of reductions in GWf, Int, E int , T and SW. While, in general, the soil columns experience comparatively slight reductions in ET, one of the most evident shifts involves a compensatory partitioning with re- Figure 14. Mean annual cycles of simulated reference (black) and tested (colored) cases for an element pair (7N, 6S) as obtained from 20-year model results. Variables include atmospheric losses (ET + S) for all evaporation, transpiration and sublimation rates, net runoff production (R), snow water equivalent (SW), vadose zone soil moisture (θ ) and groundwater storage (GW). V, VS10, VS20 and VS40 are represented by blue, green, orange and red colors, respectively. Mean annual changes ( x) have been added to each variable to compare mean monthly changes relative to each reference case.
ductions in E int and T and increases in E soil in both hillslope aspects. The degree of thinning ( V F %) appears to elicit a direct and proportional influence on the relative change of NR, T s , Int, S snow , E soil , E int , T and SW across the eight pixel pairs. A more detailed scrutiny of these trends during a typical water year is illustrated by Fig. 14 for an element pair (7N-6S), and considers the most important fluxes and reservoirs ranging from atmospheric (ET + S), surface (R, SW) and subsurface (θ , GW) components. Table 4 shows mean total annual changes across the eight element pairs (N, S) for all tested cases. Figure 14 and Table 4 indicate that larger reductions in the total atmospheric losses (ET + S) can be achieved in the north-facing slopes, particularly for the most impervious cases (e.g., VS20, VS40) and more marked during the first ET peak in March. Additionally, larger gains in runoff are achieved from the north-facing slopes especially during the winter peak and more significantly for the most impervious cases (i.e. VS20, VS40).
Regarding water reservoirs, reductions in snow water due to forest thinning are far larger for south-facing slopes where four elements (1S, 5S, 3S, 7S) evidence total depletion of Table 4. Mean annual differences between forest thinning scenarios (V, VS10, VS20, VS40) and reference case for atmospheric losses (ET + S), runoff (R), snow water equivalent (SW), vadose zone moisture (θ ) and groundwater storage (GW) across eight element pairs with contrasting (north, south) hillslope aspects.

North aspect
South aspect  48 −81.48 −81.48 −81.48 −197.54 −197.54 −197.54 −197.54 µ( θ )  snowpack between 15 and 25 days earlier than during reference conditions. The trade-offs between less snow and faster melt mechanisms are clear through the increase of element runoff and a greater recharge (GW) of the aquifer, whose groundwater table levels appear deep and sometimes disconnected from the surface channel network. The interplay of θ and GW is clear when comprehensible increases in saturated thickness lead to corresponding reductions in vadose zone water storage in the bedrock-limited soil column element.

Model assumptions and study limitations
This section explains a set of important model assumptions and limitations that help with the interpretation of the results, estimation of the scope and identification of potential lines for future work from this study. The following items are presented without an order of importance as the amount of uncertainty introduced by each of them was not quantified in a systematic fashion.
1. The model does not consider dynamic changes in vegetation physiology, re-growth and/or mortality rates. This assumption ignores the actual (probable) response of vegetation to post-treatment conditions, if thinning operations are discontinued. This would include, but is not limited to, progressive increases in basal area (and thus sapwood area), concomitant linear increase in projected leaf area index for conifers (McDowell et al., 2008) and the accompanying physiological, radiative and hydraulic responses of the overstory and understory vegetation (dePury and Farquhar, 1997;Ivanov et al., 2008;Sampson et al., 2006) being ignored. Notwithstanding, typical growth rates (woody increment) at this geographic region are of about 2 % per year, depending on the species (Worley, 1965), and so likely canopy processes would be slow to respond during the modeling period considered in this study. A misrepresentation of the vegetation evolution during post-treatment time would, most likely, result in underestimation of interception capacity and on-the-ground snow duration but overestimation of runoff rates.
2. The model does not consider gradual recovery in soil saturated hydraulic conductivity during the posttreatment condition that would, most likely, result in reduction of runoff volumes but increases in vadose zone soil moisture.
3. The uncertainty propagation from the NLDAS precipitation product to the hydrologic simulations and the lack of "ground-truth" hydrologic information (i.e. rain gauges, nested streamflow gauges, snow, evapotranspiration and soil moisture stations) hinders the entire validation process and simulation skill and constrains the comparison to only a few measuring stations of streamflow and snow. This fact seriously constrains extrapolation of results to other variables that were not verified during this modeling effort. Nonetheless, results can be fully understood relative to a base-case scenario that aimed to reproduce hydrologic conditions as real as possible.
4. Finally, the model did not analyze the effects of forest thinning in sediment and pollutants load to streams and reservoirs. Further studies should investigate the combined effects of deforestation and their subsequent shifts in water residence times from surface to groundwater reservoirs.

Summary and conclusions
This study investigated the long-term effects of simulated forest thinning for both element and basin-scale hydrologic balance and extreme discharges in a semi-arid basin of the southwestern US. We used the 4FRI forest restoration project as the context for these silvicultural operations applied to the headwaters of Tonto Creek along the Mogollon Rim, the most water productive region in Arizona. Long-term hydrologic simulations in this basin are challenging due to topographic complexity as well as the lack of ubiquitous hydrologic measurements on the terrain. In appraising the spatiotemporal water footprints of forest removal, we investigated the changes induced in the probability distribution functions that involve mean and extreme discharge events in long-term and during three distinct seasonal hydrologic conditions. The mechanisms that support this shift behavior are explored through an analysis of the inter-annual and seasonal effects on the mean and variability of hydrologic variables and the water-related outcomes induced by the occurrence of ENSO phases. Finally, an emphasis was placed on identifying the mechanism through which water transitions occur due to changes in the solar radiation, surface temperature, wind speed and water balance at the element scale in contrasting slope aspects. Our results are summarized below.
1. Forest thinning leads to a less regulated hydrologic system for mean and extreme events. A probabilistic analysis of the magnitude of recurrence of mean and extreme event conditions indicates a net increase in the annual streamflow distributions, particularly dominated by larger, consistent increases in mean and maximum events during the winter months. This shift can increase the risk of negative flood-related effects directly downstream of the treated areas. For the less impervious scenarios (V, VS10), consistent increases in runoff are not observed for the mean and higher quartiles during the dry and low flows of the pre-monsoon and monsoon seasons, leading to an even drier hydrologic condition. Consistently across seasons, impervious soils contribute to increased streamflow values.
2. Headwater forest thinning can lead to hydrologic shifts in the areas directly affected by this procedure that are reflected by anomalies in the average basin-scale integrated values. Observable basin-scale changes occur through increases in runoff (7 %) and decreases in snow-water (−16 %) and snow-covered areas (−5 %). This result is consistent with recent observations in high elevation forests (Metcalfe and Buttle, 1998;Musselman et al., 2008;Lindquist et al., 2013;Venkatamaran, 2013). Increases in soil impermeability due to removal operations exacerbate alterations, particularly in runoff volume. Climatic stressors like ENSO affect the magnitude of such re-distributions, principally through modifications in precipitation availability. For instance, El Niño appears to exacerbate runoff production while La Niña reduces snow presence due to a rainfall suppression effect.
3. At the monthly scale, forest thinning results in streamflow augmentation, particularly during the winter precipitation peak but less clearly for the monsoon season, when the most permeable scenarios instead decrease runoff yields, on average. Conversely, consistent reductions in the depth to groundwater (maximum in January), evapotranspiration (maximum in July) and snow water (maximum in April), are observed across simulated scenarios, thus lowering the historic maximum values occurring in corresponding months.

The Tonto
Creek basin presents spatially distinct responses to forest thinning characterized by local increases in runoff and generalized decreases in interception, soil moisture, evapotranspiration and snow persistence and volume, when compared to the current reference case. In terms of runoff, local increases in runoff production in heavily thinned areas and for the most impervious case (VS40) are observed. In contrast, mixed shifts in θ , but with a dominant reduction trend, are observed in heavily thinned areas, with VS40 producing the largest reduction rates. Regarding ET, except for a few increasing trends in heavily thinned transects with slightly higher surface temperature (Temp), the vast majority of the thinned area indicates decreases on ET associated with reductions in transpiration rates (T). Because impervious cases (i.e., VS40) impose increases in surface runoff production to the river network, both θ and ET decrease more drastically in this case. Forest thinning consistently reduces snow persistence and peak values in intensively thinned areas.
5. Multiple element soil column analysis indicates that gains in runoff and aquifer recharge are due to net reductions in interception, snow water equivalent and, less likely through reductions in evapotranspiration. Removal of forest canopy shading creates a nearly balanced mechanism where decreases in transpiration are compensated by increases in soil evaporation rates. The annual net radiation imbalance between north-and south-facing slopes in this north-latitudinal basin results in increased vulnerability of south-facing areas to less snow accumulation and faster melt periods by increases in surface temperature, sublimation and evaporation rates.
Despite this modeling study not considering the recovery of vegetation dynamics (e.g. re-growth) and soil hydraulic properties during the long-term simulations, the use of highly credible (Hampton et al., 2011) forest thinning projections and three additional simulation scenarios considering increases in soil imperviousness provide one set of reasonable, spatially distributed cases to identify potential effects on the mean and extreme hydrologic conditions in this semi-arid region. This situation could, specially, apply if authorities decide to maintain forest thinning operations in the long term. The results of this study are based on the use of a distributed hydrologic model that was calibrated and verified during 20 consecutive years at daily scale, using 12.5 km, 1 h resolution climate forcing from the NASA Land Cover Data Assimilation System (NLDAS; Mitchell et al., 2004) with precipitation fields locally adjusted through rain gauge data. The tuning and evaluation procedures both provided appropriate skill scores for streamflows and snow water equivalent, despite some discrepancies introduced by model forcing, initial conditions and structural errors. While calibration and validation coefficients are not optimal, model performance of-fers the possibility of quantifying changes introduced by forest thinning, independent of the model structural and parametric uncertainty, as results are primarily presented relative to model simulations made with 2006 vegetation conditions, which we adopted as current reference case.

Appendix A: Precipitation bias correction
While a global bias correction procedure (Steiner et al., 1999) provided poor rainfall adjustments, a modified local correction strategy (Seo and Breidenbach, 2002) produced much better hourly rainfall estimates due to the high spatial variability of this phenomenon. This approach uses the three closest daily ground rain gauges to correct hourly volumes of the NLDAS gridded precipitation product (R) pixels following a weighting strategy according to the following expressions: where β i = 1 if g i / r i > β t g i / r i if g i / r i < β t , (A2) r c is bias-corrected R (mm), r o is raw R at the pixel centered at µ o (mm), w i is weight given to the R-gauge pair at the ith vertex in the triangle of R-gauge pairs that encloses µ o , β i is multiplicative sample bias from the ith R-gauge pair, δ i is additive sample bias from the ith R-gauge pair, g i is gauge rainfall measurement (mm) at the ith vertex in the enclosing triangle, r i is collocating R rainfall estimate (mm) at the ith vertex in the enclosing triangle and β t is adaptable parameter that denotes the threshold for the multiplicative or additive bias.
The neighboring R-gauge pairs are identified by triangulation, which connects all available R-gauge pairs into a mesh of triangles. The weights, w i , i = 1, 2, 3, sum to unity and are inversely proportional to the distance to the neighboring R-gauge pairs in the enclosing triangle. An iterative procedure was conducted to select the best β t =1 that minimized the MSE between observed and corrected rainfall at rain gauge locations. Figure A1 illustrates the spatial distribution of precipitation for the VTS system, averaged during 21 years (1990 to 2010) as measured by (a) Thiessen polygons derived from 30 daily rain gauges, (b) raw NLDAS and (c) bias-corrected NLDAS estimations. Figure A2 shows an example scatter plot comparing daily rain gauge values (x axis) with raw and corrected NLDAS (y axis) for 1 of the 30 stations within the study region.  A set of empirical relations are used to relate remote sensing and field information to vegetation parameters and processes in our hydrologic model. Such processes account for vegetated fraction (VF), LAI, throughfall (p), and canopy storage (S) and below canopy light (Q / Q 0 ) and wind speed attenuation (U h ), in ponderosa pine forests. From historical measurements in northern Arizona across seven different ponderosa pine forest densities, 90 m resolution vegetation fraction maps were derived for pre-and posttreatment basal area conditions only (i.e. ignoring plant evolution or phenology), as reported by the small-diameter wood supply report (Hampton et al., 2011) with the following expression: where VF is the vegetation fraction (%) and BA is the measured basal area (ft 2 / Ac). LAI maps for ponderosa pine were derived following an empirical relation with basal area from field measurements in 10 study sites with different pine densities (Armstrong, 2012) through a relation that minimized residuals between observed and predicted LAI: LAI values were verified on typical ranges for ponderosa pine forests under different vegetation fraction conditions. Vegetation fraction and LAI values for non-ponderosa covered areas were extracted from the 2006 Landfire products (http://landfire.gov/) and derived from existing literature (Mendez-Barroso et al., 2013;Mitchell et al., 2004), respectively. Free throughfall coefficient (p), which accounts for the fraction of rainfall not captured by plants, and canopy capacity (S), were derived from the expressions B3 and B4 (Carlyle-Moses and Price, 2007; Mendez- Barroso et al., 2013;Pitman, 1989): The Beer-Lambert law was adopted to account for the reduction in radiative transmittance due to dense canopies (Brantley and Young, 2007;Marshall and Waring, 1986) following where K t is the optical transmission coefficient. Finally, below canopy momentum transfer by wind speed was corrected by forest density as surface rugosity factor following (Sypka and Starzak, 2013;Yi, 2008) Where U (h) is the wind speed at the height h within the canopy, in m s −1 , U H is the wind speed at the top of the canopy, in m s −1 , and H c is the top of the canopy, in m.