Carbon and ecohydrological priorities in managing woody encroachment: UAV perspective 63 years after a control treatment

Woody encroachment, including both woody species expansion and density increase, is a globally observed phenomenon that deteriorates arid and semi-arid rangeland health, biodiversity, and ecosystem services. Mechanical and chemical control treatments are commonly performed to reduce woody cover and restore ecohydrologic function. While the immediate impacts of woody control treatments are well documented in short-term studies, treatment impacts at decadal scales are not commonly studied. Using a controlled herbicide treatment from 1954 in the Sierra Ancha Experimental Forest in central Arizona, USA, we quantify woody encroachment and associated aboveground carbon accumulation in treated and untreated watersheds. Woody encroachment and aboveground carbon are estimated using high resolution multispectral images and photogrammetric data from a fixed-wing unmanned aerial vehicle (UAV). We then combine the contemporary UAV image-derived estimates with historical records from immediately before and after the treatment to consider long-term trends in woody vegetation cover, aboveground carbon, water yield, and sedimentation. Our results indicate that the treatment has had a lasting impact. More than six decades later, woody cover in two treated watersheds are still significantly lower compared to two control watersheds, even though woody cover increased in all four drainages. Aboveground woody carbon in the treated watersheds is approximately one half that accumulated in the control watersheds. The historical records indicate that herbicide treatment also increased water yield and reduced annual sedimentation. Given the sustained reduction in woody cover and aboveground woody biomass in treated watersheds, we infer that the herbicide treatment has had similarly long lasting impacts on ecohydrological function. Land managers can consider legacy impacts from control treatments to better balance carbon and ecohydrological consequences of woody encroachment and treatment activities.


