Overstory–Understory Vegetation Cover and Soil Water Content Observations in Western Juniper Woodlands: A Paired Watershed Study in Central Oregon, USA

The effects of western juniper (Juniperus occidentalis) control on understory vegetation and soil water content were studied at the watershed-scale. Seasonal differences in topsoil (12 cm) water content, as affected by vegetation structure and soil texture, were evaluated in a 96-ha untreated watershed and in a 116-ha watershed where 90% juniper was removed in 2005. A watershed-scale characterization of vegetation canopy cover and soil texture was completed to determine some of the potential driving factors influencing topsoil water content fluctuations throughout dry and wet seasons for approximately one year (2014–2015). We found greater perennial grass, annual grass, and shrub cover in the treated watershed. Forb cover was no different between watersheds, and as expected, tree canopy cover was greater in the untreated watershed. Results also show that on average, topsoil water content was 1% to 3% greater in the treated watershed. The exception was during one of the wettest months (March) evaluated, when soil water content in the untreated watershed exceeded that of the treated by <2%. It was noted that soil water content levels that accumulated in areas near valley bottoms and streams were greater in the treated watershed than in the untreated toward the end of the study in late spring. This is consistent with results obtained from a more recent study where we documented an increase in subsurface flow residence time in the treated watershed. Overall, even though average soil water content differences between watersheds were not starkly different, the fact that more herbaceous vegetation and shrub cover were found in the treated watershed led us to conclude that the long-term effects of juniper removal on soil water content redistribution throughout the landscape may be beneficial towards restoring important ecohydrologic connections in these semiarid ecosystems of central Oregon.


Introduction
The relationships between soil water content and vegetation cover are highly impacted by the ongoing shift from grassland to woodland-dominated landscapes occurring in many arid and semiarid regions worldwide. This change in vegetation-soil water dynamics has the potential to disrupt the ecological and hydrological balance of these water-limited regions [1][2][3]. This is particularly true in many semiarid landscapes across North America, from central Mexico to southwest Canada, where the significant expansion of juniper (Juniperus spp.) observed over the last two centuries is disrupting important ecohydrological functions [4]. The effects of Juniperus spp. cover on understory community structure are believed to have a domino effect on hydrologic processes [1,5]. Research in native grasslands of Oklahoma found Juniper spp. to be negatively correlated with soil water content, water storage, infiltration rates, and stream flow [6]. Further research involving juniper has emphasized its ability to utilize and influence horizontal and vertical soil water reserves throughout the intercanopy zones [7][8][9]. Encroachment tends to lead to an increase of bare ground, which significantly alters erosion and overland flow rates, depending on the amount of litter beneath the canopy patches [10][11][12].
Among the array of juniper species that have expanded across rangelands of the U.S., western juniper (Juniperus occidentalis spp. occidentalis Hook.) has been rapidly encroaching into grassland and sage steppe ecosystems of Oregon since the late 1800s and has increased from 170,000 ha in 1936 [13] to more than 1.4 million ha [14]. Western juniper densities are estimated to range from 100 to 600 individuals per hectare in the semiarid landscapes of central and eastern Oregon [15]. These numbers are a significant increase from the Euro-American settlement estimates of <5 individuals per hectare. This increase in western juniper woodlands has accelerated efforts towards the control of postsettlement populations with the intention of restoring ecologic values of rangeland resources [16]. Several studies conducted at the plot-scale show that western juniper removal could have positive results on several hydrologic processes. A study conducted by Mollnau et al. [17] found that juniper removal led to higher soil water content over the winter months, which was in part due to an increase in soil water recharge and a decrease in transpiration and interception rates. Also, western juniper removal has shown positive results towards reduced sediment yield and runoff, and increased infiltration rates and infiltration depth during plot-scale rainfall simulations [18,19]. Findings from these plot-scale research efforts improve understanding of important hydrologic functions. Yet, the underlying causal relationships that affect soil water dynamics in larger-scale juniper landscapes are not well understood. As stated by Ffolliot et al. [20], most soil water content related studies in juniper ecosystems have been conducted at the plot scale. Several authors [20][21][22] have determined the need for additional and more robust information regarding juniper encroachment and its effects on landscape-scale processes.
In this study, we aimed to enhance base knowledge on the long-term effects of western juniper control on vegetation and soil water content at the watershed scale. The objectives were to: (1) assess compositional vegetation differences between a treated (juniper removed in 2005) and an untreated watershed; and, (2) characterize topsoil water content variability across the spatial and temporal domains in both watersheds. We hypothesized that an increase in herbaceous vegetation would still be observed in the treated watershed ten years post juniper removal and that greater soil water content values would be observed in the treated watershed.

