Modelling the shrub encroachment in a grassland with a Cellular Automata Model

Arid and semi-arid grasslands of southwestern North America have changed dramatically over the last 150 years as a result of shrub encroachment, i.e. the increase in density, cover and biomass of indigenous shrubby plants in grasslands. Numerous studies have documented the expansion of shrublands in the southwestern American grasslands; in particular shrub encroachment has occurred strongly in part of the northern Chihuahuan desert since 1860. This encroachment has been simulated using an ecohydrological Cellular Automata model, CATGraSS. It is a spatially distributed model driven by spatially explicit irradiance and runs on a fine-resolution gridded domain. Plant competition is modelled by keeping track of mortality and establishment of plants; both are calculated probabilistically based on soil moisture stress. For this study CATGraSS has been improved with a stochastic fire module and a grazing function. The model has been implemented in a small area in Sevilleta National Wildlife Refuge (SNWR), characterized by two vegetation types (grass savanna and creosote bush shrub), considering as encroachment causes the fire return period increase, the grazing increase, the seed dispersal caused by animals, the role of wind direction and plant type competition. The model is able to reproduce the encroachment that has occurred in SNWR, simulating an increase of the shrub from 2% in 1860 to the current shrub percentage, 42%, and highlighting among the most influential factors the reduced fire frequency and the increased grazing intensity.


INTRODUCTION
Shrub encroachment is the increase in density, cover and biomass of indigenous woody or shrubby plants in various grasslands, especially arid and semi-arid grasslands (Van Auken 2000).The semi-arid and arid grasslands of southwestern North America have changed dramatically over the last 150 years.Although some authors have attributed the encroachment of shrub and woody plants to only one factor, most recent studies have suggested an interaction of several factors (van Auken 2000(van Auken , 2009)).Putative causes of increased shrubby plant abundance vary, and include increased grazing intensity, reduced fire frequency, other alterations in local land management practices and rising atmospheric CO 2 concentrations.A second type of factor includes nutrient levels, global climate change, spread of seed by livestock, small animal populations, and combinations of these factors (Van Auken 2000, 2009).The major cause of the encroachment of the shrub seems to be the reduction of grass biomass (fine fuel) by chronic high levels of domestic herbivores coupled to a reduction of grassland fires, which would have killed or suppressed the shrubby plants to the advantage of the grasses (van Auken 2009).The role of plant competition and the spread of seeds by introduced domestic herbivores seem to be secondary and probably modified the rate of the shrub encroachment (van Auken 2009).Bahre and Shelton (1993) showed that in the areas with historical climate change, the climate change did not influence the encroachment.
The encroachment of shrubs in North American deserts has been well documented for the Chihuahuan deserts.Therefore, we used data from the Sevilleta National Wildlife Refuge (SNWR), located in the northern Chihuahuan desert, New Mexico, to study the shrub encroachment.The SNWR shows a dramatic encroachment front of creosote bush (Larrea tridentata) into native desert grassland.Creosote bush encroachment in this grassland area has been simulated here using a new version of the Cellular Automata Tree-Grass-Shrub Simulator (CATGraSS) (Zhou et al. 2013).A fire cellular automata component has been introduced to simulate the fire effect.For each plant type, the probability of being ignited and killed by fire is a function of the fuel availability and of the vulnerability to fire.CATGraSS was also improved with grazing and seed dispersal functions, and its plant establishment algorithm has been modified.The causes considered here for the encroachment analysis are: the fire return period increase, the grazing increase, the seed dispersal caused by animals, the role of wind direction and the plant type competition.First the model was calibrated to reproduce the initial vegetation percentage in the study basin before encroachment in 1860 using a 5000-year base run simulation without encroachment causes.Then CATGraSS was used with the encroachment causes for 150 years, showing a good ability to simulate the encroachment, with a shrub increase from 2% to 42% (current percentage).