Introduction
Woody encroachment into rangelands has been commonly observed around the world (Archer et al 2017, Stevens et al 2017. In the western US, for example, pinyon-juniper has encroached into rangelands at varying rates (Sankey and Germino 2008) resulting in woody cover increases of up to 600% in the 20th century (Romme et al 2009). Encroachment can accumulate large amounts of woody biomass and associated above-and below-ground carbon (Strand et al 2008, Rau et al 2011, Sankey et al 2013, Fusco et al 2019. Recent estimates have demonstrated that woody encroachment in the western US has resulted in twice the total aboveground carbon than previously estimated (Fusco et al 2019). However, the magnitude of the carbon increase varies among ecosystems due to edaphic variables and plant species mix. Detailed carbon accounting is missing in many regions, especially in dryland ecosystems. Consequently, woody contribution to global carbon pools and fluxes is unclear.
Land managers face competing priorities in managing woody encroachment Predick 2014, Bestelmeyer et al 2018). In addition to managing for increased carbon pools and sequestration, land managers have to consider the ecohydrological and biodiversity consequences: (a) Woody encroachment increases plant water use and evapotranspiration, which subsequently reduces downslope stream water yield, an important resource in arid and semi-arid ecosystems (Ingebo and Hibbert 1974). (b) Woody encroachment reduces herbaceous cover and forage for livestock and wildlife, and increases bare ground, erosion, and connectivity of sediment and runoff, which spatially redistributes sediment and nutrients across the landscape (Ludwig et al 2005, D'Odorico et al 2010, Wang et al 2019, Williams et al 2020, Sankey et al 2021. Taken together, woody encroachment and management thereof encompass four critical elements of the ecosystem: vegetation cover, carbon accumulation, water yield, and sediment yield (figure 1).
Many encroached areas are treated with controlled burning, mechanical thinning, or herbicide applications to increase water yield, reduce sediment runoff, and improve herbaceous cover and livestock forage (Romme et al 2009). Following treatments, the decrease in shrub canopy cover and increase in bare ground cover initially result in enhanced sediment availability and increased soil erosion (Pierson et al 2015, Dukes et al 2018, Sankey et al 2021. However, herbaceous vegetation recovers over time, reduces bare ground, improves water infiltration, and decreases water runoff, soil erosion, and sedimentation (Williams et al 2019). The vegetation structure change, therefore, leads to shifts in hydrologic and sediment connectivity across spatial scales and improves ecohydrologic function (Peters et al 2020, Williams et al 2020). Due to the spatial extent of woody-encroached areas, which span hundreds of millions of hectares in USA drylands as well as all other major drylands worldwide (Knapp et al 2008, Eldridge et al 2011, the control treatment impacts have large implications for vegetation, water, sediment, and carbon budgets at the national scale (Puttock et al 2014).
Using high resolution unmanned aerial vehicle (UAV) multispectral images and photogrammetric data, we quantify and compare woody vegetation cover and aboveground carbon in treated and control watersheds known as Natural Drainages in central Arizona, USA. We then combine the contemporary UAV observations of woody and herbaceous cover with historically documented trends in water and sediment yield before and after the treatment. While most ecohydrological studies document short-term impacts, the carefully-controlled historic experiments performed 63 years ago at the Natural Drainages provide a unique opportunity to compare the long-term impacts of the competing management priorities. Our detailed UAV-derived estimates of woody and herbaceous cover compliment the historical field-based records collected throughout the 20th century (figure 1) given the comparable level of details in both estimates of vegetation cover, which subsequently impacts water and sediment yield.

Study setting and historical water and sediment records
The Natural Drainages are four adjacent watersheds in the chaparral vegetation cover type, which covers approximately one half of the National Forest lands (∼4 million ha) in Arizona. Chaparral cover type in Arizona is historically dominated by turbinella oak (Quercus turbinella), manzanita (Arctostaphylos spp), desert ceanothus (Ceanothus greggii), and sumac (Rhus spp), while perennial grass cover co-occurs at varying abundance, typically on elevations between 1300-1800 m. Much of the contemporary chaparral cover type also includes juniper (Juniperus monosperma) and pinyon (Pinus monophyla). Ranging in size from 3.7 to 7.9 ha, the Natural Drainages experimental sites were first established in 1934 within the century-old Sierra Ancha Experimental Forest (SAEF) to assess livestock grazing effects on vegetation, streamflow, and sediment (figure 1). Known as the drainages A, B, C, and D, the Natural Drainages are individual watersheds contained within well-defined ridges that separate each adjacent watershed. The Natural Drainages all occur on easterly aspects of 15%-25% slope on 1380-1515 m elevation (Ingebo and Hibbert 1974).
Streamflow, runoff, and sediment have been measured in Natural Drainages since 1935 using 90 • V-notch weirs and stage recorders in each watershed (figure 1). Perennial herbaceous cover has also been measured since 1935 in 1 m 2 quadrats, while shrubs have been monitored in 5 m 2 quadrats (Rich and Reynolds 1963). All four drainages were closed to livestock grazing between 1934 and 1938, but 80% and 40% of the grasses were grazed on drainages A and D, respectively, between 1939 and 1954. No livestock grazing has occurred in any of the drainages since 1954.
In 1954, the US Forest Service treated Natural Drainages with herbicide application to control woody encroachment and to increase water yield from streamflow owing to surface water runoff (Rich andReynolds 1963, Ingebo and. Two of the drainages, A and C, were sprayed with 6.6% 2,4-D and 2,4,5 T in diesel oil during the summer 1954, while drainages B and D were designated as control drainages (Pond 1964). In the treated drainages A and C, the basal 13 cm of all shrubs and trees were sprayed until all woody species were totally killed. No other landscape-scale disturbance or treatments have occurred since 1954 in any of the drainages.
Prior to the treatment, crown and basal intercepts of all woody species were measured in all four drainages along ten 12 m line transects in 1950 and 1954. Pre-treatment total vegetation cover ranged 32%-40% (Rich and Reynolds 1963). Among them, shrub cover varied 19%-21% (Cable 1957) (figure 1). The most abundant pre-treatment shrub species was turbinella oak at 15%-22% cover across the four drainages. Other common shrub species included two sumac species (Rhus trilobata; Rhus ovata) and desert ceanothus (Ceanothus greggii) (Pond 1964). Following the treatment, the same measurements were made in 1957and 1959(Pond 1964. Historical studies and post-treatment records, which began in 1954 and terminated in 1971 (Ingebo and Hibbert 1974), demonstrate the following ecohydrological impacts of the herbicide treatment (figure 1): (a) Treated drainages had three times more herbaceous and sub-shrub (Eriogonum wrightii; Menodora scabra; Lotus rigidus; Lotus wrightii) cover and production compared to the control drainages on all slopes and quartzite-derived soils, a common parent material in Natural Drainages.
(b) No significant differences were observed in herbaceous and sub-shrub cover and production on the diabase-derived soils, the other common parent material in Natural Drainages. (c) The herbaceous cover increase resulted in 72% decrease in annual sedimentation. (d) The two treated drainages showed a 22% increase in streamflow compared to pre-treatment baseline data. (e) The streamflow increases were attributed to the lower evapotranspiration demands by the herbaceous cover that replaced the woody vegetation.

Contemporary estimates of woody and herbaceous vegetation cover
We flew a Sensefly eBee Ag fixed-wing UAV platform (Sensefly, Lausanne, Switzerland) (figure 1) with a multispectral sensor (multiSPEC4C) in four spectral bands (spectral range center): green (550 nm), red (660 nm), red edge (735 nm), and near-infrared (790 nm). We completed three flights on 29 July 2017 at 80 m altitude above ground, which resulted in image pixel resolution of 13 cm and covered a total area of ∼30 ha (figure 2). The multispectral image was processed with a photogrammetric method known as Structure-from-Motion (  that the final orthomosaicked image from the fixedwing UAV had a root mean squared error (RMSE) of 1.8, 1.6, and 2.9 m in the X, Y, and Z dimensions, respectively. However, our comparison of the UAV image with ground-based GPS locations (<30 cm post-processing accuracy) of calibration targets indicated <1 m accuracy in the X and Y dimensions. We classified the ground versus vegetation returns in the 3D point cloud data and then used the vegetation point cloud in ENVI 5.2 software and BCAL lidar module (Harris Geospatial, Boulder, CO, USA) to calculate the maximum vegetation height as well as the vegetation height range in 13 cm raster cells, consistent with the pixel size of the multispectral data.
These vegetation height raster bands were then combined with the original four spectral bands of the multispectral image to take advantage of the woody vegetation height information in classifying the rangeland cover types (Sankey et al 2019, 2021). In addition, we calculated the normalized difference vegetation index (NDVI) using the red edge and near infrared bands and calculated the mean texture of the NDVI values as a co-occurrence measure (Haralick et al 1973). The NDVI band and its mean texture band were also added to the multispectral image to generate a final eight-band image composite.
A classification and regression tree (CART) model was used in ENVI 5.2 software to classify the final eight-band multispectral image. Our preliminary analysis indicated that the CART model produced much more accurate results than the supervised classification models in ENVI, which used the available bands weighted equally as inputs. In contrast, the CART model selects input bands weighted by their importance as input variables. The final CART model included the maximum vegetation height, green, red, red edge, and NIR bands, the NDVI band and its mean texture as predictor variables. Our training data consisted of 400-1100 corresponding image pixels at locations where field-based GPS samples were mapped for the dominant species depending on their relative abundance. The differentiallycorrected Trimble GeoXH GPS locations had <30 cm accuracies.
We classified the following dominant target cover types: pinyon-juniper (Pinus monophyla-Juniperus monosperma), manzanita (Arctostaphylos pungens), turbinella oak (Quercus turbinella), bare ground, herbaceous and other vegetation. Given the fine spatial resolution of the multispectral image, shadows associated with each large woody plant was also classified as a separate class 'shadow' . Wright buckwheat (Eriogonum spp), hollyleaf buckhorn (Rhamnus spp), yucca (Yucca spp), desert ceanothus (Ceanothus spp), and sacahuista (Nolina spp) were other less common, small sub-shrubs observed in the field. We combined these less common sub-shrubs into a single cover type 'other vegetation' . Due to their spectral and fieldobserved similarities, pinyon versus juniper could not be separately classified and were combined into a single class pinyon-juniper.
We also generated a binary classification of total woody cover versus non-woody cover. The total woody cover included pinyon-juniper, manzanita, and turbinella oak, whereas the non-woody cover class included the other vegetation, herbaceous vegetation, and bare ground classes. We then estimated total woody cover and non-woody cover in percent within 5 × 5 m cells (∼1275 pixels within each cell) in all four drainages. By counting and summarizing the total number of pixels classified as woody cover within each 5 m cell (∼3000 cells per drainage), we estimated the woody and non-woody percent cover. The mean and standard deviation of the woody cover were then calculated for each drainage and compared using a one-way analysis of variance test with Tukey's pairwise multiple comparisons in R software.

Ground-based validation data for woody cover estimates
In spring and fall 2018, we identified and mapped with a Trimble GeoXH GPS all shrub and subshrub species in the historical 5 m 2 plots (n = 24 plots) distributed across the Natural Drainages. We assessed the UAV image classification accuracies using 60-220 pixel samples per class at the corresponding locations (N = 1096). More abundant species were represented by greater numbers of samples in the accuracy assessment. We also measured the location, aboveground height, and canopy diameter of 30 individual shrubs randomly distributed throughout the Natural Drainages. Using these field-based measurements, we evaluated the accuracy in the UAV SfM-derived shrub height estimates, which were an important input for the shrub species classification and aboveground shrub biomass estimates.

Contemporary estimates of aboveground carbon
The UAV-based dominant species map delineated each of our target species at the individual crown level, whereas the UAV SfM-derived canopy height model provided estimates for each shrub. Using these outputs, we estimated total aboveground biomass (AGB) within each drainage and converted the total biomass estimates to aboveground carbon estimates. To estimate AGB of the dominant woody species, we used the following allometric equations (equations (1)- (3)) established for juniper, manzanita, and turbinella oak, respectively: where CA is individual juniper tree crown area, which has a correlation coefficient of 0.92 with AGB for juniper with a RMSE of 41.23 kg (Cunliffe et al 2020); where CA is individual manzanita shrub crown area, which has a strong correlation with AGB, but not with shrub height, and RMSE of 0.67 kg (Huff et al 2017); where CH is the individual turbinella oak canopy height and BR is the basal radius, which predict AGB with a correlation coefficient of 0.95 and a RMSE of 0.16 Log(g) (McCraw 1985). We used the mean BR value of 1.125 cm and standard deviation of 0.66 cm, which we estimated from data reported by McCraw (1985), since our remote sensing-based measurements did not include basal radius. AGB was estimated for each individual plant and then summed to determine a total AGB per species in metric tons per hectare (Mg ha −1 ) in each of the control and treated drainages. Uncertainty of AGB was estimated following the approach of Sankey et al (2013). In short, 1000 realizations of the AGB model solutions for each individual plant of each species were run in which the predictor variables (e.g. CA, CH, BR) were allowed to vary from the observed value based on a specified distribution and standard deviation. The standard deviations were Individual tree/shrub biomass was then converted to carbon (C) by multiplying the AGB by a factor of 0.5 following Strand et al (2008), Sankey et al (2013) and Cunliffe et al (2020), who estimated juniper C content using AGB. The resulting C storage map was then overlaid with the dominant species map to estimate total aboveground C storage within each drainage and to compare treated versus control drainages.

Results
The UAV image was successfully classified with an overall accuracy of 93% (table 1). The dominant woody species were all classified with high producer's accuracies, especially manzanita. When we compared the SfM-derived individual shrub height estimates to the field-based measurements, the coefficient of determination was 0.72 with a RMSE of 0.76 m (figure 3). Given this correlation, we used the SfMderived shrub heights in our UAV image classification and AGB estimates, which integrated both the individual crown area and height.
Historical studies indicate that: (a) pre-treatment total shrub cover was 19%-21% among all four drainages (Cable 1957), (b) the 1954 chemical control treatment reduced shrub cover to 0% in the two treated drainages, and (c) consequently, posttreatment historical total shrub cover was 0% and ∼20% on treated and control drainages, respectively. In comparison, our contemporary UAV-derived total woody cover was ∼11% and ∼25% in treated and control drainages, respectively (figure 4), indicating 11% and 5% increases in woody cover over the 63 year period, respectively. Total non-woody cover in the UAV image was on average 89% and 75% in treated and control drainages, respectively, with 20% and 16% herbaceous cover.
The dominant woody species in the UAV image showed varying levels of abundance. Among them, manzanita and turbinella oak covered the largest areas at 22.5% and 6% of the total area, respectively, in the Natural Drainages (figure 4). In comparison, historical studies suggest that the most abundant pretreatment shrub species was turbinella oak. It averaged 41.5% of the total woody cover (Cable 1957) and 8% of the total area and has had a relatively stable abundance over time. The historical studies also show sporadic manzanita, pinyon, and juniper were recorded (Ingebo and Hibbert 1974), which suggests that these are newly arrived woody species within the study area that have greatly expanded over the last 63 years (figure 3). When we compared the UAV-derived total woody cover among the four drainages, the treated drainages A and C had significantly lower total woody cover compared to the control drainages (p < 0.001) (figure 5). Specifically, the control drainages B and D had twice as much total woody cover compared to the treated drainages 63 years after the treatment.
The total AGB of juniper in the treated and control drainages were 5.4 and 9.1 Mg ha −1 , respectively (figure 6). The juniper biomass estimates translate to 2.7 and 4.5 Mg ha −1 in total C in treated and control drainages, respectively. The total AGB of turbinella oak was 0.02 and 0.04 Mg ha −1 in treated and control drainages, respectively, which is equivalent to total C of 0.01 and 0.02 Mg ha −1 , respectively. The total manzanita AGB was 1 and 1 Mg ha −1 in treated and control drainages, respectively (figure 6), which translate to 0.5 and 0.5 Mg ha −1 of total C, respectively. Taken together, the control drainages had accumulated much larger C at 5.1 Mg ha −1 compared to the treated drainages, which included a total C of 3.2 Mg ha −1 (figure 6). This difference is largely associated with juniper biomass difference, since the other woody species are smaller in stature and similar in distribution between the control and treated drainages, while juniper AGB comprises a large component of the total C estimates.

Woody control impacts on vegetation and aboveground carbon
This study links contemporary UAV data with historical records to create a six-decades-long monitoring dataset, which reveals three important findings. First, our results demonstrate that the herbicide treatment from 1954 has had a land-use legacy impact on the vegetation. The herbicide-treated drainages today still have significantly lower total woody cover compared to the control drainages (figure 5). Although the herbicide was only applied to turbinella oak and other woody species that were dominant 63 years ago, the impact of the treatment appears also notable on currently dominant woody species including the newly arrived pinyon pine and juniper trees. While many woodlands undergo burning and thinning treatments at large scales, they are not often continuously controlled and monitored and only a few other studies have documented legacy effects (Browning and Archer 2011). The Natural Drainages in SAEF is unique in that no other disturbance events and treatments have occurred since 1954 and provide a rare opportunity to examine the longterm impacts from the treatment. Dryland management must consider legacy land uses in making decisions regarding community structure, ecological productivity, and ecosystem services (Browning et al 2015).
Secondly, woody cover in the control drainages has also increased compared to the historical estimates from prior to the herbicide treatment. This estimate indicates woody encroachment rates in no-action management scenario and is consistent with woody encroachment rates documented in chaparral ecosystems (Gibson et al 2019) and other rangelands across the western US (Romme et al 2009). Taken together, the historical records (Pond 1964(Pond , 1971) and the UAV image-derived estimates indicate that turbinella oak has continuously had a strong and relatively stable presence at 6%-8% cover in this ecosystem for over a century, but manzanita and pinyon-juniper are newly dominating species, which are now more common and as abundant as turbinella oak across the study area. We also note that the soil parent material appears to play a major role in these observed patterns in woody species distribution (figure 4). Turbinella oak, a deep-rooted sprouter that spreads clonally, is currently distributed across much of the diabase soil parent material and fractured bedrock within Natural Drainages, whereas the seeder species of manzanita and pinyon-juniper are found on the quartzitederived soil parent material (figure 4). This is consistent with previous observations in our study region (Leonard et al 2015, Gibson et al 2019).
Third, the control drainages B and D have approximately twice as much AGB in comparison to the treated drainages. In particular, the long-living species of juniper and turbinella oak, which can live >200 years, demonstrate the substantial impacts from the treatment, while shorter-living manzanita demonstrate similar AGB in treated versus control drainages. Land managers interested in herbicide treatment should consider the long-lasting C footprints at the decadal scale. Furthermore, regional-and continental-scale US carbon estimates continue to have large uncertainties due to the incomplete C stock estimates in non-forested and dryland ecosystems (Fusco et al 2019). Localscale C stock estimates such as ours can help complete regional-scale estimates in dryland ecosystems. Although many allometric relationships have been developed for the dominant species, some of the current allometric relationships need to be further refined to improve our estimates at the species and individual canopy levels. For example, the currently available turbinella oak allometric model  . Total AGB (Mg ha −1 ) for each dominant woody species as well as their sum at the drainage-scale in treated and control drainages. AGB was first estimated for each individual plant per species and then summed to determine the total AGB in the control and treated drainages. We estimated the uncertainty in the total AGB estimate using simulations and the uncertainty is shown with the error bars.
was based on a small dataset and was very sensitive to small changes in height and basal radius over a small range of both variables, which was evidenced by the small AGB values with large uncertainty (standard deviation) in our estimates for that species ( figure 6). Furthermore, our C stock estimates do not include potential changes in aboveground herbaceous biomass and associated belowground root biomass and soil carbon, while previous studies indicate that belowground carbon pools can be much larger than aboveground carbon stock (Fusco et al 2019). Ground-based plant AGB estimates are often time-consuming and labor-intensive. We demonstrate that UAV data provide a potential alternative means to estimate C at a landscape scale (e.g. 35 ha) and monitor woody encroachment similarly to previous lidar-based studies (Sankey et al 2013) and photogrammetry studies (Cunliffe et al 2020). While airborne lidar data provide similar accuracies (Sankey et al 2017, the UAV SfM data provide a cheaper alternative to acquire 3D data over larger areas (Shin et al 2018, Sankey et al 2019. The UAVderived individual shrub crown and height estimates (figure 3) can be particularly important for the SAEF, where individual shrub canopies were permanently marked in the 1920s and have been monitored over the last century. In UAV data applications for such estimates, the accuracies and errors should be carefully considered. While the initial UAV data error reported by the Pix4D software indicated errors of up to 2.9 m in the X, Y, and Z dimensions, our comparison of the UAV data with ground-based measurements and GPS data indicated errors of <1 m in the X and Y dimensions and 0.76 m in the Z dimension. Some of the errors might be attributed to the fact that our UAV images and field-based GPS data were collected several months apart, although this time difference might have negligible impacts on the long-living species in the chaparral ecosystem. Such errors and differences have important implications for estimating and monitoring short stature woody vegetation as well as small variability in topography.

Woody control impacts on ecohydrology
The legacy impacts on woody cover and AGB have important consequences for water and sediment yield (figure 7). First, woody encroachment increases plant water use and evapotranspiration (Qiao et al 2017). In the control drainages A and C, where woody cover and AGB are twice as high compared to the treated drainages, plant water use and evapotranspiration likely continued to be greater over the last six decades. The increases in plant water use and evapotranspiration significantly reduce streamflow and groundwater recharge in woody-encroached rangelands (Wilcox et al 2008, Kormos et al 2017. The herbicide treatment at our study area altered these trends and increased water yield by as much as 22% (Ingebo and Hibbert 1974), which was also consistent with other studies in our region (Ffolliot andThorud 1974, Hibbert et al 1974). A long-term study conducted nearby, which included initial impacts from wildfire followed by periodic herbicide treatments to suppress regrowth of shrub species, demonstrated initial average increases in water yield of 300%-700% in the first few years after treatment and average of 36% increase over the 18 years following treatment (Hibbert 1982). Although we lack recent quantitative estimates of water yield, we infer that the treated watersheds within the Natural Drainages have likely have benefitted water yield for at least part of the last six decades because our UAV-based estimates indicate that: (a) woody cover and AGB in the treated drainages are still half of those at the control drainages, and (b) total herbaceous and bare ground cover in treated drainages are still greater than those in the control drainages. This trend is overall consistent with the historical records (Pond 1964, Indego and, which demonstrated that the treated drainages had lower woody cover and greater herbaceous cover and thus greater streamflow. Secondly, the herbicide treatment has likely had sustained impact on sediment yield over time. The historical records indicate that annual sediment yield decreased on average by 72% following treatment (Ingebo and Hibbert 1974). Post-treatment herbaceous cover reduces undercanopy bare ground and sediment runoff (Williams et al 2019). We infer that the treated drainages at our study site have likely continued to have lower sediment yield because: (a) our UAV-based estimates indicate that the herbaceous cover in the treated drainages continued to be high and greater than that in the control drainages, and (b) this trend is consistent with the posttreatment historical observations. In contrast, sediment yield in the control drainages likely increased over the six decades because our UAV images indicate that woody cover increased in control drainages compared to the historical vegetation cover estimates. Over time, woody encroachment results in increased bare ground under and between tree canopies and subsequent sediment runoff (Pierson et al 2015). The increasing woody cover thus increases sediment connectivity and deteriorates ecohydrologic function (Turnbull et al 2012).
Finally, the long-term impacts from woody control treatment presented here can be useful for land managers, since most observational studies are limited to shorter time periods (figure 7). Rangeland managers around the world are facing a unique challenge because woody encroachment has been shown to deteriorate rangeland health, biodiversity, and ecohydrologic resources, but it also provides opportunities for increased C accumulation (Strand et al 2008, Sankey et al 2013. Policy makers tasked to manage for increased productivity and C accumulation can expect to double the C storage over several decades by not chemically or mechanically treating woody encroachment. In contrast, our contemporary and historical records demonstrate that given the right circumstances, managers can expect greater water yield and herbaceous forage in the decades following a landscape-scale woody control treatment (Hibbert 1983). To simultaneously create both impacts, land managers might consider a mosaic of controlled and treated landscapes, where woody species are treated in favorable locations for increased water yield and forage production, while other woody patches are left to accumulate C.

Conclusion
We re-visited historical records on woody and herbaceous cover and subsequent control treatment impacts in central AZ. Combined with the contemporary UAV images, these records illustrate that woody encroachment into rangelands left uncontrolled has accumulated large aboveground C storage at decadal scales. In contrast, herbicide control has had lasting impact on woody cover and provided some ecohydrological benefits. Land managers and policy makers can compare and contrast such benefits with risks to carefully balance impacts to C and ecohydrology associated with woody encroachment. Historical and contemporary observational datasets collected at similar scales and spatial resolutions can be effectively linked to provide site-specific records to inform decision making.

Data availability statement
The data generated and/or analysed during the current study are not publicly available for legal/ethical reasons but are available from the corresponding author on reasonable request.