On Transient Semi‐Arid Ecosystem Dynamics Using Landlab: Vegetation Shifts, Topographic Refugia, and Response to Climate

Projecting how arid and semi‐arid ecosystems respond to global change requires the integration of a wide array of analytical and numerical models to address different aspects of complex ecosystems. We used the Landlab earth surface modeling toolkit (Hobley et al., 2017, https://doi.org/10.5194/esurf-5-21-2017) to couple several ecohydrologic and vegetation dynamics processes to investigate the controls of exogenous drivers (climate, topography, fires, and grazing) and endogenous grass‐fire feedback mechanisms. Aspect‐controlled ecosystems and historical woody plant encroachment (WPE) narratives in central New Mexico, USA are used to construct simulations. Modeled ecosystem response to climatic wetness (i.e., higher precipitation, lower potential evapotranspiration) on topography follows the Boyko's “geo‐ecological law of distribution.” Shrubs occupy cooler pole‐facing slopes in the dry end of their ecoclimatic range (Mean Annual Precipitation, MAP ≤ 200 mm), and shift toward warmer equator‐facing slopes as regional moisture increases (MAP > 250 mm). Trees begin to occupy pole‐facing slopes when MAP > 200 mm, and gradually move to valleys. Pole‐facing slopes increase species diversity at the landscape scale by hosting relict populations during dry periods. WPE observed in the region since the middle 1800s is predicted as a three‐phase phenomenon. Phase II, rapid expansion, requires the removal of the positive grass‐fire feedback by livestock grazing or fire suppression. Regime shifts from grassland to shrubland are marked by critical thresholds that involve grass cover remaining below 40%, shrub cover increasing to 10%–20% range, and the grass connectivity, Cg, remaining below 0.15. A critical transition to shrubland is predicted when grazing pressure is not removed before shrub cover attains 60%.

In the American Southwest, introduction of domestic grazers, fire suppression, and changes in rainfall variability are among the exogenous drivers that have accelerated the rate of shrub encroachment into grasslands since the European settlement in the 19th century (e.g., Archer et al., 2017;Van Auken, 2009). Shrub expansion is reinforced by both the loss of grass-fire feedback that has historically acted as a demographic bottleneck for expansion, and the initiation of shrub-microclimate and shrub-wind erosion feedbacks (Archer et al., 2017;D'Odorico et al., 2012). Invading shrubs increase nocturnal air temperatures through the shrub-microclimate feedback that benefits shrub survival in winters (He et al., , 2015. The shrub-wind erosion feedback promotes nutrient-rich soils under shrub canopies by trapping nutrients from aeolian sediments eroded from bare inter-shrub patches, gradually causing degraded desert states (Bowman et al., 2015;Okin et al., 2009;Turnbull et al., 2012;Wilcox et al., 2018).
As in the case of most ecosystem transitions (May, 1977;Touboul et al., 2018), grass-woodland transitions stemming from WPE can follow one of the three types: (a) continuous and reversible; (b) threshold-dependent and reversible; and (c) threshold-dependent and hysteretic (irreversible). An ecosystem transition is considered continuous and reversible when endogenous feedbacks are absent or weak enough such that exogenous drivers move the ecosystem state back and forth along similar trajectories; for example, fluctuations in vegetation cover in response to annual precipitation (Brandt et al., 2019;Moreno-de las Heras et al., 2015). Conversely, an ecosystem tends to exhibit irreversible changes when feedbacks (removed or initiated) are strong (Ratajczak, Nippert, Briggs, & Blair, 2014;Ratajczak, Nippert, & Ocheltree, 2014;Wilcox et al., 2018). In such a transition, ecosystem states do not return to previous states even when the exogenous drivers get back to levels prior to ecosystem transition. Examples include irreversible shifts from a grass to a woodland state as a result of the removal of grass-fire feedback and establishment of shrub-wind erosion feedback, which may cause land degradation. Because of its potential to cause degraded ecosystem states, being able to reliably predict WPE in water-limited ecosystems is critical for sustainable management of ecosystems under changing anthropogenic conditions. The logical question that needs to be addressed at this stage is: Are any existing numerical models suitable for modeling ecosystem shifts in water-limited ecosystems? Spatially lumped dynamical system models based on ordinary differential equations (ODE) (e.g., Lotka-Volterra type) were among the first models used to study vegetation dynamics in mixed grass and woodlands. ODE models have been formulated based on hypotheses of competition, disturbances, and predator-prey interactions to implicitly represent the landscape-scale effects of spatially explicit and complex plant-environment interactions. In particular, connections between climate and vegetation dynamics through land surface hydrological processes are oversimplified in such models, limiting their effective use in climate change research (e.g., Yatat-Djeumen et al., 2018). Therefore, ODE models have been used paradigmatically to explore mathematical properties of population dynamics, such as the emergence of thresholds and tipping points, and to sharpen science questions, rather than making management decisions (e.g., Anderies et al., 2002;D'Odorico et al., 2012;Hanan et al., 2008;May, 1977;Okin et al., 2009;Tilman, 1994;Touboul et al., 2018;Walker et al., 1981).
Spatially-explicit vegetation models can be broadly viewed in three classes: macro-scale (spatial resolution, 1-10 s of km) Terrestrial Biosphere Models (TBMs), distributed watershed ecohydrology models (10 s of m), and fine-resolution (often <10 m) vegetation pattern (self-organization) models (Fatichi et al., 2016;Fisher et al., 2014). TBMs largely rely on steady-state bioclimatic envelopes to broadly distribute plant functional types (PFT) (grass, tree, shrub) across the model domain according to climate. Once plant types are identified for each cell, TBMs simulate a collection of plant-life processes, such as physiology, phenology, carbon allocation, and senescence, coupled with vertical water and energy balance for each plant type (Fisher et al., 2014;Notaro et al., 2012). TBMs are widely used in climate change research. Watershed ecohydrology models provide spatially distributed applications of TBMs over terrain by resolving the two-way interactions between hydrologic fluxes (both vertical and lateral) and carbon and nutrient (in some models) dynamics over time scales relevant for both hydrologic and ecosystem predictions. Spatial patterns of vegetation biomass emerge from topographic controls on solar energy and redistribution of soil moisture. Spatial competition, establishment, disturbances, and mortality of species are not explicitly represented in these models (Fatichi et al., 2012;Ivanov et al., 2008;Tague & Band, 2004;Tague et al., 2013).
Organization of vegetation patterns such as fairy rings, labyrinths, and bands have been investigated with heuristic vegetation self-organization models, largely based on probabilistic state-transition rules and simple flow models. These models have been used as research tools to sharpen research questions, as a large number of transport and competition-facilitation processes conceptualized in these models have not been methodically studied in field observations. General classes of these models include: (a) reaction-advection-diffusion models (RAD) (also known as models (Turing, 1990)), in which a local function is used to express growth/decay of biomass, and spatial vegetation patterns are formed through local interactions driven by gradients of resources (water, sediment, seeds) (Klausmeier, 1999;Rietkerk et al., 2002;Saco & Moreno-de las Heras, 2013); (b) kernel-based models that represent non-local facilitation-competition interactions among plant communities (D'Odorico et al., 2006;Lefever & Lejeune, 1997;Rietkerk & Van de Koppel, 2008); (c) hybrid models that couple (a) and (b) (Thompson et al., 2008;Vincenot et al., 2016); and (d) stochastically-driven discrete cellular automaton models (CAMs) that relate state-transition probabilities to plant water stress (WS), and include rules for plant competition, grazing and fire (e.g., Caracciolo et al., 2016Caracciolo et al., , 2017Jeltsch et al., 1998;van Wijk and Rodriguez-Iturbe, 2002;Zhou et al., 2013).
Existing vegetation modeling frameworks collectively include some of the key process representations for modeling transient semiarid ecosystems. However, when taken individually, these modeling frameworks are poorly suited for handling spatial processes and are limited by long run times (e.g., TBMs), assume stationarity/linearity in spatial dynamics (e.g., RAD models), and have overly simplified process-based representations (e.g., CAMs). These constraints drastically limit the use of any of these models individually to characterize transient ecosystem response.
A flexible modeling framework that facilitates: (a) the coupling of different process representations from different models based on the needs of users; (b) rapid development and testing of new algorithms to represent transient ecosystem drivers and feedback loops; and (c) effective integration of earth surface geomorphic and ecosystems processes (e.g., vegetation-erosion feedbacks) is needed to advance the science and modeling of transient ecosystems. The Landlab open-source earth surface modeling toolkit Hobley et al., 2017) provides the fundamental building blocks needed for such a framework. In this paper, we use Landlab to address the following questions in a water-limited ecosystem in central New Mexico USA, where elevation and aspect controls on vegetation patterns have been observed at the landscape scale on pronounced topography, and active shrub-grass ecotones have been observed on flat terrain (Moreno-de las Heras et al., 2015;Yetemen et al., 2010).
1. How do exogenous drivers of topography and climate control landscape scale patterns and cover limits of vegetation in water-limited ecosystems? 2. How do interactions between grazing and endogenous grass-fire feedback influence grassland-shrubland transition? 3. Are there any predictable grass cover and connectivity thresholds that may be used as early warning indicators of rapid shrub encroachment?
In water-limited ecosystems, topographic diversity increases regional biodiversity and creates microclimatic plant refugia during dry spells, which may provide the seedling population for a quicker vegetation succession when the climate gets wetter (e.g., Cartwright et al., 2020;Hoylman et al., 2019). Therefore, with the first question, we aim to establish the relationship between vegetation cover limits in relation to water and energy availability across terrain without explicit representations of endogenous feedback loops. A frequent grassland fire regime is a major shrub disturbance that prevents shrub encroachment. The second question aims to investigate how various factors control the grass-fire feedbacks and the rate of WPE. The third question complements the second question in search of a threshold state that occurs before rapid expansion of shrubs. Identifying such thresholds would be critical for land management, as the efficacy of positive shrub feedbacks are often predicated on the reduction of grassland fire activity.
In what follows, we investigate these questions by coupling two CAMs, originally described in Nudurupati (2020), the Cellular Automaton (CA) Tree-Grass-Shrub Simulator, CATGraSS (Zhou et al., 2013), and a CA-based disturbance model described by . We selected CATGraSS (Zhou et al., 2013) because it combines basic functionality of a simplified distributed watershed ecohydrology model, such as continuous local water balance, plant growth and senescence, with probabilistic CA rules for seed dispersal, plant establishment, and mortality for multiple vegetation types, and it is computationally efficient for long-term simulations. However, CATGraSS does not represent vegetation disturbances by fire and grazing. The model of  complements CATGraSS by providing a set of endogenous rules to study grass-fire and shrub-erosion feedbacks (Wang et al., 2019).
In addressing question 1, we focus on exogenous controls of topography and climate by exploring the interactions between soil moisture-dependent local vegetation dynamics and the establishment and mortality of individual PFT (upper green blocks, Figure 1). To address questions 2 and 3, we integrate CA fire and direct removal-based grazing rules and study the role of grazing on grassland fire frequency and shrub encroachment. Models are briefly described in the next section, where we used component names (which follow CamelCase naming convention) associated with each process or a collection of processes. Further model analytical details and validation results can be found in Supporting Information S1, and in the literature (Caracciolo et al., 2014(Caracciolo et al., , 2016(Caracciolo et al., , 2017Wang et al., 2019;Zhou et al., 2013).

Soil Moisture and Vegetation Growth
Local soil moisture and vegetation dynamics are respectively represented in SoilMoisture and Vegetation components of Landlab. The SoilMoisture component tracks soil moisture using a single-layer, depth-averaged root zone water balance model for a specified PFT (tree, grass, or shrub), using an approach adapted from CATGraSS (Zhou et al., 2013). Moisture is added as pulse of individual storms and (optionally) lateral runon input [L] from upstream cells. Runoff can result from either saturation or infiltration excess runoff generation mechanisms, and it is instantaneously routed downstream in the steepest downhill direction (Collins & Bras, 2010;Yetemen, Istanbulluoglu, Flores-Cervantes, et al., 2015). Moisture loss occurs by evapotranspiration and leakage. At each iteration of the soil moisture model, water input is introduced to the soil column as an instantaneous pulse. Following this pulse, soil moisture drawdown is modeled analytically ) for a given model time step, which can be for subdaily, daily and interstorm durations. During each time step, the rate of maximum unstressed evapotranspiration (ET max ) and plant leaf area index (LAI) are assumed to be constant. LAI is used from the previous time step of the local vegetation dynamics model, while ET max is estimated for conditions within a current time step. SoilMositure component returns actual evapotranspiration, soil moisture at the beginning and end of each time step, and vegetation WS (adapted from Porporato et al. (2001)).
Vegetation component calculates Net primary productivity, NPP (g DM m −2 d −1 ) within a model time step, as a linear function of cumulative evapotranspiration and plant water use efficiency. NPP is further partitioned into aboveground live, dead, and structural biomass for each PFT. Aboveground live biomass is used to calculate plant live leaf are index (LAI l ) at the end of each model time step. Further analytical details on SoilMoisture and Vegetation components are given in Supporting Information S1. Examples of model confirmation are presented in Section 4.1 and in Supporting Information S1.
We designed two options to configure the spatial implementation of these components: fully distributed, and binned. In the distributed option soil moisture and vegetation growth are simulated in each cell of the Lanldab raster model grid individually. In the binned version, slope and aspect ranges are lumped into representative slope-aspect classes, within which representative soil moisture, local vegetation dynamics, and plant WS are simulated for each PFT and for bare soil. These variables are mapped across the watershed to the corresponding local PFT types, conditioned on slope and aspect. The binned version is computationally more efficient and can be preferred when flow routing for runoff-runon calculations is not desired. Both options can be run subdaily, daily, or interstorm time steps, by modifying the time step loop in the model driver, which is provided in the Github repository of this paper (see Acknowledgments).

Ecohydrologic Cellular Automaton Model for Spatial Vegetation Dynamics
VegCA is the Landlab component we developed for modeling soil moisture-driven spatial vegetation dynamics based CATGRaSS (Zhou et al., 2013). In VegCA, each model grid either holds a single PFT (Grass, Shrub, Tree, Shrub Seedling, or Tree Seedling) or remains empty (i.e., bare soil). The spatial pattern of PFTs evolve over time according to given probabilities for plant establishment P E , which decreases with plant WS; and mortality P M , which increases with plant WS. Plants compete to establish in bare soil cells. Shrubs can establish as new seedlings that increase the spatial frequency of shrubs, or established shrubs can expand radially by rhizomatous clonal stems (e.g., Ratajczak, Nippert, & Ocheltree, 2014). The model we propose for establishment can be used to represent both expansion mechanisms for shrubs.  (Table S1 in Supporting Information S1). The upper green block shows the componentization of CATGraSS (Zhou et al., 2013) in Landlab, used in Section 5.1. The fully coupled model is used in Section 5.3.
Mature shrubs disperse seeds to the first-ring neighbors (adjacent cells). Mature trees disperse seeds to both first and second ring neighbors (van Wijk and Rodriguez-Iturbe, 2002). The CA algorithm checks all the bare cells in the model domain, and a PFT is randomly selected for establishment based on comparison of mean live indices of plants within the first ring (for shrubs) and first and second rings (for trees) of the selected empty cell. P E is assigned based on plant live index, ɸ, of each PFT in the domain, calculated from cumulative WS within the growing season, WS: is interstorm mean plant WS [0-1], which nonlinearly quantifies a normalized deficit of mean soil moisture, , from a PFT-specific incipient stomata closure threshold, s* (Ivanov et al., 2008;Porporato et al., 2001), where s wp is soil saturation at wilting point. N is is the number of interstorm events in the year, T b (d) is the interstorm duration, and T G (d) is the length of the growing season. The mean live indices for shrubs and trees that neighbor a bare soil cell are calculated for plants in the first for shrub, SH , and both first and second ring for trees, . Mean live index for grass, , is calculated globally as an average index for all grass cells, based on the assumption of abundance of grass seeds.
Grass seeds are assumed unlimited, while mature woody plants are considered as seedling sources. Establishment of new vegetation on bare soil cells follows two mechanisms: (a) local ecohydrologic establishment that relies on soil moisture and local woody plant seed dispersal; and (b) non-local establishment that relies on landscape-scale cover fraction of woody plants (WP) and grazing rate. Here the former is described; we discussed the latter in Section 2.3. Ecohydrologic establishment uses the mean live index of a PFT surrounding a bare soil cell as a probability of establishment, with an upper limit, P E-max . Allelopathic (inhibitory) effect of shrubs on grass (Hierro & Callaway, 2003;Knipe & Herbel, 1966) is implemented by reducing the grass establishment probability by a factor, IN G . This inhibition is amplified by the number of shrub cells, N I , located in the first ring neighbors of the bare cell considered for establishment. Therefore, the aggregate water-driven local plant establishment probabilities for grass (P EW-G ), shrubs (P EW-SH ), and trees (P EW-T ) can be written as: where P E-max-G , P E-max-SH , and P E-max-T are maximum (limit) establishment probabilities of grass, shrubs, and trees respectively. For every bare cell, a PFT is randomly selected, and its corresponding P E is calculated and compared to a uniformly distributed random number ∼U(0,1). If P E is greater than the generated number, the selected PFT is established at this bare cell. The selection process of the "winner" plant type when more than one PFT has a chance of establishing at an empty cell is a challenge in CA modeling, for an alternative method, see Caracciolo et al. (2016).
Probability of mortality P M comprises three parts: drought stress (P Md ), plant aging (P Ma ), and background mortality (P Mb ) due to local disturbances (e.g., fires and grazing when not simulated explicitly). Thus, P M is defined as: Mortality probability due to drought stress is conceptualized as the difference of WS and a PFT-specific WS threshold θ PFT (van Wijk and Rodriguez-Iturbe, 2002): Mortality probability due to aging P Ma is calculated for woody plants. It begins at half the maximum age of a plant and grows linearly as plant gets older (Jeltsch et al., 1998). Total probability of mortality P M is calculated for each cell and compared to a uniform random number ∼U(0, 1). If P M is greater than the generated number, the plant is removed from the cell.
Several studies tested and confirmed the CATGraSS model for predicting local ecohydrologic states (soil moisture, biomass and LAI), evapotranspiration and landscape-scale patterns for tree, grass, and shrub PFTs in the US Southwest (Caracciolo et al., 2016;Zhou et al., 2013), Pacific Northwest (Caracciolo et al., 2017), and in a Mediterranean catchment in Sicily, IT (Caracciolo et al., 2014). The CAM can be reduced into a simplified macro-scale competition-colonization model with implicit establishment and mortality parameters when topography is neglected (e.g., Tilman, 1994), see Text S3 in Supporting Information S1 for details.

CA Rules for Disturbances and Non-Local Plant Establishment
Simple models for grazing and endogenous grassland-fire and vegetation-eolian sediment feedbacks are adapted from , who studied the role of vegetation-aeolian sediment feedbacks in the development of "islands of fertility" in shrublands. In the following, we discuss both the structure of these models and their computational implementation. Two components, ResourceRedistribution and SpatialDisturbance are developed in Landlab that update vegetation type and soil resource level. In this paper we used the SpatialDisturbance component. Details on both components are presented in Nudurupati (2020) as well as Table S1 in Supporting Information S1.
SpatialDisturbance component can simulate the initiation and spread of grassland fires as well as herbivore grazing. In arid and semiarid grasslands, natural fires are generally initiated by lightning strikes (Van Auken, 2009). The fire model starts with simulating a desired number of lightning strikes by randomly placing them across the modeled domain. Fire starts if the struck cell is vegetated and the fire vulnerability parameter of a PFT, pft_vuln [0-1] (probability of a PFT to start fire by lightning), is greater than a generated uniformly distributed random number U∼[0, 1]. The fire spreads to the vegetated neighbors if the fire susceptibility parameter of each neighboring PFT, pft_susc [0-1] (probability of a PFT to catch fire), is greater than a U∼[0,1], generated separately for each neighbor. A fire, if started or spread to a cell, converts the PFT occupying the cell to either a bare soil or a burnt PFT state. The latter is only used for aeolian erosion simulations (not discussed in this paper). A periodic boundary condition is also assumed for the spread of fires. A maximum fire area limit is included in the model as a calibration parameter, which can be used as needed. This fire model is intended for efficient long-term simulations. There are two critical limitations of our fire model. First, it is constrained to fires that are smaller than the size of the model domain. Second, it neglects fuel moisture and wind speed variables, which are used to relate the timing and spread of fires to seasonal climatology and short-term weather (e.g., Arora & Boer, 2005).
Herbivore grazing is implemented by removing grass from a cell based on the grazing pressure parameter, G, defined as the probability of converting a grass cell into a bare cell. Where G > U ∼ [0,1] grass is converted to bare soil. The influence of grazing is primarily on grass in arid and semiarid ecosystems as herbivores and introduced livestock are known to graze herbaceous vegetation (Van Auken, 2000, 2009. Woody plant (WP) establishment probability due to non-local seedling establishment, P ES , is related to the rate of grazing and natural grassland fires, using an implicit piece-wise rule varied for baseline and grazed conditions. Baseline refers to a natural grassland fire and grazing regime before the introduction of livestock grazing. We linearly relate P ES to WP cover fraction of the model domain, V WP , between a minimum and a maximum establishment probability, P 1x and P 2x , respectively; P 1x and P 2x are non-local seedling establishment probabilities that correspond to minimum and maximum WP cover fractions, V wp_Smin and V wp_Smax , respectively. x is an index that represents baseline (x = B) and grazed (x = G) conditions. α is the slope of the P ES -V WP relationship defined for the V wp_Smin -V wp_Smax range.
In this study we use V wp_Smin = 0, V wp_Smax = 0.5. For a given grazing rate of G, P 1G and P 2G are obtained from a linear approximation with defined P 1Gmax and P 2Gmax established for heavy grazing, G max .
where, P 1G and P 2G cannot be lower than P 1B and P 2B , respectively. In this paper we choose parameter values of non-local seedling establishment model consistent with  and Wang et al. (2019), who use constant seedling establishment probabilities between 0.001 (for managed grazing, G = 0.01), and 0.01 (for over grazing, G = 0.1) conditions. For baseline conditions, we use P 1B = 0.001, P 2B = 0.005 and a baseline grazing rate of G min = 0.01. This small amount of grazing is assumed to represent the impact of natural grazers. For heavy grazing (G max = 0.2), we use P 1Gmax = 0.007, P 2Gmax = 0.035.

Characterization of Vegetation Connectivity
Connectivity of grassland patches plays a commanding role in grassland-fire feedbacks and potential ecosystem shifts (Anderies et al., 2002;Wilcox et al., 2018). Types of indicators used for early warning signs of ecosystem shifts include: critical slowing-down indicators (i.e., recovery from perturbations), spatial variability indicators, and patch shape and size indicators (Burthe et al., 2016;Kefi et al., 2014). Implementations of these indicators in terrestrial ecosystems are commonly used with model-constructed fields (e.g., Guttal & Jayaprakash, 2009;Kefi et al., 2014;Thompson et al., 2021).
We introduce a landscape-scale index of vegetation connectivity, C, calculated for a selected PFT in the domain. At each vegetated cell of a certain PFT, we count the number of neighboring cells of the same PFT. Then, the number of cells in the domain that is connected to i number of neighbors, n i (i = 0,1,..8) is calculated. C is estimated as the proportion of the area that has connectivity: where, n i /N is the fraction of the area occupied by i number of connected neighboring cells of the model domain of N number of cells. We use this index to quantify the strength of grass-fire feedbacks and to search for trends that signal the onset of rapid expansion of shrubs (Phase II of WPE). Figure 7d illustrates a plot of n i /N for select years from the control simulation, and reports the calculated C for grass, C g , for each vegetation map. Note that C of a PFT, C pft , will always be smaller than the cover fraction of the same PFT in the domain, V PFT , C pft < V PFT , unless a full cover condition is observed, in which case C pft = V PFT = 1.

Study Site
Our modeling experiments are designed to broadly represent the climate and ecosystem conditions of the Sevilleta National Wildlife Refuge (SNWR) in central New Mexico, USA, which is a long-term ecological research (LTER) site. We use two domains: one with a flat topography and the other with an observed topography of a 3.3 km 2 headwater catchment located in the northwestern corner of the SNWR. The latter domain was also used by Zhou et al. (2013) (Figures 2a and 2b). This watershed has an elevation range from 1,567 m to 1,711 m, a 4.3 km west-east flowing main channel, and average hillslope gradients in the headwater valleys as high as ∼30% (Yetemen et al., 2010). A 10 m Interferometric Synthetic Aperture Radar (IFSAR) digital elevation model (DEM) is used for the modeling study (Gesch et al., 2002). Climate is characterized by a high-intensity summer monsoon season (July to September) which accounts for the ∼50% of the ∼250 mm MAP and low-intensity winter rainfalls, controlled by broad-scale frontal systems (Gutiérrez-Jurado et al., 2007).
Where topography is pronounced, ecotone boundaries and vegetation patterns are predominantly controlled by hillslope aspect, while on flat surfaces soil textural heterogeneities may play an increasing role on plant species distribution in the broader region (Buxbaum & Vanderbilt, 2007;Gosz, 1993  Plant WS is a major modeled ecohydrologic variable that reflects differences in plant rooting depth, soil moisture retention parameters (s* and s wp ), and the PFT-specific transpiration, driven by rainfall and solar radiation ( Figure 2c). In central NM grass vegetation experiences a higher overall WS across the topography than Creosotebush shrubs and Juniper trees, with minimal sensitivity to aspect and slope. Both Juniper and Creosotebush show a strong WS gradient with respect to aspect and slope. Creosotebush experiences the lowest WS of the three PFTs, and their sensitivity to aspect grows on steep slopes. Creosotebush would be the most suitable PFT for a desert environment if WS was the only constraint for plant establishment.

Model Confirmation and Design of Simulation Experiments
This section explains the choice of parameters used in Landlab, presents several examples of model confirmation, and provides details on the design of model simulation experiments used to address the research questions.
In this study, we use parameters for local SoilMoisture, Vegetation (Tables S2 and S3 in Supporting Information S1), and cellular automaton vegetation, VegCA (Table 1) components from Zhou et al. (2013), who calibrated CATGraSS at the same headwater catchment we use in this study. Their calibration did not include grazing and fire disturbances, because such disturbances are largely excluded since 1973 in SNWR. Minor modifications are made to VegCA parameters (Table S2 in Supporting Information S1) as we build explicit disturbance representations. In this study, "tree" refers to one-seed Juniper, "shrub" refers to Creosotebush, and grass refers to warm season C4 species such as black and blue grama and other various sparce surface plants.
Except for the model validation simulations where actual daily weather data is used (Supporting Information S1), in all other simulation experiments we used a single time series of seasonally varying storm depths, separated by inter-storm durations, generated using the Poisson rectangular-pulse model (e.g., Eagleson, 2005) (MAP = 250 mm, model parameters, Table S4 in Supporting Information S1). This choice is made to remove any potential consequences of using different storm sequences on vegetation patterns. Daily potential evapotranspiration (PET) for each PFT on flat surface are obtained from a sinusoidal model (parameters in Table  S5 in Supporting Information S1), fit to the Penman-Monteith estimates of daily PET for tree, grass, and shrub PFTs using daily weather data from the Deep Well weather station, located in the McKenzie Flats region of the SNWR, 30 km east of our headwater catchment study site. A ratio of sloped surface to flat surface radiation is used to obtain spatially distributed daily PET (Zhou et al., 2013). In different model experiments the model is run at a 5 m grid resolution on the DEM of the headwater catchment (Figure 2a), using both binned and distributed configurations, and on a flat domain. Unless otherwise specified, model simulations used to investigate the roles of topography and climate are initiated by randomly distributed equal proportions of each PFTs (including bare surface) and run until vegetation cover reaches a steady-state. The binned configuration is used for the watershed simulation using 108 slope and aspect bins.
In this study we verify our model with aspect-controlled vegetation patterns in our study watershed (Figure 2). Although we discuss the roles of seed dispersal strategies and grazing on the rate of WPE in relation to plant cluster size, our results currently remain theoretical, without a direct field verification of spatial plant patterns. This is partly because we did not notice strong spatial patterns other than those controlled by aspect in the vicinity of our field site. While there is evidence for plant species distribution differences controlled by soil textural differences, microclimatic and resource distribution feedbacks in the broader region (e.g., Buxbaum & Vanderbilt, 2007;D'Odorico et al., 2010), such causes are not explored in this study.

Model Confirmation
Here we provide a limited validation of Landlab by running the coupled SoilMoisture and Vegetation components (Figure 1. top block) at two sites where daily weather and various ecohydrologic observations are available (see for details in Supporting Information S1). First Landlab is run as a point model at a grassland site located near the Deep Well weather station using both daily and interstorm time steps. Time series of modeled soil moisture, evapotranspiration, and LAI are corroborated with observations ( Figure S1 in Supporting Information S1). Model runs using daily and interstom model time steps showed Nash Sutcliffe efficiency values of 0.66 and 0.81, respectively, for daily and interstorm-average soil moisture predictions against observations.
Second, the distributed model is run in the upland valley at the north-western corner of the headwater catchment study site ( At the same catchment Gutiérrez-Jurado et al. (2013) reports observed annual runoff in field plots (4 m by 2 m) ranging from 0.1% to 7% of annual precipitation. The model predicts annual runoff at the outlet of the study catchment ( Figure 2a) ∼3% (range 2.6%-3.3%) of annual precipitation. In larger watersheds within the Rio Puerco system in central New Mexico, runoffs were reported <1 up to ∼3% of annual precipitation (e.g., Molnár & Ramírez, 2001). It is also expected that growing transmission losses lead to smaller runoff ratios as watershed area gets larger (e.g., Dunkerley, 1992;Vivoni et al., 2006). Figure S1 in Supporting Information S1 shows that the model predicts seasonal ecohydrological behavior with reasonable skill for both daily and interstorm time steps. Therefore, to reduce the run time of spatially distributed simulations, we used interstorm time step for soil moisture and local vegetation dynamics, while the CA vegetation establishment and mortality model is run once a year, at the end of each growing season.

Topographic and Climatic Controls on Vegetation
Topography and climate controls on vegetation patterns (first research question) are addressed using the coupled ecohydrologic CAM (top and middle blocks of Figure 1). Results are presented in Sections 5.1 and 5.2. Fire and grazing are implicit in the constant background mortality probabilities, as our primary goal is to capture the controls of water and energy in these simulations (Table 1). First, we run the model on the DEM of the study watershed with and without the allelopathic (inhibitory) effect of shrubs on grasses (Figure 3). Although the modeled ecosystem reaches an equilibrium within roughly 2,000 years, we continued the simulations for 10,000 years to ensure ecosystem shifts do not emerge due to randomness in the representation of rainfall and plant establishment and mortality. To aid the interpretation of model results on actual topography, we perform two additional simulations on a square flat domain (500 m by 500 m), with allelopathy, one simulation for the current climate (Figures 4a-4c), and another simulation for unlimited water supply (WS = 0 for all PFTs) (Figures 4d-4f).
The role of runoff-runon processes have been studied within the context of banded and striped vegetation patterns (e.g., Saco et al., 2007;Saco & Moreno-de las Heras, 2013). Less attention has been given to the influence of runon in aspect-controlled ecosystems. Desert pavements observed on the crests of hillslopes and low infiltration capacity on EFS greatly contribute to runoff in our study site (Gutierrez-Jurado et al., 2007). Here we examine the role of runon on vegetation patterns at a smaller headwater catchment extracted from the study watershed (Figures 2a and 2b). A smaller watershed (<1 km 2 ) is used because of the computational overhead of running the distributed version of the model over several thousand years required to attain ecosystem equilibrium. We ran three simulations: a base case that represent current soil parameters; and two other simulations where bare soil infiltration capacity is reduced by half and one-tenth of the base case.
Climate sensitivity of the ecosystem is investigated using MAP between 150 and 450 mm (Section 5.2). This is achieved by scaling each storm in the generated long-term storm data for the current climate (MAP = 250), by a factor ranging from 0.6 to 1.8, for dry and wet limits, respectively ( Figure 6). This MAP range encompasses the reconstructed climate anomalies of central NM for the wetter Pleistocene and warmer and drier mid-Holocene conditions (Nudurupati, 2020), based on the work of Hall and Penner (2013). There is a negative correlation between PET and precipitation. To characterize the PET-MAP dependence, we used such a relationship from Yetemen, Istanbulluoglu, and Duvall (2015), which was developed for grassland and woodland sites in the US Midwest and Southwest. Simulations are run for the actual study watershed and a flat domain for comparison.

Grazing and Grass-Fire Feedback Controls on Shrub Encroachment
This section presents the details of model experiments used to address our second and third research questions.
Here we first discuss the selection of the fire model parameters (results in Section 5.3.1). In undisturbed natural grasslands of the Southwestern US, frequent fires, with 3-to 12-year return periods, and ∼85% plant replacement potential, have been common (e.g., Fryer & Luensmann, 2012;Schussman et al., 2006). Fire behavior in woodlands is more complex. Woodland fires typically show much longer return periods, and are highly variable in size and intensity. Broadly speaking, return intervals of replacement fires for shrublands are in the 150-300 years range; while for Pinon and Juniper woodlands, return intervals may vary from 300 to >500 years in lowlands to >1,000 years in Arizona and New Mexico Plateaus (Floyd et al., 2000;Fryer & Luensmann, 2012).
All fire models are run on a 0.25 km 2 flat terrain. Before we calibrate the fire model parameters, we first reduce the background mortality parameters of grass and woody plants by 90% and 50%, respectively, from their initial values Aspect distribution of tree cover is plotted from Lidar-based estimates reported by Zhou et al. (2013). (e, f) Time series of vegetation cover (%) of the domain, with and without shrub-grass allelopathy, respectively. Only mature shrubs and trees are included in the plots.
( Table 1). The reduced values are assumed to represent background mortality rates without fire and grazing disturbances. Each year either one or two lightning strikes are generated over the 0.25 km 2 flat domain (Orville et al., 2011). We used 0.8 for both susceptibility and vulnerability parameters for grass. These values characterize high fire susceptibility. For a 90% and connected grass cover, this gives an expected fire size of ∼70% of the model domain.
Fire susceptibility parameters for shrubs and trees are both taken ∼60% of the grass value, which we termed as moderate fire susceptibility, sh_susc = tree_susc = 0.5. A lower shrub fire susceptibility is also explored in the model, sh_susc = 0.3 (36% of gr_susc) in mixed shrub-grass simulations. Low shrub fire susceptibility represents woodlands where grassland fuels are required to spread the fire across the landscape; without grass cells, shrub fires are confined in local patches. Thus, low shrub susceptibility is deemed more realistic for the study site. Shrub and tree seedlings are assumed 20% more susceptible to fires than their corresponding mature plants ( Table 2). The fire model for trees and shrubs is calibrated individually using only the vulnerability parameters in monospecific simulations (Figures 7a and 7b). As part of this calibration, the vulnerability parameters for tree and shrub seedlings are both taken six times greater than those of the mature plants. Non-local seedling establishment is neglected in monospecific shrub and tree simulations as we relate this process to grazing.
Following the calibration of the fire model, we develop a control run on flat domain that aims to represent natural grasslands prior to the introduction of livestock grazing. Baseline grazing and baseline local and non-local seedling establishment parameters are used (Section 2.3), and the model is run for 300 years, which was sufficient to produce stable grasslands with frequent fires (Figure 7c). We then develop the following simulations to investigate shrub encroachment with the initial vegetation map obtained from the control run: 1. Effects of seedling establishment rules on shrub encroachment without fires (results in Section 5.3.2). We use local and baseline non-local plant establishment rules, without introduced herbivore grazing to establish how seedling dispersal rules control the rate and patterns of WPE ( Figure 8); 2. Effects of grazing on shrub encroachment and grass-fire feedbacks (results in Section 5.3.3). Simulations are run for both moderate and low shrub fire susceptibility by incrementing grazing rates from no grazing to heavy grazing (G max = 0.2), without and with fires. The non-local seeding rule is implemented in relation to grazing rate (Equations 9-11). The results are used to establish relationships between the time it takes for shrub dominance as a function of grazing rate, investigate the phases of shrub encroachment in relation to grass-fire feedbacks (Figures 9 and 10), and identify (if any) vegetation connectivity and cover thresholds that may be used as early warning indicators for ecosystem shift ( Figure 11);

Topographic Controls on Vegetation
Modeled vegetation maps (Figures 3a and 3c) are broadly consistent with the observed aspect-dependent vegetation patterns at the study watershed (Figures 2b and 2c). Box-whisker plots use outputs of PFT fields following ecosystem equilibrium, between model years 2,000 and 5,000, and illustrate the interannual variability of PFT cover percentages within each aspect group under ecosystem equilibrium (Figures 3b and 3d). Juniper trees grow mostly on moist PFS (N, NW, NE, slopes have ∼10%-20% tree cover on average) and to a lesser extent on east and west slopes (∼5% tree cover), consistent with the Lidar-derived tree cover for the watershed (Figures 3b and 3d). On EFS, shrub cover is slightly greater than grass cover without allelopathy (∼30% on average); with allelopathy, shrub cover (up to 40% on average) significantly dominates grasses (Figures 3b and 3d). Aspect control on vegetation fades downstream as topography becomes less pronounced, consistent with Figures 2b and 2c).
Creosotebush (shrub) is predicted as the dominant PFT at the landscape scale for the current climate (Figures 3,  4b, and 4c). Sparse grass cover is equally supported across all aspects when shrub-grass allelopathy is not modeled (Figure 3d). With allelopathy, grass cover dips up to 50% in all aspect groups. Contrary to grass and shrub cover, which show a fair amount of interannual variability (ranges up to ∼25% and ∼50% of the median cover for shrub and grass, respectively), tree cover shows much less (or almost no) interannual variability within a given aspect group. In addition, aspect distribution of tree cover is not influenced by allelopathy (Figures 3b  and 3d). The remaining simulations in this paper use shrub-grass allelopathy (IN G = 2.0, Equation 4) as it more closely reproduces the observed aspect-dependent shrub patterns in the field (Figure 2a) (see Section 3). The results presented in this section are not sensitive to the stochasticity of the stationary climate. Model converges into similar vegetation patterns when simulations are run using different generated storm time series.
In the beginning of both flat surface simulations (with and without WS), grass rapidly colonizes the empty space due to seed availability. In the water-stressed case, drought-tolerant shrubs gradually outcompete trees and grass over 500 years (Figures 4a-4c). When WS is omitted, trees dominate the ecosystem in ∼200 years as a result of their longer seed dispersal advantages over shrubs. Grass remains as secondary vegetation in both cases. For the case WS = 0, the steady-state solution of the simplified macro-scale competition-colonization model gives woody plant cover V wp = 0.59 and grass cover V g = 0.36 (see Text S3 in Supporting Information S1), which is consistent with the spatially explicit model.
The contribution of runoff-runon dynamics is not well studied in aspect-controlled ecosystems. In our experiments with varying bare soil infiltration rates at the headwater catchment (Figures 2a and 2b), watershed annual runoff proportion of precipitation grows from ∼3% (base case) to 4% and ∼7% with reduced bare soil infiltration capacity (half and one-tenth of base case of bare soil, respectively). Here we report outputs from the two endmember cases, no runon and runon with one-tenth of bare soil infiltration rate in a. As it can be seen from the maps (Figures 5a and 5b), a noticeable influence of incorporating runon contribution is that scattered trees and grass patches replace shrubs on EFS. Juniper trees benefit from runon, adding up to 15% more cover in E, SE, S, and SW slopes, where they largely occupy areas previously covered by shrubs. Grasses lose area to Juniper on W and NW facing slopes, but they too benefit from runon on other aspect classes (Figures 5c-5e). Shrubs lose ground in nearly all aspect classes. The model predicts a greater degree of variability in vegetation cover on W facing slopes when runoff-runon model is implemented. We have not further investigated the causes of this in the model. Field evidence would be necessary to investigate this model-based finding in the future.

Climatic Controls on Vegetation
Topography enables the coexistence of trees and shrubs in drier MAP ranges compared to the flat simulation ( Figure 6). A clear MAP threshold is found for the emergence and rapid response of shrubs and trees to MAP on flat domain while biome transition is more gradual on the actual topography. Sparse grass is the primary vegetation when MAP is between 150 and 200 mm in the model, largely due to the assumption of unlimited availability of grass seeds, and strong rainfall seasonality (i.e., North American Monsoon). Major spatial organization can be observed between 175 and 225 mm of MAP on PFS. As climate gets wetter, woody plants, largely shrub, first expand into PFS where Juniper trees only occur in small clusters (Figures 6c and 6d, see arrow). With MAP = 225 mm, Juniper trees mostly outcompete shrubs on PFS. Shrubs expand into the downstream valley where topography is mild, and coexist with grasses on EFS. Juniper domination expands toward the valleys with further wetting (MAP>250) and trees become the dominant cover across the landscape (MAP>300 mm) (Figures 6e and 6f).

Calibration of the Fire Model
Fire model calibration results are reported in Figure 7. Shrublands experience more frequent fires (1 fire in ∼250 years) than trees (∼1 fire in 400 years) (Figures 7a and 7b). Lower drought tolerance of trees also causes climate-driven fluctuations. In the mixed vegetation (control) simulation, a frequent grassland fire regime with ∼8-year return period emerges, and pressures WPE (Figure 7c). The connectivity distribution is right skewed with a low integrated value (C g = 0.026) after a grassland fire. It gradually gets left skewed and C g grows (C g = 0.51) as grass recovers to peak levels ( Figure 7d). Average modeled fire replacement (i.e., loss of cover) is ∼50% of the vegetation cover prior to a fire in the woodlands, and ∼80% for grassland, generally consistent with the long-term observational database of Fryer and Luensmann (2012). When the grassland simulation is rerun using only shrubs as woody plant, the outcome is nearly identical. This condition is analogous to limited native woodland cover in grasslands prior to WP encroachment in the Southwestern USA. The model output PFT field at the end of the simulation year 200 (grass cover 60%, shrub and trees<1%) is used as input for fire and grazing simulations reported next.

Effects of Seedling Establishment on Shrub Encroachment Without Fires
First, we exclude fires and grazing to examine the effects of local and non-local seedling establishment methods on WPE (Figure 8). Drought tolerant shrubs begin replacing grasses once the fire bottleneck is removed. The rate of WPE is controlled by how fast grasses die due to climatic stochasticity and shrub allelopathy, and how fast shrubs colonize empty spaces (Figure 8).
Local seed dispersal causes radial WPE, which forms interconnected structures over time (Figures 8a and 8b). This leads to a slow and gradual increase in WP cover (∼0.025% per year) within the connected grassland domain. Including the baseline global seedling establishment accelerates WPE significantly (Figure 8c). New individual plants establish and grow both laterally from parent plants and randomly with the global seeding rule. Shrubs are distributed with greater spatial frequency. This leads to a rapid decline in grass connectivity, as shrub connectivity increases more slowly. Baseline global seedling establishment yields about a factor of 5 larger WPE rate (∼0.13% per year) than the rate of WPE with local seedling dispersal (Figures 8a and 8b). However, this rate is still much lower than average reported rates in the Southwest US (∼0.5% per year) (e.g., Archer et al., 2017). Regardless of the type of seedling dispersal, shrub connectivity is smaller than grass connectivity for the same cover fraction, and cannot fully replace grass connectivity when shrub becomes the dominant PFT (Figures 8a-8c).

Effects of Grazing and Grass-Fire Feedbacks on Shrub Encroachment
Here we compare how interactions between grazing and endogenous grass-fire feedback influence grassland-shrubland transition. Shrub encroachment is exacerbated by grazing, especially with fire suppression ( Figure 9). Modeled time to shrub dominance (shrub crossover point) follows an exponential function of grazing rate, T sd = 233.27exp(-11.13 G), R 2 = 0.9632, in the absence of fires. This curve gives the lower limit for time to shrub dominance. When grassland fire is included in Landlab simulations, the delay in time to shrub dominance is controlled by the interaction of grazing rate with grass-fire feedback. Figure 10 illustrates how grazing leads to reduced fire frequency and fire area, amplifying the rate of shrub encroachment. Natural fires maintain stable grasslands with reduced cover when ≤ 0.05 (Figures 9 and 10a). Natural fires delay shrub dominance when 0.05 ≤ ≤ 0.1 , especially when shrubs have moderate fire susceptibility (Figures 10a and 10b). For G > 0.1, fires have limited ability to slow down shrub encroachment, and fire simulations collapse on the no-fire T sd -G curve.
Grazing reduces the recovery rate and peak grass cover following a fire disturbance. In the control run, grass cover usually recovers to 50%-80% range (Figure 7d). This range of recovery decreases to 40%-60% with G = 0.05 and to less than 40% with higher grazing rates ( Figure 10). Grassland will remain stable as long as the endogenous grass-fire feedback is intact and the rate of shrub establishment following each fire is less than or equal to fire-induced shrub mortality (e.g., Figure 10a, moderate fire susceptibility). Rapid expansion of shrub begins when natural grassland fire regime is replaced by less frequent and smaller fires (Figures 10b and 10c).
How is the grass-fire feedback impaired? The model provides a three-phase process for WPE, consistent with conceptual models proposed earlier based on observations (e.g., Archer et al., 2017;D'Odorico et al., 2012;Miller et al., 2000;Miller, 2005). The first phase (Phase I) is seedling establishment. New woody plants establish in the grassland at rates controlled by grass removal by grazing rate (G) and shrub fire susceptibility. The end of this phase is marked with the beginning of a period of significantly reduced fire activity (i.e., removal of a bottleneck for shrub growth). The duration of Phase I can range from several decades with high G and/or low shrub fire susceptibility (Figure 10b left and Figure 10c) to several centuries with low G and moderate shrub fire susceptibility (Figure 10a left, Figure 10b right). Decrease in grass cover to ∼40% and reduction in fire frequency, generally toward the end of Phase I, emerges as an indicator for potential future ecosystem shift to shrubland. Phase II is rapid expansion. In this phase grass fuel is not sufficient to start frequent fires. As a result, grazed grass cells are replaced at growing rates by shrubs driven by both local and grazing-enhanced non-local shrub establishment. Phase III is often referred to as canopy closure in the WPE literature (e.g., Miller, 2005;Miller et al., 2000). Here we use woody plant equilibrium, which includes both shrub canopy closure without fires (Figures 10b and 10c, left column), or newly established less-frequent shrub fire regime (Figures 10b and 10c, right column).
In Phase I, the grass connectivity, C g , fluctuates between roughly zero and 0.3, owing to frequent fires (Figure 11). At the end of Phase I, C g converges to ∼0.15, followed by its sustained drop to nearly zero, as shrubs outcompete grasses. The compensating shrub connectivity index, C s , rises rapidly to values much greater than C g of full grassland because of shrub's lower fire susceptibility that prevents frequent fires and local-seedling dispersal ability that promotes cluster growth. These three lines of evidence-loss of grass cover below %40, shrub cover in the 10%-20% range, and C g of ∼0.15-indicate potential critical threshold conditions before the onset of Phase II of the WPE in our model.

Model Comparison to Paleoecologic Evidence
Climate sensitivity of our model ( Figure 6) is broadly consistent with historical vegetation succession in the Southwestern USA (Archer et al., 2017;Hall & Penner, 2013;Holmgren et al., 2007;Van Devender, 1995). In the Southwestern USA, C 3 Juniper trees were dominant and co-exited with C 4 grasses in the wetter and cooler late Pleistocene (14 -10K yr BP) when MAP was approximately 100 mm higher than today (Hall & Penner, 2013). In our model results ecosystem is dominated by trees mixed with grasses for MAP > 350 mm, similar to late Pleistocene, except for a small fraction of shrubs, remnant from the initial condition for vegetation (i.e., model artifact) ( Figure 6).
As climate transitioned to drier conditions in the early Holocene (10-9K yr BP), similar to today's climate, C 4 grasses began to expand. The emergence and rapid establishment of C 3 shrubs as dominant plant functional type Arrow indicates overgrazing region with grazing rate of 0.1 and above  was attributed to further drying of the climate in the middle Holocene (4-5 K yr BP) and the northward migration of shrubs. At that time, MAP may have been as low as 200 mm at the field site (Nudurupati, 2020). Since then, relatively short, cooler-wetter and warmer-drier periods than the current climate followed with alternating C 4 grass and C 3 shrub dominance, respectively (Archer et al., 2017;Betancourt et al., 1990;Hall & Penner, 2013;Holmgren et al., 2007;Nudurupati, 2020;Van Devender, 1995). Model results show shrubs peak with a slightly drier MAP than today, and co-exist with grasses and trees (especially under the influence of topography) for MAP broadly in the 200-300 mm range ( Figure 6).

Water Stress and Seedling Dispersal Mediate Ecosystem Response to Climate
We found that aspect-dependent vegetation patterns under the current climate (Figure 3), as well as the ecosystem response to varying climate (Figure 6), are controlled by the interplay between plant WS and differences in seed dispersal strategies of trees and shrubs in a water-limited ecosystem. This interplay can be used as a framework to study the response of competing PFTs to climate in water-limited ecosystems.
Water stress is used for both establishment and mortality probabilities in our model. When moisture is low, the structure of the ecosystem is dictated by the order in the magnitude of WS (WS G > WS T > WS S ). On a flat domain, with the current climate, shrubs dominate the vegetation cover, as they are the most drought tolerant plants with the lowest WS. When moisture is not limited, trees outcompete shrubs due to their longer seed dispersal ranges (Figure 4). This is manifested in the form of an aspect-controlled ecotone in the topography simulation with MAP = 250 mm. Despite trees experience slightly higher WS than shrubs (∼0.15, Figure 2c) on PFS, their longer seed dispersal range gives trees a competitive advantage (Figure 3).
The model predicts a shift in the aspect distribution of woody plants as climate gets wetter. Shrubs and trees emerge on the flat domain with MAPs 200 and 250 mm, respectively. We interpret these values as representing the dry end of shrub and tree climatic limits. When climate is drier than the dry end of the shrub limit (e.g., MAP ≤ 200 mm), shrubs become limited to the cooler and wetter PFS (Figures 6c and 6d). As climate gets wetter (MAP > 200 mm) trees begin to establish and gradually become the primary vegetation type on PFS. Shrubs shift to warmer EFS slopes and flat areas (Figures 6e and 6f). Boyko (1947) observed and first termed this phenomenon the "geo-ecological law of distribution": local populations at the dry edge of a plant species climate range shift toward wetter and cooler PFS. Ackerly et al. (2020) renamed this as the hydroclimatic compensation hypothesis (HCH) and tested it in the Coast Ranges of northern California where vegetation consists of grass, shrub and woodland tree species. Ackerly et al. (2020) found that at the regional scale, populations of the same species predominantly occupy cooler PFS in hot and dry regions and shift to warmer EFS as reginal moisture increases.
Heterogeneous soil water availability induced by aspect-modulated evapotranspiration and runoff-runon processes create areas of stable ecohydrologic refugia for trees on PFS in our model. For example, the climate with MAP in the 175-250 mm range is not sufficient to support trees in the flat domain, but when given sufficient topography, trees persisted on moist PFS (Figures 6c, 6d, and 6f) and valley flats (with runon, Figure 5). At the local scale these pockets of trees are at the warm edge of their biogeographic scale, and thus can be considered sensitive hot spots to further warming (Ackerly et al., 2020). For example, in the western US, the productivity of dry forests, located on PFS and along forest ecotones to shrublands and grasslands, show greater sensitivity Figure 11. Connectivity, C, plotted for grass, shrub and total vegetation cover in 10-year increments for G = 0.085 (a, b) and G = 0.1 (c, d) with low and moderate shrub fire susceptibilities. End of woody plant encroachment phase I is marked with a vertical dashed line.
to dry climatic spells (Assal et al., 2015;Cartwright et al., 2020). However, when viewed at the regional scale, model results suggest that areas with high topographic diversity would likely host such relict tree populations during dry periods; these relict populations contribute to regional biodiversity and provide the seedling population and dispersal for a quicker succession of trees as climate gets wetter (Figures 6a and 6b) (e.g., McLaughlin et al., 2017).

Grazing and Grass-Fire Feedback Controls on Shrub Encroachment
In the absence of grazing and natural fires, our model predicts a slow linear-like increase in shrub cover (∼0.025% per year) (Figure 8a), as the initial grassland domain is filled by radial expansion of native shrub clumps (Figure 8b). Global baseline seedling makes shrub encroachment nearly 5 times faster (0.13%) than local seed dispersal, with shrubs having a greater spatial frequency and smaller cluster sizes (Figure 8c). Gradual WPE without thresholds have been reported in areas with complete fire suppression and without introduced grazers. In such cases, shrubs and trees may attain their new carrying capacities following a logistic growth model, similar to the growth curve we obtain (Figure 8c) with the inclusion of baseline non-local seedling establishment in our model (Briggs et al., 2002;Ratajczak, Nippert, Briggs, & Blair, 2014;Roques et al., 2001;Walker et al., 1981).
When natural grass-fire feedback is included in the model, a three-phase response of WPE emerges, which is consistent with the literature (e.g., Archer et al., 2017;Caracciolo et al., 2017;D'Odorico et al., 2012;Miller, 2005;Miller et al., 2000). These three phases are often termed as: seedling and early establishment (phase I), rapid expansion (phase II), and woody plant equilibrium (or canopy closure) (phase III). Grazing acts as a demographic bottleneck for grasses (e.g., Anderies et al., 2002;Wilcox et al., 2018). In our modeling experiments, the transition to Phase II is triggered with a reduction in fire frequency due to loss of grass cover and grass connectivity. Regime shift from grassland to shrubland is marked by critical thresholds that involve grass cover remaining below 40%, shrub cover increasing to 10%-20% range, and the grass connectivity, C g , transitioning from a fluctuating pattern to a decay curve with respect to time when C g < 0.15 ( Figure 11). Rates of WPE have been reported in the range of 0.1%-3.3% cover/year, that vary depending on ecoregion (Archer et al., 2017). In hot, dry regions like the Chihuahuan and Sonoran deserts, shrub encroachment rates vary from <0.1% up to 0.4% cover/year (Barger et al., 2011). In our simulations, Phase II WPE rates are in the 0.4%-0.7% cover/year range for G = 0.05-0.1 and 0.7%-0.85% cover per year range for G > 0.1 (Figures 10b and 10c).
A ∼30-year grassland manipulation experiment in a North American tallgrass prairie at the Konza Prairie Biological Station, Kansas, USA (MAP = 835 mm/y) provides direct evidence to our model-based findings. Woody plant encroachment into grasslands under grazed, ungrazed, and fire treatments were investigated in field plots (e.g., Ratajczak, Nippert, Briggs, & Blair, 2014;Ratajczak, Nippert, & Ocheltree, 2014). Shrub encroachment was observed when fire return period is greater than 3 years, first at slow rates for about 17 years (0.5% cover/ year), bringing shrub cover to 5%-10%. Rapid shrub expansion started with critical thresholds of 50%-70% grass cover and 5%-10% shrub cover (phase II) and continued for more than a decade at rates 1.4%-2.7% cover/year in both grazed and ungrazed plots. They found the increase in shrub cover correlate with a size index that quantifies the radial expansion (similar to Figures 8a and 8b) of individual shrub patches, rather than spatial frequency of shrubs. This is because the shrub species they studied reproduce by rhizomatous clonal stems.
The relationship we find between the time to shrub dominance and grazing rate (T sd -G) shows the nonlinear effects of fire-grass feedbacks (Figure 9). The cross-over from grass to shrub dominance coincides to the early stages of phase II of the WPE process ( Figure 10). The duration of Phase I greatly contributes to T sd , and it depends on the susceptibility of shrubs to fire ( Figure 10). Natural fires can delay shrub dominance beyond time scales of land management for low grazing rates. In overgrazing conditions (assumed, G > 0.1) the effect of fires greatly diminishes; the T sd -G relationship with fires collapses on that without fires; shrub dominance is observed within a time period of several decades up to 60 years. These time scales of shrub conversion are broadly consistent with those reported from observed evidence in the SNWR and Southwestern US (Archer et al., 2017;Van Auken, 2009).

Model Limitations
The Landlab ecohydrology models presented in this study are only suitable for water-limited ecosystems. Here we prioritize several key limitations for future model development and analysis: (a) Our seed dispersal rules are fairly basic. More elaborate and ecohydrologically-sound rules that involve kernel-based seed dispersal, facilitation, and competition interactions would be critical to incorporate in this model; (b) We have not explored the spatial variance and aggregation structure of PFTs, particularly grasses. A detailed analysis of grass spatial structures based on field observations and model calibration and confirmation would be critical as grass connectivity has a strong control on fire spread as well as frequency; (c) Incorporating the role of other environmental variables such as fuel moisture and wind speed on fire spread would provide a direct weather connection to fire risk; (d) Subsurface soil moisture transfer has been observed on PFS during the wet season in central New Mexico (Gutiérrez-Jurado et al., 2007. A subsurface flow transport component should be developed and its potential role on ecosystem productive should be explored.

Conclusions
Our main conclusions in response to the three questions raised in the introduction section are summarized as follows: • Modeled ecosystem response to a gradual increase in climatic wetness (i.e., higher precipitation, lower potential evapotranspiration) on topography follows the geo-ecological law of distribution (Ackerly et al., 2020;Boyko, 1947). Shrubs occupy cooler PFS at the dry end of their ecoclimatic range (MAP ≤ 200 mm), and shift to warmer EFS as regional moisture increases (MAP > 250 mm). Trees begin to occupy PFS when MAP > 200 mm, below their dry ecoclimatic range of MAP = 250 mm in the model. • Model results suggest that areas with high topographic diversity would likely host relict populations on PFS during dry periods, contributing to regional biodiversity and providing the seedling population for tree succession as climate gets wetter. • Our model illustrates a positive grass-fire feedback that prevents WPE into grasslands. WPE is predicted as a three-phase phenomenon. Phase II, rapid expansion, requires the removal of the positive grass-fire feedback by livestock grazing or fire suppression. • Regime shifts from grassland to shrubland are marked by critical thresholds that involve grass cover remaining below 40%, shrub cover increasing to 10%-20% range, and the grass connectivity, C g , remaining below 0.15. These model-based findings should be confirmed with field observations. • Relationships between time to shrub dominance as a function of grazing pressure (T sd -G) are developed with and without natural fires. Allowing natural fires during grazing can delay shrub dominance for centuries for low rates of grazing (G << 0.1).
In this paper we only explored the role of endogenous grass-fire feedback on vegetation change. Endogenous feedbacks are less constrained than exogenous drivers. Landlab lowers the barriers to implement hypotheses associated with various endogenous biophysical feedbacks. There is need to design field experiments and develop data-driven approaches that integrate data from environmental sensors and remote sensing to explore such feedbacks, parallel to model development efforts.

Data Availability Statement
The Landlab Software used for the numerical experiments presented in this manuscript is available in Hutton et al. (2020). The Python software and supporting files that include the driver scripts, plotting scripts, and data, are available in S. S. Nudurupati (2021). OAC-1147454, OAC-1450409 (Tucker, Hobley), OAC-1147519, OAC-1450338, The Oliver Fund of Tulane University (Gasparini), and EAR-1725774 (Barnhart). Landlab is a component of the Community Surface Dynamics Modeling System (CSDMS) Workbench, which is supported by NSF EAR-1831623 and EAR-2026951. Istanbulluoglu thanks Sevilleta LTER for summer faculty fellowship. We thank the editor, Prof. Peter Troch, three anonymous reviewers, and Prof. Kelly Caylor for providing constructive comments.