Description of Study Site
Our study site (43.96° Lat.; -120.34° Long.) is located 27 km northeast of Brothers, Oregon and it comprises one 116-ha watershed (treated) and one 96-ha watershed (untreated) (Figure 1). Elevation ranges from 1367 m at the outlet of the untreated watershed to 1524 m at the top of treated watershed. In the fall of 2005, juniper trees <140 years of age were cut using chain saws from the treated watershed and by the end of the summer in 2016, the boles were removed and the remaining limbs scattered. Old growth juniper trees and those that were host to wildlife were not removed [23]. Prior to juniper removal, tree canopy cover was estimated at 27% of total area. According to Fisher [24], the average percent slope for each watershed is around 25%, and the distributions of aspects is similar across both watersheds, with 35% north-facing slopes and 25% west-facing slopes. The orientation and drainage points of both watersheds are positioned in the northern portion of each watershed. Average annual precipitation (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017) at the study site is 358 mm. Most precipitation (69%) in the study area occurs as a mix of rain and snow between October and March, with intermittent rainfall events that occur during spring and summer, accounting for the reminder 31% of total precipitation [25].
Three major soil series, using the USDA official series description, Westbutte very stony loam, Madeline loam, and Simas gravelly silt loam comprise the majority of the soil types in both watersheds [24]. The Westbutte series is classified as loamy-skeletal, mixed, superactive, frigid Pachic Haploxerolls. The Madeline series is classified as clayey, smectitic, frigid Aridic Lithic Argixerolls. The Simas series is classified as fine, smectitic, mesic Vertic Palexerolls. The Westbutte and Madeline series are formed in colluvium from weathered basalt, tuff, and andesite materials and tend towards moderately shallow to deep, well drained soils. The Simas series is formed in colluvium and loess from tuffaceous sediments and tend towards very deep, well drained soils [26]. The treated watershed is primarily composed of 26% Westbutte, 48% Madeline, and 21% Simas series. The untreated watershed is composed of 50% Westbutte, 20% Madeline, and 3% Simas series [24].

Field Data Collection of Vegetation and Soil Properties
The line-point intercept sampling method, adapted from Herrick et al. [28], was used to estimate percent foliar cover, percent litter and percent bare ground represented in two individual sampling layers (i.e., top canopy and soil surface). In the summer of 2014, a total of 289 ten-meter transects were placed throughout both watersheds; 143 in the treated and 146 in the untreated (Figure 1). Transect locations were established to provide equal representation of aspect and elevation. Transects were permanently marked and aligned perpendicular to slope. The transects were established to avoid crossing ecotones between plant community types and soil types, thereby reducing risk of spatial heterogeneity effects as a result of differing abiotic and biotic factors. Vegetation points and soil surface cover were read every 1 m along the transect line. Canopy cover was recorded by species, and additional features were characterized as either herbaceous litter or woody litter (>5 mm). Soil surface measurements were of basal cover of plant species, rock (>5 mm), bedrock, moss, lichen crust, soil, embedded litter, or duff. Species functional groups were categorized as annual forb, perennial forb, annual grass, perennial grass, shrub, or tree. Data from each 10-m transect was used to estimate average percent total canopy cover by vegetative species, relative cover of each functional group, percent bare ground, and percent litter cover. Relative cover was estimated for each transect by dividing the sum of occurrences for each functional group by the sum of all occurrences.
Topsoil (12 cm) cores were collected in July 2014. Five soil cores were collected from each of the 289 transects using a soil-step probe (AMS, Inc.; American Falls, ID, USA), starting at the 2-m point. A few transects were comprised of rocky soils resulting in several unattainable cores, which reduced the number of samples from 1445 to a total of 1349. Each of the soil samples were analyzed in the lab for water content and soil texture. Sample values were calculated by averaging the results of the five core samples from each transect resulting in a total sample size of 289. Similar to soil core sampling, soil volumetric water content (θ) data was collected every two meters along each transect using a portable soil water probe (Model HydroSense II, Campbell Scientific Inc.; Logan, UT, USA), which integrates θ for the top 12 cm soil profile. The five individual θ measurements obtained were averaged for each transect, and this resulted in a total sample number of 289 samples (Treated, n = 143; Untreated, n = 146) used in the various θ analyses. In order to represent seasonal changes in θ, we collected data in July and November (dry season) of 2014, and during January, March, and May (wet season) of 2015. The duration of each data collection period was approximately three days, with the exception of the July reading, which took place over a four-week period. Data collected with the soil water content sensor were used to determine the temporal and spatial distribution of θ in both watersheds. Total precipitation for the duration of the study (July 2014 to May 2015) was 316 mm. Total rainfall for the months corresponding to θ data collection was, July (11 mm), November (64 mm), January (8 mm), March (38 mm), and May (50 mm).

