Assessing vegetation structure and ANPP dynamics in a grassland–shrubland Chihuahuan ecotone using NDVI–rainfall relationships

. Climate change and the widespread alteration of natural habitats are major drivers of vegetation change in drylands. In the Chihuahuan Desert, large areas of grasslands dominated by perennial grass species have transitioned over the last 150 years to shrublands dominated by woody species, accompanied by accelerated water and wind erosion. Multiple mechanisms drive the shrub-encroachment process, including precipitation variations, land-use change, and soil erosion–vegetation feedbacks. In this study, using a simple ecohydrological modelling framework, we show that herbaceous (grasses and forbs) and shrub vegetation in drylands have different responses to antecedent precipitation due to functional differences in plant-growth and water-use patterns. Therefore, shrub encroachment may be reﬂected in the analysis of landscape-scale vegetation–rainfall relationships. We analyse the structure and dynamics of vegetation at an 18 km 2 grassland–shrubland ecotone in the northern edge of the Chihuahuan Desert (McKenzie Flats, Sevilleta National Wildlife Refuge, NM, USA) by investigating the relationship between decade-scale (2000–2013) records of remotely sensed vegetation greenness (MODIS NDVI) and antecedent rainfall. NDVI–rainfall relationships show a high sensitivity to spatial variations on dominant vegetation types across the grassland–shrubland ecotone, and provide biophysical criteria to (a) classify landscape types as a function of the spatial distribution of dominant vegetation and to (b) decompose the NDVI signal into partial components of annual net primary production (ANPP) for herbaceous vegetation and shrubs. Analysis of remotely sensed ANPP dynamics across the study site indicates that plant growth for herbaceous vegetation is particularly synchronized with monsoonal summer rainfall. For shrubs, ANPP is better explained by winter plus summer precipitation, overlapping the monsoonal period (June–September) of rain concentration. Our results suggest that shrub encroachment was not particularly active in this Chihuahuan ecotone for the period 2000–2013. How-ever, future changes in the amount and temporal pattern of precipitation (i.e. reductions in monsoonal summer rainfall and/or increases in winter precipitation) may enhance the shrub-encroachment process, particularly in the face of ex-pected upcoming increases in aridity for desert grasslands of the southwestern USA. rates that do not support long-term vegetation dynamics for the simulated rainfall conditions.


