A statistical approach for predicting grassland degradation in disturbance-driven landscapes

The relationship between fire and long-term trends in tallgrass prairie vegetation was assessed at Fort Riley and Konza Prairie Biological Station (KPBS) in Kansas. Linear trends of surface greenness were previously estimated using BFAST and MODIS MOD13Q1 NDVI composite images from 2001 to 2010. To explain trends, fire frequency and seasonality (fire regime) was determined and each site was divided into spatial strata using administrative or management units. Generalized linear models (GLM) were used to explain trends by fire regime and/or stratification. Spatialized versions of GLMs were also computed address unexplained spatial components. Non-spatial models for FRK showed fire regime explained only 4% of trends compared to strata (7-26%). At KPBS, fire regime and spatial stratification explained 14% and 39%, respectively. At both sites, improvements in performance were minimal using both fire and strata as explanatory variables. Model spatialization resulted in a 5% improvement at FRK, but with weak spatial structure in the residuals, and was not necessary at KPBS as the existing stratification most of the spatial structure in model residuals. All models at KPBS performed better for each explanatory variable and combination tested. Fire has only a marginal effect on vegetation trends at FRK despite its widespread use as a grassland management tool to improve vegetation health, and explains much more of the trends at KPBS. Analysis of predictors from spatial models with existing stratification yielded an approach with fewer strata but similar performance and may provide insight about additional explanatory variables omitted from this analysis.


Introduction
The US Department of Defense is one of the largest land stewards in the nation [1]. Readily available and easily accessible military lands are required to support an array of training and testing activities by the armed forces of the United States. However, training activities ranging from vehicle-based maneuvers and indirect fire weapons to foot traffic are known to negatively affect environmental conditions [2]- [8].
Maintaining a land base that supports safe and realistic training operations is a significant challenge impacting the mission readiness of US military forces. A number of federal laws (e.g., Sikes Act of 1960, National Environmental Policy Act of 1969) and military regulations outline specific requirements for natural resource conservation and management plans and activities, including that sustainable use should be a programmatic goal within the US Army [9]- [11]. Given the geographic extent of training sites and their diverse settings, including deserts, tropical and temperate forests, semi-arid grasslands, and arctic landscapes, defining and agreeing on metrics to measure "sustainable use" is difficult. This difficulty is compounded when superimposing spatially and temporally unique military training disturbances on very different ecosystems.
As in temperate grasslands and savannas across the world, fire frequency and timing can have significant impacts on the spatial pattern of plant productivity, vegetation structure, and nutrient cycling in the Flint Hills ecoregion of North America [12]- [15]. For example, frequent burning typically results in prairie dominated by C4 grasses, with lower species richness and diversity, and higher average peak-season aboveground biomass. When considered with other factors such as weather and climate, topographic position, grazing and other anthropogenic disturbances, the impact of fire on grassland dynamics is even more complex [16]- [20].
In this study, promoting sustainable grasslands is considered the primary objective with sustainability being defined using results from a long-term analysis of vegetation trends derived from remotely-sensed estimates of normalized difference vegetation index (NDVI) data. The influence of fire and land use practices on grassland vegetation trends is assessed to provide land managers with site-specific feedback about progress towards their training land sustainability goals.
Spatial statistical tools are used to analyze spatiotemporal trends in vegetation and fire regime using a time series of multisource remote sensing images over the period 2001-2010. This approach has previously been used by [21] who studied the role of fire on savanna vegetation in Madagascar. As in this previous work, results from non-spatial and spatial statistical analyses with spatial stratification are compared to identify the method that best explains vegetation dynamics. The spatial generalized linear model (GLM) is a classical tool in the spatial statistics toolbox. Here, a hierarchical coupled logistic multinomial model with a spatial auto normal random component was used. However, because the spatial stratifications used in this study were not always constructed from areas with homogeneous environmental conditions, land use practices, or human impacts, a new approach was proposed and implemented to construct improved and simplified spatial stratifications that were better adapted to the response variable. This approach may benefit other research efforts in disturbed environments where the landscape exhibits similar heterogeneity. Though implemented in a basic form here, future stratifications could be defined using hierarchical spatial models that operate in a manner similar to non-supervised image classification.
This research was conducted at Fort Riley, a US Army installation located in northeastern Kansas, within the framework of the Army's Integrated Training Area Management (ITAM) program. In conjunction with installation ITAM staff, we seek to demonstrate methods for continuously monitoring, assessing and identifying trends in key indicators of military training land sustainability over appropriate timescales, rapidly analyze data, and ensure information is available when needed by military land managers. Though military lands are the focus here, these objectives are equally applicable for effective management of other grasslands and savannas in the western US, South America, Europe, Africa, and Asia.
A. Jacquin et al.

