Livestock removal increases plant cover across a heterogeneous dryland landscape on the Colorado Plateau

Livestock removal is increasingly used as a management option to mitigate the negative impacts of grazing-related disturbances on rangelands. Removal generally increases plant cover, but it is unclear when, where, and by how much plant and soil cover changes can be expected. On the Colorado Plateau, complex geology, topography, soils, and climate all interact to mediate the relationship between land cover, climate, and disturbance. In this study we used new developments in land cover mapping and analysis to assess landscape plant and bare soil cover up to 30 years after livestock removal from two grazing allotments in Capitol Reef National Park, Utah, USA. Results indicate that livestock removal increases plant cover 0.17%–0.32% per year and reduces bare soil cover 0.34%–0.41% per year, although these rates may be suppressed by warming temperatures. Soils, assessed through soil geomorphic units, played a strong but complex role in mediating land cover changes through time. These results suggest that livestock removal is an effective strategy for increasing plant cover and reducing bare soil on the Colorado Plateau, but including soil information in decision making will enhance efficiency by improving manager’s ability to prioritize management actions effectively across space and through time.


Introduction
Livestock grazing is a dominant land use in semi-arid and arid landscapes ('drylands' hereafter; Maestre et al 2016) where dry conditions preclude other uses such as crop cultivation (Thornton et al 2009, Hoover et al 2020. Grazing, like other ecological disturbances, can have significant effects on plants and soils depending on timing, intensity, and duration as well as the ecological context in which grazing occurs (Yates et al 2000, Sasaki et al 2007, Eldridge et al 2017, Staver et al 2021, Maestre et al 2022. Drylands can have especially low resilience to grazing impacts due to low growth rates, high climactic variability, and sensitive soils (Fleischner 1994, Sloat et al 2018, Souther et al 2020. Livestock removal of varying durations is increasingly used to ameliorate deleterious impacts, but previous research has found instances where removal of livestock alone does not lead to ecosystem recovery due to persistent ecosystem state change (Schlesinger et al 1990. Understanding where and when changes to livestock management may lead to desired shifts in ecosystem condition is a central priority for global land conservation (Rietkerk and van de Koppel 1997, Bestelmeyer et al 2004, 2015, Steinfeld et al 2006, Miller et al 2011, Willemen et al 2020, Durigan et al 2022. Mitigating grazing-induced land degradation requires understanding landscape context such as climate, soils, topography, management history, and their interactions (Milchunas and Lauenroth 1993, Souther et al 2020, Maestre et al 2022. In the Colorado Plateau and other drylands of the western United States, widespread and often heavy cattle and sheep grazing began in the 19th century, with changes in land use priorities and policy beginning in the mid-to late-20th century that increased oversight and management (Copeland et al 2017, Sayre 2017). Throughout this time, the landscapes of the Colorado Plateau have presented significant challenges to grazing and its management due to restrictive terrain, high intra-and inter-annual variability in precipitation and temperature, and complex soil physical properties (Hereford et al 2002. An ongoing regional 20 year drought as well as changes in land use priorities have driven significant grazing reduction and livestock removal, both of which are likely to continue under forecasted rapid regional warming (Williams et al 2022). While research in other systems has examined the relationship between plant cover and grazing in the context of landscape setting, few studies have attempted to examine their joint effects on a landscape as heterogeneous as the Colorado Plateau.
Novel gridded data products and analytical frameworks provide new opportunities for learning from past changes in land management across large spatiotemporal extents (Fick et al 2020, Allred et al 2021, Simler-Williamson and Germino 2022. In particular, Nauman et al recently developed a technique to select study reference sites by matching study areas to comparable areas of the landscape based on similar topography and soils (Nauman andDuniway 2016, Nauman et al 2017). To extend this approach to time-series of remotely sensed data, Fick et al incorporated the 'synthetic control' method common in econometric literature (Fick et al 2020. This addition constructs a 'counterfactual estimate (CE)'-i.e. a modeled covariate that estimates the value of landscape response variables in the absence of management intervention. The combination of these new tools and datasets creates powerful approaches to understand landscape management effectiveness through time. It nonetheless remains unclear how these approaches perform in studies of the often more subtle, widespread, and gradual changes associated with domestic livestock grazing on the complex drylands of the southwestern US.
The goal of this study is to assess land cover change following the permanent removal of domestic livestock from grazing allotments in Capitol Reef National Park on the Colorado Plateau. We seek to test whether the direction and magnitude of land cover changes depends on soil-geomorphic setting and how this context interacts with past grazing intensity. We also seek to evaluate the utility of the synthetic control methodology for studies of landscape change by evaluating our results both with and without a modeled CE of management intervention. To accomplish this, we fit a suite of multilevel Bayesian linear models built with gridded climate data and soil geomorphic unit (SGU) maps. Past grazing intensity is accounted for using additional Digital Eelevation Model-derived explanatory variables, while additional model covariates attempt to minimize the influence of spatial and temporal autocorrelation within our landscape response variables.

Study system
We focus our study on Capitol Reef National Park in south-central Utah, USA. Federally managed rangelands in and adjacent to Capitol Reef National Park consist of dry mesas, escarpments, and valleys and are representative of many dryland ecosystems on the Colorado Plateau . Here, livestock grazing historically occurred from approximately October to May in large valleys managed through an allotment-permit system. The establishment of Capitol Reef National Monument in 1937, its expansion in 1969, and conversion to National Park in 1971 brought policy changes which nearly eliminated livestock grazing in the park through permit expiration and buy-outs (appendix A; Copeland et al 2017, NPS 2018. Livestock were removed gradually through time, with final livestock removal occurring on two allotments ('Cathedral' and'Rock Springs') in 1999 and1989, respectively. These 'retirement' years ensured that several years of pre-and post-retirement data was available from the current satellite record . Both allotments are spatially extensive (4293 ha and 1354 ha respectively) and capture significant landscape heterogeneity in soils and topography (figure 1; Nauman et al 2019, Nauman 2021).

Land cover responses
Land cover was assessed using three remotely sensed 30 m resolution measures of plant functional group and soil cover for 1984-2020: mean annual percent of soil with no plant, rock, litter, or biological soil crust cover ('bare soil cover'); mean annual percent foliar cover of herbaceous perennial plants ('perennial herbaceous cover'); and mean percent foliar cover of all plants during the months of March, April, and May ('total spring cover'). Bare soil and perennial herbaceous cover were sourced from the Rangeland Analysis Platform version 3 in November 2021 (Allred et al 2021). LandCART 4.0 was accessed in November 2021 to calculate spring total foliar cover as the mean of AIM 2.0 indicator 'Total Foliar' for March, April, and May (Zhou et al 2020, Okin et al 2021, McCord et al 2022. Validation of Rangeland Analysis Platform data across our study area using plot data from the Northern Colorado Plateau Monitoring Network showed good (R 2 = 0.77) and modest (R 2 = 0.25) performance for bare soil and perennial herbaceous cover, respectively (n = 493; Witwicki et al 2017). Park, Utah, USA. Pixels within the allotment boundaries made up the initial pool of 'treatment' pixels used for counterfactual and inference analysis, while pixels outside the allotment made up part of the initial pool of 'reference' pixels used in counterfactual analysis. Greyed areas are part of adjacent allotments that were excluded from the study. Topography is displayed as a 45 • hillshade generated from USGS 3DEP 1/3 • DEM. Additional details on counterfactual analysis are available in Appendix C.

Explanatory variables of interest
Variables describing climate, edaphic, and livestock management factors were constructed for explanatory model building. Climate variables sourced from Daymet version 4 monthly climate summaries included the sum of monthly precipitation ('PPT'), monthly mean maximum temperature ('T max ' , • C), and monthly mean minimum temperature ('T min ' , • C; Thornton et al 2020). Because the seasonality of water availability on the Colorado Plateau does not align with the calendar year, climate variables for each study year were summarized over the 'water year' beginning October 1 of the prior calendar year and ending September 30 of the current calendar year (Hereford et al 2002). Variance in soil factors was accounted for using SGU maps (Nauman 2021). Briefly, Nauman et al initially formulated SGUs by grouping ecological site descriptions according to qualitative similarity (Caudle et al 2013. They were then refined using statistical processes to maximize their ability to explain variability in landscape production and mapped across the Upper Colorado River Basin using Random Forest classification . These soil maps capture substantial information on soil texture, topographic setting, salinity, depth, and available moisture and summarize them into groups that are amenable to further statistical analysis and interpretation.
In accordance with previous literature, we used cost-distance to water ('CDW') as a surrogate for past grazing intensity with a terrain ruggedness index as the cost surface and surface water sources as nodes (appendix B; Bailey et al 1996, 2015, Tarboton 1997, Riley et al 1999, Pringle and Landsberg 2004, Fensham and Fairfax 2008, ESRI 2019. Surface water sources used by livestock in these allotments are reliable creeks and earthen 'stock tanks' that capture surface runoff (appendix B). Areas not accessible to livestock because of terrain (ruggedness index >2.3) were excluded by assigning them an infinite cost for crossing. Ruggedness on pixels located on roads was reduced by half. To reflect uncertainty in the relationship between modeled CDW and actual grazing intensity, we normalized continuous cost-distance values as Z-scores (expressed in units of standard deviations from the within-allotment population mean, σ) and binned them into relative cost-distance categories (Raynor et al 2021): 'Very High' (⩾1σ), 'High' (<1σ and ⩾0σ), 'Low' (<0σ and ⩾−1σ), and 'Very Low' (<−1σ). We expect highest grazing intensity to occur in the 'Very Low' CDW category due to low metabolic cost for livestock movement to those areas, and lowest grazing intensity to occur in the 'Very High' cost-distance category due to high metabolic cost for livestock movement to those areas (Martina et al 2015). Water available to livestock as snow was accounted for using mean snowpack cover from November through March as a covariate sourced from Landsat Level-3 Fractional Snow Covered Area (fSCA) maps (Degen and Young 1984, Selkowitz et al 2017.

Counterfactual analysis
Counterfactual analysis conducted here built on and combined the methods of previous research (Nauman and Duniway 2016, Nauman et al 2017, Fick et al 2020Fick et al , 2022. Specifically, post-intervention estimates Table 1. Sample size and performance evaluation results for 24 multilevel Bayesian models fit to three land cover response variables on two retired grazing allotments in Capitol Reef National Park, USA. Relative model fit is quantified as log pointwise predictive density (∆ELPD), where the best-fit model is assigned a value of zero and worse-fit models have progressively more negative values. of each land cover response variable in the absence of management intervention for each study ('treatment') pixel were predicted by time series models fit to pre-intervention data. Predictions were based on observed land cover response in pools of pixels taken from currently grazed ('reference') areas outside Capitol Reef National Park. While this notation is potentially confusing, we emphasize that it is conventional within the synthetic control literature and use it here to maintain consistency with previous research . Interventions were the removal of livestock. Our 'CE' predicted landcover through time for each pixel under the assumption that livestock grazing had continued within the retired allotments as they do on pixels in our reference pool. Similar to our operational use of 'treatment' and 'reference' , the use of 'intervention' here is conventional within the synthetic control literature. Additional details are available in appendix C.

Landscape response to livestock removal
We modeled the effect of livestock removal on landscape cover as the number of years since counterfactual intervention ('YSI'). In Cathedral, the year of management intervention (YSI = 0) was defined as the first year without grazing (1999; table 1). In Rock Springs, we set the year of management intervention as the second year without grazing (1990) for analysis of bare soil and perennial herbaceous cover and the third year without grazing (1991) for total spring foliar cover. This is due to data requirements during synthetic control model training as well as limits on the temporal extent of our fractional cover products (Fick et al 2020, Allred et al 2021, Okin et al 2021. We emphasize that actual livestock removal across the landscape proceeded haphazardly in the years preceding official allotment retirement, and management records on this process are incomplete (appendix A). Nonetheless, our analysis here necessarily assumes that livestock removal was started and finished before the first year without grazing.

Model building
Data for each model was randomly under-sampled to ensure the coarser-resolution climate data was not duplicated across pixels. Furthermore, data was balanced by randomly under-sampling group strata so that a maximum of 30 samples of each unique combination of SGU, CDW, and YSI were included for analysis. Finally, PPT, T max , T min , and fSCA were centered and scaled to facilitate comparison of model coefficients (Schielzeth et al 2010). We utilize a multi-model approach to evaluate progressively more complex hypotheses while guarding against overfitting. Four multilevel Bayesian models were built for each landscape response variable in each allotment (Bürkner 2017(Bürkner , 2018: where 'y' is the land cover response; α 0 is the grand intercept; β 1 -β 4 are fixed slopes; 'PPT' is total precipitation; 'T max ' is mean annual maximum temperature; 'T min ' is mean annual minimum temperature; 'YSI' is time in years since counterfactual intervention; 'SGU' is soil geomorphic unit; 'CDW' is costdistance to water; 'fSCA' is fractional snow cover; α SGU , α CDW and α SGU:CDW are varying intercepts; β fSCA and β YSI are varying slopes. Additionally, all four models shared three covariate terms: where α CE is a varying intercept fit to our CE of cover in each pixel for each post-intervention study year, α year is a varying intercept fit to each postintervention study year to further account for interannual variability in climate, and ρ year|pixel is a first-order autocorrelation structure fit to study year grouped within pixel to account for temporal autocorrelation in land cover responses. Model (1) hypothesizes that land cover within grazing allotments is controlled by climate and soils alone and not affected by the magnitude of past grazing. Model (2) excludes soil factors and hypothesizes that only climate and the potential effects of past grazing through the interaction between mean snow cover and CDW explains land cover change. Model (3) hypothesizes that climate, soils and past grazing affect land cover and includes a varying intercept fit to the crossed levels of SGU and CDW to account for the differential sensitivity of some soil settings to past grazing. Model (4) serves as our full model, wherein we hypothesize that all previous terms are informative in addition to the interaction between time since intervention and SGU. Lastly, we chose to evaluate the CE separately because of the relative paucity of synthetic control studies in the ecology literature as well as the relative novelty of its present implementation. To accomplish this, we ran each model with and without α CE and hypothesized that its inclusion would improve model fit and confidence in parameter estimates by accounting for unmeasured sources of heterogeneity.
Model fitting utilized weakly-informative Gaussian priors and four Markov chains with 4000 total (1000 warmup) iterations (`brms`; Bürkner 2017Bürkner , 2018. Overall model fit was evaluated by calculating Bayes' correlation coefficient (R 2 ; Gelman et al 2019). Relative model fit was assessed by comparing expected log pointwise predictive density (∆ELPD) estimated using Pareto smoothed importancesampling leave-one-out cross-validation (PSIS-LOO; Piironen and Vehtari 2017. After selecting the model that best explains each landscape cover variable in our dataset, we consider fixed slopes, varying intercepts, and varying slopes significant if their estimated Bayesian posterior 90% credible interval did not include zero. All analyses utilized R version 4.1.1 (R Core Team 2022). All analysis and figure code is available upon request in R package format. Data generated during this study are available from the USGS ScienceBase-Catalog (McNellis et al 2023).

Regression model evaluation
All models converged adequately (R = 1.00 ± 0.01, supplementary results). Posterior distributions adequately matched observed data, with some bimodality in bare soil and total spring cover models and left skew in perennial herbaceous cover models. Alternative specifications that attempted to account for this skew and bimodality performed poorly. Overall fit for candidate models was modest (R 2 = 0.26-0.38, table 1), while our full model was the best fit to our data across allotments and response variables (table 1). Models without CE were not significantly different from models with CE with respect to overall performance. As such, we restrict further analysis and discussion to only our fullest model for each combination of our grazing allotments and response variables.
Model coefficient estimates (expressed in units of the land cover response per standard deviation of each explanatory variable ±90% Bayesian posterior credible interval) for gridded climate covariates were similar between allotments but differed based on response variable (figure 2). T max and T min were not significantly different from zero except for bare soil in Cathedral, where its estimated contribution was 0.91 (±0.30) and −1.01 (±0.27) for T max and T min , respectively. Precipitation was significant only for perennial herbaceous cover, where it had a coefficient of 0.33 (±0.11) in Cathedral and 0.17 (±0.10) in Rock Springs. YSI (expressed in units of land cover response per year) had a significant positive effect on total spring foliar cover (0.32 ± 0.07, 0.20 ± 0.07) and perennial herbaceous cover (0.17 ± 0.08, 0.21 ± 0.09) and a significant negative effect on bare soil (−0.41 ± 0.11, −0.34 ± 0.09) in Cathedral and Rock Springs, respectively (figure 2).

Landscape response to soil geomorphic setting
A total of 11 SGUs were analyzed across the study area, with 10 in Cathedral (no 'Finer Uplands') and 9 in Rock Springs (no 'Gypsum' or 'Saline   . Modeled interaction between cost-distance to water sources calculated based on terrain ruggedness (CDW) and water available to livestock through snow cover (fSCA). The Y-axis quantifies the change in each land cover response variable with each unit change in fractional snow cover, which is allowed to vary by cost-distance bin. Continuous CDW calculations are binned into relative CDW categories based on Z-scores calculated on the overall CDW allotment mean. High CDW implies high metabolic cost on livestock attempting to forage in that pixel, and therefore low expected grazing intensity. Low CDW implies low metabolic cost on livestock attempting to forage in that pixel, and therefore high expected grazing intensity (Martina et al 2015). Coefficients are considered significant if their 90% Bayesian credible interval (shown as error bars) does not include zero.
Springs. Perennial herbaceous cover was positively associated with 'Sandy Bottoms' in Cathedral and negatively associated with 'Gypsum' and 'Saline Hills' in Cathedral and Rock Springs, respectively. Total spring cover was positively associated with 'Deep Rocky' and 'Loamy Uplands' (both allotments) and 'Very Shallow' (Rock Springs) but negatively associated with 'Saline Hills' and 'Sandy Bottoms' (Rock Springs).

Impacts of past grazing intensity
Varying intercepts fit to CDW were significant for 3 out of 24 total combinations of allotment, land cover, and CDW (supplementary results). All three significant intercepts fit to CDW were in Cathedral: bare soil cover in 'Very High' (−14.59 ± 6.79), perennial herbaceous cover in 'Very Low' (2.41 ± 1.45), and total foliar cover in 'Very High' (9.18 ± 3.66). Varying slopes fit to fSCA for each level of CDW were significant for 20 out of 24 levels (figure 4). Slopes relating bare soil cover to fSCA were negative ([−2.12, −0.71]) while all slopes for total spring cover were positive ([0.33, 1.56]). Slopes for perennial herbaceous cover were near zero or not significantly positive ([−0.19, 0.13]).

Interactions between soils, years since intervention, and past grazing intensity
Varying slopes fit to quantify the interaction between SGU and YSI were significant for 35 out of the 57 available combinations of SGU, land cover, and allotment (figure 5; supplementary results). All SGU-YSI slopes, except 'Shallow' , were significant for at least one response variable in one allotment. 'Deep Rocky' was the only SGU that had a significant slope for all response variables in all allotments and was consistently negative for perennial herbaceous cover, negative for total spring cover, and positive for bare soil cover. 'Saline Hills' and 'Loamy Uplands' had significant slopes for all but total spring cover in Cathedral and perennial herbaceous cover in Rock Springs, respectively, but this was not consistent between allotments. 'Gypsum' was not present in Rock Springs but had a significant negative slope in all three cover responses in Cathedral. 'Sandy Uplands' and 'Very Shallow' , while present in Cathedral, only had significant slopes in Rock Springs. 'Finer Uplands' and 'Saline Bottoms' only had significant slopes for perennial herbaceous cover in Rock Springs and Cathedral, respectively.
Varying intercepts that estimated the interaction between SGU and CDW were significant for 41 out of 201 available combinations of SGU, CDW, landscape response, and allotment (figure 6). Significant interactions were generally evenly spread across responses and allotments but were concentrated in 'Loamy Uplands' and 'Deep Rocky' in Rock Springs (18/21 significant interactions) and 'Loamy Uplands' in Cathedral (8/20 significant interactions). These patterns were also apparent in the magnitude of significant interactions, although effects were somewhat stronger in bare soil cover and muted in perennial herbaceous cover.

Discussion
Livestock grazing in drylands is a complex ecosystem disturbance that interacts with climate, topography, and soils to influence the extent and composition of plant cover across landscapes (Maestre et al 2022, Milchunas et al 1988, Souther et al 2020, Raynor et al 2021. Here, we attempted to quantify the effects of livestock removal on plant and soil cover in a large, complex, and spatiotemporally-heterogenous landscape on the Colorado Plateau using new developments in land cover and soil mapping as well as spatial analysis. Overall, we found plant cover to be increasing and bare soil generally decreasing post-removal with significant interactions with temperature-but this was heavily mediated by differences in soil setting.
Our model results suggest that livestock removal on Colorado Plateau drylands resulted in a ca. 0.25% per year increase in plant cover and ca. 0.35% per year reduction in bare soil cover after accounting for significant climate and landscape heterogeneity (figure 2). Importantly, this change was detected despite significant interannual variability in climate and reduced precipitation due to the co-occurring regional 'megadrought' (Williams et al 2022). Our recovery rate is somewhat high (cf Duniway et al 2018; 0.14% yr −1 ) but is generally supported by previous work that has found that grazing in many dryland ecosystems reduces overall plant productivity and soil cover, which can be reversed with livestock removal (Lusby 1970, Valone and Sauter 2005, Yeo 2005, Rickart et al 2013, Witwicki 2020, Wolf and Mitchell 2021. Understanding landscape resilience to livestock grazing intensity is a central tenet of rangeland science. This has previously been understood by calculating the sum of the number of months each animal grazes on the landscape (Animal Use Months, Figure 6. Modeled response of landscape functional cover due to the interaction between soil classes (SGU; Nauman et al 2022) and cost-distance to water sources calculated based on terrain ruggedness. Interaction is quantified as a varying intercept (±1 S.E.) in multilevel Bayesian linear models. Cost-distance to water is summarized by binned Z-scores of a continuous cost-distance measure calculated within each retired allotment. High CDW implies high metabolic cost on livestock attempting to forage in that pixel, and therefore low expected grazing intensity. Low CDW implies low metabolic cost on livestock attempting to forage in that pixel, and therefore high expected grazing intensity (Martina et al 2015). Non-significant varying intercept terms are omitted for clarity (i.e. their Bayesian 90% credible interval included zero; supplementary results). 'AUM'; Smith et al 2017). While AUMs were measured sporadically in our study area (appendix A), our alternative measure using CDW provided key insight into the landscape response to livestock removal. Specifically, past grazing intensity terms improved model performance but trends across CDW bins were variable (figures 4 and 6). This aligns with previous literature that suggests grazing intensity, driven by topography as well as proximity to water and snow, is important for understanding grazing impacts (Fensham andFairfax 2008, Raynor et al 2021). Our results may be partly influenced by the confounding interaction of CDW and fSCA with soil water availabilitylow CDW areas tended to occur in productive run-in locations at the bottom of drainages where livestock tend to congregate. Nonetheless, our results suggest that the overall impact of grazing on a landscape may be underestimated if only areas close to water sources are considered in range management assessments. Improved measures of livestock movement across the landscape (e.g. GPS collars; Bailey et al 2021) may improve estimates of grazing intensity in future studies.
Climate had modest effects on land cover across our study area (figure 2; appendix D). Notably, warmer annual temperatures significantly increased bare soil in Cathedral, supporting other studies that show increases in temperatures on the Colorado Plateau are reducing plant cover and damaging biological soil crusts, even in the absence of livestock (Munson et al 2011, Finger-Higgens et al 2022, Phillips et al 2022. The Colorado Plateau is predicted to heat faster than any other non-polar region in North America, indicating that forage production and other ecosystem services may face significant threats under future climates (Seager and Vecchi 2010, Cook et al 2015, Munson et al 2015, Williams et al 2022. Our results found less evidence of an effect of precipitation except for a small positive effect on perennial herbaceous cover (figure 2). Prior research nonetheless suggests changes in the timing and magnitude of precipitation may have significant effects on dryland plant function, especially when compounded with increased air temperature (Schwinning et al 2008, Winkler et al 2019. While we did find evidence of post-grazing recovery, our results indicate that managers may need to carefully consider the joint effect of grazing and warming during decision making if reduction in bare soil cover is a desired intervention outcome. For example, drylands that are warming rapidly (such as the Colorado Plateau) that appear resilient now may experience a slower or even stalled recovery from grazing disturbance under a future warmer climate.
Soil systems of the Colorado Plateau are highly complex due to the interactive effects of geology, geomorphology, and pedogenesis . In this study, we accounted for this soil complexity using new developments in soil mapping and robust controls (Brungard et al 2021. Our results indicate that soils (and their interactions) are the most important factors determining post-grazing recovery on the Colorado Plateau. Importantly, the influence of soil factors on recovery from grazing depends strongly on the soils themselves. 'Shallow' soils, for example, showed no significant trends between livestock removal and land cover (figures 5 and 6; supplementary results). Part of this may be because shallow soils tend to have lower forage, higher tree cover, and exposed bedrock that may reduce their attractiveness to livestock . However, potential forage productivity did not consistently predict grazing response: productive grasslands ('Sandy Uplands') had a positive interaction between plant cover and time since livestock removal, while soils likely with less forage ('Gypsum' , 'Saline Hills') had a negative or mixed interaction (figure 5; Nauman et al 2022). The most robust interaction we found was in two productive upland soil types ('Deep Rocky' , 'Loamy Uplands') that negatively impacted plant cover (i.e. bare soil increased and plant cover decreased after livestock removal, figure 5). These were also the two soil types that interacted most strongly with past grazing intensity: CDW had a somewhat negative effect in 'Loamy Uplands' but a mixed effect in 'Deep Rocky' (figure 6). While soil type strongly influences landscape response to grazing, it is unclear what mechanisms may be driving these patterns. For example, ecological drought is complex and may be better assessed using cumulative measures of soil water deficit or measures that account for the seasonality of water inputs . Alternatively, analysis with improved maps of bare ground and perennial herbaceous cover as well as other landscape cover types (such as invasive grass and biological soil crust) may improve interpretation in particular soil settings (Weber et al 2008, Horning et al 2020. While we had modest to good agreement between modeled and measured bare ground (R 2 = 0.77) and perennial herbaceous cover (R 2 = 0.25) in our study area, we did not assess the relative performance of our remotely-sensed spring total foliar cover estimate or gridded climate variables. Lastly, use of these allotments by wildlife was not accounted for and likely increased following livestock removal (particularly for elk). Our results nonetheless point to the utility of soil mapping products, specifically those based on the concept of 'land potential' , for management or policy actions aimed at mitigating livestock impacts to rangelands. For example, livestock removal may be most impactful on 'Sandy Bottoms' and 'Saline Uplands' soil types in our study area.
Our analytical approach with the synthetic control estimates in this study diverges from past work, providing new insights into its applicability. Synthetic control econometric studies generally quantify ratios or differences between treatments and controls at the scale of US counties or even entire countries (e.g. Mao 2018, Schwarz et al 2022. We instead analyze individual pixels and used the CE as a covariate in a multilevel model, which generated new insights into the influence of CE. Including CE as a varying intercept allowed us to quantify its utility through model and parameter uncertainty, but it surprisingly did not affect model performance. This may in part be due to our inclusion of additional covariates (year and a 1st-order autocorrelation term, ρ). In particular, ρ was structured on pixel and intended to capture the strong dependence of current plant cover on prior plant cover. Monroe et al utilized a similar approach and found increased support for models that included a space-for-time control covariate, but with much smaller treatment polygons (Monroe et al 2022). Future landscape studies may find less utility in this approach than studies of smaller spatial extents, such as oil and gas well-pads or vegetation treatments (Nauman et al 2017. Livestock removal is also a subtle treatment effect and may be less amenable to synthetic control methodology, although studies using smaller livestock removal extents (e.g. exclosures) may benefit from the approach. We emphasize, however, that further validation is required to ensure the robustness of this approach in landscape studies. Further methods adapted from the synthetic control literature may also improve future studies (e.g. placebo studies; Abadie and Gardeazabal 2003). Nonetheless, the redundancy of the CE estimate here provides optimistic support for synthetic control methods in studies where robust estimation of spatial and temporal covariates may not be possible. Studies which seek to explain differences between treatment and reference areas of the landscape (rather than simply control for nuisance variance) may find significant utility in the approaches outlined here.

Conclusions
Livestock grazing is common in drylands globally, with complex impacts on landscapes. Here we leveraged new tools that can help model grazing impacts while accounting for significant variance in topography, climate, soils, and past grazing intensity. Our results suggest that livestock removal is effective at decreasing exposed bare soil and increasing plant cover, although the magnitude of response is strongly mediated by soils and geomorphology. While we found a smaller effect for climate than anticipated, key relationships with temperature suggest that increased regional temperature may limit or reverse trends due to removal of livestock. Our study suggests drylands on the Colorado Plateau appear to be currently responsive to livestock removal, but it may be a less effective management intervention under future climates.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https:// doi.org/10.5066/P9EK2PHY. Data will be available from 31 January 2023.
Additional details describing counterfactual analysis.
For each response variable within each retired grazing allotments, counterfactual analysis proceeded in 3 steps: (a) find 'treatment' and 'reference' pixels eligible for analysis, (b) build a pool of 'reference' pixels for each treatment pixel based on automated selection methods, and (c) build a time-series model for each treatment pixel and its reference pixel pool that estimates landscape responses in the absence of counterfactual intervention. In step 1, we first exclude pixels from both treatment and reference pools that were identified as roads in the U.S. Census TIGER dataset (Census Bureau 2021). We further excluded roads, livestock exclosures, and livestock trailing routes absent from public datasets that were identified using local land management expertise and internal records. Additionally, we excluded pixels that had been previously mapped as exposed bedrock with >60% confidence (Brungard et al 2021). Reference pixels were further restricted according to four criteria: pixels must occur within Bureau of Land Management grazing allotments exterior to Capitol Reef National Park, pixels must occur no more than 20 km from the borders of grazing allotments within Capitol Reef, pixels must have the same soil particle size class as the target pixel, and pixels must have salinity within 5% of the target pixel . In step 2, all reference pixels were ranked according to their Gower's similarity to each treatment pixel along 15 axes of topographic variability (see below) and select only the top 200 best-matching pixels for time series model building (table C1; Nauman andDuniway 2016, Fick et al 2022). Time series models utilized the synthetic control approach elaborated by Fick et al (2020) and predict each landscape response variable through time starting from the year of counterfactual intervention under the assumption that the intervention did not occur. In this study, interventions were the removal of livestock, and our counterfactual estimate ('CE') predicted landcover through time for each pixel under the assumption that livestock had continued to graze within the retired allotments as they do on pixels in our reference pool. We include CE as a covariate in explanatory model building to account for spatiotemporal variability in landscape cover that can be explained by topography, salinity, particle size class, and interannual stochasticity in climate.  Reef National Park, Utah, USA. Pixels within the allotment boundaries made up the initial pool of 'treatment' pixels used for counterfactual and inference analysis, while pixels outside the allotment made up the initial pool of 'reference' pixels used in counterfactual analysis. Greyed areas are (1) part of other allotments within Capitol Reef National Park that were excluded from the study, (2) excluded from SGU mapping based on the methods of Nauman and Duniway (2021), or (3) outside the 20-km radius defined as the search boundary for 'reference' pixels in counterfactual analysis. Topography is displayed as a 45 • hillshade generated from USGS 3DEP 1/3 • DEM.

Appendix D
Daymet v4 1 km gridded monthly climate data (1981-2020) across all 'treatment' pixels used in inference analysis of two retired livestock grazing allotments in Capitol Reef National Park, Utah, USA. Grey vertical bars indicate the year of livestock removal. Precipitation ('PPT') is summarized as annual sums, while maximum temperature ('Tmax') and minimum temperature ('Tmin') are summarized as annual means.