Data Analysis
A single-factor analysis of variance (ANOVA) with two-sample t-tests was conducted (Treated, n = 143; Untreated, n = 146) to evaluate the differences in percent foliar cover, litter cover, bare ground and the relative cover of each functional group throughout the watersheds. The difference between means for percent relative canopy cover of each functional group was also analyzed. An additional single-factor ANOVA was utilized to test the significance (p ≤ 0.05) of the differences between mean θ for each watershed and at each measurement period (July, November, January, March, and May).
To determine the effects of terrain indices on θ, a single-factor ANOVA was utilized to test the significance (p ≤ 0.05) of the differences between aspects (represented in the eight cardinal directions) within each individual watershed.
Linear models were used to test the main effects and interactions of the measured variables between watersheds as well as within each individual watershed (p ≤ 0.05) using RStudio statistical software (RStudio; Boston, MA, USA). We used a general linear model (hereafter, the full model) to determine the effects of the watershed treatment (watershed), measuring period (month), total canopy cover (canopy), and soil clay content (clay) on mean θ. Within the full model we included canopy × clay, watershed × canopy, and watershed × month interactions to account for potential dependencies of these factors on one another. Mean values of θ, total canopy cover, and clay for each transect were utilized for model input variables. Clay content was used as the representative textural variable for the analyzed soil cores due to its influence on water holding capacity, and therefore in θ. We then utilized reduced linear models to determine the dominant effect of the independent variables, canopy and clay, on θ across each of the five measurement periods (July, November, January, March, and May).
The program ArcMap (version 10.2.2; Redlands, CA, USA) and the geospatial interpolation method, kriging, were used to demonstrate the spatio-temporal variability of mean θ for each monitoring transect (n = 289). Ordinary kriging (OK) was chosen over the other widely accepted geospatial kriging methods since the number of θ representative data points is relatively high (~7 points ha −1 ), providing an extensive spatial representation of the catchment areas. The general approach of the OK model determines statistical and spatial relationships among measured points to produce a prediction surface of the remaining unmeasured space. It is assumed within this model that predictions are possible due to the existence of spatial correlations, where points that are close in space will display similar soil water content values. The two-step OK method first determines the variance of the points against the mean values in order to fit the model and ultimately uses those values to create the prediction surface. The model was estimated with an 8 × 8 m grid size in ArcMap, and 289 data points, which represent the mean values of each transect resulting from the 1445 total θ points. No additional terrain indices were utilized since the range in elevation, slope, and the distribution of aspects was relatively uniform. It has been documented that elaborate interpolation methods are not necessary when landscape terrain indices are constant [29,30].
Results from the soil-particle size analysis showed four main soil types exist in the two watersheds: sandy, sandy loam, loamy sand, and sandy clay loam (Table A1). Characterization of soil hydraulic parameters such as θ at permanent wilting point (θPWP ) and field capacity (θFC) is important to understand the potential of soil water availability for plant growth. Pedotransfer functions (PTF), which are empirical relationships between the soil hydraulic properties and other physical properties such as soil texture, are often used to understand soil water flow relationships [31]. However, results obtained using PTF equations are always uncertain because PTFs were developed for particular datasets and their accuracy beyond that is unknown. One way to reduce uncertainty is to test several PTF equations using the same dataset. In this study, we used pedotransfer equations developed by Bruand et al. [32], Petersen et al. [33], and by Santra et al. [34], who recently developed PTFs using soils data from arid regions. We also used the 'rosetta' PTFs developed by Schaap et al. [35], which are widely applied worldwide. We used average sand, silt, and clay content values obtained for each of the four soil types identified to estimate θPWP and θFC using these PTF equations. Then we determined available water content (AWC) from θFC-θPWP for each soil type (Table A2). We compared our measured θ values to those obtained with the PTF equations to determine if they were within the expected range of θPWP and θFC values for the different soil types present in the two watersheds.