CATGRASS MODEL
CATGraSS uses a gridded digital elevation model (DEM) in which each model cell can hold a single Plant Functional Type (PFT: shrub, shrub seedlings, or grass), or can be bare soil.The model combines the functionality of a simplified dynamic global vegetation model (DGVM), which includes the dynamics of local water balance, plant life and plant mortality, with a rulebased probabilistic cellular automata (CA) component that simulates seed dispersal and plant establishment processes (e.g.Jeltsch et al. 1996, van Wijk andRodriguez-Iturbe 2002).It is driven by daily rainfall and maximum evapotranspiration, ET max .The salient aspect of CATGraSS, among other CA models, is the spatially explicit treatment of solar radiation on the landscape.The amount of net radiation used to calculate ET max is estimated at a daily time step as a function of day of year (DOY), slope and aspect for each PFT.Local soil and plant dynamics modelled continuously at each model element include soil moisture, ET, net primary productivity (NPP), driven by ET, and its allocation to aboveground and belowground biomass pools, and plant senescence.Plant spatial processes consist of seed dispersal and probabilistic seedling establishment based on water stress.It is assumed that shrubs provide seeds to their first (8 cells) and second ring of neighbouring cells (16 cells) following Jeltsch et al. (1996) and Van Wijk and Rodriguez-Iturbe (2002).We postulate that the probability of establishment, P E , due to seedling dispersal can be related to the ecohydrological "well-being" of the seedling source in the neighbouring plants.To measure plant well-being, the plant live index, φ X , has been introduced and defined as one minus a water stress index, WS X .WS X is defined following Porporato et al. (2001) as the cumulative plant water stress normalized by the growing season length, T seas .φ X is calculated at the end of each growing season in each model cell in the first and second ring, for shrub neighbours of a bare soil cell.
The model can be forced by either observed or generated weather data (the latter for long-term simulations).Each of the model components outlined above are discussed in detail in Zhou et al. (2013).To analyse the encroachment the following new algorithms were implemented in the CATGraSS model: fire, grazing, seed dispersal caused by animals and wind effect.

Fire
A fire CA component was introduced in the model to simulate the fire effect.The vegetation probability of death caused by fire increases with the amount of the grass fuel available (Frost and Robertson 1987) and it is therefore modelled to increase with an increase of the number of cells dominated by grass.Therefore, shrub probability of being ignited and killed by fire increases with the increase of the neighbour grass cells (Jeltsch et al. 1996).We assumed a probability of fire, P F , (i.e. the reverse of fire return period, T F ) to be compared to a random number, U~(0,1), created each year.For each year, if the probability of fire is greater than the random number (P F > U~(0,1)) then the fire starts.For each cell, the plant dies (is burned) if its vulnerability to fire, V F_X , is greater than a random number, U~(0,1), created each year (V F_X > U~(0,1)).

Grazing
The grazing effect is limited to the grass.The animals eat the herbaceous vegetation reducing the grass percentage.The grazing effect is obtained as by Zhou et al. (2013), considering a constant background probability for the disturbance factor, P Mb , to be added to the annual probability of plant mortality, P M .

Seed dispersal caused by animals
The spread of seed probability (SSP) represents the probability that a seed can arrive in the bare cell from everywhere due to animals that carry the seed (birds, cows or other animals) or through animal faeces.The CA rules for plant establishment include the following steps.In each annual iteration of the algorithm, all bare soil cells are identified and a "candidate" PFT (shrub or grass) is selected to establish in each of them.Then, the plant establishment probability, P E , is calculated for the selected PFT and compared with a uniformly distributed random number ~U(0,1).If the generated number is less than P E , the selected PFT in the first step is placed in the cell; otherwise in a second step the SSP is compared with another uniformly distributed random number ~U(0,1).If the generated number is less than SSP, the selected PFT is placed in the cell, otherwise the cell is left bare for one year.

Wind effect
To consider the wind direction effect on the seed dispersal, before estimating the mean live index of mature shrub neighbours, that is the probability of establishment for the PFTs in the bare cell, the shrub plant live index, φ SH , of the cell in the same direction as the wind speed and the φ SH of the neighbour cells of the diagonal are multiplied by the wind direction factor, WD.

