Simulating the Productivity of Desert Woody Shrubs in Southwestern Texas Simulating the Productivity of Desert Woody Shrubs in Southwestern Texas

In the southwestern U.S., many rangelands have converted from native grasslands to woody shrublands dominated by creosotebush ( Larrea tridentate ) and honey mesquite ( Prosopis glandulosa ), threatening ecosystem health. Both creosotebush and mesquite have well-developed long root systems that allow them to outcompete neighboring plants. Thus, control of these two invasive shrubs is essential for revegetation in arid rangelands. Simulation models are valuable tools for describing invasive shrub growth and interaction between shrubs and other perennial grasses and for evaluating quantitative changes in ecosystem properties linked to shrub invasion and shrub control. In this study, a hybrid and multiscale modeling approach with two process-based models, ALMANAC and APEX was developed. Through ALMANAC application, plant parameters and growth cycles of creosotebush and mesquite were characterized based on field data. The devel oped shrub growth curves and parameters were subsequently used in APEX to explore productivity and range condition at a larger field scale. APEX was used to quantitatively evaluate the effect of shrub reductions on vegetation and water and soil qualities in vari - ous topological conditions. The results of this study showed that this multi modeling approach is capable of accurately predicting the impacts of shrubs on soil water resources. 80 cm sampled area. Measurements of PAR were also taken with an external sensor above the shrubs concur-rently with each below-canopy measurement. The multiple above and below readings were averaged to estimate FIPAR. FIPAR was calculated as ratio of PAR below canopy to PAR above canopy subtracted from 1.0. A subsample was harvested within each sample area for the light measurement. This subsample was brought to the laboratory for LAI estimation. In the laboratory, the subsample was weighed and then separated into green leaves, dark brown live woody material, and gray dead woody material. The leaf area was measured with a LI-3100 Area Meter (LI-COR Biosciences, Lincoln, NE, USA). LAI was calculated as leaf area of subsample (cm 2 ) divided by ground area sampled (cm 2 ), and then multiplied by the ratio of total fresh weight (g) to subsample fresh weight (g). The light extinction coefficient (k) was calculated by modified Beer’s law. The value of k was calculated as the natural log of difference between 1 and FIPAR, and then divided by LAI. For creosotebush, PAR and leaf area measurement were taken from February to March in 2016, while the measurements were taken from mesquite shrubs in April, May, and July from 1993 to 1995. gramas show low productivity in arid regions due to high water stress levels. Further studies conduct a of simulation studies to demonstrate productivity of diverse native perennial grasses in the rangelands. Identification of perennial grasses that well adapted to desert arid rangeland essential process to determine the best management strategies for these lands. modeling approach developed in this study provide a realistic decision making tool that can predict results of various rangeland management strategies, which management strategies.