Vegetation
Total Cover and Relative Canopy Cover by Functional Group No significant differences (p > 0.05) in total canopy cover (66% treated vs. 61% untreated watershed) were observed between watersheds. Greater litter cover (46%) and less bare ground (16%) cover were observed in the treated watershed ( Figure 2). Analysis of major functional groups showed the treated watershed to have higher relative cover of perennial grasses, shrubs, and annual grasses when compared to the untreated watershed ( Figure  3). There were significant (p ≤ 0.05) differences in annual and perennial grasses, tree, and shrub relative cover between the treated watershed and the untreated. There were no significant (p > 0.05) differences in forb canopy cover between the treated watershed (3%) and the untreated (2%). Juniper canopy cover was 31% in the treated watershed and <1% in the untreated (Figure 3).

Soil Water Content by Month and by Aspect
Soil volumetric water content values that were obtained from the analysis of the soil core samples collected at the beginning of the study ranged from 2% to 15% in the treated watershed and from 3% to 12% in the untreated watershed. Results from the t-test analysis showed no significant differences (p > 0.05) in the treated watershed between mean θ values obtained from the soil cores (7.8% ± 0.198%) versus mean θ values collected with the portable soil sensor (8.0% ± 0.187%) in July 2014 when the soil samples were collected. However, in the untreated watershed there were significant differences (p ≤ 0.05) between mean θ values obtained from the soil cores (8.6% ± 0.203%) and mean θ values from the portable soil sensor (7.1% ± 0.142%).
Increases in mean θ corresponding to the transition from the dry to wet season were observed in both watersheds throughout the study period (Table 1). In three (July 2014, January 2015, and May 2015) out of the five months evaluated, greater θ values were observed in the treated watershed than in the untreated. No differences in θ were observed in November 2014, and a greater θ value was obtained in the untreated watershed in March 2015. When statistically significant (p ≤ 0.05), differences in θ between the treated watershed and the untreated ranged from 1% to 3%.   Table 1) were generally higher than mean θFC (17%). Based on the results from the PTF equations, the amount of water available was relatively the same for all soil types and ranged from 7% in sandy soils to 10% in sandy clay loam soils (Table A2). Our PTF-based results are consistent with those reported by Abdallah et al. [36] who used a pressure chamber to obtain θFC (24%) and θPWP (5%) from sandy loam soil samples collected in a rangeland location near our study site.
When assessing the influence of topographical aspect on topsoil θ, we found significant differences in θ by aspect for the months of July 2014 (p < 0.001), January 2014 (p < 0.01), and May 2015 (p < 0.001) in the treated watershed, and for the months of January 2015 (p < 0.01), March 2015 (p < 0.001), and May 2015 (p < 0.001) in the untreated. In general, the treated and untreated watersheds displayed similar θ patterns during the wet season, with slightly lower mean θ values in the southern aspects. This effect was much more apparent in the untreated watershed; where mean θ in the southern and western aspects were significantly lower (p ≤ 0.05) than in the northern aspects during the wet season ( Figure 4). Soil texture by aspect was similar for both watersheds (Table A3).

Spatial Distribution of Soil Water Content
Using GIS analysis and the ordinary kriging method, we generated θ surfaces across spatial and temporal scales for both watersheds. The interpolated θ classifications ranged from 5% to 11% in July 2014; 6% to 12% in November 2014; 12% to 30% in January 2015; 15% to 32% in March 2015; and 16% to 40% in May 2015. The range of interpolated θ values was consistent with the range of measured θ values derived from the portable sensor readings for each measurement period. Figure 5 shows the seasonal changes in θ from the dry to the wet seasons in both watersheds. During the dry season months (July 2014 and November 2014), greater θ values were observed in the higher elevation hillsides in each watershed. As topsoil conditions became progressively wetter throughout winter and spring, greater θ values were observed near the stream channels and at the valley bottoms. The effects of subsurface lateral flows on θ at lower elevation or near channel sites during the wettest months (March 2015 and May 2015) were highly visible during field data collection. This was even more apparent in the treated watershed, as indicated by the θ surface image corresponding to the May 2015 records.