Study Area
Fort Riley is a 41,128 ha U.S. Army installation in northeast Kansas (39 11'N, 96 48'W), on the Kansas River, between Junction City and Manhattan within Geary, Riley and Clay counties. The installation is located within both the Flint Hill ecoregion [22] and the Central Great Plains Winter Wheat and Range land resource region [23]. Soils are productive and typically consist of clay uplands combined with loamy uplands, limy soils, and loamy lowlands. Vegetation at Fort Riley is dominated by grasslands (81%) followed by woodlands (16%) and shrublands (3%). Dominant plant species in grassland areas include big bluestem (Andropogon gerardii), Indiangrass (Sorghastrum nutans), switchgrass (Panicum virgatum), and little bluestem (Schizachyrium scoparium). These characteristic tallgrass prairie species are actively managed using prescribed fire during the fall, winter, and spring seasons. However, wildfires resulting from military training may occur during any season.
The installation is subdivided into 103 administrative units within which military training occurs. Disturbances from military training on Fort Riley include a variety of on-and off-road field maneuvers (including tracked and wheeled combat vehicle operations), mortar and artillery fire, small arms fire, and aircraft flights. For the past four decades, the majority of combat vehicle maneuvers take place in the western part of the installation.
The Konza Prairie Biological Station (KPBS) was included in this analysis to contrast vegetation change measured at Fort Riley with that at a relatively natural grassland site. KPBS (3487 ha) is owned by the Nature Conservancy and operated by the Division of Biology at Kansas State University (http://kpbs.konza.ksu.edu). One of the U.S. National Science Foundation's Long Term Ecological Research Sites, KPBS has similar vegetation, soils and climate due to its close proximity (less than 10 km) to Fort Riley.
It is assumed that any differences observed between the two sites are caused by varying anthropogenic activities and land management practices. Though not impacted by military training, vegetation at KPBS is influenced by other human actions including prescribed burning. The KPBS site is subdivided into 48 sub-watershed areas which serve as experimental units where specific combinations of fire (by both season and frequency) and grazing (including both domestic cattle and native bison) are applied. Figure 1 shows the general vegetation types and spatial stratification used by land managers for the two study areas.

Data
Data used in this analysis consists of vegetation activity change and fire regime data. All procedures used to process data for Fort Riley were also applied for KPBS. The study period is limited to 2001-2010 to match the study duration of the prerequisite vegetation trend analysis work previously completed for Fort Riley and KPBS and the lack of data needed to characterize the fire regime (high spatial resolution images) after 2010.

Vegetation Activity Change Data
A vegetation change indicator was developed to characterize long-term trends in activity over the study period ( Figure 2). The indicator was obtained using a temporal decomposition method applied to a time series of MODIS 16-day maximum value composite NDVI images (MOD13Q1) with a spatial resolution of 250 m between January 2001 and December 2010. The method, described in [25] consists of three steps: 1) NDVI time series decomposition using BFAST based on LOESS [26]; 2) analysis of the BFAST trend component using linear regression; and 3) identification of significant positive or negative trend slopes using a Student's t-test. This method has already been applied to the Fort Riley and KPBS study areas [27].