Introduction
Rangelands cover 31% of the total land base of the U.S. and occur mostly in western regions [1]. Western rangelands are mostly in arid and semi-arid regions that are subject to low and variable precipitation, high evaporative demand, nutrient poor soils, high spatial and temporal variability in plant production, and low net primary production [2]. Arid and semi-arid rangelands are susceptible to desertification as the result of cumulative threats such as extreme weather events (e.g. drought), land use change (e.g. suburbanization), inappropriate land management (e.g. livestock overgrazing), and invasion by shrubs and other woody plants [3,4]. Among these threats, plant invasions are considered as one of the most serious problems in much of the southwestern U.S. [5]. Encroachment of woody shrubs into grasslands has been commonly observed in the arid and semi-arid regions and often reported [6][7][8][9]. Encroachment can be defined as increasing density, cover and biomass of shrub and/or woody species in open canopy systems [8]. These woody shrubs are indigenous species that have increased in density or cover because of changes in climate variables (i.e. warmer and more humid conditions), land use modifications, or decreased frequency of disturbance regimes [8,10,11]. Extensive expansion of shrubs and woody plant into grasslands has caused largely irreversible changes in ecosystem function (e.g. alterations in landscape net primary production pattern and reduction plant biodiversity) accompanied by increased water erosion, runoff, and leaching. This has also resulted in decreased forage availability for domestic livestock and wildlife [8,[12][13][14][15][16][17][18].
services by increasing perennial grasses. Increase in the density of perennial grasses improve soil quality, increase plant richness, and provides forage for livestock and wildlife [13,14,17,18,38]. Range managers have employed a variety of management practices to remove existing shrubs such as fire [39], herbicide [40], and physical removal [41] (Figure 1). However, these practices have common limitations: logistical difficulties and side effects potentially harmful to habitat restoration [42]. These control efforts often target only one part of the life cycle of the invasive species [42]. Moreover, attempts at control have been largely decreasing due to increasing costs [43]. An effective control strategy for invasive shrubs should therefore address these challenges posed by high cost, logistical difficulties, high risk impacts on non-targeted species, and both invasion and vegetation dynamics related to climate change. Also, a successful control strategy should be designed to control the targeted invasive species and to predict their effectiveness under specific environmental conditions. Process-based models can be used to assist range managers in identifying best management strategies through providing various outcomes of short-and long-term western rangeland conditions responding to different land management strategies and rapid changes in climate and other physical processes [44,45]. To develop processbased model systems for assessing the impacts of creosotebush and mesquite in rangelands, it is important to understand factors that determine their distribution and abundance and how these relate to environmental factors. It is crucial to optimize their plant parameters describing growth in models.
Productivity of creosotebush clones is highly dependent on water availability [35,[47][48][49][50]. If there is sufficient water, creosotebush increases growth rate as new tillers initiate within a clone [35,51,52]. Mesquite productivity is also affected by water availability. According to Easter and Sosebee [52] and Ansley et al. [35], when irrigated, mesquite shrubs produce more foliage, have higher canopy cover, have higher transpirational water loss, and have lower rootto-shoot mass ratio than non-irrigated mesquites in western Texas. Soil type is also an important factor for creosotebush and mesquite establishment as determined by the soil nutrient availability as well as the soil physical characteristics. Soil physical characteristics are important because they influence surface infiltration and surface percolation [53]. Deeper horizons enriched by clay or calcium carbonate have deeper percolation depth and water availability, whereas fine-textured vesicular subsurface and surficial soil horizon development can limit infiltration. These soil characteristics differentially change availability of water for desert plants [53][54][55][56][57]. Landscape position also affects vegetative growth because it determines the  [46]. time interval between receipt of rain and its infiltration into the soil [58]. For example, creosotebush and mesquite do poorer on steep slopes with coarse, shallow soil [59][60][61] which have more runoff and less water available to plants [62]. Hamerlynck and McAuliffe [26] reported that branch mortality of creosotebush tended to increase on hillslopes, while no dead plants were found in alluvial sites.
Based on these results, creosotebush and mesquite growth varies with different rainfall, different soil, and different landscape position. Simulating creosotebush and mesquite growth chronological patterns in different desert rangelands and simulating the effects of control of these two invasive shrubs on vegetation and soil and water qualities are important when trying to control shrub productivity under various climate and soil conditions in the long term.
Scaling up from small-scale experiments to large scale field-based monitoring is an important step for reducing the long-term productivity of creosotebush and mesquite under various climate and soil conditions in future. Process-based models can simulate the effects of precipitation and geomorphic patterns in detail, estimating apparent contradictory effects. They can project variation in creosotebush and mesquite production across several different landscapes and climatic conditions. Such models can be used systemically and in combination of characteristics of hydrology, soil erosion, land slope, and nutrient balance, which are hard to approach theoretically or technically in field and plot experiments. Two field-based processlevel models, Agricultural Land Management Alternatives with Numerical Assessment Criteria (ALMANAC) and Agricultural Policy & Environmental eXtender (APEX), have potential to satisfy the needed characteristics in simulating creosotebush growth in desert rangelands.
The ALMANAC model is a process-oriented plant model that effectively simulates growth of a wide range of plant species [63,64]. Strength of ALMANAC is its capability to accurately simulate competition for light, nutrients, and water for several plant species [65]. APEX can be applied for whole-farm or small watershed (up to 2500 km 2 ) analyses and can evaluate plant growth and yield of plant species, with focus on soil and water quality in small-scale watersheds [66]. Both models operate on a daily time step. APEX's major components are climate, hydrology, plant growth, nutrient cycling, soil erosion, carbon cycling, and agricultural management practices [67]. This model uses the ALMANAC plant growth algorithms to predict productivity for over 100 plant species [67]. APEX calculates several surface hydrological parameters (daily runoff, plant transpiration, soil evaporation, water stress for plant growth, and lateral subsurface flow) in different climates having variable land topological characteristics [67,68]. Through APEX, the effects of control of invasive shrubs on soil quality can be calculated by the net differences in soil organic carbon (SOC) that occur with both invasive shrub control and no control sites.
Simulating plant development of evergreen desert shrubs like creosotebush requires some restructuring of the basic approach of degree days. Typically annual crops are simulated with a degree day sum from planting to physiological maturity for an annual growing season with crop specific values for the base temperature and optimum temperature [63][64][65]69]. This approach has also been applied to warm season perennial grasses with annual growing cycles for the leaf area and biomass [70][71][72][73]. Unlike creosotebush, mesquite is a perennial deciduous tree that drops its leaves each year and then resumes growth the following spring, each year possibly attaining (in the absence of environmental stress) its potential leaf area index value for that year. When applied to trees in Canada, the degree day sum is for a series of years so that the trees can develop over several years [72,74]. In those studies, the annual value for maximum leaf area index for the growing season increased each year to simulate how trees grow. The difficulties when attempting to transfer these approaches to desert evergreen shrubs are: (1) these shrubs do not lose their leaves during the winter of each season, (2) their phenological development is strongly tied to rainfall amount and patterns, and not just degree day accumulation, and (3) these shrubs can lose noticeable amounts of biomass due to tiller death during severe drought periods.
For this chapter, we used a multi-model combination approach, combining the strengths of two different models. A range of morphological characteristics of creosote and mesquite has been investigated from multiple locations. Based on field data, the ALMANAC model was used to create and optimize both the plant parameters and the growth curve. The resulting simulated biomass yields of creosotebush and mesquite were compared with the measured yields at the sampling locations. The resulting plant parameters and growth curve were subsequently incorporated into APEX to evaluate the effects of rainfall patterns and local soil and topological properties on the growth and productivity of the creosotebush within the sub-watershed scale in multiple regions of western Texas. The multi-model system can describe invasive plant growth and development interaction with environmental factors including light, temperature, soil characteristics and water availability. This is important to help understand why mesquite and creosotebush expand in rangelands. In addition, the multi-model system can quantitatively evaluate the invasive shrubs-perennial grasses competitive interactions in different environments and study the effects of control of invasive plants on soil organic matter and soil water content. This study will provide the desired outcomes in invasive plant management programs on rangelands.  [75], creosotebush morphological measurements were conducted at two sites in Pecos County (Fort Stockton 1 and 3), one site in Reeves County (Fort Stockton 2), and 10 sites in Brewster County (Alpine A, 1-9), all in Texas. Fort Stockton 1 was located in the right-of-way of Highway I-10, 91 km west of Fort Stockton. Fort Stockton 2 was also located in the right-of-way of Highway I-10, 61 km west of Fort Stockton. Fort Stockton 3 was inside Fort Stockton. Ten study sites (Alpine A, 1-9) were randomly selected within a 15 km wide distance on a large ranch 57 km south of Alpine. Alpine A was an airplane landing strip until 2005, so the creosote bushes there have been established for only 11 years.

Honey mesquite
As described by Kiniry [76], mesquite morphological measurement was conducted in the field located at the Grassland, Soil and Water Research Center near Temple in Bell County, Texas, U.S.

Soil and weather
For all study sites, elevation and soil type obtained from Web Soil Survey [77] (Table 1). Four weather stations which are closest to the study sites were selected for analysis. For Fort Stockton 2, the weather station in Balmorhea was selected, while for Fort Stockton 1 and 3, the weather station inside Fort Stockton was selected. The weather station inside Alpine was selected for Alpine 1-9. For mesquite study site, the weather station in Temple was selected. Total precipitation and maximum and minimum temperature from January, 1980 to March, 2016 were obtained from National Oceanic and Atmospheric Administration [78]. Detailed soil and weather information about these sites were described in the Kim et al. [75] and Kiniry [76]  maximum crown diameter. Total fresh weights of each shrub and a subsample were weighed immediately following harvest. The subsample was dried in a forced-air 66°C oven until dry weight was stabilized. Shrub height was measured from the ground to the top of the highest leaf. The thickest tiller which had no damage from insects and disease was collected from each shrub sample in all study sites. A total of 174 tillers, including 9 tillers for Fort Stockton 1-3 and Alpine 1-3, 15 tillers for Alpine 4-9, and 5 tillers for Alpine A, were used for measurements of radius of cross section of sampled tiller, growth ring count, and growth rate. As the growth of creosote bush in a dry year can be negligible, we assumed that no rings formed during severe drought years. Detailed information about morphological measurements is described in Kim et al. [75].

Honey mesquite
Mesquite seeds were collected at the Grassland, Soil and Water Research Center and planted in pots in greenhouse. Mesquite seedling 0.08 m tall was planted in plots on March 1992. To avoid competition with herbaceous plants, intensive hand hoeing with spraying chemical weed control were done every year. The experiment was laid out in randomized completed block design with four replication. Each replication was seven rows (5 m) wide with a length of 37 m. Fertilizer was applied in early 1994 and 1995. Each spring in 1993, 1994, and 1995, 18 trees per replication were randomly selected for measurements of plant height, stem diameters at the base and at half total height, and number of main stems. Among those 18 mesquite shrubs, 3 shrubs were randomly selected to measure aboveground and belowground biomass during 1993-1995. Detailed information about morphological measurements is described in Kiniry [76].

Intercepted light and leaf area measurements
Photosynthetically Active Radiation (PAR) measurements were taken using an AccuPAR LP-80 Ceptometer (Decagon Devices, Pullman, WA, USA) to enable calculation of fraction of PAR intercepted (FIPAR). Measurement of FIPAR was taken between 10:00 and 14:00. Multiple readings were made under the shrub canopy within an 80 cm x 80 cm sampled area. Measurements of PAR were also taken with an external sensor above the shrubs concurrently with each below-canopy measurement. The multiple above and below readings were averaged to estimate FIPAR. FIPAR was calculated as ratio of PAR below canopy to PAR above canopy subtracted from 1.0. A subsample was harvested within each sample area for the light measurement. This subsample was brought to the laboratory for LAI estimation.
In the laboratory, the subsample was weighed and then separated into green leaves, dark brown live woody material, and gray dead woody material. The leaf area was measured with a LI-3100 Area Meter (LI-COR Biosciences, Lincoln, NE, USA). LAI was calculated as leaf area of subsample (cm 2 ) divided by ground area sampled (cm 2 ), and then multiplied by the ratio of total fresh weight (g) to subsample fresh weight (g). The light extinction coefficient (k) was calculated by modified Beer's law. The value of k was calculated as the natural log of difference between 1 and FIPAR, and then divided by LAI. For creosotebush, PAR and leaf area measurement were taken from February to March in 2016, while the measurements were taken from mesquite shrubs in April, May, and July from 1993 to 1995.

Multi-model simulation of development
The field-based process-level models, ALMANAC and APEX, simulate processes of plant growth and soil water balance including light interception by leaves and dry matter production. Firstly, plant parameters were estimated based on 1) leaf area development; 2) development rate response to temperature; 3) radiation-use efficiency and physical descriptions; and 4) nitrogen and phosphorous concentrations in plant biomass ( Table 2). In addition, ALMANAC accounts for the effects of stresses such as nutrient deficiency, drought, and temperature on plant biomass and LAI [65]. Plant parameter values and plant growth curve were optimized through ALMANAC application using the field data. ALMANAC has been used to simulate a wide range of species, but not evergreen shrubs like creosotebush. Thus this study is the first attempt to simulate an evergreen shrub using ALMANAC. The developed plant parameters and plant growth curve were directly integrated into APEX model to simulate creosotebush and mesquite productions at a larger scale fields. The APEX model simulated water and soil qualities for each study site. In addition, APEX predicts the spatially distributed increase in water use by invasive creosotebush and mesquite within targeted watershed and also predicts effects of controlled and uncontrolled invasion on grass vegetation, water and soil conditions. Based on field data, two types of creosotebush can be categorized based on crown size: CB1 (crown size <9098 cm2) and CB2 (crown size >9098 cm2). CB1 is mostly composed of younger, small, conical shaped shrubs, while CB2 is mostly composed of older, larger, hemispherical shaped shrubs. The growth patterns of CB1 and CB2 are visually distinct (Figure 2a). In years with adequate water, new tillers grow within CB1, and the conical shaped CB1 becomes the hemispherical shaped CB2 (Figure 2b) [51,79]. Also, CB1 can reproduce either by seed (sexually) or clones (asexually). Biomass of creosotebush is dependent on the densities of CB1 and CB2, which vary among different topography features and climatic conditions. In the study sites, these two types of creosotebush co-exist at different densities (Figure 2).
In the ALMANAC Plants database, two separate sets of plant parameters named CB1 and CB2 were created for creosotebush ( Table 3). Since creosotebush is a treelike shrub, growth of CB1 and CB2 were simulated as tree growth. The crop category number (IDC) was set as 7 (evergreen tree crop). Most parameters for plant growth (e.g. DMLA, DLAI, DLAP1, DLAP2, HMX, CLAIYR, and EXTINC) were derived from measured values [75]. The base temperature (TG, °C), the temperature below which development ceases, and optimum temperature (TB, °C), the temperature at which development rate and growth rate were greatest, were estimated from the observed weather data from the three weather stations. According to Fisher et al. [80] (1988) and Newingham et al. [81] (2012), creosotebush grows slowly in spring, while fast vegetative growth and reproductive growth occur in summer. Therefore, TB and TG for creosotebush were determined from average temperatures in spring and summer, respectively. For both CB1 and CB2 plant database, TG was set as 12°C, while TB was set as 25°C. PPL1, lower plant density (plants 100 m −2 ) and fraction of maximum LAI at that density, and PPL2, higher plant density than PPL1 with fraction of the maximum LAI at that density, for CB1 were 10.01 and 35.03, respectively. PPL1 and PPL2 for CB2 were 2.20 and 11.85, respectively.
The ALMANAC code was modified to account for drought effects on development rate and for drought effects on plant stand. Creosotebush growth is largely affected by water availability. In addition, branch mortality of creosotebush increases as water deficit increases [26].  In the modified code, degree days (potential heat units -PHU) that drive plant development do not accumulate when water stress is less than or equal to 0.4. Water stress is defined as the ratio of the soil water available for ET divided by the water demand for the day, based on PET and LAI. Thus plant development stops under such drought stress. In addition, when water stress is less than 0.2, the potential leaf area index (DMLA) (a surrogate for plant stand density) decreases by 1%. This accounts for reduced plant stand with severe drought. When there is sufficient water for no water stress (water stress = 1), DMLA is increased by 1% to account for increased tillering. DMLA is not allowed to exceed the input potential value for the plant.
To simulate the annual growth cycles of creosotebush, the model simulates growth over 2 cycles each year: between October and April/May (winter/spring) and between April/May and October (summer) [82]. Due to the overlapping sequence of growth among years, we created a series of 2 year growth cycles for creosotebush (Figure 3). In the first year of a 2 year growth cycle, LAI slowly increases during spring, reaches to maximum LAI during summer, and maintains the maximum LAI during winter. In the second year of the 2-year growth cycle, the simulated LAI slowly increases during spring, and then rapidly increases during summer.  Twenty to fifty percent of aboveground biomass can be lost from plants every year, and the amount of this biomass loss depends upon the degree of habitat utilization by consumer organisms [83][84][85]. In addition, as the creosotebush grows older, its older branches gradually die and so the biomass may be reduced by decomposing dead older branches (Figure 2). Ludwig et al. [86] investigated creosote biomass in the growing seasons of 2 years and found that dead stem biomass was about half size of live stem biomass. According to Phillips and Comus [87], more than 60 species of insects are associated with creosotebush. Lac insects (Tachardiella larrea, a scale insect) are commonly found on stems and produce a lacquerlike substance by sucking juices out of the stem [87][88][89][90]. Termites, Gnathamitermes tubiformans, where they are abundant, have a significant impact on biomass loss of creosotebush. Termites consume mostly creosotebush leaf litter, especially older leaves, which apparently contain lower levels of antiherbivore allelochemics [91]. Johnson and Whitford [84] reported that termites annually consumed about 50% of the net primary production at a Chihuahuan Desert site.
In addition, the creosote grasshopper (Bootettix argentatus) lives on the plant and eats the small resinous leaves that creosotebush has developed to preserve water [87]. Mispagel [92] found that the grasshoppers consume from 0.8 to 1.9% of the cresosote bush's annual leaf biomass. Mammalian herbivores also consume great amount of the annual production [91]. For example, jackrabbits (Lepus californicus) eat leaves and stems of creosotebush [87]. Due to the factors listed above, though not measured directly in the study, it is expected that the annual creosotebush production is reduced between 20 and 50%. In the simulations, creosotebush biomass is reduced by removing 65% of its production without killing the plants in late October in the second year of 2-year growth cycle. Simulating the Productivity of Desert Woody Shrubs in Southwestern Texas http://dx.doi.org/10.5772/intechopen.73703 Since the plant stand densities of CB1 and CB2 varied among locations, the management parameter PLANTPOP (number of plants per 100 m 2 ) differed by location ( Table 4). Among the 13 sites sampled in the previous study [75], 12 sites were included in this study. "Alpine A" site was removed from this study due to small sample size. Values of PLANTPOP of CB1 and CB2 were determined based on the measured densities reported in Kim et al. [75]. Since creosotebush does not drop all leaves after maturity, large values of potential heat units (PHU) were assigned for CB1 and CB2 at all locations. The value of PHU should be close to the number of growing degree days for the area and should be large enough to avoid a terminating harvest operation in the simulations. The PHU values of 2000 and 5000 were used for CB1 and CB2, respectively, grown over the 2-year growth cycle at all sampling locations. A similar large value of PHU was also used in simulations of 2-year growth cycle sugarcane yields (Saccharum officinarum L.) in Hawaii [93].
After the initial year of establishment, the 2 year growth cycle was repeated within growth periods of creosotebush, and then 85% of biomass was harvested on the harvest date (February 1). Creosotebush is a slow-growing shrub that takes approximately 6 and 12 years to become CB1 and CB2, respectively [75]. Thus, the growth period or the number of years of simulation (NBYR) should be in the range of 6 and 12 years. Based on the field data, mean number of growth rings for CB1 and CB2 was used as the growth period or the NBYR for each site. Two more sets of tree plant parameters were used in the simulation: Tree1 and Tree2 (

Honey mesquite
Honey mesquite is a winter-deciduous tree or shrub. Growth of mesquite was simulated as tree growth. The plant category number (IDC) was set as 8 (deciduous tree plant). Most parameters for biomass-energy ratio (WA), and plant growth (e.g. DMLA, DLAI, DLAP1, DLAP2, HMX, CLAIYR, and EXTINC) were derived from measured values [76] ( Table 3). Other plant parameters were derived from previously published research, and the ALMANAC model's database of over 100 plant species' parameters, with minimal adjustment after comparing output with measured tree biomass data. The base temperature (TG, °C), the temperature below which development ceases and optimum temperature (TB, °C), the temperature at which development rate and growth rate were greatest were obtained from literature reviews.
Mesquite seedlings produce the highest biomass yields at 27°C [94]. Mesquite begins to leaf out least in April, increased until July [95]. So, the TG for mesquite was determined from average minimum temperature in April. For mesquite database, TG was set as 15°C, while TB was set as 27°C. Parameters FRST1and FRST2 indicate two points on the frost damage curve. Numbers before decimal are the minimum temperatures (C) and numbers after decimal are the fraction of biomass lost each day the specified minimum temperature occurs. According to Schuch and Kelly [96], mesquite can survive temperature down to −18°C, thus FRST1, was set to 18.3, while TB was set to 20.99. Ansley et al. [35] reported that the root-to-shoot mass ratio is 0.32 for mesquite grown under control (plants only obtained water from precipitation). A similar result was also observed in Kiniry [76] who reported root-total biomass ratio was 0.38. Based on these results, RTPRT1 and RTPRT2 were set as 0.

ALMANAC calibration and validation
To evaluate the plant parameters and test ALMANAC's ability to accurately simulated creosotebush biomass, simulated biomass values were compared with the measured biomass values from field measurements by estimating correlation and linear regression. Creosotebush simulated yields in 2016 were compared with measured yields across 12 sites [75]. For mesquite, measured values of LAI collected from May, June, July, and September in 1995 were reported by Kiniry et al. [76]. Simulated values of LAI in May, June, July, and September in 1995 were compared with the measured LAI by estimating correlation and linear regression.
Simulated biomass values in September 1993-1995 were compared with the measured biomass values from field measurements [76]. Relative ratio between simulated and measured dry biomass yields were calculated in each year. Non-irrigated capability class is defined as the suitability of soils for most kinds of field crops. The soils are grouped into three classes: VI, soils have severe limitations that make them generally unsuitable for cultivation and that restrict their use mainly to pasture, rangeland, forestland, or wildlife habitat; VII, soils have very severe limitations that make them unsuitable for cultivation and that restrict their use mainly to grazing, forestland, or wildlife habitat; VIII, soils and miscellaneous areas have limitations that preclude commercial plant production and that restrict their use to recreational purposes, wildlife habitat, watershed, or esthetic purposes [77]. Map of slope gradient was obtained from PSSAT [46].

APEX simulation development
United States Department of Agriculture-Natural Resource Conservation Service (USDA-NRCS) has conducted to detailed surveys of soils, geomorphology, and vegetation communities at O2 Ranch located in Brewster County, TX. Based on the report [46], creosotebush and mesquite were commonly found in the O2 Ranch. The creosotebush field study was conducted on the same ranch. Thus, the Alpine study site was used for APEX simulation ( Table 1 and Figure 4). The suitability of soils for most kinds of field crops in Brewster County are divided into three groups: Capability Class (CC) VI, VII, and VIII [77]. Most areas in O2 Ranch are classified as CC VI that contains soils having severe limitations that make them generally unsuitable for cultivation and that restrict their use mainly to pasture, rangeland, forestland, or wildlife habitat [46]. The projecting areas where vegetation simulation occurred are flatter and gently sloping between 0.5-1% (tangent multiplied by 100) slope terrain on the O2 Ranch [46]. Through APEX application, shrub productivity was simulated at all study subareas with different topographic features and climate conditions of small-scaled watershed in O2 Ranch in western Texas ( Figure 5). APEX performs these processes across channel systems to the outlet of a field through channels, subsurface flow, or ground water [46,97]. The first step in APEX model setup was to divide a small watershed or field into smaller spatial units called subareas (sub-watershed) represented with homogenous soils and topographic properties. Digital Elevation Model (DEM) was used in ArcAPEX to create the channel network in the study watershed, which was then used for describing the APEX routing scheme from one subarea to another and to the watershed outlet ( Figure 5). A DEM layer in 10-meter grids was downloaded for the study sites from USDA-NRCS: Geospatial Data Gateway [98]. The DEM was processed and reprojected using ArcGIS version 10.

2.2.3552.
A small field in Alpine had the highest elevation range (1063-1699 m) (Figure 5). In Alpine, 19 subareas were created within a watershed, but only five subareas (Subarea ID: 7,8,9,13,18) were used for field morphological collection. After subareas were delineated, the land use/ land cover, soils, and slope distributions were characterized for each subarea. The United States land use map was downloaded from the USDA-NRCS: Geospatial Data Gateway. The land use layer was processed and reprojected using ArcGIS version 10.2.2.3552. The soil data layer was imported from the U.S. SWAT2012 SSURGO soils database which is packaged and integrated with the ArcAPEX interface. The dominant soils in the five study subwatersheds in Alpine were gravelly loamy and bedrock soils. The dominant soils in the subwatersheds for Fort Stockton 1-3 were gravelly loam and loam. The historical weather data used in the ALMANAC model was reformatted and used for APEX simulation.

APEX calibration and validation
APEX was calibrated and validated with satellite image analysis from Subareas 7, 9, 11, and 18. Satellite image analyses of quantifying creosotebush canopy cover (density) in four study sites were performed using Texas Natural Resources Information System (TNRIS) (available at https://tnris.org/data-download/#!/quad/) and ImageJ (available at https://imagej.nih.gov/ij/ index.html). The TNRIS provided historical satellite images of the six sites within Subareas 7, 9, 11, and 18 between 1996 and 2016 ( Figure 5). Satellite images that were taken between January and February were used because only creosotebush has green leaves during the winter. January and February satellite images are only available in 1996 and 2015. The size of captured area Simulating the Productivity of Desert Woody Shrubs in Southwestern Texas http://dx.doi.org/10.5772/intechopen.73703 conducted to quantify plant canopy cover varies between study sites due to various topography features at the four study sites. The total areas used for quantifying plant canopy cover were between 0.73 and 5.94 km 2 . The captured satellite image was converted to 16-bit grayscale image using ImageJ (Figure 6), and the density of creosotebush was measured by quantifying the fraction of gray values over the entire image. The detail method is described in ImageJ manual (available in https://imagej.nih.gov/ij/docs/menus/image.html).
APEX simulated creosotebush yields in 1996 and 2015. All study subareas were configured with the same management inputs as in the ALMANAC model. Plant parameters for creosotebush and mesquite are described in Table 1. Plant parameters were transferred from the plant database in ALMANAC, except for PPL1, PPL2, COSD, PRY, EXTINC, and PHU. In APEX model, plant population parameter is in different units from ALMANAC. In APEX, the plant population parameter is expressed in number of plants per hectare, and PPLP1 should be larger than PPLP2 for tree crops. In addition, the extinction coefficient for calculating light interception is a fixed number, 0.65, which is the representative of crops with narrow row spacing [99]. APEX calculates PHU from time from planting to maturity (XMTU). Since the plant stand densities of CB1 and CB2 varied among locations, the management parameter PLANTPO (number of plants per hectare) differed by subareas ( Table 3). When more than one sampling site was in one subarea, average density of each CB1 and CB2 calculated from sampling sites were used. The growth period or number of years of simulation (NBYR) for each subarea was 25 years from establishment (1992-2016). The growth cycle management was followed by the 2-year growth cycle. The average ratio of simulated 2015 and 1996 yields in each subarea were compared with ratio of density of 2015 and 1996 by estimating correlation and linear regression.

APEX simulation on effects of invasive control
Each subarea is relatively homogeneous in terms of soil, land use, management, and weather.
In addition, the effects of invasive shrub control on vegetation and soil and water quality were evaluated. The invasive shrub-perennial grass competitive interactions were evaluated by simulating quantitatively changes in yield, soil organic carbon in plow depth in kg/ha (OCPD) and soil water content (SW) when only shrub grows, when only mixed perennial grasses grow, or when shrub and mix perennial grass grow together. Plant parameters for creosotebush and mesquite are described in Table 3. The same plant parameters and management for mixed perennial grasses (black grama, blue grama, and sideoat grama) were obtained from ALMANAC model that were already developed from plant data set collected in Texas [70,100]. Management for mixed perennial grasses in first simulated year consisted of fertilizer application on 1 April, planting in 50 plants m −2 for each grama on 10 April, and harvesting on 30 October. All black, blue, and sideoat gramas had 1800 PHUs. The simulation years for perennial grasses were same as shrub ages observed in the subareas. In addition, APEX quantitatively evaluated soil erosion, sediment yield, and water stress (days) when only shrubs grow, when only mixed perennial grasses grow, or when shrubs and mixed perennial grass grow together in all study subareas. Surface runoff (Q) was calculated using the modified Soil Conservation Service (SCS) curve number (CN) technique [101]. The SCS runoff CN can be adjusted by soil type, land use, land slope, soil water content and management practices. The CN and given daily rainfall value were used as inputs to compute soil erosion (RUS2) in each study subareas [99]. The average sediment concentration (CYAV) values were also calculated for all subareas through APEX simulation.

Creosotebush
Based on the growth patterns of creosotebush obtained from field measurements and the literature, a two-year growth cycle model was created. Results show that the LAI values for CB1 and CB2 gradually increase during spring, reach maximum LAIs during summer, and maintain the maximum LAIs during winter within the first year of the two-year growth cycle. In the following year, LAIs from the previous year gradually increase during spring, reach maximum LAIs during summer, decrease in late October, and slowly regrow during winter (Figure 3).
The simulated dry biomass yields of creosotebush were compared with the observed biomass yield from 12 study sites (Table 4). Various sizes and ages of creosotebush plants were found in different densities across 12 study sites. The biomass yields also varied across the sites and were mainly due to total shrub density and proportion of CB2 shrubs within the area. The greatest measured and simulated yields were observed in Alpine 8, which has 3.47 Mg ha −1 and 2.68 Mg/ha, respectively. The lowest measured and simulated yields were observed in Alpine 4, where the shrub density was 14 per 100 m 2 and only had CB1 shrubs. The r 2 between simulated and measured values for dry aboveground biomass based on the 1:1 line was 0.95.
These values indicate that the model performed well (Figure 7).

Honey mesquite
Honey mesquite is a winter-deciduous tree that drops its leaves in winter and leafs out again in late March or early April [96]. Simulated LAI values for honey mesquite gradually increase during spring, reach maximum LAIs during summer, and decrease in mid-October (Figure 8a). Each year shows a similar LAI developmental pattern, but the maximum LAI increases as tree age increases (Figure 6a). The simulated LAI values at the main development stages (May to September in 1995) were realistically simulated, with a highly significant fit (r 2 = 0.94) (Figure 8b).The simulated dry biomass yields of honey mesquite were compared with the observed biomass yields from Kiniry [76]. As LAI increased from 1993 to 1995 (Figure 8), the measured and simulated mesquite dry yields increased from 1993 to 1995 ( Table 5). The simulated dry yield production agreed well with measured yields of mesquite in 1993-1995 (Table 5). Relative ratios between simulated and measured dry biomass yields in 1993-1995 were obtained between 0.91-1.11.

APEX calibration and validation
The APEX model was validated by comparison of the simulated biomass yield pattern between 1996 and 2015 with real satellite imagery ( Table 6). According to Kim et al. [75] creosotebush size and leaf area index were highly correlated with aboveground biomass yield (both r > 0.8).
Thus, spatial patterns in canopy cover may directly reflect changes in yields since the study area is dominated by creosotebush. The APEX aboveground biomass yield changes between 1996 and 2015 agrees relatively well with creosotebush canopy cover (r 2 = 0.61) (Figure 9).
Both canopy cover and simulated above ground biomass yield increased between 1996 and 2015 in all subareas. The highest increases in yield and canopy cover were observed in Subarea 7 ( Table 6). Based on map of watershed (Figure 5), Subarea 7 was at higher altitude and may have much steeper slope areas. This may be why fewer creosotebush established in Subarea 7 in 1996. But, after establishing, creosotebush production may have exponentially increased by producing new tillers within clones in 2015. However, the canopy cover estimation (% km −2 ) and simulated biomass yield (Mg ha −1 ) has a weak relationship (r 2 = 0.19) (Figure 9). This may be because the study areas may have high topographic variation [102,103], or surface features (e.g. exposed rock and soil) which can create mixed pixels in satellite data [104]. Moreover, the size of satellite images is much smaller than simulating subareas, which can further confound the relationship with biomass.

Modeling the potential effects of control of invasive plants
Through APEX, effects of invasive shrub control on mix perennial production, soil organic, and water content were calculated when only invasive shrubs (creosotebush and mesquite) grow, when only mixed perennial grasses grow, or when they grow together in Alpine subareas. In APEX, the surface runoff was calculated to predict soil erosion and sediment concentration in channelized flows. Annual average values of these simulation results were summarized in Table 5  Source adapted from Kiniry [76]. Table 5.
The harvest year, canopy area per a tree, measured dry biomass yield (Mg ha −1 ), simulated dry biomass yield (Mg ha −1 ), and relative ratio between measured and simulated dry biomass yields of mesquite in 1993, 1994, and 1995 at Temple, TX.
Subareas 9, 11, and 18 had same soil erosion. When only perennial grasses (black-, blue-, and sideoat gramas) were planted, soil erosion increased. This is may be because perennial grasses have different root structures from creosotebush and mesquite. Vegetation roots were of substantial importance for soil reinforcement. Although perennial grasses have intense small roots that contribute more strength per unit area than the larger roots of creosotebush and mesquite [105], relative low vegetation cover, results in decreased root coverage per area, may increase soil erosion ( Table 7). The perennial grass yield was relatively low due to high water stress days ( Table 7). Under irrigation, mixtures of black-, blue-, and sideoat gramas can potentially produce 7 Mg ha −1 [100].
When only perennial grasses were planted, the sediment yield increased as soil erosion increased. Soil water content decreased when shrubs and perennial grasses grew together (Figure 10), which led to high numbers of water stress days ( Table 7). Since perennial grasses suffer water stress, productivity of perennial grasses when subjected to competition with shrubs is lower than perennial biomass yields with no competition. Among four subareas, Subarea 7 had the lowest perennial grass productivity due to high values of soil erosion and  Table 6. APEX subarea ID, subarea size, creosotebush canopy cover (%) per km 2 from satellite imagery, ratio of canopy cover between 2015 and 1996, ratio of canopy cover averaged within subarea, plant densities of CB1 and CB2, wet yield (Mg ha −1 ), and ratio of yield between 2015 and 1996. water stress days. Increased perennial grass density resulted in increased soil organic carbon stocks (Figure 10). Soil organic carbon yield also increased when perennial grasses and shrubs grew together (Figure 10).

Summary and conclusion
Overall, the modeling results reveal that the combined approach with ALMANAC and APEX is capable of accurately simulating the productivity of creosotebush and mesquite. Both models are capable of simulating variability of shrub yields depending on water availability. For example, the shrub yields from eroded soils are lower than those from uneroded lands. With developed plant parameters and growth cycle, APEX model is capable of simulating the effects of invasive shrubs on vegetation and soil and water qualities in different topological conditions. As shrub density decreases, the perennial grass richness, organic carbon yield, and water contents increase. However, the perennial grass mixtures with black-, blue-, and sideoat  Table 7. APEX subarea ID and annual averages of surface runoff, soil erosion, sediment yield, and water stress days for APEX subwatersheds (subareas) that were used in this study. gramas show low productivity in arid regions due to high water stress levels. Further studies are necessary to conduct a series of simulation studies to demonstrate productivity of diverse native perennial grasses in the rangelands. Identification of perennial grasses that are well adapted to desert arid rangeland is essential process to determine the best management strategies for these lands. This modeling approach developed in this study can provide a realistic decision making tool that can predict results of various rangeland management strategies, which will optimize management strategies.