Canopy Cover and Clay Content Effects on Soil Water Content
Full linear regression model results (data not shown) produced significant watershed × month interactions (p ≤ 0.05), and reported the canopy × clay and watershed × canopy interactions as nonsignificant (p > 0.05). Therefore, we used reduced models for each watershed and month combination to understand the source of the significance of the watershed × month interaction in the full model.
Analysis of the untreated watershed for the July 2014 soil reading showed the independent variable of canopy cover as having a significant influence on soil water content (F = 7.94, p < 0.01) and a non-significant effect of clay content (p > 0.05) ( Table 2) showed a significant effect of canopy for both the untreated (F = 2.37, p ≤ 0.01) and treated (F = 8.09, p ≤ 0.01) watersheds. Clay content in the treated watershed was also a significant (p ≤ 0.01) factor affecting mean θ values in January 2015. Interestingly, results for March 2015 showed non-significant effects of canopy for either watershed (p > 0.05), whereas clay content showed a significant effect on θ in both the untreated (F = 5.05, p ≤ 0.001) and the treated (F = 17.67, p ≤ 0.00) watersheds. The final analysis for the independent variable effects for the May 2015 reading resulted in results equivalent to the March 2015 reading. The effect of canopy cover on θ was not significantly different (p > 0.05) in both watersheds, while the effect of clay content on θ was significantly different for both the treated watershed (p ≤ 0.00) and the untreated (p ≤ 0.05) ( Table 2). Results derived from the linear models showed significant influences of total canopy cover on θ readings for July 2014, November 2014, and January 2015. A graph of one representative dry period (July 2014) and one representative wet period (March 2015) is shown in Figures 6 and 7 to illustrate the most prominent trends in mean θ by transect based on mean canopy cover for each transect. This approach was taken to help explain the relationship of θ by categorical values of total canopy cover. Figure 6 shows the decreasing trend in θ in July 2014 as total canopy cover increases in both watersheds. Results from the additional dry month (November 2014) displayed similar graphical analyses (data not shown). March 2015 resulted in an increasing trend in θ with increasing total canopy cover in both watersheds (Figure 7). The additional wet periods (January 2015 and May 2015) displayed similar graphical results (data not shown). This graphical analysis supports the results of the reduced linear model described above, which show there is a negative (p > 0.05) correlation between canopy cover and topsoil θ in the dry season months and a positive (p ≤ 0.05) correlation during the wet season months.

Discussion
The body of literature regarding large-scale understanding of ecohydrologic processes associated to the removal of juniper is still limited [21,22,37,38]. Soil water related studies in juniper woodlands have been conducted mostly at the plot-scale and there is an important need to understand soil-plant hydrological interactions at a larger scale [20]. This study provides a watershed-scale understanding of overstory-understory vegetation cover and topsoil (12 cm) water content relationships occurring in western juniper ecosystems.
Study results show greater perennial grass, annual grass, and shrub cover in the watershed where 90% of the overstory vegetation (i.e., western juniper) was removed in 2005. Tree canopy cover results (<1% treated watershed; 31% untreated watershed) are similar to those reported by Bates et al. [39] for a 25-year study conducted in southeast Oregon where they found average tree canopy cover of 0.8% in treated plots and 29.6% in untreated plots in year 13 of their study. They found that in year 25 after the treatment, juniper canopy cover was at 3.8% in the treated plots and it remained at 29.8% in the untreated plots. Given the proximity of our study site to the one reported by Bates et al. [39], and the very close difference in tree cover data observed, we believe a similar pattern of juniper density recovery would be expected in our treated watershed.
Our study results indicate that statistically significant, however marginal, positive differences on mean soil water content were generally observed in the treated watershed versus the untreated. These results (1% to 3%) are similar to those found in other juniper studies by Skau [40], and Everett and Sharrow [41], who reported soil water content differences of 1% to 2.5% and 2% to 5%, respectively. Similar soil texture conditions (Table A1) found in both watersheds may have contributed to the close soil water content conditions observed. The changes in soil water content observed during the transition from the dry to the wet season were affected by density and type of vegetation cover. Total canopy cover, and most notably juniper canopy cover, was inversely related to soil water content. Dense tree canopy cover commonly observed in highly encroached juniper landscapes, similar to that in our untreated watershed, can intercept significant amounts of precipitation; therefore, limiting the amount of water reaching the ground [25,42,43]. The role of overstory canopy cover, in regards to soil water distribution, seemed to be less important from late winter to early spring after soil water content had progressively increased throughout the wet season. Aspect was a factor influencing soil water content within each of our watersheds. In general, southern aspects receive more sunlight and become drier and warmer while northern aspects retain more moisture and are cold and humid [44]. Similar to the findings by Westerband et al. [45], who conducted a study on aspect-soil moisture relationships in juniper woodlands of New Mexico, soil water content in our two watersheds' southern aspects, and western aspect in the untreated watershed, were significantly lower than in northern and eastern aspects.
Juniper canopy cover has the potential for creating a protective barrier against solar radiation while promoting microsites and cooler soil temperatures [46]. The shading effects of juniper canopy on soil water content were more evident in the March 2015 measurements, when spring warmer temperatures dried the exposed soils in the treated watershed faster. This is consistent with a more recent study [25] we have conducted at this location, where we have been able to document greater soil water content levels occurring at topsoil (20 cm) depths under the canopy vs. the inter-canopy during early spring. This is also consistent with other studies that have documented higher soil water content under the canopy of juniper trees at certain times during the year and have attributed this condition to the interception of solar radiation interception [47,48]. Our watershed-scale results are different from those in a study conducted in western juniper in Idaho [49] where authors found no difference in topsoil water content between under canopy and inter-canopy locations.
While soil physical properties below the upper 12 cm profile evaluated in this study may be different across the landscape, these topsoil properties (i.e., moisture and texture) can provide valuable information regarding hydrological processes such as infiltration and runoff [6,48,50], which are important determinants of understory vegetation composition across the spatial domain in the study site. In arid environments, soil water content is a critical resource that largely determines vegetation community scale, structure, and diversity [51]. The spatial distribution of vegetation cover is both a cause and a consequence of soil water available [52]. The greater herbaceous vegetation and shrub cover found in the treated watershed at our study site indicates there are long-term benefits from soil water content redistribution through the landscape that can be associated with the removal of juniper over a decade ago. These reestablished ecohydrologic connections are beneficial in restoring important ecosystem functions that may be impaired by the high levels of juniper encroachment that commonly occur throughout rangeland ecosystems of the western United States.
This study adds to the body of knowledge by providing a watershed-scale understanding of the long-term effects of juniper removal on soil water and vegetation cover in a semiarid woodland ecosystem in central Oregon, in the Great Basin region of the United States. Results from this and other recent research [25,53] we have conducted at this study site show there are significant ecologic and hydrologic benefits associated with western juniper control. However, it is also evident that juniper is reestablishing at a relatively rapid pace in the treated watershed. Future work includes evaluating secondary treatment options to expand the longevity of the treatment.