Fire Regime Data
A fire regime indicator was also developed to characterize burning activities according to season and frequency by analyzing a time series of multisource remote sensing data for the period 2001-2010. The primary remote sensing data were high spatial resolution (HSR) images most appropriate for delineating burn scars at local and regional scales [28]. A total of 83 Landsat 5 TM images with a spatial resolution of 30 m were acquired and burn scars extracted using a combination of unsupervised classification of the near-infrared band with CAPI techniques [29].
Due to persistent cloud contamination during the study period, it was impossible to acquire and classify HSR images at a frequency sufficient to adequately identify burn scars. This was especially a problem in 2004-2005, To fill these gaps, 135 additional MODIS Burned Area Product (MCD45) images were used to complement the HSR time series. Produced by the MODIS Fire Team at the University of Maryland (http://modis-fire.umd.edu/pages/BurnedArea.php), the Burned Area Product provided monthly images of burned pixels at a 500 m spatial resolution. Though the spatial resolution is coarse, preventing detection of small and fragmented burns, the MCD45 images were a practical and inexpensive source of data to ensure a continuous and reliable source of fire information within the time series.
Next, monthly burned area maps were created for each year and used to calculate annual burned area maps. To ensure compatibility with the vegetation change indicator, monthly burned area maps were overlaid with a grid of the same spatial resolution. Only grid cells with a minimum burned area of 80% were retained for further analysis. Summing the annual maps over the study period yielded fire frequency (Figure 3(a)). Annual burned area maps were also created to identify when fires occurred (e.g., seasonality) (Figure 3(b)). Three seasons were defined: Fall (September-October-November), winter (December-January-February), and spring (March-April-May). In determining both fire frequency and seasonality (i.e., fire regime), it was assumed that an individual pixel could burn only once per year with the timing of that burn defining the season. With information on both fire frequency and seasonality, a fire regime indicator consisting of 13 classes was defined.

Statistical Model Development
Modeling the relationship between explanatory and response variables is a fundamental activity encountered in statistics. Multiple linear regression is used to investigate the relationship between a single response variable and several explanatory variables. Often, the response is simply a designation of one or two possible outcomes (i.e.,  a binary response). In this study, the possible outcomes are "degraded" or "non-degraded" vegetation activity and this binary response is examined in relation to fire regime and/or other anthropogenic activities represented by stratification using non-spatial and spatial generalized linear models (GLM).
The grassland vegetation activity change map was first overlaid with the fire regime data. Next, GLM models were computed on a per-pixel basis to study how well fire regime and/or stratification explained degraded vegetation activity. In the first non-spatial model, where independence between observations was assumed, a classical estimation method was used. Due to the nature of the response variable, the model was a binomial GLM with a probability of success depending on the characteristics of each observation site. The link between this probability and the linear predictor based on individual characteristics is the logit and the estimation was done by maximum likelihood [30].
Because correlation between observations was suspected, a second model accounting for spatial correlation was applied. In this model, an existing spatially random field was assumed which takes into account what is not explained in the independent model. Given this unobserved random field, the model adds to the linear predictor the value of the random field at the observation site. The maximum likelihood method, however, is now problematic because the likelihood cannot be derived analytically and must be approximated. Instead, the Bayesian Markov Chain Monte Carlo method based on Metropolis dynamics was used to estimate model parameters. The estimates are then taken as the posterior mean [31].
Both types of models (non-spatial and spatial) were evaluated with two metrics. The first was the Akaike Information Criterion (AIC) which summarizes the tradeoff between model accuracy and complexity [32] [33].
Low AIC values for a model indicate superior explanatory power. The second metric was a ratio equal to 1 minus the residual deviance divided by null deviance (1 -(residual deviance / null deviance)). If the ratio is low then, at the level of the study area, other explanatory variables not included in the model should be considered.

Definition of Tested Models and Analysis Procedures
Five models using different explanatory variables were tested and their ability to explain degraded vegetation activity assessed. The initial assumption was that vegetation degradation would be explained adequately by only the fire regime variable (model 1), or only stratification (model 2), or a combination of fire regime and stratification (model 3). However, the original stratification, especially that at Fort Riley, was not specifically designed for use in monitoring grassland vegetation dynamics. To address this, a simplified stratification more adapted to the response variable was also tested. This simplified stratification was developed by analyzing the distribution of the model predictor based only on the original stratification (model 2) as a function of the stratification units. Two additional models were then developed. One only with simplified stratification (model 4) and one with fire regime and simplified stratification (model 5).
Comparing model 1 vs. model 2 and model 2 vs. model 3 provided results sufficient to assess the effect of fire regime on grassland vegetation degradation. The comparison of model 2 to model 4 or model 5 allowed evaluation of the simplified stratification and its ability to better explain observed vegetation degradation as well as its potential for serving as a spatial framework for monitoring grassland condition in the future.

GLM Validation Procedure
Each GLM model generates a probability map for vegetation degradation. Model performance is then validated using two different reference datasets. The first is the vegetation degradation class resulting from the previously mentioned MODIS NDVI trend analysis. The second is a new HSR vegetation degradation class from a change analysis performed using two Landsat 5 TM images acquired near the beginning (08/05/2001) and end (08/30/ 2010) of the study period. This method was recommended by [34], has been proven well adapted for grassland ecosystems [35] [36], and has been previously used to validate vegetation change activity for the Fort Riley and KPBS study areas [27].
The HSR dataset consists of a raster with the same spatial resolution as the vegetation activity change map (250 m) in which each grid cell corresponds to one of two classes (degraded or non-degraded) based on the interpretation of NDVI change derived from the Landsat 5 TM images. Since the GLM model estimates only a probability of vegetation degradation, only the NDVI change class for degraded vegetation is retained during this validation step.
A contingency matrix was constructed with rows representing the reference data (Landsat 5 NDVI change class for degraded vegetation or MODIS NDVI trend class for degraded vegetation) and columns with two classes of vegetation degradation probabilities (>0.50 and 0.50).

Results
The ratio (1 residual deviance/null deviance) * 100 was used to determine how well explanatory variables in the five GLM models explained variability in the vegetation activity change indicator. Associated AIC values were also calculated to characterize the relative explanatory power of each model ( Table 1).

Impact of Fire on Vegetation Degradation
To assess the effect of fire regime on vegetation activity change, results from GLM models with only the fire regime variable (model 1) were compared with those of GLM models using only stratification (model 2). Fire alone explains very little at Fort Riley (4% -11%) versus to stratification alone (26% -31%). The marginal effect of fire on vegetation dynamics is confirmed when also considering the results of the GLM model with fire and stratification (model 3). The combination of fire regime and stratification does not improve the explanatory power of the model (27% -31%). This same applies to KPBS, except that the effect of fire on vegetation activity change is higher than Fort Riley (14% -29%).
At Fort Riley, regardless of the GLM model, using a spatial approach in the statistical analysis is needed. Spatial versions of the GLM models provide a 4% -7% improvement in performance compared to non-spatial models with the added benefit of systematically lower AIC values. The situation at KPBS is slightly different. Model spatialization is not needed if the analysis is done with the original stratification, as the percentage of vegetation activity change explained and AIC values are nearly the same with both non-spatial and spatial models. At KPBS, the original stratification already accounts for the type and distribution of explanatory variables in the models. However, when modifying the original stratification (model 4), spatialization provides a 9% improvement in model performance with a decrease in AIC values from 577 to 511.

Assessment of Simplified Stratification Model Performance
By analyzing the range of predictor values for GLM models with only stratification (model 2) for each of the original management units (103 training areas at Fort Riley and 48 sub-watersheds at KPBS), a simplified stratification was proposed. For both study areas, this simplification consists only of five strata obtained after reclassification of training areas or sub-watersheds.
Results from the GLM model with only original stratification (model 2) were compared to those of the GLM model with only simplified stratification (model 4). At Fort Riley and KPBS, simplification explains nearly the same amount of vegetation activity change but with only five strata. As already seen with the original stratification, adding the fire regime variable to simplified stratification in a GLM model (model 5) does not improve the ability to explain vegetation activity change. Therefore, the GLM model with only simplified stratification was considered most appropriate.

Validation of GLM Model with Simplified Stratification
The vegetation degradation probability maps for Fort Riley and KPBS generated by the GLM model with only simplified stratification are shown in Figure 4. Black and dark gray areas correspond to the highest probabilities of vegetation degradation as explained by the five new strata, while sites in light gray or white suggest a low probability of vegetation degradation. Evaluation of the contingency matrix helps validate the vegetation degradation probability maps from the selected GLM model ( Table 2). Agreement is high at Fort Riley (59% -76%) and moderate to relatively high at KPBS (40% -62%) based on the reference data used. These validation results suggest the new simplified stratification is capable of highlighting areas of degradation, especially at Fort Riley. Identifying concentrations, or hot spots, of vegetation degradation can focus the attention of military land managers on specific areas where anthropogenic activities are having a negative impact on vegetation activity.

Discussion
Results from this analysis provide two major advances. First, from a thematic point of view and contrary to what   was expected, fire regime did not adequately explain patterns of vegetation activity change. At both sites, based on comparison of models with only original stratification (model 2) with original stratification and fire (model 3), vegetation degradation is better explained by the original stratification. This is especially true at Fort Riley and points to disturbances other than fire as the driver of local grassland trends. Second, from a methodological point of view, the soundest approach for identifying a model that best explained vegetation activity change involved first using a model with original stratification to form a new stratification more adapted to the response variable using a spatial approach. At Fort Riley, where simplified stratification is very informative (almost 3 times more than fire regime only), the five new strata were used to characterize the probability of vegetation degradation provided by the GLM model (Figure 5(a)). Strata 4 and 5, respectively, have an average vegetation degradation probability of 0.79 and 0.75. This is significantly higher than that obtained for strata 1, 2, and 3.
The role of military training activities was also investigated. High intensity training associated with mechanized military maneuvers has been cited as the cause of increased bare soil, reduced plant cover, soil compaction, and compositional shifts in plant communities [9] [37]- [40]. A map of military training intensity was created using a combination of expert knowledge and information from related literature (Denker, pers. comm.; Hut-chinson pers. comm.; [41]). Training areas were then classified using two levels (high and low) of training intensity (Figure 6).
In strata 4 and 5, the percentage area classified as having high military training intensity is almost double that of strata 1, 2, and 3 (Figure 5(c)). Strata 4 and 5 also represent areas where vegetation degradation identified by MODIS NDVI trend analysis was most concentrated, with 80% of strata 4 and 85% of strata 5 being classified as degraded (Figure 5(b)). These same strata have the highest probabilities of vegetation degradation as predicted by the GLM model.
Though fire appears to have only a minimal impact on vegetation dynamics based on GLM model results, fire regime was also analyzed within the new simplified strata. All strata were almost entirely burned (87%) during the 2001-2010 study period at predominantly low or very low frequency (i.e., no more than one fire every three years). In strata 5, the zone with the highest level of vegetation degradation (85%) and high intensity training (43%), fire typically occurred in spring (70%). In strata 1, where vegetation degradation and training intensity were both low (0.4% and 16%, respectively), fire occurred almost equally across the spring, fall, and winter seasons (Figure 5(d)).
Spring burning in the Flint Hills ecoregion removes senescent vegetation from the previous growing season to promote growth of new grasses. Annual spring fires favor C4 vegetation over C3 species, but also cause longterm increases in bare ground [42] [43]. Native C4 grasses in Kansas are more resilient to, and recover more quickly from, damage caused by military vehicles during off-road maneuvers, while areas dominated by C3 species tend to be more sensitive training [44]. It is important to note that over a third of the annual precipitation at Fort Riley also occurs during the spring months. That spring fires are a common feature of the strata also experiencing the greatest percentage of vegetation degradation and highest training intensities suggest that the temporal correspondence of fire, rainfall, and military training are contributing to increased vegetation damage and, perhaps, soil erosion. Considering KPBS, strata 5 had an average vegetation degradation probability of 0.69% and 61% of the area is classified as degraded based on MODIS NDVI trend analysis. Strata 1 has a very low probability of degradation (0.03) and 99% of the area is classified as non-degraded. Both strata 5 (75%) and 1 (90%) were almost entirely burned during the study period with spring the most common season of burning (96% for strata 5 and 100% for strata 1). The main difference between the strata is fire frequency and other land uses. In strata 1, 56% of the area is affected by annual burns without any grazing pressure and 36% by one fire every four years with bison grazing conducted annually between May and October. In strata 5, fire frequency is variable but well balanced (54% low and 46% high frequency), but 40% of the area is under intensive grazing by cattle (northeast) and 37% is used for trail-based recreation (north central). Spring fires at KPBS seem to be a stabilizing factor for the grassland ecosystem with vegetation degradation due mainly to grazing activities combined with fire frequency and human disturbances.

Implications
There appears to be contrasting approaches to fire management at Fort Riley. In areas where training intensities are high, spring fires are perhaps used to promote grassland health where military training disturbances are significant and frequent. In other training areas, no fire regime dominates. While the literature and results here from KPBS support spring burning as a good practice from a grassland resource perspective, when combined with frequent and intense training during the wettest time of the year, it may actually be counterproductive and amplify vegetation degradation. At minimum, the influence of military training on vegetation dynamics is important and serves to minimize any potential benefit realized from a typical Flint Hills fire regime.
To confirm such conclusions, more detailed and temporally consistent training data are required. While these data do exist, they are embedded within complex US Army database management systems and are difficult to extract for use in natural resource studies. Until such a time when training data are routinely available for land managers, the new simplified strata presented here may be used as a spatial guide to modify training schedules and prevent further declines in vegetation condition where degradation has been predicted and measured. Land managers suffering from similar data deficiencies in land usage, such as off-road vehicle use in public recreation areas, may also find utility in defining custom strata for their sites.