Introduction
Land degradation is pervasive across many dryland regions, which cover over 40 % of the Earth's surface and account for about 30 % of global terrestrial net primary productivity, globally supporting about 2.5 billion inhabitants (Millennium Ecosystem Assessment, 2005).Over recent decades these dryland regions have experienced growing human and climatic pressures.The most dramatic landscape alterations resulting from these pressures are those associated with desertification, which are perceived as catastrophic and largely irreversible changes that can ultimately lead to relatively barren ecosystem states (Schlesinger et al., 1990;Okin et al., 2009).A common form of vegetation change in drylands involves the encroachment of desert shrub species into arid and semi-arid grasslands, which has already affected more than 250 million hectares worldwide throughout the US, South America, southern Africa, and Australia (D'Odorico et al., 2012;Turnbull et al., 2014).
A classic case of vegetation shift is the shrubencroachment process that has been taking place over the last 150 years in the Chihuahuan Desert in southwestern USA and northern Mexico, where large areas of grasslands dominated by C 4 perennial grass species (black grama, Bouteloua eriopoda, and blue grama, B. gracilis) have been replaced by shrublands dominated by C 3 desert shrub species (mainly creosotebush, Larrea tridentata, and honey mesquite, Prosopis glandulosa).These changes in vegetation have been accompanied by accelerated water and wind erosion (for example, Schlesinger et al., 1990;Wainwright et al., 2000;Mueller et al., 2007;Turnbull et al., 2010a;Ravi et al., 2010).A complex range of mechanisms have been suggested to explain the occurrence of this vegetation transition, including external drivers that initiate the transition, and endogenous soil erosion-vegetation feedbacks that further drive vegetation change (Turnbull et al., 2012).These internal feedbacks strongly alter the organization and distribution of both vegetation and soil resources (i.e.substrate, soil moisture, and nutrients), strengthening the vegetationchange process (Okin et al., 2009;Turnbull et al., 2010bTurnbull et al., , 2012;;Stewart et al., 2014).
The onset of the grassland-shrubland transition in the Chihuahuan Desert is thought to have started with the introduction of large numbers of domestic grazers, which may have favoured the establishment of pioneer shrubs via the creation of gaps (Buffington and Herbel, 1965;van Auken, 2000;Webb et al., 2003) and via a reduction in the frequency and intensity of natural wildfires (D'Odorico et al., 2012).Changing rainfall amount and frequency has also been invoked as one of the major external drivers of shrub encroachment, which may contribute to vegetation change by shifting competitive plant physiological advantages of grass and desert shrub species (Gao and Reynolds, 2003;Snyder and Tartowsky, 2006;Throop et al., 2012).However, there remains a lack of consensus regarding changes in rainfall in the southwestern USA over recent decades.While Petrie et al. (2014) found no significant changes in precipitation at the Sevilleta Long Term Ecological Research Site in central New Mexico, other studies have reported significant increases in both annual and winter precipitation at the Jornada Experimental Range in southern New Mexico but concurrent decreases in the size of discrete precipitation events (Wainwright, 2006;Turnbull et al., 2013).
A comprehensive understanding of how desert grasslands are responding to the present variability on both climate and land use is critical for environmental management, especially in consideration of uncertainty regarding future climate change across many dryland regions.Remote sensing of vegetation provides a valuable source of information for landscape monitoring and forecasting of vegetation change in drylands (Okin and Roberts, 2004;Pennington and Collins, 2007;Moreno-de las Heras et al., 2012).Satellite-derived chlorophyll-sensitive vegetation indices, such as the normalized difference vegetation index (NDVI), provide important information on vegetation structure (e.g.surface cover, aboveground green biomass, vegetation type) and dynamics over broad spatial domains (Anderson et al., 1993;Peters et al., 1997;Weiss et al., 2004;Pettorelli et al., 2005;Choler et al., 2010;Forzieri et al., 2011).
In drylands, where vegetation dynamics are particularly well coupled with rainfall patterns, the relationship between time series of NDVI and precipitation provides specific information on the use of water for the production and maintenance of plant biomass (Pennington and Collins, 2007;Notaro et al., 2010;Veron and Paruelo, 2010).Investigations of the relationships between NDVI and rainfall suggest that arid and semi-arid vegetation responds to antecedent (or preceding cumulative) precipitation rather than to immediate rainfall, since plant growth is affected by the history of available soil moisture (Al-Bakri and Suleiman, 2004;Schwinning and Sala, 2004;Evans and Geerken, 2004;Moreno-de las Heras et al., 2012).The length (or number of days) of antecedent rainfall that best explains the NDVI (or green biomass) dynamics of dryland vegetation (hereafter optimal length of rainfall accumulation, Olr) appears to be site-specific and strongly dependent on vegetation type (Evans and Geerken, 2004;Prasad et al., 2007;Garcia et al., 2010).Herbaceous vegetation (i.e.grass and forb life forms) and shrubs usually show important differences in the patterns of vegetation growth and water use, which mediate the responses of plant biomass to rainfall in drylands (Ogle and Reynolds, 2004;Gilad et al., 2007;Pennington and Collins, 2007;Forzieri et al., 2011;Stewart et al., 2014).Thus, the study of the relationship between the NDVI and rainfall may offer important clues for detecting broad-scale landscape changes involving grassland-shrubland transitions in arid and semi-arid landscapes.
The aim of this study is to analyse vegetation structure and dynamics at a Chihuahuan grassland-shrubland ecotone (McKenzie Flats, Sevilleta National Wildlife Refuge, New Mexico, USA).To fulfill this aim we explore the relationship between decade-scale (2000-2013) records of remotely sensed vegetation greenness (MODIS NDVI) and rainfall.Our analysis is based on a new approach that examines characteristic NDVI-rainfall relationships for dominant vegetation types (i.e.herbaceous vegetation and woody shrubs) to investigate the organization and dynamics of vegetation as a way of evaluating how the shrub-encroachment process occurs.
This paper is organized into two parts.First, we present the conceptual underpinning and theoretical basis of our study by using a simple, process-based ecohydrological model to illustrate the biophysical control of the relationship between plant biomass dynamics and antecedent rainfall for dryland herbaceous and shrub vegetation.Secondly, we empirically determine reference optimal lengths of rainfall accumula-tion (in days) for herbaceous and shrub vegetation (Olr hv and Olr s ) in a 18 km 2 Chihuahuan ecotone, and use these vegetation-type specific NDVI-rainfall metrics to (i) analyse the spatial organization and dynamics of net primary production (NPP) for herbaceous vegetation and shrubs, and to (ii) explore the impact of inter-annual variations in seasonal rainfall on the dynamics of vegetation production at the grassland-shrubland ecotone.

Theoretical basis: herbaceous and shrub plant biomass-rainfall relationships in drylands
Dryland herbaceous vegetation (i.e.grass and forb life forms) and shrubs usually exhibit important differences in the patterns of vegetation growth and water use.Herbaceous vegetation typically shows quick and intense growth pulses synchronized with major rainfall events, while the dynamics of plant biomass for shrubs are generally less variable in time (Sparrow et al., 1997;Lu et al., 2003;Garcia et al., 2010).These dissimilar growth responses are controlled biophysically by the different plant-growth and mortality rates associated with herbaceous vegetation and shrubs.While grasses and forbs are associated with high rates of plant growth and mortality, shrubs are associated with comparatively lower plant-growth and mortality rates (Ogle and Reynolds, 2004;Gilad et al., 2007).We use a simplified version of the dynamic ecohydrological model developed by Rietkerk et al. (2002) to illustrate conceptually how the vegetation-specific rates of plant growth and mortality control the relationship between the dynamics of aboveground biomass and antecedent rainfall for herbaceous vegetation and shrubs in drylands.The model consists of two interrelated differential equations: one describing the dynamics of vegetation (aboveground green biomass, B, g m −2 ) and the other describing soil-moisture dynamics (soil-water availability, W , mm).
Changes in plant biomass are controlled by plant growth and mortality: where plant growth is a saturation function of soil-moisture availability, and is determined by the maximum specific plant-growth rate (g max , day −1 ), the permanent wilting point or minimum availability of soil moisture for vegetation growth (W 0 , mm), and a half saturation constant (k w , mm).Plant senescence (biomass loss) is controlled by a plantspecific mortality coefficient (m, day −1 ).Soil-water dynamics are controlled by rainfall infiltration, plant transpiration, and soil-moisture loss due to both deep drainage and direct evaporation: where water infiltration is modelled as a saturation function of plant biomass, characterized by the minimum proportion of rainfall infiltration in the absence of vegetation (i 0 , dimensionless), a half saturation constant (k i , g m −2 ), and daily precipitation (P , mm day −1 ).Plant transpiration is controlled by plant growth, and is modulated by a plant-waterconsumption coefficient (c, L g −1 ).Finally, water losses to both deep drainage and direct evaporation are modelled as a linear function of soil-water availability, with a rate r w (day −1 ).A Maple 9.5 (Maplesoft, Waterloo, Canada) code for this model is available for download as a supplement to this article (Code 1).Two sets of plant-growth and mortality coefficients were applied to this model to simulate vegetation dynamics for a herbaceous species (g max = 0.32 day −1 , m = 0.05 day −1 ) and a shrub (g max = 0.12 day −1 , m = 0.03 day −1 ), following criteria established in previous studies (Ogle and Reynolds, 2004;Gilad et al., 2007).Plant-biomass dynamics for these two vegetation types (Fig. 1a) were modelled using a northern Chihuahuan Desert 15-year daily precipitation series obtained at the Sevilleta National Wildlife Refuge (Sevilleta LTER, http://sev.lternet.edu/data/sev-1;mean annual rainfall 238 mm) and a set of parameters obtained from literature suited to dryland environments: (Rietkerk et al., 2002;Gilad et al., 2007;Saco and Moreno-de las Heras, 2013).
Using this model, we explored the strength of the plant biomass-precipitation relationship as a function of the length of rainfall accumulation (Fig. 1b).We have applied Pearson's R correlation between the simulated plant biomass for both the herbaceous and the shrub species and antecedent rainfall series using various lengths of rainfall accumulation, i.e. for any time t i in the plant biomass series, the rainfall in the preceding day (t i−1 ), the cumulative rainfall in the two preceding days (t i−1:i−2 ), in the three preceding days (t i−1:i−3 ), and so on.Modelling results show that the plant biomass-rainfall correlation is maximized at 52 days of cumulative rainfall for the simulated herbaceous species (Olr hv = 52 days) and is maximized at 104 days of cumulative rainfall for modelled shrub species (Olr s = 104 days; Fig. 1b).This result indicates that the simulated herbaceous species responds to short-term (∼2 months) antecedent rainfall for the production of plant biomass, while the simulated shrub species responds to a longer period of antecedent precipitation to support plant dynamics.Here, ARain hv and ARain s are defined as the antecedent rainfall series that optimize those vegetation-type specific relationships (i.e.time series of precedent rainfall with accumulation lengths Olr hv for herbaceous vegetation and Olr s for shrubs, Fig. 1a).Further analysis using a range of plausible values for the plantmortality and maximum plant-growth coefficients (Fig. 1c) indicates that Olr increases largely by reducing the characteristic plant-mortality and growth rates of vegetation, and therefore suggests a strong influence on vegetation type.Sen-sitivity analysis of Olr to other model parameters (Supplement Fig. S1 in the online supporting information of this study) indicates that W 0 , k w , k i , and c have negligible effects on simulated Olr values.Reductions on bare soil infiltration (i 0 ) and increases on water loss by direct evaporation and/or deep drainage (r w ) can impact Olr hv and Olr s values, ultimately amplifying the differences we obtained between vegetation types.Other factors not explicitly considered in our model, such as differences in root structure, may also reinforce herbaceous and shrub differences in timescale plant responses to antecedent precipitation (Reynolds et al., 2004;Collins et al., 2014).
The simple model presented in this study provides a good starting point for addressing general differences in plant responses to antecedent precipitation for different vegetation types in drylands.Overall, our modelling results illustrate conceptually the distinct dependence of the relationship between plant biomass and antecedent precipitation on vegetation type, particularly when comparing the dynamics of dryland herbaceous and shrub vegetation.
In the following part of this study, we empirically determine and use metrics of reference vegetation-type specific relationships between aboveground green biomass and antecedent rainfall (i.e.optimal Olr hv and Olr s lengths, and corresponding ARain hv and ARain s series) to explore the spatial organization and NPP dynamics of herbaceous and shrub vegetation at a semi-arid grassland-shrubland ecotone.

Study area
This study is conducted in the Sevilleta National Wildlife Refuge (SNWR), central New Mexico, USA, the location of the Sevilleta Long Term Ecological Research (LTER) site.The SNWR is located in the northern edge of the Chihuahuan Desert, and is a transition zone between four major biomes: the Chihuahuan Desert, the Great Plains grasslands, the Colorado Plateau steppe, and the Mogollon coniferous woodland (Fig. 2a).Livestock grazing has been excluded from the SNWR since 1973, following 40 years of rangeland use.Due to the biome-transition nature of the SNWR, minor variations in environmental conditions and/or human use can result in large changes in vegetation composition and distribution at the refuge (Turnbull et al., 2010b).Analysis of aerial photographs and soil-carbon isotopes indicates that the extent of desert shrublands has considerably increased over the grasslands in regions of the SNWR over the last 80 years (Gosz, 1992;Turnbull et al., 2008).
Our study area is an 18 km 2 grassland-shrubland ecotone within the McKenzie Flats, an area of gently sloping terrain on the eastern side of the SNWR (Fig. 2b).This study area extends over two LTER core sites established in 1999 (Fig. 2c): a desert shrubland (Creosotebush SEV Figure 1.Simulated dryland biomass-rainfall relationships for herbaceous and shrub vegetation: (a) modelled biomass dynamics for a herbaceous (green) and a shrub (red) species, (b) strength of the biomass-precipitation relationship (Pearson's R correlation) using different lengths of rainfall accumulation for the simulated herbaceous and shrub species (values above the dotted grey line are significant at P < 0.05), and (c) optimal rainfall accumulation length (Olr) as a function of the plant-growth and mortality rates.ARain hv and ARain s lines in panel (a) represent the antecedent rainfall series that best correlate with the simulated series of herbaceous and shrub biomass, respectively (i.e.time series of precedent rainfall with rainfall accumulation lengths Olr hv for herbaceous vegetation and Olr s for shrubs).The green and red dots in panel (c) indicate optimal rainfall accumulation lengths obtained for the simulated herbaceous (Olr hv , 52 days) and shrub (Olr s , 104 days) species, respectively.The (grey) "vegetation extinction" area in panel (c) reflects combined values of plant-growth and mortality rates that do not support long-term vegetation dynamics for the simulated rainfall conditions.LTER core site) dominated by creosotebush, and a grassland (Black Grama SEV LTER core site) dominated by black grama.The central and northeastern parts of the study area are mixed black and blue grama (Bouteloua eriopoda and B. gracilis, respectively) grasslands.The abundance of creosotebush (Larrea tridentata) in the grasslands is generally low, although smaller shrubs and succulents (e.g.Gutierrezia sarothrae, Ephedra torreyana, Yucca glauca, Opuntia phaeacantha) can be common.The abundance of perennial grass species decreases to the southern and southwestern parts of the study area, where creosotebush stands are widely distributed with in general low (although variable in time) amounts of annual forbs and grasses.Soils are Turney sandy loams (Soil Survey Staff, 2010) with about 60 % sand and 20 % silt content (Muldavin et al., 2008;Turnbull et al., 2010b).The climate is semi-arid, with mean annual precipitation of ∼ 240 mm that is made up of 57 % falling in the form of high-intensity convective thunderstorms during the summer monsoon (June-September) and the remainder being received as low-intensity frontal rainfall and snow (October-May).The mean annual daily temperature is 14 • C, with a winter average of 6 • C and a summer average of 24 • C. Daily air temperature rises over 10 • C in the beginning of April, leading to the onset of the yearly cycles of vegetation growth (Weiss et al., 2004).Vegetation growth in the study area generally peaks between July and September, coinciding with the summer monsoon (Muldavin et al., 2008).

Vegetation measurements (remotely sensed and ground-based) and rainfall data
We use temporal series of NDVI as a proxy of aboveground green biomass in our study area.NDVI is a remotely sensed chlorophyll-sensitive vegetation index that correlates with www.biogeosciences.net/12/2907/2015/Biogeosciences, 12, 2907-2925, 2015 green biomass in semi-arid environments (Anderson et al., 1993;Huete et al., 2002;Veron and Paruelo, 2010).Differences in soil background brightness can generate important uncertainties in relating NDVI levels to dryland vegetation, especially when vegetation cover is low and soil type is heterogeneous in space (Okin et al., 2001).Despite these uncertainties, multiple studies have demonstrated the usefulness of NDVI for examining primary production and vegetation structure in arid and semi-arid ecosystems (for example, Weiss et al., 2004;Choler et al., 2010;Moreno-de las Heras et al., 2012), and particularly in Chihuahuan landscapes with sparse vegetation (30-50 % cover) similar to those included in this study (Peters and Eve, 1995;Peters et al., 1997;Pennington and Collins, 2007;Notaro et al., 2010).
We compiled decade-scale (2000-2013) series of NDVI with a 16-day compositing period from the MODIS Terra satellite (MOD13Q1 product, collection 5, approx.250 m resolution).We used the NASA Reverb search tool (NASA EOSDIS, http://reverb.echo.nasa.gov/) to download the corresponding MODIS tiles.The data were re-projected to UTM WGS84 and further resampled to fit our 18 km 2 study area (335 pixels; 231.5 m pixel resolution after re-projection to UTM coordinates).We checked the reliability layer of the acquired MODIS products and discarded those NDVI values that did not have the highest quality flag value (less than 1 % of data).Missing values were interpolated using a second-order polynomial.To reduce inherent noise, the NDVI time series were then filtered by applying a Savitzky-Golay smoothing algorithm, as recommended by Choler et al. (2010).
To validate remote-sensing analysis of the spatial distribution of vegetation types, the dominance of herbaceous vegetation, shrubs, perennial grass, forbs, and creosotebush plants was recorded at a set of 27 control points (Fig. 2c) using the point-intercept method (Godin-Alvarez et al., 2009).Vegetation presence/absence of the aforementioned vegetation types was recorded every metre using a 2 cm diameter, 1.2 m tall, metal rod pointer along five 50 m long linear transects that were laid at each control point at random directions (without overlapping).Dominance was determined as the relative abundance of a particular vegetation type in relation to the total amount of vegetated points found per linear transect.
Reference information on aboveground net primary production (NPP) was obtained from a pre-existing, decadescale (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011) data set (Sevilleta LTER, http://sev.lternet.edu/data/sev-182).This data set was recorded in a set of 10 sampling webs distributed within the Black Grama and Creosotebush SEV LTER core sites (five webs per core site, Fig. 2c).Each sampling web consisted of four 25 m 2 square sub-plots located in each cardinal direction around the perimeter of a 200 m diameter circle, with four 1 m 2 quadrats spatially distributed in the internal corners of the 25 m 2 subplots.A detailed description of the methods that were applied for the development of the SEV LTER field NPP data set can be found in Muldavin et al. (2008).Briefly, species-specific plant standing biomass was estimated three times per year (in February-March, May-June, and September-October) using allometric equations, and NPP was calculated seasonally for spring (the difference in plant biomass form March-May), summer (from June to September), and autumn/winter (from October to February).For this study, we have used lumped records of annual net primary production (ANPP) for herbaceous vegetation and shrubs that were spatially averaged at the core site scale.ANPP for each yearly cycle of vegetation growth has been calculated as the sum of the seasonal NPP records (i.e.spring + summer + autumn/winter).
Daily rainfall information for this study was obtained from an automated meteorological station located in the study site (the Five Points weather station, SEV LTER, Fig. 2c; Sevilleta LTER, http://sev.lternet.edu/data/sev-1).The meteorological station is equipped with a rain gauge that records precipitation on a 1 min basis during periods of rain.

Reference NDVI-rainfall metrics for herbaceous vegetation and shrubs
We explored reference NDVI-rainfall relationships for herbaceous vegetation and shrubs in the Black Grama and Creosotebush SEV LTER core sites (where vegetation is dominantly herbaceous and shrub, respectively) using the 2000-2013 NDVI time series (averaged from five MODIS pixels in each site, covering a total of 1200 m 2 per site).
Pearson's correlations between NDVI and antecedent precipitation series were calculated for the two sites using various lengths of rainfall accumulation (1-300 days).Optimal lengths of rainfall accumulation for herbaceous vegetation and shrubs (Olr hv and Olr s , respectively) were then determined as the length of rainfall accumulation (in days) of the antecedent precipitation series that maximized the correlations between NDVI and rainfall in the black grama-and the creosotebush-dominated core sites, respectively.Growth of non-dominant herbaceous vegetation in arid shrublands can make the detection of shrub-specific NDVI-rainfall metrics (i.e.Olr s ) difficult due to the emergence of secondary Olr hv values, particularly in wet years with strong herbaceous production (Moreno-de las Heras et al., 2012).We applied detailed analysis of the NDVI-rainfall relationships in the core sites for each annual cycle of vegetation growth to facilitate discrimination of the Olr hv and Olr s metrics.Our approach assumes linearity between rainfall and both NDVI values and green biomass, which has been broadly demonstrated to occur for dryland vegetation (Evans and Geerken, 2004;Choler et al., 2010;Notaro et al., 2010;Veron and Paruelo, 2010;Moreno-de las Heras et al., 2012) and particularly in our grassland-shrubland desert ecotone (Pennington and Collins, 2007;Muldavin et al., 2008).
The optimal antecedent rainfall series determined in the core sites for herbaceous vegetation (ARain hs , with Olr hv length of rainfall accumulation) and shrubs (ARain s , with Olr s rainfall accumulation length) were further used in our 18 km 2 ecotone to classify landscape types and to decompose local NDVI signals into greenness components for herbaceous and shrub vegetation.

Spatial distribution of vegetation types and landscape classification
We applied analysis of the relationship between local series of NDVI and the reference ARain hv and ARain s antecedent rainfall series to determine the spatial distribution of dominant vegetation and classify landscape types over our 18 km 2 ecotone study area.This analysis builds on the assumption that spatial variations in the NDVI-rainfall relationship reflect spatial differences in the dominance of vegetation types.We assume that areas dominated by herbaceous vegetation (or shrubs) will show a strong NDVI-rainfall relationship for the herbaceous-characteristic ARain hv (or the shrub-characteristic ARain s ) antecedent rainfall series along the study period.
The strength of the relationship between NDVI and rainfall (quantified using Pearson's R correlation between NDVI and antecedent precipitation) was calculated for every MODIS pixel in the study area using the reference ARain hv and ARain s antecedent rainfall series.Correlation values were determined for each cycle of vegetation growth (April-March) in 2000-2013.In order to reduce data dimensionality, we applied principal component analysis (PCA) using the calculated correlation coefficients as variables for analysis (28 variables resulting from the two vegetation-specific antecedent rainfall series and the 14 growing cycles).We studied further the relationship between the main PCA factors and ground-based dominance of vegetation types using the reference vegetation distribution data set (27 control points).Finally, we used the empirical relationships between vegetation dominance and the main PCA factors to discriminate differentiated landscape types across the study area: grass-dominated (GD), grass-transition (GT), shrubtransition (ST), and shrub-dominated (SD) landscapes.

NDVI decomposition and transformation into herbaceous and shrub ANPP components
Time series of NDVI at any specific location reflects additive contributions of background soil and the herbaceous and woody shrub components of vegetation (C bs ,C hv , and C s , respectively) for that particular site (Lu et al., 2003): Montandon and Small (2008) carried out in situ measurements of field spectra convolved by the MODIS bands to determine the background soil contribution to NDVI in the SNWR.They obtained a soil NDVI value of 0.12 for Turney sandy loam soils, which are broadly distributed across the McKenzie Flats.Analysis of the local MODIS NDVI time series revealed that this soil background reference value broadly matches the minimum NDVI values for our study area.Application of reference soil values in NDVI decomposition and normalization methodologies provides an efficient standardization approach for characterizing the background soil baseline, particularly in areas with homogeneous soils (Carlson and Ripley, 1997;Roderick et al., 1999;Lu et al., 2003;Choler et al., 2010).Soil background NDVI may change with soil-moisture content (Okin et al., 2001).
Although this effect can be especially important for dark organic-rich soils, soil-moisture variations have shown little impact in desert-type bright sandy and sandy-loam soils, like those represented in the study area (Huete et al., 1985).Therefore, a constant value of 0.12 was applied to subtract the background soil baseline (C bs ) from the NDVI time series, obtaining a new set of soil-free series (NDVI O ): We applied the reference herbaceous-and shrubcharacteristic antecedent rainfall series, ARain hv and ARain s , to partition single time series of soil-free NDVI (NDVI O ) into separate contributions for herbaceous vegetation (C hv ) and woody shrubs (C s ) across our study area.This approach is based on the assumption that the primary determinant of the dynamics of both NDVI and green biomass in Chihuahuan landscapes is the rainfall pattern (Huenneke et al., 2002;Weiss et al., 2004;Muldavin et al., 2008;Pennington and Collins, 2007;Notaro et al., 2010;Forzieri et al., 2011), and therefore the partial contributions of herbaceous vegetation and shrubs to NDVI can be estimated as a function of their characteristic dependency on antecedent rainfall.In other words, we assume that C hv and C s for any t i are proportional to ARain hv and ARain s .The NDVI components for herbaceous vegetation and shrubs were partitioned using the following two-step NDVI-decomposition procedure (Maple 9.5 code for analysis provided as online supporting material of this article; Code 2).First, we applied first-order least-squares optimization of the relationship between soil-free NDVI (NDVI O ) and the vegetation-type specific antecedent rainfall series (ARain hv and ARain s for herbaceous vegetation and shrub, respectively): where h and s represent vegetation-type specific rainfall-NDVI conversion coefficients for the herbaceous and shrub components.
Secondly, we used the determined coefficients h and s to calculate the weights of C hv and C s on the time series (i.e. the predicted percentage contribution of each vegetation type over the predicted totals for any t i ).Seasonal variations in other environmental factors (e.g.temperature, day length) may influence NDVI dynamics for Chihuahuan vegetation, shaping the responses of vegetation to precipitation (Weiss et al., 2004;Notaro et al., 2010) The 2000-2013 time series of NDVI were decomposed into separate contributions of herbaceous vegetation and shrubs for the Black Grama and Cresotebush SEV LTER core sites.We used the reference 2000-2011 field NPP data set to study the relationship between the decomposed NDVI time series and ground-based estimates of herbaceous and shrub NPP for the core sites.The sum of the herbaceous and the shrub NDVI components ( NDVI veg.type ) were calculated for each growing cycle of vegetation (April-March).We further determined the relationships between field ANPP estimates of herbaceous and shrub vegetation and NDVI veg.type .Finally, we applied the signaldecomposition procedure to every single NDVI time series of the 335 MODIS pixels contained within our study area.The established core site NDVI-ANPP relationships were used to estimate herbaceous and shrub ANPP across the 18 km 2 study site.

Spatiotemporal dynamics of vegetation production and impact of seasonal precipitation on herbaceous and shrub ANPP
We used the remotely sensed ANPP estimations and landscape-type classification (GD, grass-dominated; GT, grass-transition; ST, shrub-transition; and SD, shrubdominated landscapes) to analyse the spatiotemporal dynamics of ANPP along our study grassland-shrubland ecotone, applying repeated-measures ANOVA with time as within subjects factor and landscape type as between subjects factor.Departures from sphericity were corrected using the Greenhouse-Geisser F-ratio method for repeated-measures ANOVA (Girden, 1992).The 2000-2013 activity of the shrub-encroachment phenomenon for the established landscape types (GD, GT, ST, and SD) was explored applying Pearson's R correlation between shrub contribution to total ANPP and time.
We used three different seasonal precipitation metrics to analyse the impact of inter-annual variations in seasonal precipitation on the production of herbaceous and shrub vegetation at our ecotone: (i) preceding non-monsoonal rainfall (Rain PNM , from October to May) that takes place before the summer peak of vegetation growth, (ii) summer monsoonal precipitation (Rain SM , from June to September), and (iii) late non-monsoonal rainfall (Rain LNM , from October to March) that takes place at the end of the annual cy-cles of vegetation growth.The effects of seasonal precipitation on herbaceous and shrub ANPP for the established landscape types (grass-dominated, grass-transition, shrubtransition, and shrub-dominated landscapes) were explored by applying Pearson's R correlation.Effect significance and size was determined using a general linear model (GLM) that includes the different sources of seasonal precipitation (Rain PNM , Rain SM , and Rain LNM ) as covariates, landscape type (LT) as a factor, and the interaction terms between landscape type and seasonal precipitation (LT : Rain PNM , LT : Rain SM , and LT : Rain LNM ).

Patterns of greenness and reference NDVI-rainfall metrics in the core sites
Inter-and intra-annual variations of NDVI show similar patterns of vegetation greenness for both the Black Grama and the Creosotebush core sites (Fig. 3a).The signal generally peaks slightly in spring (May) and strongly in summer (July-September).The lowest NDVI values are observed between February and April.Summer peaks in NDVI values are, however, less marked in the Creosotebush core site.In addition, the NDVI signal for the creosotebush-dominated site generally shows a autumn (October-November) peak that is especially important during particular growing cycles (2000-2001, 2001-2002, 2004-2005, 2007-2008, 2009-2010).
Correlations between NDVI and antecedent precipitation using rainfall-accumulation lengths of 1-300 days indicate that an optimal short-term cumulative rainfall period of 57 days best explains the NDVI variations for the dominant herbaceous vegetation of the grassland site (ARain hv antecedent rainfall series, with Olr hv accumulation length; Fig. 3; see also Supplement Fig. S2 in the online supporting information of this study for details on the annual cycles of vegetation growth).For the Creosotebush core site (with dominant shrub vegetation and subordinate forbs and grasses), the short-term, 57-day antecedent rainfall series ARain hv also has an important impact on the strength of the NDVI-rainfall relationship, particularly for three consecutive growing cycles with strong summer precipitation (2006-2007, 2007-2008, and 2008-2009; summer precipitation for the period is 40 % above the long-term mean).However, the NDVI-rainfall correlation in this shrub-dominated site generally peaks using a much longer optimal cumulative rainfall period of nearly 145 days (ARain s series, with Olr s length).

Spatial distribution of vegetation types and landscape classification
PCA of the NDVI-rainfall correlation coefficients (per growing cycle) for the reference 57-and 145-day antecedent rainfall series (i.e.ARain hv and ARain s with Olr hv and Olr s rainfall accumulation lengths, respectively, for all MODIS pix- els contained within our study area) shows that PCA factor 1 (about 40 % of total data variance) reflects a landscape gradient that discriminates the two reference responses of vegetation greenness to antecedent rainfall (Fig. 4a and  b).The correlation between the NDVI and the short-term antecedent rainfall series ARain hv increases to the negative side of factor 1 (particularly for growing cycles 2001-2002, 2002-2003, 2005-2006, and 2012-2013), while the correlation with the 145-day antecedent rainfall series (ARain s ) increases to the positive side of the this factor (particularly for cycles 2000-2001, 2002-2003, 2005-2006, and 2006-2007;Fig. 4b).Analysis of the relationship between PCA factor 1 and vegetation dominance for the ground-based set of control points indicates that this landscape gradient is explained by the field distribution of dominant vegetation types since the dominance of herbaceous vegetation and shrubs increases to the negative and positive side of PCA factor 1, respectively (R 2 approx.0.90, Fig. 4c).Four different landscape types (GD, GT, ST, and SD) are defined in the 18 km 2 study area as determined by the spatial projection of the relationship between PCA factor 1 and field dominance of herbaceous and shrub vegetation (Fig. 4c and  d).SD, ST, and GT landscapes are distributed in the southwestern part of the study site, while GD landscapes are located in the central and northeastern parts of the area (Fig. 4d  and e).

NDVI transformation into herbaceous and shrub ANPP components
Temporal decomposition of NDVI into partial herbaceous and shrub vegetation components results in very different outputs for the reference Black Grama and Creosotebush core sites (Fig. 5a).The herbaceous component (which is derived from the relationship between NDVI and the reference 57-day antecedent rainfall series, ARain hv ) prevails in the grass-dominated reference site, while the shrub component (which is a function of the reference 145-day antecedent rainfall series, ARain s ) comprises the leading NDVI fraction in the shrub-dominated reference site.The annual sums of herbaceous and shrub NDVI components for the reference core sites show a strong linear agreement (R 2 ≥ 0.65; P < 0.001) with ground-based measurements of ANPP (Fig. 5b), while the remote-sensing ANPP estimations yield a root-mean-square error of 26 g m −2 (NRMSE 12 %, Fig. 5c).
Spatial projection of the reference NDVI-ANPP relationships across the 18 km 2 study area displays a contrasted distribution of mean 2000-2013 ANPP for herbaceous and shrub vegetation (Fig. 5d and e).Herbaceous ANPP is mainly distributed in the central and northeastern parts of the study site, contributing to > 80 % of total ANPP.Conversely, shrub ANPP is concentrated in the southwestern edge of the study area.

ANPP spatiotemporal dynamics and impact of seasonal precipitation on herbaceous and shrub primary production
Remotely sensed estimations of ANPP are significantly impacted by landscape type (F 3, 334 = 48.6,P < 0.01), with grass-dominated sites supporting in general higher levels of vegetation production (Fig. 6a).However, landscape-type effects are variable in time (landscape type × time interaction: F 14, 1515 = 57.2,P < 0.01).Year-to-year variability of ANPP is particularly large for the grass-dominated sites, which show higher levels of ANPP than the transition and shrub-dominated landscapes for highly productive years (Fig. 6a).For growing cycles with low primary production there are no significant ANPP differences or the differences are reversed, with shrub-dominated sites showing higher production than grass-dominated sites (e.g.2003-2004 cycle, Fig. 6a).Analysis of the temporal evolution of shrub contribution to total ANPP over 2000-2013 reflects significant (although very weak) positive correlations with time for the grass-and shrub transition landscapes (Fig. 6b).The same analysis at the individual pixel level, however, does not show any significant correlations between shrub contribution to total ANPP and time.
Exploratory analysis of the influence of seasonal precipitation on remotely sensed estimations of ANPP indicates different responses for herbaceous and shrub vegetation (Fig. 7).Herbaceous ANPP strongly correlates with monsoonal summer precipitation for all landscape types (Fig. 7a).The slope of the relationship between herbaceous ANPP and monsoonal summer (June-September) precipitation decreases for the shrub-transition and shrub-dominated landscapes.Conversely, shrub ANPP strongly correlates with both preceding non-monsoonal (October-May) and monsoonal summer (June-September) precipitation for all landscape types (Fig. 7b).
General linear model results confirm the exploratory observations of the relationships between remotely sensed estimations of ANPP and seasonal precipitation (Table 1).Model results identify both monsoonal summer precipitation (Rain SM ) and the interaction between Rain SM and landscape type as the most important contributors (effect size, η 2 > 10 %; P < 0.001) to the total variance comprised in ANPP data for herbaceous vegetation.Similarly, non-monsoonal summer precipitation (Rain PNM ) and monsoonal summer precipitation (Rain SM ) are identified as the leading contributors to shrub ANPP.

Vegetation-growth patterns and reference NDVI-rainfall metrics for herbaceous and shrub vegetation
Analysis of time series of NDVI provides important information on the dynamics of vegetation growth in drylands (Peters et al., 1997;Holm et al., 2003;Weiss et al., 2004;Choler et al., 2010).NDVI trends in the grass-dominated site show strong peaks centred in the summer season (Fig. 3a), which agrees with both field and remotely sensed observations of the dynamics of aboveground biomass for desert grasslands dominated by Bouteloua eriopoda and B. gracilis in the area (Peters and Eve, 1995;Huenneke et al., 2002;Muldavin et al., 2008;Notaro et al., 2010).For the shrub-dominated site, summer peaks in the NDVI signal are smaller, and for particular years both spring and late-autumn peaks can exceed summer greenness.Accordingly, the timing of plant growth for Larrea tridentata (which dominates the reference shrubland site) has been shown to vary from year to year, since this species has the ability to shift the temporal patterns of vegetation growth to take advantage of changes in resource availability (Fisher et al., 1988;Reynolds et al., 1999;Weiss et al., 2004;Muldavin et al., 2008).
The analysis of the relationships between NDVI and precipitation provides further insights into plant water-use patterns and, hence, on vegetation function and structure (Pennington and Collins, 2007;Veron and Paruelo, 2010;Notaro et al., 2010;Garcia et al., 2010;Forzieri et al., 2011;Moreno-de las Heras et al., 2012).Temporal trends in NDVI for the reference grass-and shrub-dominated SEV LTER sites are explained by antecedent (or preceding cumulative) rainfall amounts, reflecting the coupling of the history of plant-available soil moisture with vegetation growth (Fig. 3).Correlations between NDVI and precipitation indicate that plant-growth pulses for the grass-dominated site are associated with short-term antecedent rainfall (ARain hv series; 57 days optimal length, Olr hv ).For the shrub-dominated landscape, vegetation greenness shows a strong association with longer-term antecedent precipitation (ARain s series; optimal length 145 days, Olr s ), although, importantly, NDVI dynamics for this site also correlate with the 57-day cumulative rainfall series.Previous work on the analysis of NDVIrainfall relationships found similar variations in the length of the antecedent rainfall series that best explain the dynamics of vegetation greenness, suggesting that such differences result from site variations in dominant vegetation (Evans and Geerken, 2004;Prasad et al., 2007;Garcia et al., 2010).
Given the strong relationship between time-integrated NDVI values and ground-based ANPP estimations for our site (Fig. 5b), our herbaceous and shrub exploratory modelling results provide a biophysical explanation for the range of variations found in the NDVI-rainfall relationships (Fig. 1).The length of the cumulative precipitation series that optimizes the relationship between plant biomass and antecedent rainfall (Olr) appears to be a function of the characteristic water-use and plant-growth pattern of dryland vegetation, that are largely influenced by the plant-growth and mortality rates of vegetation (Fig. 1c).Vegetation growth and water use strongly differ for herbaceous and shrub life forms in drylands (Sparrow, 1997;Ogle and Reynolds, 2004;Gilad et al., 2007;Garcia et al., 2010), in which case plant biomass dynamics respond to short-term and long-term antecedent precipitation, respectively (Fig. 1a-b).Olr variations in the reference SEV LTER core sites may, therefore, be expressed as a function of the dominant vegetation types (Fig. 3): the strong and quick responses of greenness to short-term precipitation (ARain hv ) in the grass-dominated Black Grama core site characterize herbaceous growth for the area, while the slow responses of NDVI to mediumterm precipitation (ARain s ) in the shrub-dominated Cresotebush core site define the characteristic pattern of vegetation growth for shrubs in the ecotone.The high correlation be-tween ARain hv and NDVI values in the shrub-dominated Creosotebush core site (Fig. 3b) can be explained by the growth of non-dominant herbaceous vegetation (mainly annual forbs), which can be especially important during wet years (Muldavin et al., 2008;Baez et al., 2012) (Evans and Geerken, 2004;Munkhtsetseg et al., 2007;Garcia et al., 2010;Moreno-de las Heras et al., 2012).31.9 3 0.000 2.1 Abbreviations: ANPP r.sensing , remotely sensed annual net primary production; Rain PNM (Oct-May) , preceding non-monsoonal rainfall; Rain SM (Jun-Sep) , monsoonal summer rainfall; Rain LNM (Oct-Mar) , late non-monsoonal rainfall; LT, landscape type; ":", interaction terms; η 2 , eta squared (effect size).Notes: η 2 values in bold are > 10 % (effects that contribute more than 10 % to the total variance comprised in ANPP r.sensing ).

Spatial distribution and net primary production of herbaceous vegetation and shrubs
Our results indicate that the relationship between temporal series of remotely sensed NDVI and antecedent precipitation is highly sensitive to spatial differences in dominant vegetation (Fig. 4).The main PCA factor (explaining about 40 % variance in data) extracted using the annual NDVI responses (i.e. the Pearson's R coefficients) to the reference 57-and 145-day characteristic antecedent rainfall series (ARain hv and ARain s series, respectively) accurately discriminates the behaviour of herbaceous and shrub vegetation for the 18 km 2 study area (Fig. 4b-c), hence providing a robust approach for classifying landscapes as a function of the dominance of vegetation types using coarse-grained remotely sensed data (Fig. 4d).This parsimonious approach offers a practical alternative to other more complex remote-sensing methodologies for the analysis of the spatial distribution of vegetation types in mixed systems, such as spectral mixture analysis (SMA; Smith et al., 1990), which may be difficult to apply in this Chihuahuan case study since both the mixed nature and fine-grained distribution of vegetation in the area (patches of grass and shrubs are typically < 1 m 2 and 0.5-5 m 2 , re-spectively; Turnbull et al., 2010b) can impose serious drawbacks on the detection of reference spectral signatures for pure herbaceous and shrub vegetation using coarse-grained MODIS data.Implementing SMA-based approaches for the analysis of vegetation distribution and landscape classification in drylands using medium-and coarse-grained data is very challenging since it requires significant amounts of ancillary data (e.g.laboratory-based or field multi-date spectra for vegetation types) to solve data uncertainties generated by surface heterogeneity, which is often not feasible (Somers et al. 2011).
The relationships of vegetation greenness to ARain hv and ARain s also provide criteria for decomposing and transforming the NDVI signal into structural components of primary production for this study.Lu et al. (2003) applied seasonal trend decomposition to partition NDVI into (cyclic) herbaceous and (trend) woody vegetation in Australia.They assumed a long-term weak phenological wave and a strong annual response for determining the shrub and herbaceous components of vegetation, respectively.Our approach relies on the use of differences in biophysical properties of herbaceous and shrub vegetation related to the coupling between vegetation growth and precipitation for decomposing the NDVI signal, rather than apparent differences in the seasonality of vegetation greenness alone.As expected, signal decomposition outcomes indicate that the herbaceous component of the NDVI leads the temporal trends for the grass-dominated reference Black Grama core site, while the shrub component largely dominates the NDVI signal for the Creosotebush core site (Fig. 5a).
Although affected by data dispersion, the annual sums of decomposed NDVI strongly agree with field estimations of ANPP for herbaceous and shrub vegetation (R 2 ≥ 0.65, Fig. 5b), resulting in a small root-mean-square error for our remote-sensing ANPP estimates (26 g m −2 , NRMSE 12 %, Fig. 5c) that is within the lower limit of reported errors by other NDVI decomposition studies (for example, Roderick et al., 1999;DeFries et al., 2000, Hansen et al., 2002;Lu et al., 2003;with NRMSE ranging 10-17 %).Other dryland studies have found important levels of data dispersion when relating fine-grained field ANPP to coarse-scale NDVI values (Lu et al., 2003;Holm et al., 2003;Pennington and Collins, 2007;Veron and Paruelo, 2010).Major sources of data dispersion for this study are most likely associated with the high spatial variability of ANPP in the analysed systems.For instance, field estimations have shown that ANPP for both grass-and shrub-dominated Chihuahuan landscapes are affected by important levels of spatial variability, primarily due to the patchiness of vegetation cover (Huenneke et al., 2002;Muldavin et al., 2008).

Spatiotemporal dynamics of ANPP and impact of seasonal precipitation on herbaceous and shrub primary production
Cross-scale interactions between vegetation composition, individual plant characteristics, and climatic drivers (e.g.variations in precipitation amount and seasonality) have an important role in determining primary production patterns in arid and semi-arid ecosystems (Peters, 2002;Snyder and Tartowsky, 2006;Pennington and Collins, 2007;Notaro et al., 2010;Baez et al., 2013).Analysis of the spatiotemporal dynamics of ANPP in our ecotone indicates that grassdominated sites, although very importantly affected by yearto-year variability, generally support higher primary production than transition and shrub-dominated landscapes, particularly for wet years with high ANPP levels (Fig. 6a).This result is consistent with other shrub-encroachment studies which have found associations between shrub proliferation and ANPP reductions in dry North American grasslands (Huenneke et al., 2002;Knapp et al., 2008).
Our results suggest that primary production is differently controlled by seasonal precipitation for herbaceous and shrub vegetation across the 18 km 2 Chihuahuan Desert ecotone (Fig. 7, Table 1).Monsoonal summer precipitation (June-September) controls ANPP for herbaceous vegetation (Fig. 7a), while ANPP for shrubs is better explained by the preceding year's non-monsoonal (October-May) plus the summer monsoonal precipitation in the present year (Fig. 7b).Accordingly, field observations of ANPP for Chihuahuan landscapes found that grassland primary production is particularly coupled with monsoonal rainfall, while desert shrublands appear to be less dependent on summer precipitation (Fisher et al., 1988;Reynolds et al., 1999;Huenneke et al., 2002;Muldavin et al., 2008;Throop et al., 2012).
Differences in the distribution of rainfall types, soilmoisture dynamics, and rooting habits of dominant plant species may explain the variable impact of seasonal precipitation on herbaceous and shrub ANPP for the studied Chihuahuan landscapes.Monsoonal summer precipitation (July-September, approx.60 % annual precipitation) generally takes place in the form of high-intensity thunderstorms that infiltrate shallow soil depths (top 15-35 cm) (Snyder and Tartowsky, 2006).Summer soil-water resources for plant production are ephemeral and strongly affected by evapotranspiration, which typically reduces soil moisture to prestorm background levels in 4-7 days after rainfall (Turnbull et al., 2010a).C 4 grasses (Bouteloua eriopoda and B. gracilis), which dominate herbaceous vegetation in the analysed ecotone, concentrate active roots in the top 30 cm of the soil and intensively exploit ephemeral summer soil moisture for plant growth (Peters, 2002;Muldavin et al., 2008).Preferential spatial redistribution of runoff to grass patches following summer storms further enhances plant production for black and blue grama (Wainwright et al., 2000;Pockman and Small, 2010;Turnbull et al., 2010b).
Non-monsoonal precipitation (about 40 % of annual precipitation, primarily from November to February) typically falls in the form of long-duration low-intensity frontal rainfall that often percolates to deep soil layers (Snyder and Tartowsky, 2006).Larrea tridentata, the dominant C 3 shrub in the studied ecotone, has a bimodal rooting behaviour that facilitates the use of both shallow and deep soil moisture for plant production (Fisher et al., 1988;Reynolds et al., 1999;Ogle and Reynolds, 2004).Deep creosotebush roots (70-150 cm depth) may acquire winter-derived soil-water resources that are unavailable to grass species, while active roots near the surface (20-40 cm depth) may serve to access summer-derived shallow soil moisture for plant growth (Gibbens and Lenz, 2001).The observed reduction in summer rain-use efficiency of herbaceous vegetation for the shrub-transition and shrub-dominated landscapes (i.e.variations on the slope of the relationship between herbaceous ANPP and summer precipitation, Fig. 7a) suggests competitive effects of creosotebush for the use of shallow water sources, probably associated with the large spatial extent of near-surface active roots (the radial spread of which typically ranges between 2 and 6 m; Gibbens and Lenz, 2001).Alternative, landscape changes induced by shrub encroachment (i.e.increased runoff and erosion) may reduce the ability of grass patches to capitalize on horizontal redistribution of runoff for plant growth after summer storms (Wainwright et al., 2000;Turnbull et al., 2012;Stewart et al. 2014).
Conceptual and mechanistic models of vegetation change suggest that vegetation composition in arid and semi-arid landscapes is likely to be highly sensitive to climate change, and point at variations in the amount and distribution of precipitation as a major driver of shrub encroachment into desert grasslands (Peters, 2002;Gao and Reynolds, 2003;Snyder and Tartowsky, 2006).Overall our results agree with those findings and suggest that changes in the amount and temporal pattern of precipitation comprising reductions in monsoonal summer rainfall and/or increases in winter precipitation may enhance the encroachment of creosotebush into desert grasslands dominated by black and blue grama.Analysis of long-term rainfall series indicates that winter precipitation has increased during the past century in the northern Chihuahuan Desert, particularly since 1950, probably associated with the more frequent occurrence of El Niño-Southern Oscillation events for that period (Dahm and Moore, 1994;Wainwright, 2006).This pattern of precipitation change may be responsible, at least in part, for past increases in woody shrub abundance over desert grasslands in the southwestern USA (Brown et al., 1997;Snyder and Tartowsky, 2006;Webb et al., 2003).Our results suggest that shrub encroachment was not particularly active in the studied ecotone for the period 2000-2013 (Fig. 6b).Accordingly, Allen et al. (2008), in a recent study on creosotebush plant architecture and age structure, indicated that the most important pulses of shrub encroachment for this area took place between 1950 and 1970.Precise estimation of shrub cover applying seg-mentation methods in time series of high-resolution imagery could help to accurately determine the intensity of the shrubencroachment phenomenon under the present variability in precipitation for our grassland-shrubland ecotone.
Climate-change projections for the area suggest a general picture of increased aridity in the next 100 years, with increased evaporation due to higher summer temperatures, as well as increased drought frequency (Christensen and Konikicharla, 2013).The capacity of L. tridentata to switch between different soil-water sources (i.e.summer-derived ephemeral shallow soil moisture and more stable deep soilwater reserves derived from winter rainfall) and adapt the timing of vegetation growth to take advantage of changes in resource availability makes this C 3 shrub less susceptible to predicted increases in aridity than C 4 grasses that are strongly dependent on summer precipitation (Reynolds et al., 1999;Throop et al., 2012;Baez et al., 2013).Current increases in atmospheric CO 2 concentrations may also contribute to reduce the competitiveness of C 4 grasses for the use of soil-water resources against C 3 desert shrubs (Polley et al., 2002).Remaining desert grasslands in the southwestern USA may, therefore, be increasingly susceptible to shrub encroachment under the present context of changes in climate and human activities.

Conclusions
In this study we applied a new analytical methodology for the study of the organization and dynamics of vegetation at a grassland-shrubland Chihuahuan ecotone with variable abundance of grasses (primarily Bouteloua eriopoda and B. gracilis) and shrubs (mainly Larrea tridentata), based on the exploration of the relationship between time series of remotely sensed vegetation greenness (NDVI) and precipitation.Our results indicate that the characteristics of the NDVI-rainfall relationships are highly dependent on differences in patterns of water use and plant growth of vegetation types.In fact, NDVI-rainfall relationships show a high sensitivity to spatial variations on dominant vegetation types across the grassland-shrubland ecotone, and provide biophysically based criteria to study the spatial distribution and dynamics of net primary production (NPP) for herbaceous and shrub vegetation.The analysis of the relationship between NDVI and precipitation therefore offers a powerful methodology for the study of broad-scale vegetation shifts comprising large changes in the dominance of vegetation types in drylands using coarse-grained remotely sensed data, and could be used to target areas for more detailed analysis and/or the application of mitigation measures.
Analysis of remotely sensed NPP dynamics at the grassland-shrubland ecotone reflects a variable performance of dominant vegetation types.Herbaceous production is synchronized with monsoonal summer rainfall, while shrub NPP shows a flexible response to both summer and winter pre-cipitation.Overall our results suggest that changes in the amount and temporal pattern of precipitation (i.e.reductions in summer precipitation and/or increases in winter rainfall) may intensify the shrub-encroachment process in the studied desert grasslands of the southwestern USA, particularly in the face of predicted general increases in aridity and drought frequency for the area.
The Supplement related to this article is available online at doi:10.5194/bg-12-2907-2015-supplement.

Figure 2 .
Figure 2. Study area: (a) location of the Sevilleta National Wildlife Refuge (SNWR) and distribution of major New Mexico biomes, (b) regional location of the study area (McKenzie Flats, SNWR), and (c) detailed location of the study site (18 km 2 area) and general view of the reference SEV LTER Black Grama (right) and Creosotebush (left) core sites.Map (a) follows the Sevilleta LTER classification of New Mexico biomes (Sevilleta LTER, http://sev.lternet.edu/content/new-mexico-biomes-created-sevlter).Source for background image in panels b) and (c): 2009 National Aerial Imagery Program (USDA Farm Service Agency).

Figure 3 .
Figure 3. Reference NDVI-rainfall relationships at the SEV LTER Black Grama and Creosotebush core sites: (a) 2000-2013 MODIS NDVI time series for the core sites and (b) strength of the NDVI-rainfall relationship (Pearson's R correlation) for the core sites using different lengths of rainfall accumulation (maximum correlations, R max , for the annual cycles of vegetation growth are shown together with the 2000-2013 mean trend; detailed correlograms for each growing cycle can be found in Supplement Fig. S1).R values above the dotted grey line are significant at P < 0.05.ARain hv and ARain s lines in panel (a) represent the antecedent rainfall series that best correlate with the NDVI series for the Black Grama and Creosotebush Core sites (i.e.time series of precedent rainfall with rainfall accumulation lengths Olr hv and Olr s , respectively).Reference Olr hv and Olr s values in panel (b) represent the optimal rainfall accumulation lengths for herbaceous vegetation (57 days) and shrubs (145 days), respectively.

Figure 4 .
Figure 4. Principal component analysis (PCA) of the NDVI-rainfall correlation coefficients for the herbaceous-and shrub-specific antecedent rainfall series ARain hv and ARain s (57-and 145-day cumulative rainfall series, respectively) and resulting landscape-type classification across the 18 km 2 study area: (a) PCA projection of cases (MODIS pixels), (b) PCA projection of variables (per growing cycle NDVIantecedent rainfall correlation scores), (c) landscape-type classification (GD, grass-dominated; GT, grass-transition; ST, shrub-transition; and SD, shrub-dominated landscapes) as a function of the relationship between PCA factor 1 and field-estimated vegetation dominance for a reference set of 27 control points, (d) spatial distribution of landscape types in the study area, and (e) general view and characteristics of the landscape types.MODIS pixel locations for the ground control points are highlighted in panel (a).Vector labels in panel (b) indicate the dates of the yearly cycles of vegetation growth (April-March).Source for background image in panel (d): 2009 National Aerial Imagery Program (USDA Farm Service Agency).

Figure 5 .
Figure 5. NDVI decomposition and transformation into partial annual net primary production (ANPP) components for herbaceous and shrub vegetation: (a) decomposed NDVI time series of herbaceous and shrub vegetation for the reference SEV LTER Black Grama and Creosotebush core sites, (b) relationships between field ANPP and the (per growing cycle) annual integrals of herbaceous and shrub NDVI components for the SEV LTER core sites, (c) remotely sensed ANPP estimates against field ANPP determinations (root-mean-square error, RMSE, and normalized error, NRMSE, of the estimates are shown within the plot), (d) remotely sensed ANPP estimations of herbaceous and shrub vegetation (mean for the 2000-2013 series), and (e) herbaceous and shrub contribution to total ANPP (mean for the 2000-2013 series) across the 18 km 2 study area.

Figure 6 .
Figure 6.Spatiotemporal dynamics of remotely sensed ANPP: (a) ANPP differences between landscape types (grass-dominated, grasstransition, shrub-transition, and shrub-dominated landscapes) over 2000-2013 and (b) 2000-2013 temporal variations of the shrub contribution to total ANPP for the different landscape types (Pearson's R correlations of shrub ANPP contributions with time).Different letters in panel (a) for each cycle of vegetation growth indicate significant differences between landscape types at P < 0.05 (tested using repeatedmeasures ANOVA and post hoc Tukey HSD tests).Dotted lines in panel (b) represent weak (R < 0.40) correlations.Displayed correlations are significant at P < 0.05.Numbers in plot (b) indicate correlation coefficients.
. Similarly, Moreno-de las Heras et al. (2012) found, in dry open shrublands of central Australia (Olr s values about 220 days), the emergence of secondary Olr hv metrics on the study of local NDVI-rainfall relationships (approx.85-day antecedent rainfall length) caused by the growth of non-dominant herbaceous vegetation.Overall, Olr values determined for herbaceous and shrub vegetation in this work are in agreement with the range of characteristic antecedent rainfall series reported in other studies to best describe green biomass dynamics for arid and semi-arid grasslands (1-3 months) and woody shrublands (4-8 months)

net/12/2907/2015/ Biogeosciences, 12, 2907-2925, 2015 observed
seasonality of the original NDVI time series in the decomposed signals for herbaceous and shrub vegetation, the predicted weights (or percentage contributions) of the fitted vegetation components were reassigned to the NDVI levels of the original time series, obtaining the final NDVI components for herbaceous vegetation and shrubs (C hv , and C s , respectively).Computed soil background baseline (C bs ) plus the partitioned NDVI components for herbaceous vegetation (C hv ) and shrubs (C s ) totals the original NDVI levels of the temporal series for any point in time and space.
. In order to preserve the www.biogeosciences.

Table 1 .
Main effects and interactions of seasonal precipitation (preceding non-monsoonal rainfall, October-May; monsoonal summer rainfall, June-September; late non-monsoonal rainfall, October-March), and landscape type (four levels: grass-dominated, grass-transition, shrub-transition, and shrub-dominated landscapes) on remote-sensing-estimated annual (per growing cycle, April-March) net primary production for herbaceous vegetation and shrubs.