Conclusions
Western juniper control provides long-term benefits to the landscape including greater herbaceous vegetation and sagebrush cover that can be both a result and a cause of greater soil water availability. This study provides important information regarding watershed-scale characterization of seasonal soil water variability and overstory-understory vegetation cover in juniper-dominated ecosystems.
Author Contributions: Grace Ray and Carlos G. Ochoa developed the study design and conducted field data collection and analyses. Tim Deboodt and Ricardo Mata-Gonzalez provided expert knowledge used in data analysis and interpretation. All co-authors contributed to the writing of the manuscript.

Acknowledgments:
The authors gratefully acknowledge the continuous support of the Hatfield High Desert Ranch, the U.S. Department of Interior Bureau of Land Management-Prineville Office, and the OSU's Extension Service, in this research effort. Our thanks go to John Buckhouse, Michael Fisher, and Mike Borman who had provided critical guidance throughout this long-term study journey. In addition, we want to thank the multiple graduate and undergraduate students from Oregon State University, and volunteers, who participated in various field data collection activities related to the results here presented. This study was funded in part by the Oregon Beef Council, Oregon Watershed Enhancement Board, USDA-NIFA, and the Oregon Agricultural Experiment Station.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix
Most of the soil samples collected were categorized as sandy loam (66.9%), followed by sandy clay loam (25.9%), then by loamy sand (6.5%), and by sandy (0.7%) soil types. Table A1 shows mean sand, silt, and clay content values Table A2 shows mean θ values at permanent wilting point (θPWP), field capacity (θFC), and available water content (AWC) for each soil type identified at our study site.  Table A2. Mean θ at permanent wilting point (θPwp) and field capacity (θfc) for the four soil types identified in both watersheds. Mean θfc and θPwp were calculated with pedotransfer equations based on mean soil particle distribution for each soil type. Available water content (AWC) is θFC-θPWP.  [35] <1 5 5 In general, soil texture were similar across aspect for both watersheds (Table A3). Particle size distribution in the treated watershed ranged from 72% to 76% for sand, 9% to 11% for silt, and 15% to 18% for clay content. In the untreated watershed, particle distribution ranged from 71% to 75% for sand, 9% to 12% for silt, and 16% to 18% for clay.