SITE DESCRIPTION
To simulate the shrub encroachment, we selected a site in the SNWR, located in the northern Chihuahuan desert, approximately 32 km north of Socorro, New Mexico.More than 5% of the annual precipitation falls during the North American Monsoon (MAP ≈ 250 mm), and mean monthly temperatures range between 2.5°C in January and 25°C in July.Climatic data have been registered from Deep Well Weather Station Site (DWWSS) maintained by the Sevilleta Long Term Ecological Research (LTER), from 1990 to 2008.In the SNWR area the predominant wind speed direction is from southeast to northwest (source: http://www.eldoradocountyweather.com).The land cover map of the SNWR, overlain by the boundary of the study site and the weather station considered are shown in Fig. 1.
The study site is an alluvial fan deposit of the Sierra Ladrones formation, sandy loam in texture, overlain by a gravelly desert pavement.A Sevilleta 10 m Interferometric Synthetic Aperture Radar (IFSAR) DEM (source: http://sev.lternet.edu/)was used for the modelling study.This area is located in the elevation range of 1558 m to 1634 m, the slope angle ranges from flat surfaces to as high as 45.9° on hillslopes, with an average of 3.9°.The study site has an area of 7.34 km 2 .In Fig. 1 the actual basin vegetation distribution is shown in the SNWR Land Cover map (28.5 m resolution).The grass includes Black grama, Blue grama and Galleta (C4), while the modelled shrub is Larrea tridentata with a current coverage equal to 42%.

Model parameters
All the PFTs are defined at a 5 m grid cell resolution, which is approximately the size of a mature tree.Vegetation parameters used here for simulating local water balance and plant dynamics come from the work of Zhou et al. (2013).

Assumption for the analysis
The model has been calibrated to reproduce the initial vegetation percentage in the study area before encroachment in 1860 with a 5000-year simulation without encroachment causes (fire return period T F constant and equal to 10 years and no grazing).The model has been forced by daily rainfall and potential evapotranspiration (PET).The precipitation time series has been generated with a simplified stochastic Poisson process calibrated with the historical observed precipitation.The PET annual cycle has been obtained from a stationary cosine function fitted to mean daily values of PET calculated from the daily Penman-Monteith equation using the DWWSS data.The model was initially run with a random vegetation distribution, characterized by the same probability of assignments for all PFTs and bare soil.The shrub distribution in 1860 is not known in pattern and percentage: the shrub percentage was very low according to Knapp et al. (2008) and in this study a value of 2% was obtained from calibration at the end of the 5000-year simulation.The time series of percent coverage of PFTs in the modelled study area over the 5000 years is shown in Fig. 2. The vegetation percentage in 1860 obtained at the end of the 5000th year is: shrub 2%, grass 78%, bare soil 20%.
For the encroachment analysis, starting from 19 years of observed weather data (from 1990 to 2008) of the DWWSS, 150 years of weather data were obtained using the Advanced Weather GENerator model, AWE-GEN, (Fatichi et al. 2011) and were used in input to the CATGraSS model for the encroachment simulation from 1860 to 2010.Stationary precipitation and temperature were considered because there was no statistically significant variation of these variables in the past century.

Encroachment factors parameters used in the CATGraSS
The causes for the encroachment that were considered in this case study are summarized in Fig. 3: the grazing increase, the fire return period increase, the seed dispersal caused by animals, the role of wind direction and the plant type competition.
Van Auken (2000Auken ( , 2009) ) demonstrated that grazing increased from 1860 to 2010 causing a reduction of the grass and an increase of the shrub.In fact, in 1860 grazing was introduced in this area and it progressively increased from 1860 to 2010.According to Van Auken (2000, 2009), P Mb , which controls the background probability for grazing factor, is supposed to linearly increase from 0.02 to 0.10.Parmenter (2008) stated that the fire return period, T F , increased from 10 years in 1860 to 100 years in 2010.This increase is caused by the reduction of the fuel (grass) eaten by herbivores and is one of the possible causes of the shrub increase because by increasing the fire return period the shrub settle capacity increased.In fact, a fire return period equal to 10 years is able to maintain the semi-arid grasslands of the pre-1900s and control shrub cover.In CATGraSS we modelled the fire effect by considering a T F linear increase from 10 years to 100 years and fixing the vulnerability to fire, V F_X , equal to 0.8 for grass and to 0.11 for shrub after a sensitivity analysis.
With regard to the seed spread, Van Auken (2000, 2009) stated that the introduction of the animals in the study area from 1860 caused the possibility that a seed can arrive in the bare cell from everywhere.To take into account this aspect within the model framework, the SSP has been introduced in the model and its value has been obtained from calibration as equal to 0.005.
Finally the wind causes a directionality of the seed dispersal.This effect has been taken into account by multiplying the shrub plant live index, φ SH , of the cells in the same direction as the wind speed for WD equal to 2 and the neighbour cells of the wind direction for a WD equal to 1.5.

RESULTS AND DISCUSSION
Before simulating the encroachment using the encroachment factors, the base run simulation, i.e. the base simulation without encroachment factors was done.In this simulation there is no increase of the shrub but the shrub percentage decreases from 2 to 0.65% and the grass percentage increases from 78 to 82.43%.Then the single encroachment factors were added to see the effects of each factor on explaining the rate of ecosystem change.We have combined Fire Return Period (F), Grazing (GR), Wind (W) and Seed Dispersal caused by animals (SD), first using a single cause then combining the causes in pairs, three by three, and all the causes together for a total of 16 factor combinations.In Table 1 the final shrub (SH) and grass (G) percentages for each factor combination are shown.The best factor combinations are highlighted in bold italic.The fire frequency reduction and the increased grazing intensity have the greatest influence on the encroachment.In fact, fire causes the greatest increase of the shrub (from 0.65 to 10.32%) and grazing causes the greatest reduction of grass (from 82.43 to 72.23%).Using only a single cause, F has provided the greatest influence.Combining two factors, F and SD showed the greatest influence while combining three factors, the greatest influence is provided by F, SD and GR.
The final vegetation distribution after the encroachment using all the factors is shown in Fig. 4(a), while the time series of percent coverage of PFTs in the modelled catchment over 150 years is shown in Fig. 4(b).Combining all the factors, F-G-W-SD, the model is able to simulate  the encroachment, with an increase of the shrub from 2% to 42%, similar to the current vegetation percentage in Fig. 1, at the same time, simulating the shrub encroachment from south to north.

CONCLUSIONS
In this study, we have modelled the past and future dynamic of the shrub encroachment in Sevilleta National Wildlife Refuge, New Mexico, by improving the ecohydrological Cellular Automata model, CATGraSS.The causes that have been considered for shrub encroachment in the case study are: the fire return period increase, the grazing increase, the seed dispersal caused by animals, the role of wind direction, the shrub-grass inhibition effect (i.e.allelopathy) and the plant type competition.The model is able to simulate the encroachment, simulating an increase of the shrub from 2% to 42% (actual vegetation percentage) highlighting among the most influential factors the reduced fire frequency that causes the greatest increase of shrub, and the increased grazing intensity that causes the greatest reduction of grass, in accordance with Van Auken (2000, 2009).Therefore, the model could be used for the assessment of encroachment in future decades in the study area.

Fig. 1
Fig. 1 Study site in central New Mexico: the land cover map of the SNWR, overlain by the watershed boundary of the study site and the Deep Well meteorological station.

Fig. 2
Fig. 2 Time series of percent coverage of PFTs in the modelled catchment over the 5000 years.

Fig. 3
Fig. 3 Conceptual model illustrating the causes of the encroachment.

Fig. 4
Fig. 4 (a) Final vegetation distribution combining all the factors; (b) time series of percent coverage of PFTs in the modelled catchment over 150 years.

Table 1
Final vegetation percentage for each factors combination (bold type indicates the most important factors/combinations of factors.